Spatial and temporal heterogeneity of the density of Borrelia burgdorferi-infected Ixodes ricinus ticks across a landscape: A 5-year study in southern England.

The density of Borrelia burgdorferi-infected Ixodes ricinus nymphs (DIN) was investigated during 2013-2017 across a Lyme disease-endemic landscape in southern England. The density of nymphs (DON), nymph infection prevalence (NIP), and DIN varied across five different natural habitats, with the highest DIN in woodland edge and high biodiversity woodlands. DIN was significantly lower in scrub grassland compared to the woodland edge, with low DON and no evidence of infection in ticks in non-scrub grassland. Over the 5 years, DON, NIP and DIN were comparable within habitats, except in 2014, with NIP varying three-fold and DIN significantly lower compared to 2015-2017. Borrelia garinii was most common, with bird-associated Borrelia (B. garinii/valaisiana) accounting for ~70% of all typed sequences. Borrelia burgdorferi sensu stricto was more common than B. afzelii. Borrelia afzelii was more common in scrub grassland than woodland and absent in some years. The possible impact of scrub on grazed grassland, management of ecotonal woodland margins with public access, and the possible role of birds/gamebirds impacting NIP are discussed. Mean NIP was 7.6%, highlighting the potential risk posed by B. burgdorferi in this endemic area. There is a need for continued research to understand its complex ecology and identify strategies for minimizing risk to public health, through habitat/game management and public awareness.

INTRODUCTION

Over recent years there has been an increasing number of human cases of Lyme disease in England (Tulloch et al., 2019). The causative agents are the spirochetes Borrelia burgdorferi sensu lato. These bacteria cause vector‐borne wildlife zoonoses, cycling between wildlife, through transmission by ixodid ticks, principally the deer/sheep tick, Ixodes ricinus (Linnaeus, 1758) (Acari: Ixodidae). Various genospecies of B. burgdorferi s.l. occur in the UK, predominantly B. afzelii, B. garinii, B. valaisiana and B. burgdorferi sensu stricto (Bettridge et al., 2013; Cull et al., 2021; Hansford et al., 2015; Hansford et al., 2021; James et al., 2013; Millins et al., 2016). In England, B. afzelii cycles between ticks and small mammals, such as wood mice (Apodemus sylvaticus), yellow‐necked mice (Apodemus flavicollis) and bank voles (Myodes glareolus) (Hanincová et al., 2003; Hoodless et al., 1998; Humair et al., 1993), with B. garinii and B. valaisiana cycling between ticks and birds, such as woodland birds like blackbird (Turdus merula) and robin (Erithacus rubecula) (Hubálek et al., 1996), and game birds such as pheasant (Phasianus colchicus) (Craine et al., 1997; Hoodless et al., 1998, 2003). Other mammals such as red squirrels (Sciurus vulgaris) and grey squirrels (Sciurus carolinensis) further contribute to transmission cycles (Craine et al., 1995; Millins et al., 2015). Livestock and deer act as key tick hosts, but are not directly involved in the transmission of Borrelia (Jaenson & Talleklint, 1992). Humans are incidental hosts, as they become infected via a bite from an infected tick, but do not contribute to the natural transmission cycle.

Our ability to understand the public health risks from Borrelia infected ticks is largely dependent upon understanding the transmission cycle of the spirochete within nature, the biotic and abiotic factors that impact it, as well as human exposure to ticks. For the latter, spending time in tick habitat, either in the countryside or in urban and peri‐urban areas, is a major risk factor for acquiring tick bites (Cull et al., 2018; Hansford et al., 2017), some of which can result in Borrelia transmission. The variability in nymphal tick density (‘density of nymphs’; DON) and Borrelia infection prevalence in nymphal ticks (‘nymphal infection prevalence’; NIP) within different habitats are less well understood. These factors are critical in determining the public health outcomes, i.e. whether a tick bite could potentially lead to infection and the possible clinical implications thereafter. DON and NIP begin to provide an assessment of the potential hazard of the risk of Borrelia infection. These values can be combined to calculate the ‘density of infected nymphs’ (DIN) in the landscape and to observe how this varies over time and between habitats, potentially offering a better overall assessment of risk.

It seems likely that Borrelia infection within ticks, and the dominance of one genospecies over another may vary in space and time across a landscape, depending on the varying abundance of different animal hosts, and their relative roles as tick hosts and reservoirs of infection. The population densities of these mammals and birds are in turn determined by biotic and abiotic factors, such as habitat, weather and climate, and various predator/prey interactions (Medlock et al., 2018; Medlock & Leach, 2015). There are also anthropogenic factors that might influence these transmission cycles, such as land and game management (Hoodless et al., 1998). It is these factors that also impact the density of ticks; whether they are infected with Borrelia or not.

It cannot be assumed that all ‘tick suitable habitats’ support the same density of ticks, or that all tick habitats will contain ticks with the same Borrelia infection prevalence, nor with the same proportion of B. burgdorferi genospecies. The mosaic habitat structure of a landscape, its varying connectivity and the variability of wild animals in space and time are likely to lead to variation in both DIN and the relative proportions of genospecies infection. However, it is possible that there are patterns that can be identified and possibly extrapolated. For example, does Borrelia infection in ticks vary between the interior and edge of woodland and does this vary between different types of woodland? Does Borrelia infection in ticks vary between woodland and grazed grassland habitats and does the presence of scrub in grazed grassland affect tick density and Borrelia infection? Do tick infections and the density of infected ticks vary between years at the same site? Indeed, does this variability fluctuate consistently across a range of habitats within a region as a result of changes in regional weather?

This study investigates DON, NIP and DIN across several habitats at a regional scale in southern England to assess whether there is any heterogeneity of Borrelia infected I. ricinus nymphs at a landscape scale. This study builds on a recent snap‐shot survey on variation in NIP in recreational sites across England and Wales during 2014–2019 (Cull et al., 2021), by focussing on one of those specific locations (predominantly the Cranborne Chase Area of Outstanding Natural Beauty, and northern part of New Forest National Park) in south Wiltshire during 2013–2017, which was identified as an area with the highest NIP. This area is reported to have a high incidence of Lyme disease (Tulloch et al., 2019) and surrounds the city of Salisbury, where a high prevalence of Borrelia infected ticks have been reported (Hansford et al., 2017). Developing an understanding of spatial and temporal heterogeneity of Borrelia infection in I. ricinus is the first step to further research that investigates the ecological drivers for Borrelia infection in ticks in southern England.

METHODS

Ecological surveys

An area of South Wiltshire (40 km × 25 km) bounded by the southern edge of the military training ground of Salisbury Plain to the north, and the counties of Dorset and Hampshire to the south, east and west were chosen as a study site. The area is dominated by Cretaceous and Jurassic bedrocks, with predominantly calcareous grasslands and woodlands in the south, Cretaceous greensand woodlands to the north, Jurassic clays to the west and areas of Tertiary clay beds to the south‐east (Barron, 1976). Within this area, all publicly accessible locations in the countryside were appraised as possible field locations for regular tick collections. A 5 km grid was placed over the study area with two accessible sites chosen in each complete grid square and one site in a partial grid square.

Thirty‐seven publicly accessible sites were surveyed during the study (Appendix S1). Each site was classified according to one of the following five habitats: (a) high biodiversity woodland (HB), (b) low biodiversity woodland (LB), (c) woodland edge (WE), (d) calcareous scrub grassland (CS) or (e) calcareous grassland (CG). Arable and pasture farmland were not surveyed due to low suitability for I. ricinus. Sampling in CS was targeted around the scrub. Sampling in CS and CG was completed away from the woodland edge habitat.

To categorize woodland as either high or low biodiversity, each woodland was appraised according to its woodland structure (canopy, understorey, shrub, herb and ground layers) in line with the Dansereau index (Dansereau, 1951). Those woodlands that were predominantly coniferous or a monoculture (e.g., Fagus sylvatica) were categorized as low biodiversity. Those with varied canopy and understorey, with areas of hazel coppice and ancient woodland flora, were categorized as high biodiversity. Calcareous grassland was categorized by the presence/absence of hawthorn (Crataegus monogyna) scrub.

The field sites included 10 high biodiversity woodlands, 7 low biodiversity woodlands, 8 woodland edge habitats, 8 calcareous scrub grasslands, and 4 non‐scrub calcareous grasslands. Each site was visited once per year, during spring (late April/early May) during 2013–2017. Nymph and adult ticks were collected by flagging a 1 m2 piece of cloth over the vegetation (Milne, 1943). At each survey site, approximately 30–70 5 m2 transects were sampled. For efficiency, fewer transects were carried out at sites with higher tick density (e.g., over 50 nymphs collected by transect 30). Surveys were carried out on dry days between 10 a.m. and 5 p.m. Ticks were collected and stored at −80°C until morphological identification to species level (Hillyard, 1996) and pathogen analysis could be carried out. Larvae were not collected.

Pathogen screening

Following identification, a subset of I. ricinus nymphs was sorted into individual tubes for processing and Borrelia screening. Where possible, 50 ticks were tested for Borrelia infection from each site for each time interval. Where less than 50 ticks were collected, all ticks were tested. DNA extracts were prepared by alkaline hydrolysis of ticks in NH4OH (Sigma Aldrich, St. Louis, MO, USA) following Hansford et al. (2015): 100 μl of NH4OH (1 M) were added to each tick and then tubes were placed on a heating block at 100°C for 25 min. Tubes were then briefly centrifuged and heated again at 100°C for 15 min with tube lids open to evaporate the ammonium until approximately 50 μl remained. Negative extraction controls (i.e., tubes containing no ticks) were included in the process. Lysates were then stored at 4°C overnight before qPCR.

DNA extracts were tested for the presence of Borrelia DNA using a modified pan‐Borrelia qPCR assay against the Borrelia 16S rRNA gene (Parola et al., 2011) using primers 5′‐AGCCTTTAAAGCTTCGCTTGTAG‐3′ (forward), 5′‐GCCTCCCGTAGGAGTCTGG‐3′ (reverse) and probe 5′‐CCGGCCTGAGAGGGTGAAC GG‐BHQ1‐3′. PCR plates were prepared on a QIAgility system (QIAgen) using 10 μl TaqManTM Fast Universal PCR Master Mix (Applied Biosystems, Waltham, MA, USA), 1 μl primer/probe mix (final concentration 900 nM each primer and 250 nM probe), 4 μl PCR grade H2O, and 5 μl sample per well. PCR grade H2O was used as a negative control and B. burgdorferi s.s. DNA is used as a positive control. The PCR was performed on a QuantStudio 7 Flex real‐time PCR system (Applied Biosystems), using a program consisting of an initial 20 s at 95°C followed by 40 cycles of 3 s at 95°C and 30 s at 60°C.

Positive samples were typed by sequencing of the 5S‐23S rRNA intergenic spacer (Alekseev et al., 2001) to determine B. burgdorferi genospecies (primers: forward 5′‐GAGTTCGCGGGAGAGTAGGTTATTGCC‐3′, reverse 5′‐TCAGGGTACTTAGATGGTTCACTTCC‐3′). Reactions were prepared in 50 μl volumes, containing 5 μl 10× PCR reaction buffer (‐MgCl2), 1 μl 10 mM dNTPs, 1.5 μl 50 mM MgCl2, 2 μl each primer (from 10 μM stock), 0.2 μl Platinum Taq DNA polymerase (Invitrogen), 33.3 μl PCR grade H2O, and 5 μl sample. Reactions were carried out in a Labtech G‐Storm Thermocycler with the following cycling conditions: 5 min at 94°C, followed by 10 cycles of 94°C for 10 s, 70°C for 30 s (lowering by 1°C each cycle) and 72°C for 30 s, then 40 cycles of 94°C for 20 s, 60°C for 30 s and 72°C for 30 s, with a final extension of 72°C for 7 min. PCR products were purified using Agencourt AMPure XP beads (Beckman Coulter, Brea, CA, USA) according to the manufacturer's protocol, and then sent to the PHE Genomic Services and Development Unit (Colindale, London) for Sanger sequencing. Borrelia burgdorferi genospecies was determined by inputting trimmed consensus sequences (using DNAStar Lasergene software, Madison, WI, USA) into the NCBI BLAST (Basic Local Alignment Search Tool) algorithm and genospecies were assigned if they had >97% identity. To confirm results, sequences were aligned using the MUSCLE algorithm and a phylogenetic tree was built using maximum likelihood (DNAStar Lasergene). Borrelia garinii and the closely related, newly reported Borrelia bavariensis were commonly reported matches for many samples, as they cannot be genetically discriminated via the 5S‐23S rRNA intergenic spacer region. Due to the dominance of B. garinii previously reported in England, they are geographically highly unlikely to be from B. bavariensis and therefore reads matching either B. garinii or B. bavariensis were assigned as B. garinii.

Statistical analyses

All data analysis was conducted in R version 4.0.2 (R Core Development Team, 2020). Two sites were surveyed under light rain in 2014 (Great Durnford and Lower Woodford) and were excluded from analyses focussing on DON (nymphs per 100 m2). The number of I. ricinus nymphs collected per 5 m transect was used as the response variable in a Generalized Linear Mixed Effect Model (GLMM), in the lme4 package (Bates et al., 2015), using a Poisson distribution. Habitat type (HB, LB, WE, CS, CG), year (2013–2017) and the interaction between habitat type and year were included as fixed effects. An observation level random effect, to account for overdispersion (Harrison, 2014) and site names were included as random terms. Explanatory variables were assessed based on the lowest AICc value (Brewer et al., 2016) using the dredge function from the MuMIn package (Barton, 2019).

NIP was modelled using a binomial GLMM with a logit link. The response variable was the proportion of nymphs infected and habitat type, year and the interaction term between habitat type and year were included as fixed terms. An observation level random effect, to account for overdispersion (Harrison, 2015) and site names were included as random terms. The dredge function was used for model selection and the model with the lowest AICc was kept. While the aim was to test 50 nymphs per site for Borrelia infection, nymph density was too low in some sites. Using Daniels (1999) to calculate the sample size needed to accurately estimate the prevalence and based on an average prevalence of Borrelia of 3.9%–8.5% (Rauter & Hartung, 2005; Strnad et al., 2017), sites with <35 nymphs tested were excluded from the analysis.

The number of infected nymphs (DIN, infected nymphs per 100 m2) for each site was used as the response variable in a GLMM using a Poisson distribution with an offset of distance surveyed on a log scale. Habitat type, year and the interaction between habitat type and year were included as fixed effects. An observation level random effect, to account for overdispersion (Harrison, 2014) and site names were included as random terms. Explanatory variables were assessed based on the lowest AICc value. When appropriate, Tukey HSD post‐hoc comparisons were conducted to assess pairwise differentiation among levels of categorical variables. For all analyses, significant p‐values appear in the figures and the effect size is represented by confidence intervals.

RESULTS

Over the five spring surveys 2013–2017, 14,727 I. ricinus ticks (845♂, 799♀, 13,083 nymphs) were collected, with approximately 3000–3400 nymph and adults collected each year, except in 2014, where only ~1750 were collected (Table 1). The proportion of nymphs to adult I. ricinus was ~12.5:1 in 2013 and 2015 and ~7.5:1 for the other 3 years (2014, 2016, 2017), although the adult density was twice as high in 2016 and 2017 compared to 2014. For each habitat, location and survey interval, DON (nymphs per 100 m2) and DIN (infected nymphs per 100 m2) were calculated (Table 1, Appendix S1).

TABLE 1

Average density of questing Ixodes ricinus nymphs (DON; [n/100m2]), I. ricinus nymphal infection prevalence (NIP) with Borrelia sensu lato and density of Borrelia‐infected I. ricinus nymphs (DIN; [n/100m2]) for each of the five habitats surveyed in South Wiltshire across 5 years (2013–2017)

Year	Habitat	DON (n/100 m2 ± SD)	NIP (%) [95% CI]	DIN (n/100 m2 ± SD)
2013	Calcareous grassland (CG)	0.1 ± 1.3	0	0
	Calcareous scrub grassland (CS)	18.7 ± 45.7	6.1 [2.7–9.4]	1.7 ± 3.0
	High biodiversity woodland (HB)	68.3 ± 119.2	5.8 [3.6–7.9]	4.4 ± 2.7
	Low biodiversity woodland (LB)	20.3 ± 40.4	4.2 [1.1–7.2]	0.8 ± 1.0
	Woodland edge (WE)	51.5 ± 86.9	6.9 [4.1–9.7]	3.5 ± 3.1
# ticks	3027 (125♂ 113♀ 2789 n a )	33.8 ± 77.9
NIP	67/1137; 5.9% [95% CI: 4.5–7.3]
2014	Calcareous grassland (CG)	0.1 ± 1.3	0	0
	Calcareous scrub grassland (CS)	8.6 ± 21.0	3.0 [0.4–5.6]	0.4 ± 0.5
	High biodiversity woodland (HB)	25.8 ± 45.4	4.7 [2.6–5.6]	1.8 ± 2.6
	Low biodiversity woodland (LB)	12.8 ± 23.3	1.0 [0–2.3]	0.2 ± 0.3
	Woodland edge (WE)	30.8 ± 52.2	4.5 [2.3–6.7]	1.1 ± 2.1
# ticks	1890 (137♂ 128♀ 1625n)	16.9 ± 36.9
NIP	40/1081; 3.7% [95% CI: 2.6–4.8]
2015	Calcareous grassland (CG)	0.2 ± 2.2	0	0
	Calcareous scrub grassland (CS)	15.5 ± 33.3	6.1 [3.0–9.2]	1.3 ± 2.0
	High biodiversity woodland (HB)	68.6 ± 88.5	5.9 [3.7–8.0]	3.9 ± 2.8
	Low biodiversity woodland (LB)	24.5 ± 36.7	8.0 [4.9–11.1]	1.3 ± 1.5
	Woodland edge (WE)	74.7 ± 99.4	11.1 [8.0–14.2]	7.9 ± 7.0
# ticks	3342 (151♂ 127♀ 3064n)	36.7 ± 69.0
NIP	108/1381; 7.8% [95% CI: 6.4–9.2]
2016	Calcareous grassland (CG)	0.2 ± 1.8	0	0
	Calcareous scrub grassland (CS)	18.1 ± 47.6	5.0 [2.1–7.9]	1.4 ± 2.2
	High biodiversity woodland (HB)	62.6 ± 94.4	11.2 [8.4–14.0]	6.7 ± 4.9
	Low biodiversity woodland (LB)	25.1 ± 38.2	8.9 [5.8–12.0]	2.4 ± 2.5
	Woodland edge (WE)	60.6 ± 87.2	12.8 [9.6–16.1]	8.3 ± 7.8
# ticks	3359 (220♂, 231♀ 2908n)	34.5 ± 69.2
NIP	144/1416; 10.2% [95% CI: 8.6–11.7]
2017	Calcareous grassland (CG)	2.0 ± 8.1	0	0
	Calcareous scrub grassland (CS)	18.4 ± 37.6	3.1 [1.1–5.2]	0.5 ± 0.6
	High biodiversity woodland (HB)	49.7 ± 63.6	8.9 [6.3–11.4]	3.9 ± 2.8
	Low biodiversity woodland (LB)	29.2 ± 38.9	12.9 [9.4–16.4]	3.6 ± 2.8
	Woodland edge (WE)	46.9 ± 83.2	11.1 [7.8–14.4]	5.5 ± 5.3
# ticks	3109 (212♂ 200♀ 2697n)	30.5 ± 55.6
NIP	134/1487; 9.0% [95% CI: 7.6–10.5]
All years	Calcareous grassland (CG)	0.5 ± 4.0	0	0
	Calcareous scrub grassland (CS)	15.7 ± 38.0	4.6 [3.4–5.9]	1.1 ± 1.9
	High biodiversity woodland (HB)	53.5 ± 85.9	7.4 [6.3–8.5]	4.1 ± 3.5
	Low biodiversity woodland (LB)	22.3 ± 36.2	7.9 [6.5–9.4]	1.7 ± 2.2
	Woodland edge (WE)	51.5 ± 83.2	9.5 [8.2–10.9]	5.5 ± 6.0
# ticks NIP	14,727 (845♂ 799♀ 13,083n 493/6502; 7.6% [95% CI: 6.9–8.2]

Abbreviations: CG, calcareous grassland; CI, confidence interval; CS, calcareous scrub grassland; HB, high biodiversity woodland; LB, low biodiversity woodland; SD, standard deviation; WE, woodland edge.

n: nymphs.

Density of questing nymphs (DON)

The density of nymphs (DON) varied between habitats both within a year and between years (Figure 1, Table 1), and the selected model included the interaction term between year and habitat. DON was consistently very low (0.5) in CG in all years while the two habitats with consistently the highest DON were HB woodland (53.5) and WE (51.5), with lower DON reported in LB (22.3) and CS (15.7) (Table 1 includes standard deviation). Interestingly DON seemed to be lower in 2014 across all locations (overall DON in 2014 was 16.9, compared to 30.5–36.7 in all other years).

FIGURE 1

Density of questing Ixodes ricinus nymphs (DON) by site, habitat and year in South Wiltshire. Horizontal black lines represent the medians and red lines the means. Sites are as follows; 1: Bentley, 2: Farley, 3: Grovely east, 4: Grovely west, 5: Gutch Common, 6: Pitton, 7: Stockton, 8: Sutton Mandeville, 9: West Knoyle, 10: Whiteparish, 11: Charlton AS, 12: Chilmark, 13: Cholderton, 14: Great Durnford, 15: Hindon, 16: Newtown, 17: Redlynch, 18: Woodford, 19: Deanland, 20: Donheads, 21: East Knoyle, 22: Hamptworth, 23: Langley, 24: Wilton, 25: Winterslow, 26: Bowerchalke, 27: Swallowcliffe, 28: Chicklade, 29: Coombe Bissett, 30: Fovant, 31: Pepperbox, 32: Win green, 33: Wylye, 34: Alvediston, 35: Burcombe, 36: Figsbury, 37: Yarnbury. HB, high biodiversity; LB, low biodiversity

By considering a mean DON of 80 as a ‘high DON threshold’, between six and nine of the 37 sites exceeded this threshold each year, except in 2014 when only two sites reported mean DON >80. Most of the field sites reporting high DON were in HB (6/10 reported a mean DON >80 in at least one year) and WE (5/8), compared to LB (1/7), CS (1/8) or CG (0/4). The highest reported DON in any location was 200 (SD: 179.8) in an HB woodland in 2016, where DON exceeded the threshold every year (Appendix S1). For the 10 HB locations, DON appeared to decline from 2013 to 2014, then increase again for 2015 before declining over the next 2 years (Figure 2). The only exceptions are for two sites where the highest mean DON was observed in 2016. For all WE sites, a similar trend is observed over the 5 years. There is one exception, which reported 2015 as the lowest year (Figure 2).

FIGURE 2

Inter‐annual variation in density of Ixodes ricinus nymphs (DON) within and between sites in South Wiltshire in (a) high biodiversity woodland (HB) and (b) woodland edge (WE) habitats. Raw data is presented here

Predicted DON is shown in Figure 3 and Table 2 and p‐values for each pairwise comparison between habitat type and years are available in Appendix S2. In all 5 years, DON was significantly lower in CG compared to all other habitats (p < 0.05 for each pairwise comparison). In four of the years (excluding 2013) DON was significantly lower in CS compared to WE, significantly lower in CS compared to HB in three of the years (2015–2017) and significantly lower in CS compared to LB in 2017 (p < 0.05 for each pairwise comparison). In terms of temporal variation within each habitat, DON in CG remained constantly low throughout the study period (p > 0.05). Regarding all other habitats, DON was lower in 2014 compared to other years (p < 0.05 for each pairwise comparison). The highest DONs were recorded in WE and HB during 2015 and 2016 (Figure 3).

FIGURE 3

Predicted density of Ixodes ricinus nymph (DON) in South Wiltshire habitats in 2013–2017: Calcareous grassland (green), calcareous scrub grassland (grey), woodland edge (blue), low biodiversity woodland (red) and high biodiversity woodland (orange). Error bars represent 95% confidence interval. Different numbers (1–3) show significant difference among habitat type within a same year and different letters (a–c) show significant difference among year for each habitat type based on Tukey HSD tests (p < 0.05)

TABLE 2

Predicted density of Ixodes ricinus nymphs (DON) based on GLMM output and 95% confidence intervals by habitat and year

	HB DON [95% CI]	LB DON [95% CI]	WE DON [95% CI]	CS DON [95% CI]	CG DON [95% CI]
2013	37.0 [20.4–67.0]	10.5 [5.1–21.4]	25.5 [13.1–49.5]	6.8 [3.4–13.3]	0 [0–0.3]
2014	13.3 [7.4–24.2]	7.1 [3.5–14.6]	19.1 [9.7–37.3]	3.2 [1.6–6.2]	0 [0–0.32]
2015	37.9 [20.9–68.5]	15.5 [7.6–31.6]	45.2 [23.3–87.6]	6.2 [3.2–12.2]	0.1 [0–0.4]
2016	37.9 [20.9–68.5]	14.7 [7.2–29.9]	38.6 [19.9–74.9]	7.1 [3.6–13.8]	0.1 [0–0.4]
2017	30.4 [16.8–55.0]	17.3 [8.5–35.1]	23.4 [12.0–45.4]	7.1 [3.6–13.8]	0.7 [0.2–2.1]

Abbreviations: CG, calcareous grassland; CI, confidence interval; CS, calcareous scrub grassland; HB, high biodiversity woodland; LB, low biodiversity woodland; WE, woodland edge.

Borrelia infection prevalence in questing nymphs (

Out of 6502 nymphs tested, 493 were positive, giving a global NIP of 7.6%. Borrelia NIP for each location is reported in Table 3. Twenty of the 37 sites reported NIP >10% in at least 1 year where the sample size exceeded 30 nymphs (n ≥ 30), with one WE site exceeding 10% every year (12%–18%). NIP >10% (n ≥ 30) was found at 4 sites in 2013, 3 in 2014, 9 in 2015, 10 in 2016 and 12 in 2017. NIP >10% (n ≥ 30) occurred most in WE habitats (15 times), compared with HB (12), LB (8), CS (2) and CG (none). NIP ≥20% was reported on eight occasions (4 occasions in both WE and LB sites) (Table 3).

TABLE 3

Nymphal Borrelia‐infection prevalence (NIP) in Ixodes ricinus and sample sizes for all locations in South Wiltshire, presented for all 5 years (2013–2017)

Site	2013	2014	2015	2016	2017	Total
HB 1	6 (3/50)	0 (0/50)	2 (1/50)	0 (0/50)	2 (1/42)	2 (5/242)
HB 2	14 (7/50)	4 (2/51)	12 (6/50)	16 (8/50)	6 (3/50)	10 (26/251)
HB 3	4 (2/50)	4 (2/50)	6 (3/50)	14 (7/50)	8 (4/50)	7 (18/250)
HB 4	6 (3/50)	2 (1/50)	4 (2/50)	6 (3/50)	2 (1/50)	4 (10/250)
HB 5	0 (0/40)	6 (2/32)	2 (1/50)	30 (15/50)	18 (9/50)	11 (27/250)
HB 6	2 (1/50)	12(6/50)	8 (4/50)	8 (4/50)	10 (5/50)	8 (20/250)
HB 7	6 (3/50)	0 (0/12)	8 (4/50)	2 (1/50)	12 (6/50)	6 (14/250)
HB 8	4 (2/50)	6 (2/32)	6 (3/50)	14 (7/50)	10 (5/50)	8 (19/232)
HB 9	18 (2/11)	0 (0/2)	18 (2/11)	22 (5/23)	19 (6/32)	19 (15/79)
HB 10	6 (3/50)	6 (3/50)	2 (1/50)	6 (3/50)	2 (1/50)	4 (11/250)
WE 11	6 (3/50)	10 (5/50)	10 (5/50)	8 (4/50)	8 (4/50)	8 (21/250)
WE 12	4 (2/50)	0 (0/50)	2 (1/49)	12 (6/50)	8 (4/50)	5 (13/250)
WE 13	2 (1/50)	5 (2/38)	14 (7/50)	6 (3/47)	3 (1/40)	6 (14/225)
WE 14	14 (7/50)	12 (6/50)	18 (9/50)	18 (9/50)	15 (7/47)	15 (38/247)
WE 15	5 (2/39)	6 (2/36)	30 (15/50)	24 (12/50)	29 (14/48)	20 (45/223)
WE 16	20 (1/5)	0 (0/9)	5 (2/38)	28 (14/50)	0 (0/16)	14 (17/118)
WE 17	4 (1/25)	0 (0/50)	6 (3/50)	4 (2/50)	4 (2/50)	4 (8/225)
WE 18	10 (5/50)	0 (0/50)	2 (1/50)	2 (1/50)	14 (7/50)	6 (14/250)
LB 19	2(1/50)	2 (1/50)	0 (0/50)	0 (0/50)	4 (2/50)	2 (4/250)
LB 20	0 (0/5)	0 (0/6)	10 (5/50)	9 (4/46)	12 (6/50)	10 (15/157)
LB 21	33 (1/3)	0 (0/11)	22 (11/50)	22 (11/50)	20 (10/50)	13 (33/250)
LB 22	0 (0/50)	0 (0/50)	0 (0/50)	6 (3/50)	6 (3/50)	2 (6/250)
LB 23	NA	0 (0/17)	4 (1/27)	3 (1/29)	8 (4/50)	5 (6/123)
LB 24	8 (4/50)	2 (1/50)	6 (3/50)	4 (2/50)	12 (6/50)	6 (16/250)
LB 25	10 (1/10)	0 (0/19)	17 (4/23)	16 (8/50)	28 (14/50)	18 (27/152)
CS 26	0 (0/3)	NA	0 (0/2)	0(0/1)	0 (0/3)	0 (0/7)
CS 27	6 (3/50)	4 (2/50)	8 (4/49)	4 (2/50)	8 (4/50)	6 (15/249)
CS 28	8 (1/12)	0 (0/1)	12 (2/17)	2 (1/45)	0 (0/23)	4 (4/98)
CS 29	6 (1/16)	0 (0/2)	0 (0/13)	0 (0/2)	0 (0/11)	2 (1/44)
CS 30	0 (0/2)	6 (1/17)	14 (7/50)	0 (0/8)	0 (0/50)	12 (15/127)
CS 31	0 (0/50)	4 (2/50)	2 (1/49)	8 (4/50)	4 (2/50)	4 (9/249)
CS 32	7 (1/15)	0 (0/19)	0 (0/50)	4 (2/50)	0 (0/50)	2 (3/184)
CS 33	12 (6/50)	0 (0/26)	NA	15 (2/13)	6 (3/50)	8 (11/139)
CG 34	0 (0/1)	NA	NA	0 (0/1)	0 (0/3)	0 (0/5)
CG 35	NA	NA	0 (0/1)	NA	0 (0/22)	0 (0/23)
CG 36	NA	0 (0/1)	0 (0/2)	0 (0/1)	NA	0 (0/4)
CG 37	NA	NA	NA	NA	NA

Note: All locations with nymph infection prevalence of 10% or above and where >30 were tested are in bold.

Abbreviations: CG, calcareous grassland; CS, calcareous scrub grassland; HB, high biodiversity woodland; LB, low biodiversity woodland; NA, no data available as no ticks were collected; WE, woodland edge.

If we consider variation in NIP over the 5 years (Figure 4), there appears to be a similar trend to that reported earlier for DON, with the lowest infection prevalence in 2014 (3.7%) and the highest in 2015–2017 (7.8%–10.2%) (see Table 1 including confidence intervals [CI]). NIP was significantly lower in 2014 compared to 2016 (p = 0.003) and 2017 (p = 0.009) (Figure 5).

FIGURE 4

Variation in the prevalence of Borrelia in questing Ixodes ricinus (NIP) during 2013–2017: In calcareous grassland (green), calcareous scrub grassland (grey), woodland edge (blue), low biodiversity woodland (red) and high biodiversity woodland (orange) habitats. Error bars represent 95% confidence interval

FIGURE 5

Variation in the predicted prevalence of Borrelia in questing Ixodes ricinus nymphs (NIP) in South Wiltshire across years. Error bars represent 95% confidence interval. Statistical significance from Tukey tests: **p < 0.01

For HB, 7.4% of nymphs were infected with Borrelia. This varied from 4.7% in 2014 to 11.2% in 2016. For WE, 9.5% of nymphs were infected, varying from 4.5% in 2014 to 3 years of 11.1%–12.8% in 2015–2017. For LB, 7.9% of nymphs were infected, with infection prevalence as low as 1.0% in 2014, increasing to 8.0%–12.9% in 2015–2017. For CS, 4.6% of nymphs were infected, varying between 5.0%–6.1% in 2013, 2015 and 2016, with the lowest infection prevalence in 2014 (3.0%). No nymphs were infected with Borrelia in CG habitat (Table 1 includes CI).

Of the 493 positive nymphs, 33 were infected with B. afzelii (6.7%), 132 with B. garinii (26.8%), 64 with B. valaisiana (13.0%), 58 with B. burgdorferi s.s. (11.8%) and 206 could not be identified to species level (41.8%). Borrelia garinii is most common with the highest NIP (2.0% [95% CI: 1.7–2.4]), followed by B. valaisiana (1.0% [0.7–1.2]), B. burgdorferi s.s. (0.9% [0.7–1.1]) and B. afzelii (0.5% [0.3–0.7]). For HB, WE and LB the proportions of genospecies varies, with B. garinii accounting for 22.4% (including untyped sequences) in HB, 23.5% in WE, and 41.1% in LB (Table 4, Figure 6a). However, for B. valaisiana the association was partly inverse with 15.8% of infected ticks in HB, 11.2% in WE and 13.1% in LB (Table 4, Figure 6a). For the remaining typed sequences, 4.7%–7.9% of sequences from infected nymphs in the three woodland habitats were infected with B. afzelii (Table 4, Figure 6a), but NIP was higher in CS. Here, the proportion of B. afzelii infected ticks accounted for 13.7% with a lower proportion than B. garinii (21.6%, but a comparable proportion of B. valaisiana [9.8%]). Borrelia burgdorferi s.s. was detected in four of the habitats (CS: 13.7%, HB: 12.1%, LB: 8.4%, WE: 12.9%) (Table 4, Figure 6a). All sequences were uploaded to GENBANK (Accession numbers OL848157‐OL848443).

TABLE 4

Borrelia‐infection prevalence in nymphal Ixodes ricinus and 95% confidence intervals by habitat in South Wiltshire for all typed and untyped genospecies of Borrelia

	B. afzelii (%) [95% CI]	B. garinii (%) [95% CI]	B. valaisiana (%) [95% CI]	B. burgdorferi. sensu stricto (%) [95% CI]	Untyped Borrelia (%) [95% CI]
CG	0	0	0	0	0
CS	0.6 [0.2–1.1]	1.0 [0.4–1.6]	0.5 [0.1–0.9]	0.6 [0.2–1.1]	1.9 [1.1–2.7]
HB	0.6 [0.3–0.9]	1.7 [1.1–2.2]	1.2 [0.7–1.6]	0.9 [0.5–1.3]	3.1 [2.4–3.8]
LB	0.4 [0–7]	3.3 [2.3–4.2]	1.0 [0.5–1.6]	0.7 [0.2–1.1]	2.6 [1.8–3.5]
WE	0.4 [0.1–0.8]	2.2 [1.6–2.9]	1.1 [0.6–1.5]	1.2 [0.7–1.7]	4.5 [3.6–5.5]

Note: Identities of >97% were considered for matches. Clustering on the phylogenetic tree was used to finalize sequences.

Abbreviations: CG, calcareous grassland; CI, confidence interval; CS, calcareous scrub grassland; HB, high biodiversity woodland; LB, low biodiversity woodland; WE, woodland edge.

FIGURE 6

Proportion of questing Ixodes ricinus nymphs (%) in South Wiltshire infected with Borrelia afzelii, B. garinii, B. valaisiana, B. burgdorferi sensu stricto and untyped Borrelia depending on habitat type (a) and year (b). HB, high biodiversity woodland; LB, low biodiversity woodland. Identities of >97% were considered for matches. Clustering on the phylogenetic tree was used to finalize sequences

If we compare relative proportions of genospecies over the years (Figure 6b), B. garinii was most common in 4 years (24%–31.5% vs. 10% in 2014) and B. valaisiana most common in 2013–2014 (20%–20.9% vs. 9%–13% in 2015–2017). For the years 2015–2016, B. burgdorferi s.s. accounted for 18.8%–25.9% of positive ticks, but this was much lower in the other years (0%–5%). Finally, B. afzelii accounted for a small proportion of infected ticks in 3 years (4.6%–6.9%) except in 2013 (14.9%), with no evidence of B. afzelii in 2014.

Density of infected nymphs (DIN)

Calculating DIN as a mean of each site by habitat, the highest mean DIN was 5.5 in WE, followed by 4.1 in HB, 1.7 in LB and 1.1 in CS (Table 1, includes standard deviations). DIN in all CG sites was zero. DIN ≥10 was reported on 11 occasions at eight sites during 2015–2017 (Appendix S1, Figure 7) throughout the study period, including at three HB sites and five WE sites. The highest DIN measurements were reported in the same WE site as 22.7 and 19.4 (Appendix S1). DIN ≥5 was reported at nine sites in 2013 (5 HB, 3 WE, 1 CS), three in 2014 (2 HB, 1 WE), eight in 2015 (3 HB, 4 WE, 1 CS), 12 in 2016 (6 HB, 4 WE, 2 LB) and seven in 2017 (2 HB, 3 WE, 2 LB).

FIGURE 7

Density of Borrelia‐infected Ixodes ricinus nymphs (DIN; size of pie chart) and genospecies dominance at field sites in South Wiltshire: (CG [calcareous grassland] habitats in light green, CS [calcareous scrub grassland] in grey, HB [high biodiversity woodland] in orange, LB [low biodiversity woodland] in red and WE [woodland edge] in blue). Empty circle represents sites where DIN is zero. Site numbers are: 1: Alvediston, 2: Bentley, 3: Bowerchalke, 4: Burcombe, 5: Charlton all saints, 6: Chicklade, 7: Chilmark, 8: Cholderton, 9: Coombe Bissett, 10: Deanland, 11: Donheads, 12: East Knoyle, 13: Farley, 14: Figsbury, 15: Fovant, 16: Great Dunford, 17: Grovely east, 18: Grovely west, 19: Gutch Common, 20: Hamptworth, 21: Hindon, 22: Langley, 23: Newtown, 24: Pepperbox, 25: Pitton, 26: Redlynch, 27: Stockton, 28: Sutton Mandeville, 29: Swallowcliffe, 30: West Knoyle, 31: Whiteparish, 32: Wilton, 33: Win green, 34: Winterslow, 35: Woodford, 36: Wylye, 37: Yarnbury

The selected model to explain variations in DIN included year and habitat as predictors. Predicted DIN was significantly lower in 2014 (0.68, 95% CI: 0.4–1.2) compared to 2015 (1.5, 95% CI: 0.9–2.3, p = 0.04), 2016 (2.0, 95% CI: 1.3–3.1, p < 0.001) and 2017 (1.7, 95% CI: 1.1–2.7, p = 0.005) but similar to 2013 (1.2, 95% CI: 0.7–2.0, p > 0.05), while there was no difference between the other years (Figure 8a). Regarding the effect of habitat type on DIN, modelled DIN was significantly lower in CS (0.7, 95% CI: 0.4–1.3) compared to WE (2.0, 95% CI: 1.3–3.2, p = 0.006) but similar across to the other habitats (HB: 1.6, 95% CI: 1.0–2.5; LB: 1.3, 95% CI: 0.8–2.1; p > 0.05) (Figure 8b).

FIGURE 8

Predicted variation in the density of Ixodes ricinus nymphs infected (DIN) with Borrelia in South Wiltshire across years (a) and habitat types (b). Error bars represent 95% confidence interval. Statistical significance from Tukey tests: *p < 0.05, **p < 0.01, ***p < 0.001. DIN, density of infected nymphs; HB, high biodiversity; LB, low biodiversity

DISCUSSION

In the study area of south Wiltshire, a typical, south‐west England endemic focus for both high densities of I. ricinus and human cases of Lyme disease, DON, NIP and DIN varied in spring between habitats and between years (Figure 7).

From the five typical lowland England habitats studied, calcareous grassland, in the absence of scrub, was found to be largely unsuitable for questing I. ricinus (typically DON <2), even in the presence of livestock. However, as calcareous grassland becomes invaded by scrub (usually by hawthorn), to create scrub grassland, the density of questing nymphs increases (DON = 8–18). This has important habitat management implications for ticks. Scrub management on grazed grassland can be an important management strategy for minimizing the survival of I. ricinus and thus reducing tick‐borne disease in humans and other animals. The small number of I. ricinus collected from scrub‐free grassland and tested for Borrelia sp. were all negative, compared to 4.6% of questing nymphs in calcareous scrub that were Borrelia positive. Scrub encroachment can be an issue for the habitat management of calcareous grassland. A grassland without scrub encroachment is usually more desirable for rare calcareous flora, however for some grassland invertebrates, the scrub can be favourable (Kirby, 2013). Scrub provides an increased surface area for ecotonal habitat that supports the microclimate for off‐host tick survival, as well as shelter belts for a range of wildlife and livestock (Medlock et al., 2008). However, nymphal densities in calcareous scrub grassland were rarely high, with mean DON significantly lower than in high biodiversity woodland and woodland edge habitats in three of the years, with very high DON (>80) only reported on one occasion in calcareous scrub grassland habitat over the 5 years.

The three woodland habitats showed variation in DON between years (significantly lower in 2014), with overall annual mean densities in HB and WE habitat higher compared to LB habitat, with the presence of high nymphal densities recorded in >50% of HB and WE sites, but only one out of seven LB sites, with the highest DON in WE. This high I. ricinus nymphal density in ecotonal/woodland edge habitat contrasts with the mosaic fragment ecotones found in scrub grassland but supports previous findings highlighting the importance of the ecotone along woodland paths/rides (i.e., linear grassy trackway) (Medlock et al., 2012) and agri‐environment field margins (Medlock et al., 2020). Woodland edge acts as an interface habitat for animals and provides movement corridors, where nearby shelter is provided for animals browsing in open countryside (Medlock et al., 2020). Guidance to the public on tick exposure highlights the importance of woodland, perhaps it should also highlight woodland edges (e.g., where footpaths are located) as having equal risk. The importance of woodland edge as a tick habitat may also have implications for livestock health, for example in the management of bovine babesiosis, particularly in areas where cattle are grazed adjacent to woodland. Although exposure to I. ricinus is largely seasonal, careful management and avoidance (e.g., use of temporary fences) of these edge habitats could be advisable to minimize risk.

The presence and abundance of questing ticks are important indicators of risk but alone do not directly relate to disease hazard, but it is its combination with NIP that defines the hazard. A recent snap‐shot Borrelia study in Wales showed that NIP may be extremely low (i.e., 0%–1%; Cull et al., 2021), highlighting that habitat and host diversity drive infection prevalence. There are at least four notable genospecies of B. burgdorferi s.l. currently circulating in British I. ricinus, each with its own transmission cycle, generally cycling between mammals and ticks (B. afzelii), birds and ticks (B. garinii, B. valaisiana), or possibly birds/mammals and ticks (B. burgdorferi s.s.). The relative abundance of fauna, the densities of possible dilution hosts, like deer and livestock, and the density of nymphs will invariably impact nymph infection prevalence both spatially within and between habitats and between years.

In this study, global NIP across all locations and years was 7.6%, which is considerably higher than in other UK studies 3.3%–3.9% in Northern England and Wales (Bettridge et al., 2013; Hall, 2018; Hansford et al., 2015) and 1.7%–5.6% in Scotland (Gandy, 2020; James et al., 2013; Millins et al., 2016); but consistent with that reported for some of the study sites in southern England (Cull et al., 2021). This 5‐year study also showed considerable variation between years, with NIP sometimes 2–3 times lower than a ‘normal’ year (3.7% in 2014 cf. 7.8%–10.2%), with low NIP mirroring the significantly low DON.

Infection prevalence also varied between habitats, with no evidence of infected ticks in calcareous grassland, highlighting a possible protective effect of livestock grazed scrub‐free grassland. The importance of woodland edge as a Borrelia risk habitat is evidenced by the fact that infection prevalence above 10% was common in woodland/edge habitats, with one particular WE site recording NIP >10% every year (sometimes >30%).

Previous studies on Borrelia infected ticks have shown a variation in the dominance of Borrelia genospecies across the UK, with B. afzelii and B. burgdorferi s.s. more common in Scotland and Northern England (Bettridge et al., 2013; James et al., 2013; Millins et al., 2016), and with B. garinii dominating in southern England (Cull et al., 2021; Layzell et al., 2018). Prior to this study, and one in Richmond Park (Hansford et al., 2021), there had been very few reports of B. burgdorferi s.s. in England (Hoodless et al., 1998). Finding sequences commonly attributable to B. burgdorferi s.s. may have important implications for the clinical manifestations of Lyme disease in England, as in North America, where it is the dominant genospecies, it is the cause of Lyme disease‐associated arthritis.

The dominance of B. garinii/valaisiana in this study further supports recent work (Cull et al., 2021). These genospecies are associated with birds, both woodland and game birds (Hoodless et al., 1998), with some evidence of a possible role for grey squirrels (Millins et al., 2015). If we consider these two genospecies alone, B. garinii and B. valaisiana accounted for 3.4% prevalence compared with B. burgdorferi s.s. and B. afzelii accounting for 1.4%.

The study area is a mixture of arable and pasture (mainly sheep/cattle) farmland, interconnected with woodland and hedgerows. All study locations provide a range of habitats for a number of bird species, including woodland birds that have been implicated elsewhere in circulating bird‐associated Borrelia and feeding juvenile life stages of I. ricinus (Hubalek et al., 1996; Humair et al., 1993). Perhaps notably though, in contrast to other parts of the European range of I. ricinus, large parts of southern England, and particularly the study area, have extensive rearing and release of non‐native pheasant (Tapper, 1999). Previous research highlighted the role of released gamebirds in feeding (and also mopping up from the environment) a large number of larval/nymphal I. ricinus, as well as amplifying existing cycles of bird‐associated Borrelia (Hoodless et al., 1998). Pheasants have been reported with median infestations of 23.5 larvae per bird in June and 26–59 nymphs per bird from April–August (Hoodless et al., 1998). In areas where NIP has been previously reported as 0%, infection prevalence in ticks collected from pheasants has been as high as 22%, suggesting a role for pheasant in infecting ticks (Craine et al., 1997). This role as a competent amplifying host for B. garinii/valaisiana was further supported (Kurtenbach et al., 1998) with evidence that B. afzelii does not survive the pheasant blood meal, which may explain the high infection prevalence and dominance for B. garinii/valaisiana over B. afzelii. Further work to investigate the role of pheasant in Borrelia transmission is now needed, particularly given that 47 million pheasants are released each year in Britain (Aebischer, 2019) and that they contribute significantly to bird biomass (Blackburn & Gaston, 2018). Comparing infection prevalence of bird‐associated Borrelia between locations with high, low and no pheasant release densities would help to further investigate the potential impact of pheasant release on Borrelia transmission cycles.

Interpreting Borrelia infection prevalence with small sample sizes can be misleading unless it is combined with the DON to calculate DIN, as high infection prevalence in an area with low tick density is not always indicative of high risk. Here we only tested a maximum of 50 ticks at each sampling interval, although this is consistent with a recent study (Cull et al., 2021) using the same methodology, and was also devised based upon the logistics of testing ticks from 37 sites over 5 years. Across four of the years (2013, 2015, 2016, 2017) DIN was not significantly different, although notably higher in 2016, with DIN significantly higher in 3 years (2015–2017) compared to 2014. DIN in woodland edge was significantly higher than in calcareous scrub grassland while it was similar across the other habitats, with a density of B. garinii/valaisiana infected nymphs higher in the three woodland habitats compared to calcareous scrub grassland. Increased biodiversity might be expected to provide greater equilibrium to parasite densities, but that does not appear to be the case here. Our measure of ‘biodiversity’ was generated as a measure of woodland structure using the Dansereau index, and not a measure of faunal diversity. Eight sites recorded DIN >10; five were woodland edge sites and three high biodiversity woodland, with an average of one infected nymph, found every 19.1 (WE) and 24.7 m (HB). This contrasts with every 62.5 (LB) and 96 (CS) m for the other two habitats. This confirms that Borrelia infected ticks in this South Wiltshire landscape are extremely common, and although there was no measure of human exposure, it may explain why Wiltshire has one of the highest incidence rates for Lyme disease, particularly when you consider some locations had a DIN >22, equating to one infected nymph every 4.4 m.

One outstanding question is the unusual year of 2014 where DON, NIP and DIN were significantly lower than in other years. We can only speculate that low nymph densities in 2014 may be due to: (a) fewer larval ticks having been fed the previous year due to low tick host densities, (b) fewer ticks overwintering, or (c) just fewer ticks active in 2014 possibly due to weather conditions. It may be that poor winter survival of nymphs led to a much reduced nymphal activity in the spring of 2014. Indeed, during January/February 2014, southern England recorded >200% of the normal monthly rainfall (Met Office, 2021). Perhaps the heavy rains affected nymph survival, although the mean daily temperatures in spring were normal for activity. This impact on DON is also replicated in NIP, suggesting that there must have been some effect from the previous year(s) to lead to lower NIP. Any reduced larval feeding in 2013 would have led to lower tick densities in 2014, but also an effect on their blood meal as infection prevalence was also significantly lower. The sequencing data reveals that infection prevalence with B. afzelii was absent in 2014 compared to other years, and so perhaps there was something that also impacted small mammals in 2013, a reported mast year (i.e., a year when trees produce a bumper crop of fruits and nuts) (Bogdziewicz et al., 2020; Hanley et al., 2019). The previous year, 2012, was one of the poorest years for tree fruiting for 10 years. Perhaps the low mast year of 2012 led to a low density of small mammals in 2013, which accounted for fewer larval ticks being fed in 2013 by small mammals and lower densities of nymphs, lower densities of infected nymphs and lower proportions of B. afzelii infected nymphs in 2014.

It does appear though that the tick‐Borrelia transmission cycle is susceptible to collapse in the face of a combination of ecological and climatic events, whether that be an impact of weather on masting (and subsequent poor fruiting events), a crash in tick‐host populations, or extreme weather events that affect tick survival and activity. Whatever the reason, this apparent low NIP and DIN may have affected clinical cases of Lyme disease. If we consider human cases from 2005 to 2017, the only year where Lyme disease cases declined from the previous year was 2014. If we compare incidence rates (per 100,000 population) over these 5 years for England and Wales, a similar trend to the DIN trend is observable: 1.64 in 2013; 1.49 in 2014; 1.83 in 2015; 1.94 in 2016; 2.70 in 2017 (Tulloch et al., 2019). Aside from ecological factors, the apparent reduction in Lyme disease incidence could also have been influenced by human behaviour, or other factors not considered in this study.

This study provides evidence that nymph density, Borrelia prevalence and the density of Borrelia‐infected nymphs vary spatially and temporally within and between habitats, thus presenting variation in the risk to public health. The findings are specific to this particular area of southern England and may not be reflective of other parts of the country, where habitats and tick hosts differ in their distribution and abundance. However, such evidence on the complexities of the spatial and temporal variation of Borrelia infected ticks can be used to inform risk awareness and risk management. Further studies that consider these ecological datasets alongside human tick bite exposure and subsequent Lyme disease transmission are needed to better translate these metrics into measurable public health risk factors. If it is accepted that higher DIN results in increased Lyme disease transmission, routine monitoring of DIN and targeting high DIN areas for intervention (e.g., habitat modification, host control or public health awareness) could have a significant impact on public health.

CONFLICT OF INTEREST

The authors declare no potential conflict of interest.

AUTHOR CONTRIBUTION

JM and AV conceived the project. JM, AV, KH, LM, BC, EG conducted the field work, MC, KH, BC, LM conducted the pathogen testing, SG, EG performed statistical analysis, SP supported sequencing analysis. JM led the writing of the manuscript, all authors contributed to the manuscript. SG prepared all the figures.

Appendix S1 The density of Ixodes ricinus nymphs (DON) (Table 1) and density of Borrelia‐infected Ixodes ricinus nymphs (DIN) (Table 2) for each site in South Wiltshire over the 5 years (2013–2017). DON>80 and DIN >5 are in bold. HB = high biodiversity woodland, LB = low biodiversity woodland, WE = woodland edge, CS = calcareous scrub grassland, CG = calcareous grassland

Click here for additional data file.

Appendix S2 P‐values based on post hoc Tukey HSD tests for each pairwise comparison of the density of questing nymphs (DON) between habitats for each year (Table 1) and between years for each habitat (Table 2) for the. HB: High biodiversity woodland and LB: Low biodiversity woodland

Click here for additional data file.