
==== Front
Curr Res Parasitol Vector Borne Dis
Curr Res Parasitol Vector Borne Dis
Current Research in Parasitology & Vector-borne Diseases
2667-114X
Elsevier

S2667-114X(24)00040-2
10.1016/j.crpvbd.2024.100209
100209
Research Article
Spatial and temporal variation of five different pathogens and symbionts in Ixodes ricinus nymphs in the Netherlands
Hartemink Nienke nienke.hartemink@wur.nl
a⁎
Gort Gerrit a
Krawczyk Aleksandra I. b
Fonville Manoj c
van Vliet Arnold J.H. d
Takken Willem e
Sprong Hein hein.sprong@rivm.nl
ce⁎⁎
a Biometris, Wageningen University & Research, Droevendaalsesteeg 1, 6708 PB, Wageningen, the Netherlands
b Jacob Blaustein Center for Scientific Cooperation, The Jacob Blaustein Institutes for Desert Research, Ben-Gurion University of the Negev, Midreshet Ben-Gurion, Israel
c Centre for Infectious Disease Control, National Institute for Public Health and the Environment, Antonie van Leeuwenhoeklaan 9, 3720 BA, Bilthoven, Bilthoven, the Netherlands
d Earth Systems and Global Change Group, Wageningen University & Research, Droevendaalsesteeg 3a, 6708 PB, Wageningen, the Netherlands
e Laboratory of Entomology, Wageningen University & Research, Droevendaalsesteeg 1, 6708 PB, Wageningen, the Netherlands
⁎ Corresponding author. Biometris, Wageningen University & Research, Droevendaalsesteeg 1, 6708 PB, Wageningen, the Netherlands. nienke.hartemink@wur.nl
⁎⁎ Corresponding author. Centre for Infectious Disease Control, National Institute for Public Health and the Environment, Antonie van Leeuwenhoeklaan 9, 3720 BA, Bilthoven, Bilthoven, the Netherlands. hein.sprong@rivm.nl
21 8 2024
2024
21 8 2024
6 10020927 5 2024
25 7 2024
16 8 2024
© 2024 The Authors
2024
https://creativecommons.org/licenses/by-nc-nd/4.0/ This is an open access article under the CC BY-NC-ND license (http://creativecommons.org/licenses/by-nc-nd/4.0/).
The incidence of diseases caused by pathogens transmitted by the tick Ixodes ricinus vary over time and space through incompletely understood mechanisms. An important determinant of the disease risk is the density of infected ticks, which is the infection prevalence times the density of questing ticks. We therefore investigated the spatial and temporal variation of four pathogens and one of the most abundant symbionts in Ixodes ricinus in questing nymphs over four years of monthly collections in 12 locations in the Netherlands. The infection prevalence of all microbes showed markedly different patterns with significant spatial variation for Borrelia burgdorferi (s.l.), Neoehrlichia mikurensis, Rickettsia helvetica, and Midichloria mitochondrii, significant seasonal variation of B. burgdorferi (s.l.), N. mikurensis, and M. mitochondrii and a significant interannual variation of R. helvetica. Despite its ubiquitous presence, no spatio-temporal variation was observed for the infection prevalence of B. miyamotoi. The variation in infection prevalence was generally smaller than the variation in the density of nymphs, which fluctuated substantially both seasonally and between locations. This means that the variation in the densities of infected nymphs for all pathogens was mostly the result of the variation in densities of nymphs. We also investigated whether there were positive or negative associations between the symbionts, and more specifically whether ticks infected with vertically transmitted symbionts like M. mitochondrii and R. helvetica, have a higher prevalence of horizontally transmitted symbionts, such as B. burgdorferi (s.l.) and N. mikurensis. We indeed found a clear positive association between M. mitochondrii and B. burgdorferi (s.l.). The positive association between R. helvetica and B. burgdorferi (s.l.) was less clear and was only shown in two locations. Additionally, we found a clear positive association between B. burgdorferi (s.l.) and N. mikurensis, which are both transmitted by rodents. Our longitudinal study indicated strong between-location variation, some seasonal patterns and hardly any differences between years for most symbionts. Positive associations between symbionts were observed, suggesting that infection with a (vertically transmitted) symbiont may influence the probability of infection with other symbionts, or that there is a common underlying mechanism (e.g. feeding on rodents).

Graphical abstract

Image 1

Highlights

• The infection prevalence of five Ixodes ricinus symbionts in the Netherlands varies.

• Spatial variation (between locations) was observed for four symbionts.

• Midichloria mitochondrii and Borrelia burgdorferi (s.l.) prevalence was highest in summer.

• Densities of infected nymphs are higher in summer, due to higher nymph densities.

• M. mitochondrii and B. burgdorferi (s.l.) were positively associated at tick level.

Keywords

Acarological risk
Lyme borreliosis
Neoehrlichiosis
Relapsing fever borreliosis
Spotted fever rickettsiosis
Prevalence
Seasonality
==== Body
pmc1 Introduction

The incidence of tick-borne diseases (TBDs) has been predicted to increase in Europe due to a combination of human behaviour, and climatic and environmental changes (Medlock et al., 2013; Köhler et al., 2023). Annually, tens of thousands of people in Europe acquire Lyme borreliosis (LB) and tick-borne encephalitis (Dobler et al., 2020; Stark et al., 2023). In contrast, the incidences of other TBDs in Europe are currently unknown (Azagi et al., 2020; Hoornstra et al., 2022). These TBDs are caused by a plethora of (potential) tick-borne pathogens (TBPs) transmitted by Ixodes ricinus, including Borrelia burgdorferi (sensu lato), Anaplasma phagocytophilum, Babesia divergens, Babesia microti, Babesia venatorum, Borrelia miyamotoi, Neoehrlichia mikurensis, Rickettsia helvetica, Rickettsia monacensis, and Spiroplasma ixodetis (Sprong et al., 2018). Previous studies have shown that all these (potential) TBPs are present in I. ricinus in the Netherlands (Sprong et al., 2009; Coipan et al., 2013; Esser et al., 2022; Krawczyk et al., 2022a).

With annually approximately one and a half million tick bites in the Netherlands alone (Hofhuis et al., 2015) and roughly one third of the ticks being infected with at least one (potential) TBP other than B. burgdorferi (s.l.), numerous individuals are exposed to TBPs each year (Jahfari et al., 2016). Except for Lyme borreliosis, tick-borne encephalitis and a few cases of B. miyamotoi disease (Hoornstra et al., 2022) and one case of human granulocytic anaplasmosis in 1999 (van Dobbenburgh et al., 1999), other I. ricinus-borne diseases in humans have not been reported in the Netherlands. Lack of awareness, case definitions, laboratory diagnostics, as well as a non-characteristic clinical presentation are among the reasons why these other TBDs often go undiagnosed. Another explanation is that the risk of infection with these (potential) TBP is very low or far more heterogeneously distributed than B. burgdorferi (s.l.), the causative agents of Lyme borreliosis. Indeed, B. burgdorferi (s.l.) has been found in all locations sampled for ticks in the Netherlands (Tijsse-Klasen et al., 2010; Coipan et al., 2013; Hartemink et al., 2021; Krawczyk et al., 2022a). The risk of acquiring a TBD depends on many different factors, of which the most important are the density of questing I. ricinus ticks infected with one or more TBP, and more specifically the density of infected nymphs (DIN), since nymphs are responsible for the majority of human bites (Hofhuis et al., 2017) and the level of human exposure to ticks (Lambin et al., 2010; Zeimes et al., 2014). The DIN is in turn determined by the density of nymphs (DON) and the nymphal infection prevalence (NIP). The DON can here be considered as a proxy for the probability of acquiring a tick bite, whereas the NIP can be regarded as an estimate of the probability that the bite is from an infected tick. The NIP is thus an important determinant of the risk of acquiring a TBD.

Insights in the spatial and temporal variation of this prevalence for the different pathogens would be useful to gauge where and when the risk of getting infected is the highest. While the spatio-temporal variation in prevalence in the Netherlands has been studied for B. burgdorferi (s.l.) (Hartemink et al., 2021) to some extent, not much is known about spatio-temporal variation in several other (potentially pathogenic) bacteria that depend on I. ricinus for their survival and propagation, i.e. B. miyamotoi, N. mikurensis, R. helvetica, and M. mitochondrii. Whether or not they cause disease in humans or not, these microorganisms can be surely considered bacterial symbionts due to the intimate relationship with their I. ricinus host. For convenience, we therefore call them ‘symbionts’ here, irrespective of their pathogenicity for humans. In this paper, we aim to determine the spatial and temporal variation in the infection prevalence of these five I. ricinus symbionts. The symbionts vary in their transmission routes. Rickettsia helvetica and M. mitochondrii are found in all tick life stages and are thought to be mainly vertically transmitted (Krawczyk et al., 2022b), while B. burgdorferi (s.l.), and N. mikurensis are horizontally transmitted, between I. ricinus and their hosts (mainly rodents). For B. miyamotoi, both horizontal and vertical transmission have been reported (Burri et al., 2007; Rollend et al., 2013). The extent of influence of abiotic and biotic factors on the prevalence of these microorganisms likely varies depending on the type of symbiont-tick relationship. Vertically transmitted symbionts are believed to play a role in the fitness or reproduction of the ticks, and therefore, all factors influencing ticks may affect the degree of transovarial and transstadial transmission of symbionts (Corbin et al., 2017). The prevalence of horizontally transmitted symbionts has been shown to be primarily affected by the presence and abundance of vertebrate reservoir hosts (Gilbert et al., 2012; Hofmeester et al., 2017; Takumi et al., 2019). In addition, symbionts, which are transmitted from infected nymphs to uninfected larvae via vertebrates, may be influenced by climatic conditions that dictate the seasonal synchrony of tick life stages. Finally, the symbionts might directly or indirectly affect each other and either increase or decrease the prevalence of co-infections (Vautrin and Vavre, 2009). Whether they do is poorly investigated and the subject of this study.

2 Materials and methods

2.1 Data collection

The nymphs used in this study were systematically collected in 12 sites in the Netherlands within the framework of the Nature’s Calendar project, for which details are provided elsewhere (Gassner et al., 2011; Hartemink et al., 2019, 2021; Takken et al., 2017a). In total 15,090 questing I. ricinus nymphs were collected between January 2012 and December 2015. The 12 sites (with abbreviations used in the tables and figures) are Bilthoven (BI), Dronten (DR), Ede (ED), Gieten (GI), Hoog Baarlo (HB), Kwade Hoek KH), Montferland (MO), Schiermonnikoog (SC), Twiske (TW), Vaals (VA), Veldhoven (VE) and Wassenaar (WA) (Fig. 1). Of these 15,090 collected nymphs, 8230 were tested for presence of symbionts. Given the low numbers of tick captures in the winter months (December–February) and the fact that this hampers the estimation of prevalence, we excluded the observations from these months from the analysis. As shown in Supplementary file 1: Table S1, the effect on the number of observations was very small. The analysis in this paper is thus performed on the 8104 tested nymphs out of a total of 14,948 nymphs that were captured in the months March-November of the years 2012–2015. Questing ticks were collected every month at each site by dragging a white cotton cloth (1 m2) over two transects of 100 m2 each, for a total area of 200 m2 per location (Hartemink et al., 2021). The cloth was inspected every 25 m, and ticks were collected using forceps. All samples were placed in individual containers with 70% ethanol, counted, and stored at −20 °C until DNA extraction. The ixodid species and life stage identification were performed at the Entomology group of Wageningen University.Fig. 1 Overview of the locations of the sampling sites in the Netherlands.

Fig. 1

2.2 DNA extraction and bacteria detection

Collected nymphs were processed individually. DNA from a subset of nymphs per month per location was extracted with ammonium hydroxide as previously described (Wielinga et al., 2006). The lysates were stored at 4 °C prior to molecular analysis. Samples were analysed with different (multiplex) real-time PCRs, based on various target genes from the microorganism of interest, namely B. burgdorferi (s.l.) (Heylen et al., 2013), B. miyamotoi (Hovius et al., 2013), N. mikurensis (Jahfari et al., 2012), R. helvetica (de Bruin et al., 2015), and M. mitochondrii (Garcia-Vozmediano et al., 2022). All qPCRs were carried out on a LightCycler 480 (Roche Diagnostics Nederland B.V, Almere, The Netherlands) in a final volume of 20 μl with iQ multiplex Powermix, 3 μl of sample (∼5% of tick lysate), primers with end concentration of 0.2 μM, and probes. Positive plasmid controls and negative water controls were used on every plate tested. Cycling conditions included an initial activation of the iTaq DNA polymerase at 95 °C for 5 min, followed by 60 cycles of a 5 s denaturation at 95 °C followed by a 35 s annealing-extension step at 60 °C (ramp rate 2.2 °C s−1 and a single point measurement at 60 °C) and a final cooling cycle of 37 °C for 20 s. To minimize contamination and false-positive samples, the DNA extraction, PCR mix preparation, sample addition, and qPCR analyses were performed in separated air-locked dedicated laboratories.

2.3 Statistical analysis

The analysis was performed in the statistical program R (R Core Team, 2023). We investigated the spatial and temporal variation in nymphal infection prevalence (NIP, for which results are in the main text), as well as in the DON and DIN (with results in the Supplementary file 1). The NIP was calculated as the number of positive nymphs divided by the number of tested nymphs. The DIN was calculated by multiplying the infection prevalence (NIP) and the density of nymphs (DON). The DIN and DON were expressed as the number of nymphs per 100 m2.

We used generalized linear mixed models (GLMMs) to analyse how the NIP, DON and DIN depend on the year, month and location. GLMMs were used because this type of model can handle non-normally distributed data such as count data as well as (multiple) random effects. We used the glmmTMB package (Brooks et al., 2017), which handles betabinomial distributions as well as crossed and nested random effects. For the NIP we used the model that assumes a beta-binomial distribution, and a logit link function. We included the variables Location, Year and Month as fixed effect. Locations within Year were included as random effect (1 | Year: Location), to account for the fact that observations from the same year-location combination are not independent. We tested for each fixed factor, Location, Year and Month, whether they had any predictive power in the model, by applying a Likelihood Ratio Test (LRT) on the single term deletions (using the drop1() function). Estimated marginal means were calculated for each month, year and location, using the emmeans package (Lenth, 2020). Estimated marginal means should be interpreted as predicted averages over all factors (based upon the fitted model); e.g. the estimated marginal mean NIP for a specific month is averaged over all years and all locations. Given that there were only four years, we did not test for (linear) trends over the years by running the GLMM models with year as numerical fixed factor. The DON and DIN were analysed using a similar GLMM model, with a negative binomial distribution, and a log link function.

We assessed whether the observed spatial and temporal variation in the NIP could be explained by climatic factors. To this end, weather data were obtained from the data center of the Royal Netherlands Meteorological Institute (KNMI). Daily values for minimum and maximum temperature and relative humidity, based on interpolated data from 34 automatic weather stations, were extracted for each of the 12 locations. The mean temperature was calculated from the minimum and maximum temperature. From the mean temperature and relative humidity, we calculated the saturation deficit according to the definition by Randolph and Storey (1999). We tested whether using the mean temperature (Temperature) and saturation deficit (SD) in the months before the sampling as explanatory variables, improved the model fit. We did this by comparing several models: (1) a null model, with only a random effect for Location within Year; (2) a model with Year, Month and Location; (3) a model with Temperature and SD in the month before the sampling took place; (4) a model with Year, Month, Location and Temperature and SD in the month before the sampling took place; (5) a model with Temperature and SD in the two months before the sampling took place; and (6) a model with Year, Month, Location and Temperature and SD in the two months before the sampling took place. All six models had a random effect for Location within Year.

2.4 Co-infection analysis

Lastly, we analysed whether there were any positive or negative correlations between the symbionts, that is, whether ticks infected with one symbiont were more (or less) likely to be infected with another symbiont.

Co-infection of ticks by pairs of symbionts was tested by fitting generalized linear mixed models to the binary absence/presence records for the two symbionts per tick, using PROC GLIMMIX of the SAS program (SAS Institute Inc, 2017). A binomial distribution (with n = 1) and logit link function was taken. To correct for possible confounders fixed effects of Species, Location, Year and Month combinations and two-way interactions of Species with Location, Year and Month were introduced into the models, and the two binary observations (for the two symbionts) per tick were allowed to correlate. The correlation was tested to be zero, using the conditional independence test as available in the software.

We also calculated the odds ratios for each combination of symbionts, based on the outcomes of logistic regression models where the binomial presence/absence of a symbiont in a tick is predicted using the binary presence/absence score of another symbiont as explanatory variable (with Year, Month and Location as fixed effects, and Location within Year as random effect). The odds ratio for a pair of symbionts, let us for example refer to them as symbiont A and B, represents the ratio of the odds of infection with A given infection with B and the odds of infection with A when B is not detected.

3 Results

The raw number of captured and tested nymphs and the infection percentage for each of the locations are shown in Table 1. Out of 14,948 captured nymphs, a total number of 8104 nymphs was tested. The overall infection percentage was highest for M. mitochondrii (69.7%) and lowest for B. miyamotoi (2.4%). The overall infection prevalences of N. mikurensis, B. burgdorferi (s.l.) and R. helvetica were 7.4%, 12.9% and 19.7%, respectively.Table 1 Raw data for infection percentage per location.

Table 1Location	Total no. of nymphs	No. of tested nymphs	Number of positive nymphs (%)	
M. mitochondrii	B. miyamotoi	N. mikurensis	B. burgdorferi	R. helvetica	
BI	408	352	237 (67.3)	11 (3.1)	29 (8.2)	24 (6.8)	37 (10.5)	
DR	392	360	263 (73.1)	6 (1.7)	18 (5.0)	83 (23.1)	56 (15.6)	
ED	983	721	464 (64.4)	17 (2.4)	62 (8.6)	97 (13.5)	48 (6.7)	
GI	2803	1195	921 (77.1)	21 (1.8)	133 (11.1)	142 (11.9)	227 (19.0)	
HB	1436	816	651 (79.8)	20 (2.5)	30 (3.7)	35 (4.3)	62 (7.6)	
KH	84	73	66 (90.4)	1 (1.4)	9 (12.3)	20 (27.4)	22 (30.1)	
MO	1328	781	560 (71.7)	15 (1.9)	61 (7.8)	107 (13.7)	121 (15.5)	
SC	159	155	87 (56.1)	1 (0.6)	26 (16.8)	25 (16.1)	31 (20.0)	
TW	973	620	385 (62.1)	15 (2.4)	13 (2.1)	52 (8.4)	120 (19.4)	
VA	550	447	294 (65.8)	7 (1.6)	7 (1.6)	60 (13.4)	34 (7.6)	
VE	2899	1251	898 (71.8)	36 (2.9)	126 (10.1)	156 (12.5)	123 (9.8)	
WA	2933	1333	819 (61.4)	42 (3.2)	85 (6.4)	245 (18.4)	716 (53.7)	
Total	14,948	8104	5645 (69.7)	192 (2.4)	599 (7.4)	1046 (12.9)	1597 (19.7)	
Notes: For each location, the total number of nymphs captured, the number of nymphs tested, and the number and percentage (in brackets) of nymphs tested positive for each symbiont are presented.

3.1 Spatial variation in NIP

We modelled the effect of the Month, Year, and Location on the infection prevalence of the nymphs. Except for B. miyamotoi, all symbionts showed significant spatial variation, with Location being a significant variable in the model, as tested by the LRT on the single-term deletions (Table 2). Note that the estimated prevalences obtained from this model are adjusted for Month and Year, and are hence slightly different from the raw data presented above.Table 2 Significance of the Likelihood Ratio Test in the single term deletions.

Table 2	Spatial variation	P-value of LRT test for “Location”	Seasonal variation	P-value of LRT test for “Month”	Interannual variation	P-value of LRT test for “Year”	
M. mitochondrii	Yes	3.37 × 10−6∗∗∗	Yes	4.36 × 10−4∗∗∗	No	7.3 × 10−1	
R. helvetica	Yes	3.18 × 10−14∗∗∗	No	1.2 × 10−1	Yes	1.84 × 10−3∗∗	
B. miyamotoi	No	3.27 × 10−1	No	5.35 × 10−1	No	9.49 × 10−1	
N. mikurensis	Yes	8.02 × 10−5∗∗∗	Yes	1.719 × 10−2∗	No	1.19 × 10−1	
B. burgdorferi (s.l.)	Yesa	5.26 × 10−6∗∗∗a	Yes	3.379 × 10−3∗∗	No	6.47 × 10−1	
Note: ***P < 0.0005, **P < 0.005, *P < 0.05.

a Since the model for B. burgdorferi (s.l.) had convergence issues, the result for Location had to be obtained using Location as random effect, instead of Location with Year.

The estimated prevalence for each of the locations is shown in Fig. 2. The prevalence of B. burgdorferi (s.l.) was highly variable between locations, with very high prevalences in Kwade Hoek and Dronten and low prevalences in Bilthoven and Hoog Baarlo (Fig. 2A). The mean estimated B. miyamotoi prevalence (Fig. 2A) was very low, between 0.6% (in Schiermonnikoog) and 3.2% (in Bilthoven). For M. mitochondrii, the prevalence was high in all locations, ranging between 53% (in Schiermonnikoog) and 88% (in Kwade Hoek) (Fig. 2B). The prevalence of N. mikurensis varied ranged between 1.3% (in Vaals) and 13% (in Schiermonnikoog) (Fig. 2C). The prevalence of R. helvetica was clearly highest in Wassenaar (53%) and ranged between 6% and 28% in the other locations (Fig. 2D).Fig. 2 Estimated marginal mean prevalence per location for B. burgdorferi (s.l.) and B. miyamotoi (A), M. mitochondrii (B), N. mikurensis (C) and R. helvetica (D). The white dots represent the estimated marginal mean prevalence and the bars the asymmetric 95% confidence intervals. Locations are ordered based on the mean estimated B. burgdorferi (s.l.) prevalence (high to low). Note the different scales on the Y-axes.

Fig. 2

3.2 Seasonal variation in NIP

For both B. burgdorferi (s.l.) and M. mitochondrii, we observed a clear seasonal pattern, with an increase in the estimated mean nymphal infection prevalence (NIP) from spring to summer, and a slight decrease from summer to autumn (Fig. 3A and B). The seasonal pattern for N. mikurensis was different; the NIP first increased between March and June, then decreased until August and then increased again until November (Fig. 3C). For B. burgdorferi (s.l.), M. mitochondrii and N. mikurensis, the variable Month was significant, as tested by a LRT on the single term deletions (Table 2), which indicates that Month has predictive power in the model. For R. helvetica, the mean prevalence increased slightly towards late summer (August) and then decreased again (Fig. 3D), but Month was not a significant variable in the model (Table 2). For B. miyamotoi, there is no pattern visible in the prevalence and Month was not significant.Fig. 3 Estimated marginal mean prevalences per month for B. burgdorferi (s.l.) and B. miyamotoi (A), M. mitochondrii (B), N. mikurensis (C), and R. helvetica (D). The white dots represent the estimated marginal mean prevalence and the bars the asymmetric 95% confidence intervals. Note the different scales on the Y-axes.

Fig. 3

3.3 Interannual variation in NIP

For most symbionts, no substantial changes in NIP occurred in the four-year period of this study (Fig. 4A–C). Only for R. helvetica, the prevalence increased slightly over the years (Fig. 4D). The four-year period of the study is too short to test for trends over time (by including year as numerical fixed factor), but we did test whether including the factor Year improved the model using the LRT. This test showed that only for R. helvetica, Year was a significant variable in the model, while this was not the case for the other four symbionts (Table 2).Fig. 4 Estimated mean prevalences for the four different years for B. burgdorferi (s.l.) and B. miyamotoi (A), M. mitochondrii (B), N. mikurensis (C), and R. helvetica (D). The white dots represent the estimated marginal mean prevalence and the bars the asymmetric 95% confidence intervals. Only for R. helvetica (panel D), a significant difference between the years was observed, with increasing prevalences over the years.

Fig. 4

3.4 Climatic factors

We investigated whether the observed spatial and temporal variation could be explained by climatic factors, such as the temperature or saturation deficit (SD) in the months before the sampling, by comparing six different models.

As judged by the AIC values of the models (Supplementary file 1: Table S2), the best-fitting model for most symbionts was model 2, with only Year, Month and Location as explanatory variables. Adding climatic factors (models 4 and 6) or replacing Year, Month and Location by climatic factors (models 3 and 5) resulted in similar or higher AIC values and did not improve the fit.

3.5 The effect of density of nymphs and the infection prevalence on the DIN

The risk for humans to be bitten by an infected tick depends, besides human behaviour (activities, protection, etc.), on the density of infected ticks, and more specifically on the density of infected nymphs (DIN), as tick bites on humans are predominantly done by nymphs (Mather et al., 1996; Eisen et al., 2010; Diuk-Wasser et al., 2012). The DIN is therefore also referred to as the acarological risk of exposure to tick-borne pathogens (Eisen et al., 2010; Pepin et al., 2012; Takken et al., 2017b). While the results of our study show that there is variation in the prevalence (NIP) over space and time, most of the variation in acarological risk (DIN) is due to variation in the DON. The estimated value for DON is shown in Supplementary file 1: Fig. S1 for each location, month and year. In Supplementary file 1: Figs. S2–S6, the NIP (panels A, B and C) and DIN (panels E, F and G) are shown for each pathogen.

For most symbionts, the values for the DIN (panels E, F and G in Supplementary file 1: Figs. S2–S6), had very similar patterns to the patterns in the DON (Supplementary file 1: Fig. S1, panels A, B and C). For R. helvetica, we observed that the high values for DON and NIP cause a high DIN in Wassenaar, while in other locations (Kwade Hoek, Schiermonnikoog and Twiske), high values of the NIP are counteracted by the lower DON values and the overall DIN was not much higher than for locations with lower NIP values (Supplementary file 1: Fig. S6, panels A and D)

This means that as a rule-of-thumb, we can say that for most pathogens, the risk is highest in areas and at times with many questing nymphs present, and that the variation in prevalence (NIP) is less than the variation in the DON.

3.6 Analysis of co-infections

The results of the analysis of the correlations between the pairs of symbionts are shown in Table 3. Correlation estimates above zero indicate a positive association (higher infection prevalence for one symbiont when the other symbiont is present). Estimates below zero indicate a negative association (lower prevalence when the other symbiont is present). The results for the Chi-square test (χ2-value and P-value) are also shown. The last column shows the odds ratios and their confidence intervals. An odds ratio above one indicates a positive association and below one a negative association.Table 3 Results of analysis of associations between symbionts; correlation estimate (CE), Chi-square (χ2-value) test and odds ratios (with 95% confidence intervals).

Table 3Combination of symbionts	CE	χ2-value	P-value	Odds ratio (95% CI)	
B. burgdorferi (s.l.)	N. mikurensis	0.131	139.20	<0.0001	3.1 (2.5–3.8)	
B. burgdorferi (s.l.)	M. mitochondrii	0.085	58.09	<0.0001	1.9 (1.6–2.2)	
B. burgdorferi (s.l.)	R. helvetica	0.039	12.26	0.0005	1.4 (1.2–1.6)	
B. burgdorferi (s.l.)	B. miyamotoi	0.021	3.58	0.0583	1.4 (0.9–2.1)	
N. mikurensis	B. miyamotoi	0.020	3.27	0.0706	2.0 (1.2–3.6)	
N. mikurensis	M. mitochondrii	−0.019	2.84	0.0921	0.8 (0.7–1.0)	
R. helvetica	B. miyamotoi	0.015	1.84	0.1753	1.1 (0.8–1.6)	
M. mitochondrii	B. miyamotoi	−0.011	1.00	0.3163	0.8 (0.6–1.2)	
N. mikurensis	R. helvetica	0.006	0.31	0.5766	1.1 (0.9–1.4)	
R. helvetica	M. mitochondrii	0.004	0.14	0.7034	1.0 (0.9–1.1)	

Overall, we found several significant positive associations between symbionts (Table 3). The strongest and most significant correlations were found between B. burgdorferi (s.l.) and N. mikurensis (P < 0.0001) and between B. burgdorferi (s.l.) and M. mitochondrii (P < 0.0001). A positive correlation was also shown between B. burgdorferi (s.l.) and R. helvetica (P = 0.0005).

The other combinations of symbionts did not show significant associations when using a level of significance of 0.05 (α = 0.05). Some associations had P-values between 0.05 and 0.10: B. miyamotoi and B. burgdorferi (s.l.) (P = 0.0583); B. miyamotoi and N. mikurensis (P = 0.0706) and N. mikurensis and M. mitochondrii (P = 0.0921).

4 Discussion

4.1 Spatial variation and overall prevalences

In this study, based on monthly observations from March till November for four consecutive years in 12 locations in the Netherlands, we explored the spatio-temporal variation in infection prevalence for five symbionts in I. ricinus nymphs. The findings for B. burgdorferi (s.l.) have been described and discussed previously (Hartemink et al., 2021) and here are presented for comparison only. The overall observed prevalence was highest for M. mitochondrii (70%). The overall prevalences for R. helvetica (20%), B. burgdorferi (s.l.) (13%) and N. mikurensis (7%) were much lower, and the lowest prevalence was observed for B. miyamotoi (2%).

While all five symbionts were found in all 12 locations in this study, there were marked differences in the level of spatial variation. Both M. mitochondrii and B. miyamotoi were distributed rather evenly over the locations. In contrast, the prevalence of R. helvetica varied substantially between locations; the prevalence in Wassenaar (53%) was approximately eight times as high as that in Hoog Baarlo, Ede and Vaals (6–7%). The spatial variations in the prevalence of these symbionts are generally in line with previous findings from cross-sectional studies in the Netherlands (Coipan et al., 2013; Cochez et al., 2015; Silaghi et al., 2016; Garcia-Vozmediano et al., 2022; Krawczyk et al., 2022a).

High prevalences (between 60% and 80%) of M. mitochondrii in I. ricinus nymphs have recently been reported in locations in the Netherlands, Belgium, Sweden and Italy (Garcia-Vozmediano et al., 2022). In the same paper, lower prevalences were reported for the UK and Germany (45% and 36%, respectively). The high prevalence of M. mitochondrii may not be surprising, since it is thought to provide fitness advantages for the tick. For instance, it has been shown that M. mitochondrii-free I. ricinus females, which were generated by an antibiotic treatment, produce larvae that are less successful during blood-feeding (Guizzo et al., 2023). The cause for the statistically significant differences in infection prevalence between locations is not known. Possibly, the male-female ratio in local populations plays a role, as the M. mitochondrii symbiont appears to be ubiquitous in the females of I. ricinus and less prevalent in males across its distribution; one study observed a significantly higher prevalence in females (100 %) than in males (44 %) (Sassera et al., 2006). This suggests that a higher prevalence could be an indication of a higher proportion of the nymphs being female, but this is just a hypothesis that needs further research.

The overall prevalence for N. mikurensis of 7% is in line with the average prevalence of around 6–8% reported in a review based on studies from 17 different European countries (Silaghi et al., 2016). In the UK, N. mikurensis has not been reported so far (Silaghi et al., 2016; Olsthoorn et al., 2021). Findings from different sites in Central Europe (Slovakia, Austria and Czech Republic) ranged from 1.1% to 23.5% (Derdáková et al., 2014), indicating that there can be substantial variation between sites.

The findings for R. helvetica were in also accordance with previous studies in the Netherlands (Krawczyk et al., 2022a). For comparison, the prevalence of this symbiont in other European countries was found to be much lower and ranged from 0.04% in Scotland to 9.3% in Germany (Olsthoorn et al., 2021). In one study conducted in southern Scandinavia, the individual nymph prevalence (calculated from pool prevalences) was estimated to range from 0% to 22% depending on the location (Kjær et al., 2020), indicating spatial variation.

The low prevalence of B. miyamotoi in this study is in accordance with the results of a meta-analysis, where the average prevalence of B. miyamotoi in I. ricinus ticks in Europe, based on 84,688 ticks, was 1.0% (Hoornstra et al., 2022). Studies in countries nearby the Netherlands found infection prevalences of 2% in Germany (Crowder et al., 2014) and an even lower prevalence in a study in England where only 2 out of 349 nymphs collected during 2013 and 2014 tested positive for B. miyamotoi (Hansford et al., 2017). It is an intriguing question how this symbiont is still able to persist at all locations, despite the very low prevalence (and hence the low numbers of infected ticks). While its transmission is partly vertical and partly horizontal, both routes seem to rely on small numbers of infections: the transstadial retention rate is low (Lynn et al., 2019) and the short duration of infection in rodents implies that the chances for a tick to bite an infected rodent are low. The latter may be compensated by a high transmission probability per bite (from host to tick). The fact that the prevalence of B. miyamotoi is so low, also means that the numbers of observed infected ticks are very low per location. Due to these low numbers, any spatial variation that would be present, would not be statistically significant.

4.2 Seasonal patterns

The results of our analysis show a striking seasonal pattern with highest infection prevalences in summer for both B. burgdorferi (s.l.) and M. mitochondrii. For B. burgdorferi (s.l.), this pattern had been reported before (for the same locations, in the periods 2006–2011 (Takken et al., 2017b) and 2009–2016 (Hartemink et al., 2021)). In the latter paper, the seasonal pattern was speculated to be linked to seasonality in host availability (e.g. rodent populations are strongly seasonal, with lower numbers in winter) or seasonality in the B. burgdorferi (s.l.) prevalence in hosts (Hartemink et al., 2021). However, our results show that the vertically transmitted symbiont M. mitochondrii, which is not transmitted to or from hosts as far as we know, shows a similar pattern. This suggests that there may be a yet unknown environmental or (micro)climatic factor that affects M. mitochondrii and B. burgdorferi (s.l.) in a similar way. Alternatively, a factor underlying the pattern in M. mitochondrii may indirectly be responsible for the pattern in B. burgdorferi (s.l.), given the positive association between the two symbionts. An association between M. mitochondrii infection and increased tick weight has been found (Bakker et al., 2023), and this higher weight could in turn be associated with the opportunity to take larger blood meals and having a higher chance of acquiring a B. burgdorferi (s.l.) infection from an infected host. However, these associations need to be studied further, both in the laboratory and in the field. As far as we know, there are no other studies on the seasonal variation of M. mitochondrii.

We can speculate that while it may be advantageous for the tick to carry certain symbionts, it may also have fitness costs. For example, M. mitochondrii has been suggested to have a positive effect on the blood meal size that ticks can take (Guizzo et al., 2023), although we have to keep in mind that this result is based on poor blood-feeding of offspring from ticks that have been treated with antibiotics, which maybe have other effects as well. Ticks may also benefit from Borrelia spp. infection by an increased lifespan (more fat and more resistance to desiccation) and by an increased questing period (less need to move to the litter zone to rehydrate), which enhances tick chances to find a host and to subsequently transmit the pathogens (Herrmann and Gern, 2015). However, there are likely also costs that come with infection, as symbionts may require energy for maintenance. How these fitness benefits and fitness costs play out at different life stages and different times of the year is unknown, but if we would like to unravel the seasonality of infection prevalence this would deserve attention.

The observed pattern in B. burgdorferi (s.l.) is in concordance with earlier findings in southern Germany, where nymphal infection prevalence was found to rise from 9.3% in spring to 16.1% in summer (Fingerle et al., 1994), but not with other findings, from e.g. Sweden and Norway, where prevalences in late spring and early summer were found to be higher than in late summer and autumn. More information from other locations would be useful to see how spatially heterogeneous the seasonality patterns are.

For R. helvetica, our study found the highest prevalence in August. A study in Thuringia, Germany, also found the highest prevalence in July and August (Hildebrandt et al., 2010) for Rickettsia spp., but there are other studies that found higher prevalences in May than in summer (Kantsø et al., 2010).

For N. mikurensis the seasonal pattern consisted of a slightly higher prevalence in June than in the late summer. A previous study in a low-prevalence region in Norway found no seasonal patterns and another study in the Netherlands found the highest prevalence in autumn (Coipan et al., 2013).

The low prevalences of B. miyamotoi make it hard to find any significant patterns or associations.

4.3 Co-infection

Symbionts might directly or indirectly affect each other and either increase or decrease the prevalence of co-infections. This is especially interesting if vertically-transmitted symbionts (M. mitochondrii and R. helvetica, and partially B. miyamotoi) might influence the chance of acquiring infection with the horizontally transmitted symbionts like B. burgdorferi (s.l.) and N. mikurensis.

Indeed, we found a significant positive correlation between B. burgdorferi (s.l.) and M. mitochondrii (P < 0.0001), which is in line with previous research (Krawczyk et al., 2022b). Since M. mitochondrii is transmitted vertically and transstadially, it is assumed to be present in the tick throughout the tick’s lifetime. This means that its presence is unlikely to be affected by B. burgdorferi (s.l.), which is acquired by a larva during its first blood meal on an infected host. It is much more likely that it is the other way around; the presence of M. mitochondrii might increase the probability that a tick acquires an infection with B. burgdorferi (s.l.). The fitness effects described above, such as improved blood-feeding of M. mitochondrii-infected larva (Guizzo et al., 2023), may play a role here, as larger blood meals might lead to higher chances of acquiring B. burgdorferi (s.l.) during the feeding.

We found a very strong positive correlation between N. mikurensis and B. burgdorferi (s.l.). Since B. burgdorferi (s.l.) infections in I. ricinus in the Netherlands mainly consist of B. afzelii (Hartemink et al., 2021), which is associated with rodents, this does not come as a surprise, as both symbionts are mainly transmitted between rodents and ticks. Ticks that carry one of these symbionts are very likely to have fed on a rodent at the larval stage. There is a probability that the rodent was co-infected by both N. mikurensis and B. afzelii, which could lead to simultaneous transmission of both symbionts to the tick. Our results are in line with the results of an analysis of co-infections in nymphs collected in 19 sites (different than in this study) in the Netherlands using the Hmsc package, that revealed positive associations between B. burgdorferi (s.l.) and N. mikurensis (Krawczyk et al., 2022b). It is also in accordance with earlier findings, in southern Sweden, where 2.1% of the host-seeking ticks was co-infected with both N. mikurensis and B. afzelii, which is significantly more than expected under random co-occurrence. As the infection intensity (number of bacterial cells) of N. mikurensis was not affected by co-infection with B. afzelii, and vice versa, it was hypothesized that the two symbiont co-occur due to simultaneous transmission of these two tick-borne pathogens between ticks and rodents (Andersson et al., 2013). A study in Sweden also showed a higher-than-expected prevalence of co-infection of N. mikurensis and B. afzelii in bank voles (Andersson et al., 2014).

A positive correlation was also shown between B. burgdorferi (s.l.) and R. helvetica (P = 0.0005). When we applied a simple χ2 test to the separate locations (not corrected for Month or Year), it turned out that the significant positive association was found in only two specific locations, Wassenaar and Veldhoven. While the prevalence of R. helvetica was very high in Wassenaar, it was rather low in Veldhoven and why the association was observed in only these specific locations is not clear yet. It could be just due to chance, as the P-value is not that low, especially taking into account that we are not correcting for multiple testing here. A study in northern Germany found no association between B. burgdorferi (s.l.) and Rickettsia spp. infections (Knoll et al., 2021). Elucidating the roles of R. helvetica in tick physiology might help understand its spatial variation and its interactions with other tick-borne pathogens and symbionts.

For B. miyamotoi, we found non-significant positive associations with B. burgdorferi (s.l.) (P = 0.0583) and with N. mikurensis (P = 0.0706). This could be linked to the fact that B. miyamotoi is (partially) transmitted horizontally between ticks and rodents (Hoornstra et al., 2022), similar to the association between N. mikurensis and B. burgdorferi (s.l.). The low infection grades for B. miyamotoi and the resulting low numbers of observations of co-infections with B. miyamotoi, may explain the lack of significant results for this symbiont.

For N. mikurensis and M. mitochondrii, a non-significant negative association was observed (P = 0.0921). This is somewhat surprising, since both are positively associated with B. burgdorferi (s.l.), and one would expect the same for the combination of the two. Also, another study, also in the Netherlands, but on different sites than in this study, found a positive association between M. mitochondrii and N. mikurensis (Krawczyk et al., 2022b). The other combinations of symbionts did not show any significant associations or even trends towards associations.

In short, we observed significant spatial variation between the sites in this study for all symbionts except B. miyamotoi and seasonal variation for all symbionts except B. miyamotoi and R. helvetica. Given the differences between the sites in terms of vegetation, (unmeasured) host composition and climatic factors, this may not be surprising. We were unable to attribute the variation to short-term climatic factors, indicating that a more intricate web of interacting factors may determine the prevalences at the local scale. These interaction factors may include the presence of other symbionts, as suggested by the observed associations. Since (long-term) climatic factors also partly determine vegetation and other biotic factors, location and climatic factors will always be correlated and disentangling the effect of (short-term) climatic factors and location proves to be difficult, perhaps even more so since the life span of ticks and the time between infection and onward transmission can be so long, that short-term effects cannot be observed. When we zoom out and look at the prevalences at a larger scale (averaged over time and the different locations), a picture emerges where each symbiont has a certain range within which the prevalences lie: very high for M. mitochondrii, very low for B. miyamotoi, and in between for the other symbionts, suggesting that symbiont-specific factors like transmission routes, host species, duration of bacteraemia, etc., play an important role in determining the transmission dynamics. For the risk for humans, it is also important to realize that the density of infected nymphs is determined by both the nymphal infection prevalence and the density of nymphs, with the latter being the more variable in time and space. The variation in risk is thus to a large extent determined by the variation in abundance of nymphs.

5 Conclusions

Our findings indicate strong between-location variation, some seasonal patterns and hardly any differences between years for most symbionts. Positive associations between symbionts were observed, suggesting that infection with a (vertically transmitted) symbiont may influence the probability of infection with other symbionts, or that there is a common underlying mechanism (e.g. feeding on rodents). Studying the effect of climatic conditions on the spatial and temporal variation of ticks and tick-borne pathogens require longer longitudinal studies, complemented with laboratory measurements to identify underlying mechanisms on tick survival in the presence of various combinations of tick symbionts.

Funding

This study was financially supported by the 10.13039/501100002999 Dutch Ministry of Health, Welfare, and Sports (VWS) . The funder had no role in the study design and interpretation, or the decision to submit the work for publication.

Ethical approval

Not applicable.

CRediT authorship contribution statement

Nienke Hartemink: Methodology, Formal analysis, Investigation, Data curation, Writing – original draft, Writing – review & editing, Visualization. Gerrit Gort: Methodology, Investigation, Writing – review & editing. Aleksandra I. Krawczyk: Formal analysis, Investigation, Data curation, Writing – review & editing. Manoj Fonville: Methodology, Formal analysis, Writing – review & editing. Arnold J.H. van Vliet: Conceptualization, Investigation, Resources, Writing – review & editing, Project administration. Willem Takken: Conceptualization, Investigation, Resources, Writing – review & editing, Supervision, Project administration. Hein Sprong: Conceptualization, Investigation, Writing – review & editing, Supervision, Project administration, Funding acquisition.

Declaration of competing interests

The authors declare that they have no known competing financial interests or personal relationships that could have appeared to influence the work reported in this paper.

Appendix A Supplementary data

The following is the Supplementary data to this article:

Supplementary file 1. Supplementary Table S1. Effect of removing the data from the winter months on the numbers of tested nymphs, and on the number of nymphs tested positive for each of the symbionts. Supplementary Table S2. Comparison of the fit for models with and without climatic factors. The models with the lowest AIC value are marked in green. Supplementary Table S3. Primers and probe used for quantitative multiplex real-time PCR assay. Supplementary Fig. S1: the estimated values for DON for each location (panel A), month (panel B) and year (panel C). Supplementary Fig. S2: the estimated values for NIP (panel A, B and C) and DIN (panel D, E and F) for each location, month and year for Borrelia burgdorferi (s.l.). Supplementary Fig. S3: the estimated values for NIP (panel A, B and C) and DIN (panel D, E and F) for each location, month and year for B. miyamotoi. Supplementary Fig. S4: the estimated values for NIP (panel A, B and C) and DIN (panel D, E and F) for each location, month and year for M. mitochondrii. Supplementary Fig. S5: the estimated values for NIP (panel A, B and C) and DIN (panel D, E and F) for each location, month and year for N. mikurensis. Supplementary Fig. S6: the estimated values for NIP (panel A, B and C) and DIN (panel D, E and F) for each location, month and year for R. helvetica.Multimedia component 1

Multimedia component 1

Data availability

The data supporting the conclusions of this article are included within the article and its supplementary file.

Acknowledgements

The authors are highly appreciative of the monthly tick collections done by the group of volunteers in the Nature Calendar project. We thank the managers of several nature reserves for permission to collect ticks in often closed areas, which has greatly facilitated the work presented here. We would like to thank the reviewers for their thoughtful comments and highly appreciate their efforts toward improving our manuscript.

Appendix A Supplementary data to this article can be found online at https://doi.org/10.1016/j.crpvbd.2024.100209.
==== Refs
References

Andersson M. Bartkova S. Lindestad O. Råberg L. Co-infection with “Candidatus Neoehrlichia mikurensis” and Borrelia afzelii in Ixodes ricinus ticks in southern Sweden Vector Borne Zoonot. Dis. 13 2013 438 442 10.1089/vbz.2012.1118
Andersson M. Scherman K. Råberg L. Infection dynamics of the tick-borne pathogen “Candidatus Neoehrlichia mikurensis” and co-infections with Borreliaafzelii in bank voles in southern Sweden Appl. Environ. Microbiol. 80 2014 1645 1649 10.1128/AEM.03469-13 24375128
Azagi T. Hoornstra D. Kremer K. Hovius J.W.R. Sprong H. Evaluation of disease causality of rare Ixodes ricinus-borne infections in Europe Pathogens 9 2020 150 10.3390/pathogens9020150 32102367
Bakker J.W. Begemann H.L.M. Fonville M. Esser H.J. de Boer W.F. Sprong H. Koenraadt C.J.M. Differential associations of horizontally and vertically transmitted symbionts on Ixodes ricinus behaviour and physiology Parasites Vectors 16 2023 443 10.1186/s13071-023-06025-3 38017525
Brooks M.E. Kristensen K. van Benthem K.J. Magnusson A. Berg C.W. Nielsen A. glmmTMB balances speed and flexibility among packages for zero-inflated generalized linear mixed modeling R J 9 2017 378 10.32614/RJ-2017-066
Burri C. Cadenas F.M. Douet V. Moret J. Gern L. Ixodes ricinus density and infection prevalence of Borrelia burgdorferi sensu lato along a north-facing altitudinal gradient in the Rhône Valley (Switzerland) Vector Borne Zoonotic Dis. 7 2007 50 58 10.1089/vbz.2006.0569 17417957
Cochez C. Heyman P. Heylen D. Fonville M. Hengeveld P. Takken W. The presence of Borrelia miyamotoi, a relapsing fever spirochaete, in questing Ixodes ricinus in Belgium and in the Netherlands Zoonoses Public Health 62 2015 331 333 10.1111/zph.12154 25212814
Coipan E.C. Jahfari S. Fonville M. Maassen C. van der Giessen J. Takken W. Spatiotemporal dynamics of emerging pathogens in questing Ixodes ricinus Front. Cell. Infect. Microbiol. 3 2013 36 10.3389/fcimb.2013.00036 23908971
Corbin C. Heyworth E.R. Ferrari J. Hurst G.D.D. Heritable symbionts in a world of varying temperature Heredity 118 2017 10 20 10.1038/hdy.2016.71 27703153
Crowder C.D. Carolan H.E. Rounds M.A. Honig V. Mothes B. Haag H. Prevalence of Borrelia miyamotoi in Ixodes ticks in Europe and the United States Emerg. Infect. Dis. 20 2014 1678 1682 10.3201/eid2010.131583 25280366
de Bruin A. van Leeuwen A.D. Jahfari S. Takken W. Földvári M. Dremmel L. Vertical transmission of Bartonellaschoenbuchensis in Lipoptena cervi Parasites Vectors 8 2015 176 10.1186/s13071-015-0764-y 25889985
Derdáková M. Václav R. Pangrácova-Blaňárová L. Selyemová D. Koči J. Walder G. Špitalská E. Candidatus Neoehrlichia mikurensis and its co-circulation with Anaplasma phagocytophilum in Ixodes ricinus ticks across ecologically different habitats of Central Europe Parasites Vectors 7 2014 160 10.1186/1756-3305-7-160 24693971
Diuk-Wasser M.A. Hoen A.G. Cislo P. Brinkerhoff R. Hamer S.A. Rowland M. Human risk of infection with Borrelia burgdorferi, the Lyme disease agent, in eastern United States Am. J. Trop. Med. Hyg. 86 2012 320 327 10.4269/ajtmh.2012.11-0395 22302869
Dobler P.D.G. Erber D.W. Bröker D.M. Schmitt P.D.H.-J. The TBE Book 3rd ed. 2020 Global Health Press Pte Ltd
Eisen R.J. Eisen L. Girard Y.A. Fedorova N. Mun J. Slikas B. A spatially-explicit model of acarological risk of exposure to Borrelia burgdorferi-infected Ixodes pacificus nymphs in northwestern California based on woodland type, temperature, and water vapor Ticks Tick Borne Dis. 1 2010 35 43 10.1016/j.ttbdis.2009.12.002 20532183
Esser H.J. Lim S.M. de Vries A. Sprong H. Dekker D.J. Pascoe E.L. Continued circulation of tick-borne encephalitis virus variants and detection of novel transmission foci, The Netherlands Emerg. Infect. Dis. 28 2022 2416 2424 10.3201/eid2812.220552 36288572
Fingerle V. Bergmeister H. Liegl G. Vanek E. Wilske B. Prevalence of Borrelia burgdorferi sensu lato in Ixodes ricinus in southern Germany J. Spiroch. Tick Dis. 1 1994 41 45
Garcia-Vozmediano A. Tomassone L. Fonville M. Bertolotti L. Heylen D. Fabri N.D. The genetic diversity of Rickettsiella symbionts in Ixodes ricinus throughout Europe Microb. Ecol. 84 2022 613 626 10.1007/s00248-021-01869-7 34580739
Gassner F. Van Vliet A.J.H. Burgers S.L.G.E. Jacobs F. Verbaarschot P. Hovius E.K.E. Geographic and temporal variations in population dynamics of Ixodes ricinus and associated Borrelia infections in The Netherlands Vector Borne Zoonotic Dis. 11 2011 523 532 10.1089/vbz.2010.0026 21083369
Gilbert L. Maffey G.L. Ramsay S.L. Hester A.J. The effect of deer management on the abundance of Ixodes ricinus in Scotland Ecol. Appl. 22 2012 658 667 10.1890/11-0458.1 22611862
Guizzo M.G. Hatalová T. Frantová H. Zurek L. Kopáček P. Perner J. Ixodes ricinus ticks have a functional association with Midichloria mitochondrii Front. Cell. Infect. Microbiol. 12 2023 1081666 https://doi:10.3389/fcimb.2022.1081666
Hansford K.M. Fonville M. Gillingham E.L. Coipan E.C. Pietzsch M.E. Krawczyk A.I. Ticks and Borrelia in urban and peri-urban green space habitats in a city in southern England Ticks Tick Borne Dis. 8 2017 353 361 10.1016/j.ttbdis.2016.12.009 28089123
Hartemink N. Van Vliet A. Sprong H. Jacobs F. Garcia-Martí I. Zurita-Milla R. Takken W. Temporal-spatial variation in questing tick activity in The Netherlands: The effect of climatic and habitat factors Vector Borne Zoonotic Dis. 19 2019 494 505 10.1089/vbz.2018.2369 30810501
Hartemink N. van Vliet A.J.H. Gort G. Gassner F. Jacobs F. Fonville M. Seasonal patterns and spatial variation of Borrelia burgdorferi (sensu lato) infections in Ixodes ricinus in The Netherlands Parasites Vectors 14 2021 121 10.1186/s13071-021-04607-7 33627166
Herrmann C. Gern L. Search for blood or water is influenced by Borrelia burgdorferi in Ixodes ricinus Parasites Vectors 8 2015 6 10.1186/s13071-014-0526-2 25560984
Heylen D. Adriaensen F. Van Dongen S. Sprong H. Matthysen E. Ecological factors that determine Ixodes ricinus tick burdens in the great tit (Parus major), an avian reservoir of Borrelia burgdorferi s.l Int. J. Parasitol. 43 2013 603 611 10.1016/j.ijpara.2013.02.007 23597868
Hildebrandt A. Krämer A. Sachse S. Straube E. Detection of Rickettsia spp. and Anaplasma phagocytophilum in Ixodes ricinus ticks in a region of Middle Germany (Thuringia) Ticks Tick Borne Dis. 1 2010 52 56 10.1016/j.ttbdis.2009.11.005 21771511
Hofhuis A. Harms M. van den Wijngaard C. Sprong H. van Pelt W. Continuing increase of tick bites and Lyme disease between 1994 and 2009 Ticks Tick Borne Dis. 6 2015 69 74 10.1016/j.ttbdis.2014.09.006 25448421
Hofhuis A. Kassteele J. van de Sprong H. Wijngaard C.C. van den Harms M.G. Fonville M. Predicting the risk of Lyme borreliosis after a tick bite, using a structural equation model PLoS One 12 2017 e0181807 10.1371/journal.pone.0181807
Hofmeester T.R. Sprong H. Jansen P.A. Prins H.H.T. van Wieren S.E. Deer presence rather than abundance determines the population density of the sheep tick, Ixodes ricinus, in Dutch forests Parasites Vectors 10 2017 433 10.1186/s13071-017-2370-7 28927432
Hoornstra D. Azagi T. van Eck J.A. Wagemakers A. Koetsveld J. Spijker R. Prevalence and clinical manifestation of Borrelia miyamotoi in Ixodes ticks and humans in the northern hemisphere: A systematic review and meta-analysis Lancet Microbe 3 2022 e772 e786 10.1016/S2666-5247(22)00157-4 36113496
Hovius J.W.R. De Wever B. Sohne M. Brouwer M.C. Coumou J. Wagemakers A. A case of meningoencephalitis by the relapsing fever spirochaete Borrelia miyamotoi in Europe Lancet 382 2013 658 10.1016/S0140-6736(13)61644-X 23953389
Jahfari S. Fonville M. Hengeveld P. Reusken C. Scholte E.J. Takken W. Prevalence of Neoehrlichia mikurensis in ticks and rodents from North-West Europe Parasites Vectors 5 2012 74 10.1186/1756-3305-5-74 22515314
Jahfari S. Hofhuis A. Fonville M. van der Giessen J. van Pelt W. Sprong H. Molecular detection of tick-borne pathogens in humans with tick bites and erythema migrans, in The Netherlands PLoS Negl. Trop. Dis. 10 2016 e0005042 10.1371/journal.pntd.0005042
Kantsø B. Bo Svendsen C. Moestrup Jensen P. Vennestrøm J. Krogfelt K.A. Seasonal and habitat variation in the prevalence of Rickettsia helvetica in Ixodes ricinus ticks from Denmark Ticks Tick Borne Dis. 1 2010 101 103 10.1016/j.ttbdis.2010.01.004 21771515
Kjær L.J. Klitgaard K. Soleng A. Edgar K.S. Lindstedt H.E.H. Paulsen K.M. Spatial patterns of pathogen prevalence in questing Ixodes ricinus nymphs in southern Scandinavia, 2016 Sci. Rep. 10 2020 19376 10.1038/s41598-020-76334-5
Knoll S. Springer A. Hauck D. Schunack B. Pachnicke S. Fingerle V. Strube C. Distribution of Borrelia burgdorferi s.l. and Borrelia miyamotoi in Ixodes tick populations in northern Germany, co-infections with Rickettsiales and assessment of potential influencing factors Med. Vet. Entomol. 35 2021 595 606 10.1111/mve.12537 34180074
Köhler C.F. Holding M.L. Sprong H. Jansen P.A. Esser H.J. Biodiversity in the Lyme-light: Ecological restoration and tick-borne diseases in Europe Trends Parasitol. 39 2023 373 385 10.1016/j.pt.2023.02.005 36890021
Krawczyk A. Röttjers L. Fonville M. Takumi K. Takken W. Faust K. Sprong H. Quantitative microbial population study reveals geographical differences in bacterial symbionts of Ixodes ricinus Microbiome 10 2022 120 10.1186/s40168-022-01276-1 35927748
Krawczyk A.I. Röttjers S. Coimbra-Dores M.J. Heylen D. Fonville M. Takken W. Tick microbial associations at the crossroad of horizontal and vertical transmission pathways Parasites Vectors 15 2022 380 10.1186/s13071-022-05519-w 36271430
Lambin E.F. Tran A. Vanwambeke S.O. Linard C. Soti V. Pathogenic landscapes: Interactions between land, people, disease vectors, and their animal hosts Int. J. Health Geogr. 9 2010 54 20979609
Lenth R. emmeans: Estimated marginal means, aka least-squares means R package version 1.4.6 https://CRAN.R-project.org/package=emmeans 2020
Lynn G.E. Breuner N.E. Eisen L. Hojgaard A. Replogle A.J. Eisen R.J. An immunocompromised mouse model to infect Ixodes scapularis ticks with the relapsing fever spirochete, Borrelia miyamotoi Ticks Tick Borne Dis. 10 2019 352 359 10.1016/j.ttbdis.2018.11.017 30503357
Mather T.N. Nicholson M.C. Donnelly E.F. Matyas B.T. Entomologic index for human risk of Lyme disease Am. J. Epidemiol. 144 1996 1066 1069 10.1093/oxfordjournals.aje.a008879 8942438
Medlock J.M. Hansford K.M. Bormane A. Derdakova M. Estrada-Peña A. George J.-C. Driving forces for changes in geographical distribution of Ixodes ricinus ticks in Europe Parasites Vectors 6 2013 1 10.1186/1756-3305-6-1 23281838
Olsthoorn F. Sprong H. Fonville M. Rocchi M. Medlock J. Gilbert L. Ghazoul J. Occurrence of tick-borne pathogens in questing Ixodes ricinus ticks from Wester Ross, Northwest Scotland Parasites Vectors 14 2021 430 10.1186/s13071-021-04946-5 34446082
Pepin K.M. Eisen R.J. Mead P.S. Piesman J. Fish D. Hoen A.G. Geographic variation in the relationship between human Lyme disease incidence and density of infected host-seeking Ixodes scapularis nymphs in the eastern United States Am. J. Trop. Med. Hyg. 86 2012 1062 1071 10.4269/ajtmh.2012.11-0630 22665620
R Core Team R: A language and environment for statistical computing 2023 R Foundation for Statistical Computing Vienna https://www.R-project.org/
Randolph S. Storey K. Impact of microclimate on immature tick-rodent host interactions (Acari: Ixodidae): Implications for parasite transmission J. Med. Entomol. 36 1999 741 748 10.1093/jmedent/36.6.741 10593075
Rollend L. Fish D. Childs J.E. Transovarial transmission of Borrelia spirochetes by Ixodes scapularis: A summary of the literature and recent observations Ticks Tick Borne Dis. 4 2013 46 51 10.1016/j.ttbdis.2012.06.008 23238242
SAS Institute Inc SAS/STAT® 14.3 User's Guide: High-Performance Procedures 2017 SAS Institute Inc Cary, NC, USA https://support.sas.com/documentation/onlinedoc/stat/143/stathpug.pdf
Sassera D. Beninati T. Bandi C. Bouman E.A.P. Sacchi L. Fabbi M. Lo N. ‘Candidatus Midichloria mitochondrii’, an endosymbiont of the tick Ixodes ricinus with a unique intramitochondrial lifestyle Int. J. Syst. Evol. Microbiol. 56 2006 2535 2540 10.1099/ijs.0.64386-0 17082386
Silaghi C. Beck R. Oteo J.A. Pfeffer M. Sprong H. Neoehrlichiosis: An emerging tick-borne zoonosis caused by Candidatus Neoehrlichia mikurensis Exp. Appl. Acarol. 68 2016 279 297 10.1007/s10493-015-9935-y 26081117
Sprong H. Azagi T. Hoornstra D. Nijhof A.M. Knorr S. Baarsma M.E. Hovius J.W. Control of Lyme borreliosis and other Ixodes ricinus-borne diseases Parasites Vectors 11 2018 145 10.1186/s13071-018-2744-5 29510749
Sprong H. Wielinga P.R. Fonville M. Reusken C. Brandenburg A.H. Borgsteede F. Ixodes ricinus ticks are reservoir hosts for Rickettsia helvetica and potentially carry flea-borne Rickettsia species Parasites Vectors 2 2009 41 10.1186/1756-3305-2-41 19732416
Stark J.H. Pilz A. Jodar L. Moïsi J.C. The epidemiology of Lyme borreliosis in Europe: An updated review on a growing public health issue Vector Borne Zoonotic Dis. 23 2023 139 141 10.1089/vbz.2022.0068 37071398
Takken W. van Vliet A.J.H. Verhulst N.O. Jacobs F.H.H. Gassner F. Hartemink N. Acarological risk of Borrelia burgdorferi sensu lato infections across space and time in The Netherlands Vector Borne Zoonotic Dis. 17 2017 99 107 10.1089/vbz.2015.1933 27893309
Takken W. van Vliet A.J.H. Verhulst N.O. Jacobs F.H.H. Gassner F. Hartemink N. Acarological risk of Borrelia burgdorferi sensu lato infections across space and time in The Netherlands Vector Borne Zoonotic Dis. 17 2017 99 107 10.1089/vbz.2015.1933 27893309
Takumi K. Sprong H. Hofmeester T.R. Impact of vertebrate communities on Ixodes ricinus-borne disease risk in forest areas Parasites Vectors 12 2019 434 10.1186/s13071-019-3700-8 31492171
Tijsse-Klasen E. Fonville M. Reimerink J.H. Spitzen-van der Sluijs A. Sprong H. Role of sand lizards in the ecology of Lyme and other tick-borne diseases in The Netherlands Parasites Vectors 3 2010 42 10.1186/1756-3305-3-42 20470386
van Dobbenburgh A. van Dam A.P. Fikrig E. Human granulocytic ehrlichiosis in western Europe N. Engl. J. Med. 340 1999 1214 1216 10.1056/NEJM199904153401517
Vautrin E. Vavre F. Interactions between vertically transmitted symbionts: Cooperation or conflict? Trends Microbiol. 17 2009 95 99 10.1016/j.tim.2008.12.002 19230673
Wielinga P.R. Gaasenbeek C. Fonville M. De Boer A. De Vries A. Dimmers W. Longitudinal analysis of tick densities and Borrelia, Anaplasma, and Ehrlichia infections of Ixodes ricinus ticks in different habitat areas in The Netherlands Appl. Environ. Microbiol. 72 2006 7594 7601 10.1128/AEM.01851-06 17028227
Zeimes C.B. Olsson G.E. Hjertqvist M. Vanwambeke S.O. Shaping zoonosis risk: Landscape ecology vs. landscape attractiveness for people, the case of tick-borne encephalitis in Sweden Parasites Vectors 7 2014 370 10.1186/1756-3305-7-370 25128197
