Temporal dynamics of the tick Ixodes ricinus in northern Europe : epidemiological implications

This is an electronic reprint of the original article. This reprint may differ from the original in pagination and typographic detail. Author(s): Title: Year: Version: Please cite the original version: All material supplied via JYX is protected by copyright and other intellectual property rights, and duplication or sale of all or part of any of the repository collections is not permitted, except that material may be duplicated by you for your research use or educational purposes in electronic or print form. You must obtain permission for any other use. Electronic or print copies may not be offered, whether for sale or otherwise to anyone who is not an authorised user. Temporal dynamics of the tick Ixodes ricinus in northern Europe : epidemiological implications Cayol, Claire; Koskela, Esa; Mappes, Tapio; Siukkola, Anja; Kallio, Eva Cayol, C., Koskela, E., Mappes, T., Siukkola, A., & Kallio, E. (2017). Temporal dynamics of the tick Ixodes ricinus in northern Europe : epidemiological implications. Parasites and Vectors, 10(1), Article 166. https://doi.org/10.1186/s13071-017-2112-x 2017 RESEARCH Open Access Temporal dynamics of the tick Ixodes ricinus in northern Europe: epidemiological implications Claire Cayol * , Esa Koskela, Tapio Mappes, Anja Siukkola and Eva R. Kallio Abstract Background: Tick-borne pathogens pose an increasing threat to human and veterinary health across the northern hemisphere. While the seasonal activity of ticks is largely determined by climatic conditions, host-population dynamics are also likely to affect tick abundance. Consequently, abundance fluctuations of rodents in northern Europe are expected to be translated into tick dynamics, and can hence potentially affect the circulation of tick-borne pathogens. We quantified and explained the temporal dynamics of the tick Ixodes ricinus in the northernmost part of its European geographical range, by estimating (i) abundance in vegetation and (ii) infestation load in the most common rodent species in the study area, the bank vole Myodes glareolus. Results: Ixodes ricinus nymphs and adult females, the life stages responsible for the most of tick bites in humans, peaked in May-June and August-September. Larvae and nymphs were simultaneously active in June and abundance of questing larvae and nymphs in the vegetation showed a positive association with bank vole abundance. Moreover, infesting larvae and nymphs were aggregated on bank voles, and the infestation of bank voles with I. ricinus larvae and nymphs was positively associated with bank vole abundance. Conclusion: Our results indicate early summer and early autumn as periods of increased risk for humans to encounter I. ricinus ticks in boreal urban forests and suggest a 2 years life-cycle for I. ricinus with two cohorts of ticks during the same year. Moreover, we identified a simultaneous activity of larvae and nymphs which allows co-feeding on the rodent host, which in turn supports the transmission of several important zoonotic tick-borne pathogens. Finally, we showed that a high density of the rodent host may enhance the risk that ticks and, potentially, tick-borne pathogens pose to human health. Keywords: Ixodes ricinus, Rodent host, Seasonality, Public health, Population dynamics Background Tick-borne pathogens are a growing burden for European public health policies [1–3]. The current observed increase in tick-borne disease incidence in Europe may be explained by the geographical expansion of Ixodes ricinus, the growing share of space between humans and wild animals, and the improvement of diagnostics tools [4–6]. The epidemiology of tick-borne zoonoses, such as Lyme borreliosis, anaplasmosis or tick-borne encephalitis (TBE), depends on tick abundance and population dynamics, infection prevalence within the tick population, and land use that may affect human exposure to ticks [7, 8]. In order to predict the risks that tick-borne diseases pose to humans, an assessment of factors underlying the temporal variation of tick abundance is necessary. The abundance of I. ricinus varies in time and space and is highly dependent on environmental conditions, including habitat quality, host availability, and abiotic conditions [9–12]. In northern Fennoscandia, at the northernmost part of the European range of I. ricinus, abiotic conditions undergo extreme seasonal variation; there are 145 to 160 days of snow cover with short daylengths, during which ticks are not active. This is followed by a quick elevation in temperature leading to a short summer with long day-lengths [13]. In these * Correspondence: [email protected] Department of Biological and Environmental Science, University of Jyväskylä, P.O. Box 35, FI-40014 Jyväskylä, Finland © The Author(s). 2017 Open Access This article is distributed under the terms of the Creative Commons Attribution 4.0 International License (http://creativecommons.org/licenses/by/4.0/), which permits unrestricted use, distribution, and reproduction in any medium, provided you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons license, and indicate if changes were made. The Creative Commons Public Domain Dedication waiver (http://creativecommons.org/publicdomain/zero/1.0/) applies to the data made available in this article, unless otherwise stated. Cayol et al. Parasites & Vectors (2017) 10:166 DOI 10.1186/s13071-017-2112-x conditions, I. ricinus activity is likely to show distinctive seasonal patterns, which have not been characterized to date (but see [14] for southwest Finland). Ixodes ricinus is dependent on vertebrate hosts to complete its life-cycle. Larvae typically feed on small vertebrates, such as rodents; nymphs, the more common biting stage for humans, parasitize mostly medium-sized mammals; and adults feed mainly on large hosts, such as deer [15, 16]. The population dynamics of ticks and rodents are expected to be linked: some studies have indicated delayed density dependence of questing nymphs on rodent abundance, suggesting that high rodent abundance provides augmented opportunities for successful larvae feeding and nymph development [17, 18]. The bank vole (Myodes glareolus) is a common rodent species throughout Europe [19]; this species is commonly infested by immature I. ricinus [16, 20, 21]. In northern Europe, vole population abundance shows both seasonality, driven by seasonal breeding, and multiannual density fluctuations shaped by predation, food availability and food quality [22–24]. These seasonal and multiannual density fluctuations are likely to be translated into the dynamics of ticks, and consequently, into the epidemiology of tick-borne pathogens. To date, there are few studies that have investigated the association between the dynamics of cyclic small rodents and ticks [25]. The bank vole is also an important reservoir host for many tick-borne pathogens, such as Borrelia afzelii, tickborne encephalitis virus (TBEV) and Babesia microti [26, 27]. Typically, tick larvae acquire infections from an infected rodent host that has become infected while feeding infected nymph(s) [10]. Alternatively, larvae acquire infections via simultaneous feeding with infected nymphs without systemic infection of the host [28, 29]. Infectivity is transstadially maintained in the tick to the following life stage [30]. Here, we present results from a 4 years of longitudinal bank vole monitoring and tick sampling in central Finland at the northernmost part of the European range of I. ricinus, where abiotic conditions undergo extreme seasonal variation. Our primary aim is to characterize temporal dynamics and quantify the importance of host related factors and abiotic conditions on temporal dynamics of I. ricinus. We also aim to identify seasonal patterns that are relevant for tickborne pathogen circulation, with the ultimate goal of providing information concerning the risk of tickborne diseases in our study area. Methods Study area Sampling took place monthly from May to October in 2012 −2015 in four periurban forests in the Jyväskylä area in Central Finland: (Kylmänoro (62°13′36.220″, 25° 45′1.739″); Jyskänlaakso (62°13′55.398″, 25°49′ 34.269″); Hämeenlahti (62°12′40.119″, 25°47′11.052″); and Sippulanniemi (62°11′9.019″, 25°44′58.147″) [31]. One trapping period within a month will be referred to as “session”in the following paragraphs. Forests were dominated by Scots pine (Pinus sylvestris) and silver birch (Betula pendula) or by spruce (Picea abies). The herbaceous stratum was typically composed of Vaccinium myrtillus,V. vitis-idaea,Maianthemum bifolium,Linnaea borealis and Oxalis acetosella. Tick dragging Monthly tick dragging was performed during or within a few days of the vole trapping, using a 1 × 1 m cotton flannel flag sewed to a wooden rod [12]. The fabric was randomly dragged over the vegetation for 300 −500 m per site around the rodent trapping transects and checked every 20−25 m for ticks, which were removed with tweezers and stored in alcohol at -20 °C. No dragging was performed during rain. In October 2014, due to early snow cover, dragging was not performed. Due to the duration, coverage and interval of the dragging (less than 30 min, 300–500 m 2 once a month in each site) it is unlikely that the flag dragging affected the overall tick population abundance and it should not have interfered with the ticks parasitizing rodents in the area. Vole trapping and tick infestation on voles As the active tick population consists in parasitizing, questing and resting ticks, sampling targeted questing ticks and parasitizing ticks on their rodent host. This latter buffers the effects of microclimate changes and rodent sampling, in particular, also buffers the effect of the patchy distribution of larvae [32]. Vole trapping was carried out with two lines of 10 Ugglan Special multiple-capture live traps (Grahnab Company, Sweden), positioned 10 −15 m apart, located near to rodent burrows. Traps were prebaited for 1−3 nights with sunflower seeds (Helianthus annuus), after which traps were set with sunflower seeds (for food) and a piece of potato (for water) for two consecutive nights. Wood shavings were provided as bedding in wet or cold weather. Traps were checked once per day and trapped voles were handled and sampled before release close to their capture site. Bycatch of species other than voles, as well as recapture of the same individual during the same session, were released immediately on site. All trapped voles were marked individually with electronic identification chips (microchip Trovan Unique™), which were injected subcutaneously at their first capture. During each capture, voles were identified, body mass was measured as a proxy for age (as in e.g. [33]), and sex and reproductive condition were recorded. The presence Cayol et al. Parasites & Vectors (2017) 10:166 Page 2 of 11 of fleas was recorded and all voles were examined for ticks, with special attention to the area around the ears and face. All ticks were removed with tweezers and stored in alcohol at -20 °C until further identification. All ticks - both those removed from rodents and those collected from vegetation - were identified to species level and life stage under a dissection microscope using morphological identification keys [34–36]. Species identification of seven ticks identified as I. ricinus and three as I. trianguliceps was further confirmed with PCR following a method described elsewhere [37]. Briefly, PCR targeted the mitochondrial 16S rRNA gene and the amplicons obtained were successfully sequenced for eight of the ten ticks. Thereafter sequence identity was determined by BLAST search against the NCBI Nucleotide database and the obtained sequences confirmed our morphological tick identification. We assessed the overall bank vole population abundance by computing the overall minimum number of voles alive (MNA) at a given trapping session (t) as follows: total number of individuals caught at a given trapping session (t) summed with the total number of individuals marked when caught during subsequent sessions, but not caught at (t) [38]. We trapped 658 bank voles, an average of 1.53 times (range 1–6), for a total of 1007observations for which all variables described above were available. The minimum number of voles alive per session varied from 5 (in May 2013) to 120 individuals (in September 2014). Three other rodent species were bycaught, consisting of 52 observations of yellow-necked mouse (Apodemus flavicollis), one observation of field vole (Microtus agrestis), and two observations of house mouse (Mus musculus) (Additional file 1: Figure S3). Statistical analysis Ticks in vegetation We characterized the temporal activity of I. ricinus in the vegetation (i.e. collected by flagging), by examining tick questing activity separately for each life stage, i.e. larvae, nymphs, adults (males and females), in relation to the following variables: year (2012 −2015), month (May-October), estimated bank vole abundance per given session (MNA), abundance of other life stages present during the same session (number of ticks/ 100 m 2 ), and the abundance of previous tick life stages collected during the previous session (for larvae: adult, for nymph: larvae, for adults: nymphs). To further identify the effect of current climatic conditions on tick activity, we computed the mean daily saturation deficit (SatDef, in millimetres of mercury) during tick flagging days, based on daily average humidity (in percent) and daily average temperature (in °C) [9, 39, 40] recorded at the meteorological station of Nenäinniemi in Jyväskylä, located 0.72−3.7 km from the study sites (http:// www.jyv-weather.info/index.php) (Additional file 1: Figure S1). SatDef was used as an explanatory variable rather than month, with which it showed collinearity. Thus, the second set of models included SatDef and its second-degree polynomial term SatDef 2 , MNA, current and previous tick abundances as described above. Furthermore, the abundance of nymphs and females pooled together was also modelled with two sets of models: the first one included vole abundance, month and year and the second one included year, vole abundance, SatDef and SatDef 2 . Models were fitted using generalized linear mixed models (GLMM) with a negative binomial error distribution (with log-link function) and site was included as a random effect to control for potential pseudoreplication [41]. To take into account the variation in the distance flags were dragged, an offset term (log(distance flagged)/100) was introduced in the models. The model selection (provided in Additional file 1: Table S2) was an automated selection process starting from the full model and based on AICc (Akaike Information Criteria corrected for small sample size [42]), using dredge function in R software. We kept the most parsimonious model that lay within 2AICc difference from the best model fitted [42] (Additional file 1: Tables S1 and S2). Ticks infesting voles Ixodes ricinus infestation load on bank voles was examined separately for larvae and nymphs. We assessed whether tick infestation showed seasonality and/or between year variation and whether it was affected by individual host characteristics or by concomitant parasitism (by other tick stages, other tick species or fleas). For that purpose, we fitted a GLMM with a negative binomial error distribution to test the fixed effects of month, year, vole sex, body mass (centred value) and its second order polynomial term, presence of fleas, presence of other life stages of I. ricinus and I. trianguliceps, body mass*vole abundance (MNA) interaction term and body mass*sex interaction term. ‘Trapping site’and ‘vole individual nested in the trapping site’were included as random effects in the models. Model selection was performed as described before except that we utilized the function drop1 in R software (Additional file 1: Tables S4 and S5). All statistical analyses were performed with R version 3.2.3 (2015, The R Foundation for Statistical Computing), and using the packages stats (http://www.R-project.org/), MASS (https://cran.r-project.org/web/packages/MASS/ index.html), glmmADMB (http://glmmadmb.r-forge.rproject.org/) and MuMIn (https://cran.r-project.org/ web/packages/MuMIn/index.html). Cayol et al. Parasites & Vectors (2017) 10:166 Page 3 of 11 Results Ticks in vegetation We sampled and identified 943 I. ricinus larvae, 867 nymphs, 239 (adult) females and 294 males from the vegetation. The mean abundance of I. ricinus per session and per area varied from 0 to 22.7 ticks/100 m 2 when considering all tick life stages and from 0 to 6.25 ticks/ 100 m 2 when taking into account only female adults and nymphs. Overall, the density of questing ticks collected from vegetation was 7.1/100 m 2 (Additional file 1: Figure S2). The ratio between I. ricinus larvae, nymphs and adults was 3.5:3.3:2.0. In addition, one I. trianguliceps nymph was identified. Models revealed unimodal questing patterns for larvae, which were mostly found in June. Conversely, a bimodal questing pattern was found for nymphs as well as nymph and females modelled together, with the highest abundances found in May-June and September, which therefore appears as the higher risk period for tick bites on humans. Questing adults (males and females) were more abundant in May-June and August-September, and their abundance varied with year, with the highest abundance found in 2015 (Table 1, Fig. 1). For any given month, the abundances of questing I. ricinus larvae and nymphs showed positive associations with vole abundance (Tables 1, 2; Fig. 2; Additional file 1: Figure S4). For each addition of one individual to the bank vole population, an increase of larvae abundance by approximately 3% and of nymph abundance by 1% was predicted (Fig. 2; Additional file 1: Figure S4). The abundance of questing larvae was positively associated with saturation deficit, while the abundance of questing adults showed a non-linear relationship with saturation deficit. The number of adults found in vegetation was positively associated with saturation deficit until an optimal value (3.16 mm Hg), after which the abundance of adult ticks was negatively affected by any further increase in saturation deficit. Nymph abundance was not associated with saturation deficit in the best model selected (Table 2). We found a positive relationship between the number of questing larvae and the abundance of adults observed in the vegetation one session before. Nymph abundance increased with adult abundance during the same flagging session whereas adult abundance showed a negative relationship with nymph abundance during the previous session (Table 2). Ticks infesting voles From bank voles, two tick species were identified: I. trianguliceps, the vole tick and I. ricinus. The proportion of infestation with either of these tick species was 75.8%. The ratio of I. ricinus larvae to nymphs found feeding on bank voles was 13:1. The total number of ticks Table 1 Selected best model for the abundance of tick questing in the vegetation with estimated coefficients (in log scale), explained by vole abundance, month (May taken as reference) and year (2012 as reference) Y = Larva abundance Estimate (SE) z-value P-value Intercept -1.424 (0.667) -2.13 0.033 June 1.891 (0.742) 2.55 0.011 July 0.003 (0.778) 0.00 0.997 August -1.100 (0.974) -1.13 0.259 September -0.987 (1.035) -0.95 0.340 October -1.707 (0.979) -1.74 0.081 Vole abundance 0.028 (0.013) 2.14 0.032 Random effect: site σ 2 = 0.46 (SD = 0.68) Negative binomial dispersion parameter 0.38 (SE = 0.07) Y = Nymph abundance Estimate (SE) z-value P-value Intercept 0.193 (0.375) 0.52 0.607 June -0.401 (0.275) -1.46 0.145 July -1.571 (0.313) -5.02 <0.005 August -1.312 (0.351) -3.74 <0.005 September -0.628 (0.390) -1.61 0.107 October -2.730 (0.404) -6.76 <0.005 Vole abundance 0.013 (0.004) 2.99 0.003 Random effect: site σ 2 = 0.38 (SD = 0.62) Negative binomial dispersion parameter 2.88 (SE = 0.70) Y = Adult (male + female) abundance Estimate (SE) z-value P-value Intercept -0.766 (0.393) -1.95 0.051 June -0.203 (0.226) -0.90 0.368 July -0.600 (0.238) -2.52 0.012 August 0.395 (0.215) 1.84 0.066 September 0.279 (0.214) 1.30 0.192 October -1.082 (0.292) -3.71 <0.005 2013 0.288 (0.196) 1.47 0.142 2014 0.306 (0.202) 1.51 0.131 2015 0.923 (0.191) 4.82 <0.005 Random effect: site σ 2 = 0.39 (SD = 0.63) Negative binomial dispersion parameter 7.48 (SE = 2.81) Y = Female + Nymph abundance Estimate (SE) z-value P-value Intercept 0.514 (0.350) 1.47 0.142 June -0.422 (0.227) -1.86 0.063 July -1.509 (0.258) -5.84 <0.005 August -1.021 (0.282) -3.62 <0.005 September -0.599 (0.324) -1.85 0.064 October -2.430 (0.323) -7.53 <0.005 Vole abundance 0.012 (0.004) 3.04 0.002 Random effect: site σ 2 = 0.37 (SD = 0.60) Negative binomial dispersion parameter 4.42 (SE = 1.10) σ 2 is the variance attributable to random effect. Number of observations: Total = 88; Site = 4 Abbreviations:SD standard deviation, SE standard error Cayol et al. Parasites & Vectors (2017) 10:166 Page 4 of 11 sampled from voles was 3564, out of which 14 ticks could not be identified due to poor condition (Additional file 1: Table S3 and Figure S5). Models revealed a clear seasonal pattern in the infestation burden of I. ricinus larvae on bank voles (Table 3, Fig. 4): larval infestation underwent seasonality, with a peak in June and a trough in August-October. The highest infestation level was in 2013 and the lowest in 2014. In addition, bank vole infestation load with I. ricinus nymphs underwent seasonal variation, with a peak in May, but was stable between years (Table 4, Fig. 3). For any given month, nymph infestation on voles was positively associated with bank vole abundance (Table 4). Similarly, larval infestation level increased with bank vole abundance, but the increase was more pronounced among female bank voles than among males (Table 3, Fig. 4). Moreover, infestation with larvae was positively associated with the amount of questing larvae observed in the environment (Table 3), whereas the bank vole infestation load with nymphs was not associated with the amount of questing nymph (i.e. the abundance of questing nymphs was not selected in the best model, Additional file 1: Table S5). Tick infestation intensity on a host varied with individual characteristics such as age, sex and co-infestation. The oldest male bank voles (i.e. those with highest body mass) were the most intensely infested with larvae (Table 3). Moreover, bank vole infestation load with I. ricinus larvae was positively associated with co-infesting I. trianguliceps females and nymphs and I. ricinus nymphs (Table 3), whereas the infestation with I. ricinus nymphs increased with the presence of I. trianguliceps larvae and females (Table 4). In addition, the infestation load with nymphs showed a non-linear relationship with body mass: infestation load increased until voles reached 32.4 g, whereupon any further increase in body mass led to a reduction of the infestationburden(Table4). Discussion In this study, we characterized the temporal dynamics of I. ricinus by assessing its abundance in the vegetation and its infestation load in one of its main host in the northernmost part of its range. We focused on bank voles, which parasitism provides insightful information concerning the local immature tick communities. Fig. 1 Predicted number ± standard error (SE) of larvae, nymphs and adults in 100 m 2 of vegetation per month. Predictions are based on GLMM shown in Table 1 Table 2 Selected best model for the abundance of ticks questing in the vegetation with estimated coefficients (in log scale), explained by the vole abundance, the amount of ticks in other stages in vegetation during the previous session and/or during the current session, and the saturation deficit (SatDef) and its second degree polynomial term (SatDef 2 ) Y = Larva abundance Estimate (SE) z-value P-value Intercept -5.426 (1.002) -5.41 <0.005 Vole abundance 0.029 (0.009) 3.11 0.002 Amount of adult ticks during the previous session 1.007 (0.308) 3.27 0.001 SatDef 0.969 (0.192) 5.03 <0.005 Random effect: site σ 2 = 4.59e -06 (SD = 0.002) Negative binomial dispersion parameter 0.34 (SE = 0.06) Y = Nymph abundance Estimate (SE) z-value P-value Intercept -0.279 (0.233) -1.20 0.232 Amount of adult ticks during the same session 0.381 (0.167) 2.28 0.023 Random effect: site σ 2 = 0.098 (SD = 0.31) Negative binomial dispersion parameter 1.12 (SE = 0.21) Y = Adult (male + female) abundance Estimate (SE) z-value P-value Intercept -1.294 (0.461) -2.81 0.005 SatDef 0.621 (0.239) 2.60 0.009 SatDef 2 -0.098 (0.037) -2.63 0.009 Amount of nymph during the same session 0.222 (0.075) 2.97 0.003 Amount of nymph during the previous session -0.135 (0.071) -1.91 0.056 Random effect: site σ 2 = 0.33 (SD = 0.57) Negative binomial dispersion parameter 3.04 (SE = 0.81) Y = Female + nymph abundance Estimate (SE) z-value P-value Intercept 0.011 (0.316) 0.04 0.971 Vole abundance 0.006 (0.003) 1.98 0.048 Random effect: site σ 2 = 0.27 (SD = 0.52) Negative binomial dispersion parameter 1.53 (SE = 0.29) σ 2 is the variance attributable to random effect. Number of observations: Total = 88; Site = 4 Abbreviations:SD standard deviation, SE standard error Cayol et al. Parasites & Vectors (2017) 10:166 Page 5 of 11 Moreover, we identified risk periods - when humans are likely to encounter tick bites - in boreal forests and seasonal patterns that might be relevant for tick-borne pathogen circulation. Tick seasonality We identified that the highest tick abundance was in early summer (May-June) and early autumn (AugustSeptember), which are consequently the periods of increased risk for humans to encounter I. ricinus ticks in boreal forests. The same pattern of bimodal questing activity was previously found in southern Finland for nymph and adult ticks in coniferous and deciduous forests, whereas larvae showed a bimodal occurrence with a larger peak in September than in June [14]. Overall, two types of tick questing activity patterns have been described in Europe: in highly seasonal climates, such as those in central Europe, a bimodal questing activity with early spring and autumn peaks has been described for all life stages of I. ricinus [43]. However, in milder climates, with less climatic variation between seasons, only one peak of activity was observed for all life stages; in either spring or early summer [43]. In the present study, nymphs and adults showed bimodal activity, whereas larvae showed a unimodal activity pattern. This unimodal activity pattern could arise from egg production during the preceding year, the product of which overwintered as eggs or as larvae [44] or from egg production during the same spring. It could be argued that the inclusion of a year*month interaction term in the model would have captured between year seasonal variations suggested by the raw data (Additional file 1: Figure S1), and would have revealed both unimodal and bimodal activity patterns for larvae. However, data from a longer time series would be needed in order to clarify this point. Our results seem to indicate the coexistence of two age cohorts of ticks during the same year. Larvae detected on bank voles and in the vegetation in early summer become nymphs in autumn, which can exhibit two different behaviors: immediate questing behavior in autumn; or activity postponed until the next spring after a behavioral diapause [45–47]. In our study, the largest peak of nymphal activity was observed in spring, Fig. 2 Predicted number ± standard error (SE) of larvae and pooled nymphs and females in 100 m 2 of vegetation explained by vole abundance. Predictions are based on GLMM shown in Table 2 Table 3 Selected best model for I. ricinus larvae infestation load on an individual bank vole with estimated coefficients (in log scale) explained by month (from May to October, with May as a reference), year (from 2012 to 2015, with 2012 as a reference), sex (female as a reference), body mass in grams (centred values), presence of I. trianguliceps females and nymphs, presence of I. ricinus nymphs, vole abundance during the same session, questing larvae in vegetation during the same session, the interaction between centred body mass and sex and the interaction between sex and vole abundance. We defined site and individual nested in site as nested random structure Estimate (SE) z-value P-value Intercept -0.923 (0.318) -2.91 0.004 June 0.477 (0.243) 1.96 0.050 July -0.691 (0.277) -2.49 0.013 August -0.900 (0.342) -2.63 0.009 September -1.734 (0.413) -4.20 <0.005 October -2.768 (0.376) -7.36 <0.005 2013 0.720 (0.150) 4.79 <0.005 2014 -0.688 (0.275) -2.50 0.012 2015 -0.248 (0.169) -1.47 0.142 Male 0.996 (0.219) 4.55 <0.005 Body mass 0.020 (0.010) 2.03 0.043 Presence of I. trianguliceps female 0.402 (0.154) 2.61 0.009 Presence of I. trianguliceps nymphs 0.202 (0.101) 2.00 0.046 Presence of I. ricinus nymphs 0.526 (0.132) 3.97 <0.005 Vole abundance 0.033 (0.005) 6.43 <0.005 Amount of questing larvae during the same session 0.027 (0.009) 2.86 0.004 Interaction: Sex(Male)*Body mass 0.048 (0.016) 3.02 0.003 Interaction: Sex (Male)*Vole abundance -0.009 (0.003) -3.20 0.001 Random effects Site σ 2 = 0.06 (SD = 0.25) Individual nested in site σ 2 = 0.22 (SD = 0.47) Negative binomial dispersion parameter 1.70 (SE = 0.24) σ 2 is the variance attributable to random effect. Number of observations: Total = 1007; Site = 4, Site:Individual = 658 Abbreviations:SD standard deviation, SE standard error Cayol et al. Parasites & Vectors (2017) 10:166 Page 6 of 11 suggesting that the second behavior was predominant [45]. In addition, we found a peak in questing adults 2 years after the largest larval infestation, indicating a probable 2-year period between larvae blood meal and adults, and suggesting a 2 to 3 year life-cycle from egg to adult for I. ricinus in our study area. Furthermore, in an additional model, the amount of questing nymphs was explained, amongst others explanatory variables, by the total amount of larvae that fed on bank vole the year before (GLMM negative binomial: estimate (±SE) = 0.005 ± 0.002, P= 0.0071, see Additional file 1: Tables S6 and S7). This model confirmed firstly, that the variation in bank vole larval infestation was translated into nymph abundance and secondly, a 1 year delayed relationship between bank vole larval infestation and questing nymphs. We observed an effect of saturation deficit on larval and adult questing behavior, but not on nymphs, as described in other studies [40, 45]. Ticks respond to microclimate, but climatic variations measured in this study presumably reflect only roughly microclimatic variations and could explain the lack of association between saturation deficit and nymph activity found in our study. On the other hand, nymphs might also be acclimatized to local conditions and therefore their questing behavior may vary compared to nymphs studied in other locations [48, 49]. This is further supported by the optimal saturation deficit value of 3.16 mm Hg over which the adult questing activity decreased, when an optimum of 4.4 mm Hg has been previously noted elsewhere [9]. Table 4 Selected best model for I. ricinus nymph infestation load on an individual bank vole with estimated coefficients (in log scale) explained by month (from May to October, with May as reference), sex (female as reference), presence of I. trianguliceps larvae and females and presence of I. ricinus larvae, centered body mass and its squared value. We defined site and individual nested in site as nested random structure Estimate (SE) z-value P-value Intercept -2.994 (0.617) -4.86 <0.005 June -1.325 (0.385) -3.44 <0.005 July -1.360 (0.429) -3.17 0.002 August -2.103 (0.518) -4.06 <0.005 September -3.043 (0.643) -4.73 <0.005 October -2.956 (0.732) -4.04 <0.005 Male 1.787 (0.298) 6.00 <0.005 Body mass 0.219 (0.036) 6.14 <0.005 Body mass 2 -0.009 (0.003) -2.80 0.005 Presence of I. trianguliceps larvae 0.709 (0.247) 2.87 0.004 Presence of I. trianguliceps female 1.012 (0.318) 3.18 0.002 Vole abundance 0.014 (0.006) 2.30 0.021 Random effects Site σ 2 = 0.75 (SD = 0.87) Individual nested in site σ 2 = 0.01 (SD = 0.09) Negative binomial dispersion parameter 1.00 (SE = 0.46) σ 2 is the variance attributable to random effect. Number of observations: Total = 1,007; Site = 4; Site:Individual = 658 Abbreviations:SD standard deviation, SE standard error Fig. 3 Predicted number ± standard error (SE) of larvae and nymphs on a vole per month. Predictions are based on GLMM shown in Tables 3 and 4 Fig. 4 Predicted number ± standard error (SE) of larvae per bank vole (male and female separated) by vole abundance. Predictions are based on GLMM shown in Table 3 Cayol et al. Parasites & Vectors (2017) 10:166 Page 7 of 11 Our ratio larvae:nymph:adult was 3.5:3.3:2.0 and was therefore different from theoretical biological expectation (100:10:2) [30], indicating a possible underestimation of nymphs, and particularly of larvae, which show an important patchiness in distribution. Indeed, the blanket dragging technique is limited by variability in sampling efficiency given the nature of the substrate, the wind speed during sampling, and the height, type and growth stage of the vegetation [50, 51]. Moreover, the total tick population is not accessible by flagging given that diapausing ticks, parasitizing ticks, quiescent ticks or rehydrating individual are not questing in vegetation. Associating bank vole screening to the blanket dragging provided a broader view of immature ticks’population by beneficiating from the buffer effect that the host offers against the larvae patchiness and the drop of activity in case of unfavourable microclimate [32]. The two approaches are complementary and had led to similar results supporting further the idea that bank voles play an important role as host for immature ticks in the area. Dynamics of tick in immature stages and bank vole population are related Questing and parasitizing abundances of larvae and nymphs showed positive associations with bank vole abundance during a given session. Regarding nymphs, this positive relationship might arise from better engorgement success for larvae in high bank vole abundance. However, regarding larvae, this correlation does not imply a causative relationship since bank voles do not contribute to larvae production, which relies on large mammal availability [52]. Consequently, this positive relationship between the abundance of larvae and bank voles might reflect large-mammal density variations or might reveal a functional response: larvae may increase questing behavior in response to increased chemical signals produced by large bank vole populations [53]. This hypothesis requires further attention and needs to be experimentally quantified. Additionally, abundance of other species known to host adult stages needs to be quantified. The largest burden of nymph parasitism in voles was observed in May, whereas peaks of questing nymphs in the vegetation were observed in May-June and September. In May, vole populations are mainly composed of overwintered sexually active adults; highly mobile males exhibit large home ranges in their search for receptive females [54]. Therefore, the probability of encountering questing nymphs present at a low level in the recovering spring vegetation is increased [55]. In September, bank vole contact rate with nymphs might be lower due to taller vegetation, which allows nymphs to quest higher on plants, where they can contact larger mammal hosts [47]. Moreover, the bank vole develops an acquired resistance to ticks, leading to a significant reduction of infestation success after the first infestation [56, 57], which could lead to poor infestation success during the second nymph peak in September. However, our data (Table 3) provide little support for this hypothesis as regards larvae infestation that increases with animal weight, which is used here as a proxy for age, when a decrease in the relationship was expected under acquired immunity hypothesis. As a consequence, the main period for larval and nymph co-infestation on bank voles is in early summer. The epidemiological consequences of these co-infestations are discussed below. Concerning larvae, we identified an infestationpeakinJune,whichisinaccordancewiththe peak of larvae questing activity and in accordance with previous surveys [58]. Male bank voles were more commonly infested with nymphs than females and the infestation increased with bank vole abundance. This sex-specific infestation load has been described previously [59] and may not only be due to the immunosuppressing role of testosterone [60–62], but also to sex-specific behavioral differences, e.g. in home range sizes [63]. Surprisingly, we found larvae infestation differs with population density; females carried more larvae at high population density, whereas males carried more nymphs at any population density. A different use of vertical space by bank vole males and females in high population densities can be hypothesized, leading males to come into contact with more nymphs that quest higher in vegetation, whereas females, which exhibit aggressive defensive behavior against intruders during the reproductive season [64], would stay close to the ground, i.e. at larvae level. More attention should be paid to the use of vertical space by bank voles in order to clarify the potential role of vertical space use causing differences between individuals in their tick infestation load. Our data show a concomitant early summer questing activity between larvae and nymphs, and a parasitic aggregation between larvae and nymphs of I. ricinus on bank voles,which are relevant from an epidemiological point of view. The simultaneous activity of larvae and potentially infected nymphs occurs when rapidly rising temperatures in spring allow the simultaneous emergence of larvae and nymphs from overwintering diapause. In these conditions, pathogen transmission from infected nymphs to susceptible larvae can occur via simultaneous feeding on the same host, even without systemic infection of the host. This co-feeding transmission pathway is important for several zoonotic tickborne pathogens, especially those with short-lived or non-systemic infections in the rodent host, such as Anaplasma phagocytophilum or tick-borne encephalitis virus (TBEV), respectively [65–68]. Cayol et al. Parasites & Vectors (2017) 10:166 Page 8 of 11