Sand lizards ( Lacerta agilis ) decrease nymphal infection prevalence for tick-borne pathogens Borrelia burgdorferi sensu lato and Anaplasma phagocytophilum in a coastal dune ecosystem

1. Understanding which factors determine tick-borne disease hazard can contribute to effective disease control. In Europe, the hazard of the pathogens Borrelia burgdorferi s.l. and Anaplasma phagocytophilum is determined by local tick densities (mainly Ixodes ricinus ) and the reservoir competence of the host species community. Sand lizards ( Lacerta agilis ) are common hosts for larvae and nymphs


| INTRODUC TI ON
A key element in preventing tick-borne diseases is understanding which factors determine hazard to humans (Swei et al., 2020).In Europe, Ixodes ricinus ticks are important vectors of various pathogens of veterinary and public health concern (Heyman et al., 2010).
Two pathogens that are frequently found in I. ricinus are bacteria of the Borrelia burgdorferi species complex and Anaplasma phagocytophilum.An infection with B. burgdorferi sensu lato causes Lyme borreliosis in humans, the most frequently transmitted tick-borne disease in both Europe and North America.Infection with A. phagocytophilum may cause human granulocytic anaplasmosis but is rarely diagnosed or reported for various reasons, among others the presence of non-characteristic "flu-like" symptoms and a lack of effective laboratory diagnostic tools (Matei et al., 2019).A. phagocytophilum is also of veterinary importance, causing tick-borne fever in both wild and domestic ruminants and in dogs.In Europe, the transmission dynamics of both pathogens relies on horizontal transmission between I. ricinus ticks and wildlife reservoir species (Hauck et al., 2020).An important eco-epidemiological question is how changes in wildlife community composition affect the densities of infected ticks, and hence tick-borne disease hazard.
As humans are most frequently infected by the nymphal stage of I. ricinus, the local infection hazard is quantified as the density of infected nymphs (DIN).DIN is the product of two factors: the density of questing nymphs (DON) and the nymphal infection prevalence (NIP).As I. ricinus ticks need vertebrate blood in all three life stages (Földvári, 2016), DON depends on the densities of suitable host species.The immature stages of I. ricinus feed on a wide variety of species, including small to medium-sized mammals, reptiles and ground-foraging birds.Adult stages almost exclusively feed and propagate on larger mammals, such as roe deer.The presence and density of these species affects DON (Hofmeester et al., 2017).I. ricinus spends more than 90% of its life cycle in the vegetation.To maintain homeostasis and prevent dehydration during the long offhost periods, favourable microclimatic conditions are crucial for tick survival and activity.As a consequence, I. ricinus is restricted to forests, shrub and grassland with dense vegetation, where humidity at ground level is >80%.(Estrada-Peña et al., 2013).
While the availability of host species is essential for tick survival and, thus, affects DON, the relative densities of these species within the host community play a key role in pathogen transmission dynamics, and hence NIP.Host species differ considerably in their reservoir competence, that is, the ability of a species to effectively maintain and transmit a pathogen to uninfected ticks.Realized reservoir competence also considers the species-specific average individual tick burden, an indication of the number of ticks successfully feeding on the host (Brunner & Ostfeld, 2008).Competent reservoir species for B. burgdorferi s.l. are rodents and ground-foraging birds (Hofmeester et al., 2016).In contrast, the main reservoir hosts for A. phagocytophilum in Europe are ungulates and medium-sized mammals like the European hedgehog Erinaceus europaeus (Fabri et al., 2021;Silaghi et al., 2012).The relative contribution of these host species to feeding ticks determines the infection prevalence.Local variation in host species assemblage composition is, therefore, likely to influence tick-borne disease hazard (Krawczyk et al., 2020).
While high relative densities of competent reservoir hosts increase NIP in questing I. ricinus, incompetent species can have the opposite effect.Here, special attention should be directed at the role of lizards, as the majority of lizard species appear to be refractory to infection with most B. burgdorferi genospecies, the Mediterranean variant B. lusitaniae being the exception (Lane, 1990;Majláthová et al., 2006;Tijsse-Klasen et al., 2010).Where B. lusitaniae is absent, such as Central and Western Europe, lizards may function as dilution hosts, so that NIP should be lower in species assemblages where lizards feed a large proportion of larvae.Several studies have been conducted to investigate this effect, focusing on different geographical scales and lizard species.Some studies have found a negative correlation between lizard density and NIP, whereas others have found no correlation (Mannelli et al., 2012;Ostfeld & Keesing, 2000;Tijsse-Klasen et al., 2010;Wright et al., 1998).It remains unclear what the eco-epidemiological role of lizards is in the transmission of tick-borne pathogens.The contradictory nature of previous findings may be due to differences in host species community that were not accounted for.Because the role of lizards as hosts for I. ricinus larvae is always relative to the role of other host species in the same area, it is essential to consider the composition of the local host species assemblage as a factor in the analysis.
Here, we quantify the role of sand lizards (Lacerta agilis) as a potential dilution host for B. burgdorferi s.l. and A. phagocytophilum in a heterogeneous landscape, while taking into account local community competence (CC).L. agilis is an important host for larvae and a public health perspective, this underlines the importance of preserving early successional habitat, as encroaching shrubs are associated with higher tick-borne disease hazard, and vegetation removal might be a solution to reduce hazard in coastal dunes.The high degree of spatial heterogeneity in the abundance of tickborne pathogens also poses opportunities to manage recreational activities to limit human exposure to tick-borne diseases.
community competence, dilution host, Europe, habitat suitability, Ixodes ricinus, Lyme borreliosis, reservoir, vector-borne diseases nymphs of I. ricinus in central and Western Europe, while being an incompetent host for both B. burgdorferi s.l. and A. phagocytophilum (Tijsse-Klasen et al., 2010).We sampled questing I. ricinus from the vegetation across six different habitats along a successional gradient in a coastal dune ecosystem associated with varying densities of L. agilis.Per habitat type, we calculated CC based on the estimated density and tick burden of L. agilis and other relevant host species.
Furthermore, we estimated DIN for each habitat type based on tick density and infection prevalence to show the spatial variation in tickborne disease hazard in a mosaic landscape.To accurately determine local tick burden of L. agilis and confirm its reservoir incompetence, we tested feeding ticks collected from local L. agilis for our focal pathogens.To our knowledge, this study is the first to assess the role of L. agilis as a dilution host in the context of local CC for B. burgdorferi s.l. in the field and the first to link local prevalence of A. phagocytophilum to the abundance of lizards.This is a key step to further understand how these pathogens persist in local wildlife assemblages and which biotic factors determine tick-borne disease hazard.and open shrub (sparse shrub and herbs).See Appendix S1 for a detailed description of the habitat types.

| Collection of questing Ixodes ricinus from vegetation
In total, we collected questing nymphs from 83 sites (forest:  (Gray & Lohan, 1982).At 5 m intervals, nymphs and adults attached to the fabric were collected.Larvae were counted but not collected.Because the transects did not always generate sufficient ticks to estimate NIP, we continued sampling around the transects until at least 100 individual nymphs were collected or until 2 h of blanket-dragging had passed.This additional sampling did not count towards our estimates for DON.

| Collection of ticks from Lacerta agilis
We regularly searched for L. agilis in suitable habitat throughout the entire study area at and around all transects.Individual lizards were captured either by hand or with a noose, taking care to handle the animal cautiously to prevent autotomy of the tail.Adults were sexed by visual inspection, for subadults we followed the method of Eplanova & Roitberg, 2015.Lizards were weighed after tick removal, using a precision scale (On Balance Tuff-Weigh 100 × 0.01 g).Larvae and nymphs were removed using forceps and stored individually in 0.2 mL Eppendorf PCR tubes.All lizards were released immediately after the procedure.All methods involving animals were approved by and carried out in accordance with the guidelines and regulations of permit AVD261002016604 issued by the ethical committee of RAVON and permit FF/75A/2016/015, issued by the Netherlands Enterprise Agency.

| Determining species assemblage
In addition to L. agilis, 10 other species (five mammals and five birds) were included in the analysis as relevant hosts for I. ricinus larvae (Table 1).We combined information from literature and local species monitoring datasets to estimate the relative densities of host species per habitat type (Table S2.1).For all species except L. agilis, we retrieved values for larval burden and pathogen prevalence in feeding larvae from a recent metanalysis of the literature by Fabri et al. (2022) and Hofmeester et al. (2016).See Appendix S2 for more details.

| Calculation of community competence and realized community competence
To estimate CC, we followed the framework of Hofmeester et al. (2016), which is an adaptation of the equations in Mather et al. (1989).Relative importance (RI i ) of each host species in feeding I. ricinus larvae was calculated using Equation (1); where D i is the population density, LB i the mean individual larval burden and ∑ n j D j LB j the total amount of larvae fed by the host community.CC is calculated using Equation (2).Multiplying the reservoir competence of each species (RC i ) with LB i and D i yields the total number of larvae infected by each species.CC is the sum of all infected larvae divided by the total number of fed larvae.
Local CC was calculated for each of the six habitat types, representing the proportion of fed larvae that become infected given the local species assemblage.We then calculated realised community competence (RCC; Equation 3), which reflects the absolute number of fed larvae that become infected and is theoretically comparable with DIN.Division by 10,000 is necessary to scale RCC from ind./km 2 to ind./100 m 2 to be able to compare it with measured DIN.

| Data analysis
We screened individual ticks for B. burgdorferi s.l. and A. phagocytophilum.Positive samples for B. burgdorferi s.l. were subsequently analysed to genotype, samples positive for A. phagocytophilum were analysed to ecotype (Jahfari et al., 2014).Ticks collected from L. agilis were also tested for Borrelia miyamotoi, Neoehrlichia mikuriensis, Babesia microti and Spiroplasma ixodetis.
See Appendix S2 for a description of the molecular protocols used.
All statistical analyses were performed using R version 4.1.2(R Core Team, 2021).We used QGIS for spatial data analysis (QGIS.org, 2021).Hourly measurements of saturation deficit (SD) were obtained from the Royal Netherlands Meteorological Institute (KNMI, 2022a).We averaged data of the two weather stations closest to the study area: Wijk aan Zee (~18 km) and Voorschoten (~24 km).
We used a negative-binomial generalized linear model (GLM) with a log-link to test the effects of habitat type, day of year (to account for seasonality), weight and sex on individual tick burden of L. agilis (larvae and nymphs combined).Similarly, a negative-binomial GLM with a log-link was used to test for an effect of habitat type and possible confounding effects of SD and day of year on DON (n/100m 2 ) at each sampling site.To test whether NIP in questing ticks was associated with CC, we used binomial generalized linear mixed modelling (GLMM) to analyse the log-likelihood of nymphs testing positive for either pathogen.Transect was included as a random effect.Models including possible confounding effects of day of year and the CC of surrounding habitat (100-m radius) on NIP are included in Appendix S4.We used packages MASS for GLMs (Venables & Ripley, 2002) and lme4 for GLMMs (Bates et al., 2015).

| Pathogens detected in ticks feeding on Lacerta agilis
We collected 1038 larvae and 324 nymphs feeding on L. agilis, of which 1031 larvae and all nymphs were successfully analysed.
Only two larvae and one nymph were infected with B. burgdorferi s.l., each collected from a different individual.A. phagocytophilum was detected in five larvae and two nymphs, collected from four different lizards.Results for other pathogens are included in Table S3.2.

| Relative importance of host species
The relative importance of L. agilis in feeding larvae was largest in grassland and open shrub (Figure 2a).In both habitat types, they accounted for more than 80% of all fed larvae.In less suitable habitat types for L. agilis, such as dense shrub and forest, rodents accounted for the highest proportion of fed larvae.Other vertebrates were of minor importance in comparison with L. agilis, A. sylvaticus, M. glareolus and E. erinaceus.

Reservoir competence (n) b
Borrelia burgdorferi s.l.Of all 320 I. ricinus nymphs positive for B. burgdorferi s.l., we could identify 174 samples to genotype.The most prominent genotype found was Borrelia afzelii (75%) followed by B. spielmanni (14%;

Reptiles
Table S4.4).For A. phagocytophilum, 136 out of 162 positive samples were identified as ecotype I. Two samples were identified as ecotype II, and 24 samples could not be classified as either ecotype (Table S4.5).
Pathogens detected in adult ticks are included in Table S4.6.

| Density of infected nymphs
Bracken had the highest density of nymphs infected with B. burgdorferi s.l.(Figure 4c), almost twice the density found in forest.
Lowest DIN was measured in grassland, open shrub and open forest.
Density of nymphs infected with A. phagocytophilum was highest in bracken and forest.DIN in the other habitat types was much lower.
For both pathogens, the DON had a stronger effect on DIN compared with NIP, since the variation in DON between habitat types was higher (Figure 5).Measured DIN was much lower than the pre- dicted RCC in open shrub and grassland, whereas it was much higher than RCC in bracken and forest habitats (Figure 4c,d).Thus, the contrasts between open and more densely vegetated habitat types were more pronounced than expected.

| DISCUSS ION
L. agilis abundance in the study area strongly reduced estimated CC.
Indeed, measured NIP was lowest in habitat types where L. agilis fed a large proportion of I. ricinus larvae.These findings suggest that L.
agilis abundance plays a significant role in reducing the prevalence of both pathogens in questing nymphs.We found that overall, DON most strongly influenced DIN, which clearly shows the importance of abiotic conditions in determining tick-borne disease hazard.

| Mean tick burden and pathogens found in ticks collected from Lacerta agilis
Tick burden of L. agilis was considerably higher than previous estimates, indicating that L. agilis is indeed a relevant host for I. ricinus in the study area (Richter & Matuschka, 2006;Tijsse-Klasen et al., 2010).Mean tick burden was higher in individuals captured in habitat types where DON was also high, suggesting that L. agilis encounters more ticks where they are more abundant, in accordance with previous work on tick burdens of P. major (Heylen et al., 2013).
Participation of adult L. agilis in reproductive behaviour may also lead to increased tick burden, as shown by Wieczorek et al. (2020).
This may explain why males have a higher tick burden than females.
During the mating season, males are generally more mobile, defending a territory and actively guarding females (Olsson, 1993).
Additionally, higher testosterone levels may lead to compromised immune function, known as the immunocompetence-handicap (Hughes & Randolph, 2001;Olsson et al., 2000).As expected, only a fraction of the larvae feeding on L. agilis was infected with B. burgdorferi s.l. or A. phagocytophilum.While we cannot rule out that the lizard infected these larvae, it is more likely that the infection stems from previous (incomplete) feeding on an infected host or from cofeeding, which may occur sporadically (Voordouw, 2015).

| Nymphal density and the effects of weather and season
We found that the DON in the study area increased along the successional gradient, likely a direct result of microclimatic conditions.We suspect that host availability played a minor role in determining DON, due to the high abundance of fallow deer, which are important propagation hosts (Hofmeester et al., 2017).DON was by far highest in bracken.Due to the limited spatial extent of this habitat type in the study area, we were only able to sample three bracken sites.Nonetheless, we decided to include this habitat type because it is relevant to site management.The fern profits from atmospheric deposition of nitrogen, while being avoided by herbivores due to its toxicity (Ouden, 2000).However, vegetation structure may also have biased measured DON, as drag sampling success depends on the contact between cloth and vegetation (Tack et al., 2011), which may have led to an underestimation of DON in thorny H. rhamnoides shrub.Besides habitat type, weather conditions and season strongly influence tick survival and activity (Randolph & Storey, 1999).As expected, DON decreased from the start of the field season until the end (Gassner et al., 2011).Here, it must be noted that our data represent a snapshot in time, and differences between years could not be accounted for.In the Netherlands, 2021 was characterized by an unusually cold spring and wet summer, while the previous 3 years were exceptionally dry (KNMI, 2022b).Climate change scenarios predict irregular weather events to increase in frequency, underlining the necessity to conduct multi-season field studies to acquire comparable datasets (IPCC, 2022).

| NIP and the role of Lacerta agilis in the context of community competence
Our estimation of CC was largely determined by the relative densities of the four main host species.Other host species were of comparatively little relevance, in agreement with the findings of Hofmeester et al., 2016.By deriving mean tick burden and infection prevalence of the host assemblage from literature, we accept the possibility that these estimated values differ in our study area due to local conditions that we did not account for.For example, it is likely that, as observed here in L. agilis, tick burdens of other host species is lower in habitat associated with lower DON (Heylen et al., 2013).This might explain why we consistently overestimate CC for both pathogens.Even with some variation in the actual contribution of host species to the CC, we expect our results to be robust (Table S2.1; Figure S2.1).
The fact that differences in NIP could be observed on a habitat level despite the high degree of spatial heterogeneity indicates that NIP is determined at a fine spatial resolution.CC of surrounding habitat had no effect on NIP, which leads us to conclude that spillover from adjacent habitat patches is limited.A possible follow-up would be to measure the relative abundance and larval burden of relevant host species as well as NIP at a finer spatial resolution, so that also the indirect effects of spatial heterogeneity and the host community structure on DON and NIP can be estimated.In any case, the prevalence of both pathogens was lower in habitats where L. agilis fed the majority of I. ricinus larvae.This strongly suggests that L. agilis acted as the main dilution host for both B. burgdorferi s.l.(5) and density of infected nymphs for A. phagocytophilum (6).and A. phagocytophilum.While L. agilis inhabits a wide range of open habitats, the coastal dune ecosystem may be uniquely suited to observe this effect in situ, as rodent densities are low in the early successional habitats where L. agilis is abundant.In other ecosystems, where rodents and L. agilis occur sympatrically, differences in NIP might be more difficult to detect.
In epidemiological research, it is often stated that high biodiversity reduces CC for pathogens that rely on a natural reservoir (Randolph & Dobson, 2012).However, most authors emphasize that species composition and relative abundances of specific species predict NIP more accurately (LoGiudice et al., 2008).Here, we found large variation in NIP in a relatively small host species community.NIP then depends on the competence of the most prominent host species.Our results present an example of how the dominance of an incompetent reservoir host can significantly lower NIP.Here, the increase in NIP followed the successional gradient, which can be explained by shifts in host species assemblage.In early successional habitat types, most I. ricinus larvae feed on an incompetent host (L.agilis), whereas in late successional habitat types, competent hosts like rodents feed most of the larvae.In the United States, Ginsberg et al. (2021) and Ostfeld and Keesing (2000) have shown that the overall increase in Lyme disease prevalence at higher latitudes may be based on the same mechanism of shifts in tick-host associations from incompetent host-dominated to competent host-dominated habitat.After all, the successional gradient in coastal dunes is not unlike the latitudinal gradient in climate and vegetation.

| Density of infected ticks and implications for hazard
Although our study highlights an important role for L. agilis as dilution host, our findings also indicate that DON is a more important factor than NIP in determining DIN.Given that DON is largely determined by abiotic conditions and habitat structure, this implies that the off-host environment has a larger impact on tick-borne disease hazard than the reservoir competence of the host community (Gandy et al., 2021).The differences between expected DIN (i.e.RCC) and observed DIN clearly illustrate this.In habitat types suitable for I. ricinus (forest, bracken), measured DIN far exceeded RCC, as infected hosts continuously produce infected nymphs, and abiotic conditions are optimal for tick survival and activity.In habitat types unsuitable for I. ricinus (open shrub, grassland), the relationship was reversed, and measured DIN was much lower than the host community could theoretically produce, most likely because survival and activity of I. ricinus are limited.

| CON CLUS ION
The current study shows how habitat heterogeneity shapes the local host species assemblage and thereby the CC for tick-borne diseases.
Incompetent hosts (L.agilis) are most prominent in early successional stages but are replaced by competent hosts (rodents) in late successional stages, resulting in an increase in infection prevalence along the successional gradient.Thus, L. agilis acts as a dilution host in early successional habitat.However, tick density plays an even more important role than infection prevalence in determining disease hazard.Early successional habitat is characterized by low tick density, as microclimatic conditions are unfavourable for I. ricinus.The conservation and restoration of nutrient-poor habitat types in coastal dunes and similar ecosystems can, therefore, lower human infection hazard for tick-borne pathogens both directly (through lowering nymphal density) and indirectly (through lowering CC, resulting in lower infection prevalence).In the present study, the highest DIN was measured in bracken, making this the most relevant habitat type for hazard reduction.As bracken is invasive in Dutch coastal dunes, active vegetation management by mowing and removal is an option for reducing hazard to the public.Additionally, closing or constructing walking paths in-or outside of hazardous areas may reduce risk of tick bites for humans.The hazard maps created in the scope of this study are an example of how the results of ecological research may serve as a tool for the management of natural areas.

2. 1 |
Field data collection 2.1.1 | Study area Fieldwork was conducted between May and September 2021 in the Amsterdam Water Supply Dunes, located at the Dutch North Sea coast (Figure 1).Permission to enter and collect samples was granted by Waternet, on behalf of the municipality of Amsterdam.Most of the area consists of fixed grey dunes, partly covered with Hippophae rhamnoides shrubland and dune forests.The host species community for I. ricinus is relatively limited.Dama dama are the only ungulate species in the area, as roe deer (Capreolus capreolus) have disappeared almost entirely.The only rodent species found in significant densities are wood mice (Apodemus sylvaticus) and bank voles (Myodes glareolus).Lacerta agilis is the only reptile species present.Six main habitat types can be distinguished: forest (closed canopy), dense shrub (with high amount of undergrowth), bracken (Pteridium aquilinum vegetation outside of forests), open forest (sparse canopy), grassland (herbaceous vegetation with occasional trees and shrubs) dense shrub: n = 15, bracken: n = 3, open forest: n = 14, open shrub: n = 18, grassland: n = 10).Sampling points were selected using a stratified random sampling approach.Location and sequence of sampling points were assigned using a GIS random sampling tool (QGIS.org,2021).We maintained a minimum distance of 100 m between sampling points and 20 m to the nearest edge of the habitat type to account for possible edge effects.Tick collection only took place on days without rain and air temperatures above 10°C (Perret et al., 2000).At each site, DON was determined by dragging a 1 m 2 white cloth across a transect of F I G U R E 1 Reference map of the study area.Points indicate the placement of sampling points within the area.13652664,2023, 6, Downloaded from https://besjournals.onlinelibrary.wiley.com/doi/10.1111/1365-2664.14379by CochraneAustria, Wiley Online Library on [21/08/2023].See the Terms and Conditions (https://onlinelibrary.wiley.com/terms-and-conditions)on Wiley Online Library for rules of use; OA articles are governed by the applicable Creative Commons License 100 m

F
Mean and standard deviation of the density of Ixodes ricinus nymphs (n/100 m 2 ) per habitat type.Original measurements are added as grey dots for reference.

F
Upper bar plots showing nymphal infection prevalence of (a) Borrelia burgdorferi s.l. and (b) Anaplasma phagocytophilum.The values for community competence (CC) represent the result of our calculations.Points representing transects with a sample size lower than two were removed from the plot to improve readability (n = 4).Lower bar plots show measured density of infected nymphs for (c) B. burgdorferi s.l. and (d) A. phagocytophilum.The realized community competence (RCC) represents the expected density of infected larvae produced by the host community per 100 m 2 .F I G U R E 5 Maps showing the different habitat types in the study area (1), the mean density of questing nymphs (2), the nymphal infection prevalence (NIP) for Borrelia burgdorferi s.l.(3), NIP for Anaplasma phagocytophilum (4), the density of infected nymphs for B. burgdorferi s.l.