Factors affecting patterns of tick parasitism on forest rodents in tick-borne encephalitis risk areas, Germany.

Identifying factors affecting individual vector burdens is essential for understanding infectious disease systems. Drawing upon data of a rodent monitoring programme conducted in nine different forest patches in southern Hesse, Germany, we developed models which predict tick (Ixodes spp. and Dermacentor spp.) burdens on two rodent species Apodemus flavicollis and Myodes glareolus. Models for the two rodent species were broadly similar but differed in some aspects. Patterns of Ixodes spp. burdens were influenced by extrinsic factors such as season, unexplained spatial variation (both species), relative humidity and vegetation cover (A. flavicollis). We found support for the 'body mass' (tick burdens increase with body mass/age) and for the 'dilution' hypothesis (tick burdens decline with increasing rodent densities) and little support for the 'sex-bias' hypothesis (both species). Surprisingly, roe deer densities were not correlated with larvae counts on rodents. Factors influencing the mean burden did not significantly explain the observed dispersion of tick counts. Co-feeding aggregations, which are essential for tick-borne disease transmission, were mainly found in A. flavicollis of high body mass trapped in areas with fast increase in spring temperatures. Locally, Dermacentor spp. appears to be an important parasite on A. flavicollis and M. glareolus. Dermacentor spp. was rather confined to areas with higher average temperatures during the vegetation period. Nymphs of Dermacentor spp. mainly fed on M. glareolus and were seldom found on A. flavicollis. Whereas Ixodes spp. is the dominant tick genus in woodlands of our study area, the distribution and epidemiological role of Dermacentor spp. should be monitored closely.

Introduction

Rodents are important hosts for the immature stages of hard ticks, and when taking a blood meal, ticks may transmit a range of tick-borne disease agents of medical and veterinary significance. In central Europe, ticks of the Ixodes ricinus complex can be infected with and subsequently transmit pathogens such as Borrelia burgdorferi spirochaetes, gram-negative bacteria of the family Anaplasmataceae and Rickettsiaceae and tick-borne encephalitis virus (Kurtenbach et al. 2002; Labuda and Nuttall 2004; Parola et al. 2005). Ticks of the genus Dermacentor, especially D. reticulatus, appear to expand their range in Germany (Dautel et al. 2006). Ticks of this genus are competent vectors of B. burgdorferi, rickettsia bacteria and tick-borne encephalitis virus and may also transmit protozoan piroplasms (Kahl et al. 1992; Randolph et al. 1996; Jongejan and Uilenberg 2004).

A key element for quantifying transmission rates of these pathogens is the vector burden on host individuals. Specifically, not only the mean per capita burden but also the level of aggregation within the population affects the intrinsic growth rate R 0 of a pathogen (Woolhouse et al. 1997; Hartemink et al. 2008). Tick-host systems and indeed most host–parasite systems are characterised by heterogeneities with respect to the probability of hosts being exposed to parasites and in turn to spread these among the population (Shaw et al. 1998; Lloyd-Smith et al. 2005). Parasite distributions within their host populations are generally most adequately described by a negative binomial distribution, characterised by the mean (μ) and a dispersion parameter (σ) (Shaw et al. 1998). The level of parasite aggregation is of crucial importance because the intrinsic growth rate of a pathogen R 0 increases with the degree of parasite aggregation (Woolhouse et al. 1997). The observed level of aggregation is often in accordance with the ‘80/20’ rule, where a fraction of approximately 20% of the population is responsible for approximately 80% of the disease transmission (Woolhouse et al. 1997; Perkins et al. 2003). Identifying characteristics of this rather small proportion of hosts which are responsible for the majority of disease spread is a central task for parasitological research and for designing effective control mechanisms (Woolhouse et al. 1997; Perkins et al. 2003; Lloyd-Smith et al. 2005). For the two tick-borne disease agents of major medical importance, B. burgdorferi, and especially for tick-borne encephalitis virus, co-feeding transmission (i.e. the pathogen transmission from infected nymphs to naïve larvae, all feeding on the same host and not involving amplification within the host) is of major importance for the maintenance of the pathogen in the tick population (Jones et al. 1987; Randolph et al. 1996; Hartemink et al. 2008). Thus, identification of factors affecting (1) mean tick burden, (2) level of aggregation and (3) the simultaneous presence of larval and nymphal ticks on a single host is central for understanding transmission dynamics of tick-borne diseases.

Several studies have investigated the patterns of tick parasitism in forests of Europe. Most studies largely focussed on single or few factors affecting individual mean tick burden in rodent populations such as seasonality (Radda 1968; Radda et al. 1969; L’Hostis et al. 1996; Randolph et al. 1999), extrinsic spatial factors such as habitat type or structure (Boyard et al. 2008; Paziewska et al. 2010), microclimate (Randolph and Storey 1999) or intrinsic features such as rodent species, sex, age and body mass (Randolph 1975; Nilsson and Lundqvist 1978; Matuschka et al. 1991; Humair et al. 1993; Perkins et al. 2003; Harrison et al. 2010). Recently, the abundance of larger ungulates such as roe deer (Capreolus capreolus) which are key hosts for adult ticks (Vor et al. 2010; Kiffner et al. 2010) has been hypothesised to contribute to increased tick densities (Gilbert 2010) and consequently to increased tick-borne disease incidence in humans (Linard et al. 2007; Rizzoli et al. 2009). However, effects of variable deer densities on individual tick burdens of rodents have rarely been tested (but see Harrison et al. 2010). Furthermore, high rodent densities might ‘dilute’ individual burdens since the available ticks might be spread across many hosts (Schmidt et al. 1999).

Given the wide range of factors potentially affecting individual tick burden, simultaneous testing of these variables is needed. Even if sufficient data are available, conventional statistical approaches are often inappropriate as they usually assume a constant aggregation level (Shaw et al. 1998; but see Brunner and Ostfeld 2008). In order to analyse a data set of tick burdens (Ixodes spp. and Dermacentor spp.) on yellow-necked mice (Apodemus flavicollis) and bank vole (Myodes glareolus), the two dominant rodent species in central European woodlands, we adopted a flexible modelling approach, general linear models for location, scale and shape (GAMLSS, Stasinopoulos and Rigby 2007). This statistical framework allows modelling of the mean and the dispersion parameter of a negative binomial distribution as a function of explanatory variables. In this framework, we test whether (1) external abiotic (e.g. season) factors influenced individual tick burdens; (2) external biotic factors, such as relative humidity, temperature, vegetation cover, roe deer and rodent density affected individual tick loads; and (3) intrinsic factors such as sex, age or body mass were correlated with tick burdens. Additionally, we tested whether the dispersion parameter varied with those parameters that influenced mean burdens.

In order to identify variables of categories (1)–(3) which might predict co-feeding, we used a statistical framework with binomial structure following Perkins et al. (2003). Additionally, we tested whether co-feeding of larvae and nymphs were correlated with the increase in spring temperature relative to the mean temperature of January. Fast temperature increase in spring time is thought to be the main driver for the seasonal synchrony of larvae and nymphs and thus for co-feeding aggregations (Randolph and Sumilo 2007).

Materials and methods

Study sites

We selected nine different forest patches (mean size 1,150 ha, range: 520–1,710 ha) in three forest districts (Beerfelden, Dieburg, Lampertheim) in the southern part of Hesse, Germany. The forest districts were located in counties (Bergstraße and Odenwaldkreis, Darmstadt-Dieburg and Bergstraße, respectively) defined as risk areas for TBEV (Robert Koch-Institute 2007).

Small mammal trapping

We conducted repeated rodent trapping in the first weeks of September 2007; May, July and August 2008; and May, July and August 2009. In each forest patch, we established two trapping grids on randomly selected intersections of a superimposed 1 km × 1 km grid. Each of the 18 trapping grids consisted of 36 Sherman live traps placed systematically in a 50 m × 50 m square (10-m inter-trap distance). We used fresh apple parts to bait rodents and placed hay into the traps to provide nesting material. We operated trapping grids for four consecutive nights and controlled the traps every morning. We transferred caught rodents in a plastic bag and identified them to species level based on morphological traits. We released non-target animals (e.g. Soricidae) immediately and euthanized rodents with CO2. The trapping and euthanizing protocol was authorised by the responsible authority. Overall, we caught 270 rodents, whereas bank voles M. glareolus (143 individuals) and yellow-necked mice A. flavicollis (106 individuals) combined represented 92% of all captures. Other species [A. sylvaticus (n = 15), Microtus agrestis (n = 4), M. arvalis (n = 1), and Mus musculus (n = 1)] were captured infrequently and were not considered for the statistical analyses. For each rodent, we assessed basic biometric characteristics such as sex (male/female), age (sub-adult, adult) and body mass. We carefully screened each rodent for ticks by combing the fur and by intensively searching the ears, head, throat, toes and tail. All detected ticks were removed using forceps and transferred in sterile tubes and stored at −80°C. Ticks were determined to genus (Dermacentor spp. (Koch 1844) and Ixodes spp. (Latreille 1795)) and stage (larvae and nymphs, no adults were found on rodents). Sample individuals determined to species level belonged to I. ricinus (Linneaus 1758) and Dermacentor reticulatus (Fabricius 1794). Since a rodent removal protocol was necessary for further virus screening (TBEV and Hantavirus) of the rodents (see, e.g. Ulrich et al. 2009), we approximated rodent density as number of individuals per 100 corrected trap nights (rodents per 100CTN). Trapping effort was corrected for closed traps without captures or captures of non-target animals. We calculated three different density indices: a mice density index (A. flavicollis, Apodemus sylvaticus and M. musculus captures per 100CTN), a vole density index (M. glareolus, M. agrestis and Microtus arvalis captures per 100CTN) and a rodent density index (all rodent species combined per 100CTN). During each trapping session, we visually estimated the percentage of vegetation cover in the herb layer in four categories: 0–24%, 25–49%, 50–74% and 75–100%.

Climatic data

At each trapping grid, we placed a weather data logger (Thermo/Hygro Button 23, Maxim Integrated Products, Inc., Sunnyvale, USA) at a tree trunk near the forest floor and with minimal exposition to solar radiation to record relative humidity and temperature. Since data loggers were not operated for the entire study period and frequently failed to store data, it was not possible to relate relative humidity and temperature to the grid- and time-specific individual tick burden. Hence, we calculated average temperature and relative humidity for the vegetation period 2009 (1 March–30 September) for each trapping grid. According to Randolph and Sumilo (2007), we estimated the spring temperature increase from February to April 2009 corrected by the mean temperature of January 2009.

Roe deer density estimation

We estimated densities of roe deer using line transect methodology (Buckland et al. 2001) and analysed the data with the software package Distance 5 Release 2 (Thomas et al. 2010). In early March 2008 and 2009, we drove a fixed circuit (mean length ± SD 18.3 ± 3.3 km) in each forest area. We repeated each circuit in one of the following nights. We counted roe deer with three persons; one person driving the car slowly (6–12 km h−1) and screening for animals on the transect line and two persons scanning both sides of the transect line with handheld spotlights. We measured sighting distances with a laser rangefinder and sighting angles with a compass. Considering that we used forest roads and hence that transects were not distributed randomly, our estimates should be regarded as density indices. However, these indices allow comparisons of roe deer densities among different forest areas and years. Because the numbers of roe deer sightings forest per area per year were low (mean 15.8 ± 6.4 SD), we pooled roe deer sightings according to the predominant terrain of the forest area. Based on Akaike’s information criterion (AIC) values, these pooled detection functions indicated a better fit than forest-area-specific detection functions. We discarded the largest 5% of the distances and used half-normal key function with cosine series expansion to fit the detection functions. Using these stratum-specific detection functions and the size-bias regression method to estimate cluster size, we estimated area- and year-specific roe deer densities. Because mean roe deer densities in the nine forest patches remained remarkably stable between 2008 and 2009 (Kendalls tau 0.93, p < 0.001, n = 9), we used the 2008 estimate also for the year 2007.

Modelling approach

For predicting host-species- and tick-genus-specific models of individual larval burdens, we ran several general additive models for location, scale and shape, defining the distribution as negative binomial (NBD type I) (Shaw et al. 1998). The modelling procedure was performed with the ‘gamlss’ package (Stasinopoulos and Rigby 2007) implemented in R (R Development Core Team 2005). Similar to Brunner and Ostfeld (2008), we used a stepwise forward model selection procedure. We started with the most basic extrinsic factors (seasonality, forest district) potentially influencing mean (μ) larval tick burdens. Then, always selecting the model with the lowest Akaike’s information criterion corrected for sample size (AICc), we tested whether inclusion of further extrinsic variables (climatic variables, vegetation cover in the shrub layer, roe deer density and rodent density) improved the models (Burnham and Anderson 2002). Further on, we tested whether intrinsic individual characteristics of rodents (sex, age and body mass) improved the model fit. Based on the selected model explaining the mean larval burden, we tested whether addition of a variable dispersion parameter (σ) enhanced the model fit. Since the level of aggregation is likely to be correlated with μ (Shaw et al. 1998), we tested whether σ was affected by those variables explaining μ. We used logistic regression, to test which factors affected presence of nymphal ticks on rodent individuals. Analogous to the larval models, we used a stepwise forward model selection procedure.

Results

Larval burden

Larval Ixodes spp. burdens on A. flavicollis (mean = 19, range = 0–129, SD = 22) were on average higher (Mann–Whitney U test, z = −8.96, p < 0.001) than on M. glareolus (mean = 6, range = 0–86, SD = 14). Almost all (98%, 104/106) A. flavicollis individuals were parasitized with at least one Ixodes spp. larvae while the larval Ixodes spp. prevalence in M. glareolus was 68% (98/144). Larval ticks were highly aggregated. In A. flavicollis, 20% of the most infested individuals harboured 56% of the entire larval Ixodes spp., and in M. glareolus the same proportion fed 81% Ixodes spp. larvae.

For predicting larval Ixodes spp. burdens on A. flavicollis, the model selection procedure provided most support for model K (Table 1). This model suggested that larval burdens were influenced by sampling month, with mean larval burdens being highest in July and lower in May and September (Fig. 1). Mean larval burdens were higher in the forest district Beerfelden compared to the forest districts of Dieburg and Lampertheim and slightly declined with increasing relative humidity during the vegetation period. Individual burdens were also associated with vegetation cover in the shrub layer whereas stands with ≥25% vegetation cover were associated with higher tick burdens compared to stands with <25% vegetation cover. There was statistical support that mean larval burden of A. flavicollis decreased with increasing rodent density. Further on, adult A. flavicollis showed higher larval burdens than sub-adult conspecifics. There was no support for including a variable dispersion parameter. Models which included dispersion parameters either as a function of month, forest district, relative humidity, vegetation cover, indexed rodent density or host age had poorer fits (based on AICC, models not shown) compared to model K without a variable dispersion parameter. Model K suggested the use of a constant dispersion parameter (−0.67, SE ± 0.14).

Table 1

Support for models explaining mean larval tick burdens (Ixodes spp.) on yellow-necked mice (A. flavicollis) and bank vole (M. glareolus)

Model	Apodemus flavicollis					Myodes glareolus
	Model letter	P	AICC	Δi	w i	Model letter	P	AICC	Δi	w i
Individual characteristics
μ is a function of age, sex and body mass	Q	15	804.72	3.11	0.09	Q	11	677.17	1.05	0.16
μ is a function of sex and body mass	P	14	808.59	6.97	0.01	P	10	678.43	2.31	0.08
μ is a function of sex and age	O	14	804.06	2.44	0.12	O	10	677.27	1.15	0.15
μ is a function of age and body mass	N	14	802.08	0.46	0.32	N	10	676.12	0.00	0.27
μ is a function of sex	M	13	815.93	14.32	0.00	M	9	689.77	13.65	0.00
μ is a function of body mass	L	13	805.94	4.33	0.05	L	9	677.71	1.59	0.12
μ is a function of age	K	13	801.61	0.00	0.41	K	9	676.55	0.43	0.22
Best model from below	J	12	813.60	11.99	0.00	I	8	691.02	14.90	0.00
Rodent density
μ is a function of rodent density	J	12	813.60	11.99	0.00	J	8	696.27	20.15	0.00
μ is a function of vole density	I	12	815.00	13.39	0.00	I	8	691.02	14.90	0.00
μ is a function of mice density	H	12	813.77	12.16	0.00	H	8	699.34	23.22	0.00
Best model from below	F	11	813.68	12.07	0.00	B	7	697.12	21.00	0.00
Roe deer density
μ is a function of roe deer density	G	12	814.51	12.90	0.00	G	8	698.11	21.99	0.00
Best model from below	F	11	813.68	12.07	0.00	B	7	697.12	21.00	0.00
Vegetation
μ is a function of vegetation cover	F	11	813.68	12.07	0.00	F	10	700.33	24.21	0.00
Best model from below	C	8	814.94	13.33	0.00	B	7	697.12	21.00	0.00
Climatic factors
μ is a function of relative humidity and temperature	E	9	817.33	15.71	0.00	E	9	701.68	24.56	0.00
μ is a function of temperature	D	8	817.42	15.81	0.00	D	8	698.66	22.54	0.00
μ is a function of relative humidity	C	8	814.94	13.33	0.00	C	8	699.36	23.24	0.00
Best model from below	B	7	815.67	14.06	0.00	B	7	697.12	21.00	0.00
Spatial factors
μ is a function of forest district	B	7	815.67	14.06	0.00	B	7	697.12	21.00	0.00
Best model from below	A	5	819.54	17.93	0.00	A	5	714.21	38.09	0.00
Seasonal dynamics
μ is a function of season	A	5	819.54	17.93	0.00	A	5	714.37	37.62	0.00

Models are grouped from bottom (simple models) to top (more complex). We first addressed the fundamental extrinsic factors season and forest district and then included further extrinsic factors (climatic factors, vegetation, roe deer density and rodent density) and intrinsic factors (age, body mass and sex). In each case, the best model from the set below was chosen based on minimum AICc values. P indicates the number of parameters used for fitting each model, Δ is the difference in AICc and w is the AICc weight based on all models (Burnham and Anderson 2002). The models with most support are in italics

Fig. 1

The effect of a month, b forest district, c relative humidity during the vegetation period, d vegetation cover in the shrub layer, e indexed rodent density and f host age on the mean larval (Ixodes spp.) burden on yellow-necked mice (A. flavicollis) as predicted by model K (Table 1). Dashed lines indicate standard errors

In order to explain variation in larval Ixodes spp. burdens on M. glareolus, we found most support for model N (Table 1). Similar to the A. flavicollis–Ixodes spp. larvae model, this model suggested including the sampling month and the forest district as explanatory variables (Fig. 2). Individual burden declined with increasing indexed vole densities. Among the intrinsic factors, host age (adult > sub-adult) and body mass were positively correlated with mean larval burdens. Again, models with variable dispersion parameters performed worse than model N with a constant (0.21, SE ± 0.17) dispersion parameter.

Fig. 2

The effect of a month, b forest district, c indexed vole density, d host age and e host body mass on the mean larval (Ixodes spp.) burden on bank vole (M. glareolus) as predicted by model N (Table 1). Dashed lines indicate standard errors

Larval Dermacentor spp. burdens were not statistically different (Mann–Whitney U, z = −0.97, p = 0.333) between A. flavicollis (mean = 3, range = 0–88, SD = 12) and M. glareolus (mean = 3, range = 0–100, SD = 13). Prevalence of larval Dermacentor spp. was slightly higher in A. flavicollis (20%, 21/106) than in M. glareolus (15%, 21/144). The larval Dermacentor spp. was highly aggregated whereas the 20% of the most infested A. flavicollis and M. glareolus individuals fed 100% of the counted Dermacentor spp. larvae.

Models for explaining larval Dermacentor spp. burden on A. flavicollis indicated no significant effect of sampling month, and thus models did not include this variable. Mean larval Dermacentor spp. burden on A. flavicollis was best explained by model M (Table 2). Once more, this model indicated an effect of forest district (Fig. 3). Further on, the model suggested that Dermacentor spp. burden were positively correlated with the average temperature during the vegetation period and that individual burdens were higher in male A. flavicollis compared to female conspecifics. Models incorporating a variable dispersion parameter (either modelled as a function of forest district, temperature or host sex) had higher AICC values, and we thus favoured the model with a constant dispersion parameter (2.30e + 00, SE ± 2.80e–01).

Table 2

Support for models explaining mean larval tick burdens (Dermacentor spp.) on yellow-necked mice (A. flavicollis) and bank vole (M. glareolus)

Model	Apodemus flavicollis					Myodes glareolus
	Model letter	P	AICC	Δi	w i	Model letter	P	AICC	Δi	w i
Individual characteristics
Μ is a function of age, sex and body mass	Q	8	242.49	3.93	0.05	–b
Μ is a function of sex and body mass	P	8	240.86	2.30	0.10	–b
Μ is a function of sex and age	O	8	240.27	1.71	0.14	–b
Μ is a function of age and body mass	N	8	246.00	7.44	0.01	–b
Μ is a function of sex	M	7	238.56	0.00	0.32	–b
Μ is a function of body mass	L	7	243.94	5.38	0.02	F	7	245.52	2.15	0.19
Μ is a function of age	K	7	243.70	5.14	0.02	E	7	245.21	1.84	0.22
Best model from below	D	6	242.33	3.77	0.05	D	6	243.37	0.00	0.56
Rodent density
Μ is a function of rodent density	J	7	244.51	5.95	0.02	–b
Μ is a function of vole density	I	7	243.47	4.91	0.03	–b
Μ is a function of mice density	H	7	242.42	3.86	0.05
Best model from below	D	6	242.33	3.77	0.05	D	6	243.37	0.00	0.56
Roe deer density
Μ is a function of roe deer density	G	7	243.97	5.41	0.02
Best model from below	D	6	242.33	3.77	0.05	D	6	243.37	0.00	0.56
Vegetation
Μ is a function of vegetation cover	F	9	242.59	4.03	0.04
Best model from below	D	6	242.33	3.77	0.05	D	6	243.37	0.00	0.56
Climatic factors
Μ is a function of relative humidity and temperature	E	7	244.40	5.84	0.02	–b
Μ is a function of temperature	D	6	242.33	3.77	0.05	D	6	243.37	0.00	0.56
Μ is a function of relative humidity	C	6	243.43	4.87	0.03	C	6	252.23	8.86	0.01
Best model from below	B	5	242.45	3.89	0.05	B	5	250.36	6.99	0.02
Spatial factors
Μ is a function of forest district	B	5	242.45	3.89	0.05	B	5	250.36	6.99	0.02
Best model from below	–/–a					–/–a
Seasonal dynamics
Μ is a function of season	A	5	242.35	3.79	0.05	A	5	256.98	13.61	0.00

Models are grouped from bottom (simple models) to top (more complex). We first addressed the fundamental extrinsic factors season and forest district and then included further extrinsic factors (climatic factors, vegetation, roe deer density and rodent density) and intrinsic factors (age, body mass and sex). In each case, the best model from the set below was chosen based on minimum AICc values. P indicates the number of parameters used for fitting each model, Δ is the difference in AICc and w is the AICc weight based on all models (Burnham and Anderson 2002). The models with most support are in italics

aThe effect of season was not significant (p > 0.10) and hence was not included in further models

bRedundant combination of variables

Fig. 3

The effect of a forest district, b average temperature during the vegetation period and c host sex on the mean larval (Dermacentor spp.) burden on yellow-necked mice (A. flavicollis) as predicted by model M (Table 2). Dashed lines indicate standard errors

For predicting larval Dermacentor spp. burdens on M. glareolus, model D was selected (Table 2). Several model combinations were not possible due to redundant factor combinations. The selected model suggested similar effects (forest district and average temperature) as the larval Dermacentor spp. model for A. flavicollis, except that it did not include host sex as an explanatory variable (Fig. 4). Also, for this model, a constant dispersion parameter was considered (1.92, SE ± 0.29).

Fig. 4

The effect of a forest district and b average temperature during the vegetation period on the mean larval (Dermacentor spp.) burden on bank vole (M. glareolus) as predicted by model D (Table 2). Dashed lines indicate standard errors. Note the different scales on the y-axes

Nymphal burden

Ixodes spp. nymphs were found on 16% (17/106) of captured A. flavicollis with a maximum of eight nymphs on one individual. Prevalence of nymphs on M. glareolus was 10% (15/144) whereas one individual was infested by 25 nymphs.

A. flavicollis individuals parasitized by Ixodes spp. nymphs had higher average larval Ixodes spp. burdens (mean = 47 ± SE 9) compared to individuals without nymphs (13 ± 1). Individuals with an Ixodes spp. nymph also fed more Dermacentor spp. larvae (12 ± 6) than individuals without a nymph (1 ± 5).

We observed similar patterns in M. glareolus: Individuals parasitized by an Ixodes spp. nymph fed on average more Ixodes spp. larvae (22 ± 6) than conspecifics not feeding nymphs (4 ± 1). The same individuals, however, fed similar numbers of Dermacentor spp. larvae (2 ± 1) compared to individuals without Ixodes spp. nymphs (3 ± 1).

The model selection approach suggested that model P (Table 3) best explained prevalence of Ixodes spp. nymphs on A. flavicollis. The logistic regression model provided support that presence of Ixodes spp. nymphs on A. flavicollis was associated with sampling month (July > May, p = 0.025, September > May, p = 0.71) and forest district (Lampertheim < Beerfelden, p = 0.0048; Dieburg < Beerfelden, p = 0.18) and that it was positively correlated with spring warming rate (coefficient = 1.33, p = 0.031) and body mass (coefficient = 0.14, p = 0.04) of individual yellow-necked mice.

Table 3

Support for models explaining the presence (¤) of nymphal ticks (Ixodes spp.) on yellow-necked mice (A. flavicollis) and bank vole (M. glareolus)

Model	Apodemus flavicollis					Myodes glareolus
	Model letter	P	AICC	Δi	w i	Model letter	P	AICC	Δi	w i
Individual characteristics
¤ is a function of age, sex and body mass	U	9	83.18	2.86	0.06	U	8	92.96	4.38	0.01
¤ is a function of sex and body mass	T	8	80.85	0.53	0.18	T	7	91.31	2.73	0.03
¤ is a function of sex and age	S	8	83.44	3.12	0.05	S	7	90.85	2.27	0.03
¤ is a function of age and body mass	R	8	82.55	2.23	0.08	R	7	90.76	2.18	0.03
¤ is a function of sex	Q	7	82.70	2.38	0.07	Q	6	90.46	1.88	0.04
¤ is a function of body mass	P	7	80.32	0.00	0.24	P	6	89.20	0.62	0.08
¤ is a function of age	O	7	83.23	2.91	0.06	O	6	88.69	0.11	0.10
Best model from below	C	6	83.12	2.80	0.06	M	5	88.58	0.00	0.10
Rodent density
¤ is a function of rodent density	N	7	85.09	4.77	0.02	N	5	89.23	0.65	0.07
¤ is a function of vole density	M	7	85.17	4.85	0.02	M	5	88.58	0.00	0.10
¤ is a function of mice density	L	7	85.22	4.90	0.02	L	5	90.66	2.08	0.04
Best model from below	C	6	83.12	2.80	0.06	E	4	88.81	0.22	0.09
Roe deer density
¤ is a function of roe deer density	K	7	85.30	4.98	0.02	K	5	90.82	2.24	0.03
Best model from below	C	6	83.12	2.80	0.06	E	4	88.81	0.22	0.09
Vegetation
¤ is a function of vegetation cover	J	9	89.53	9.21	0.00	J	7	92.23	3.65	0.02
Best model from below	C	6	83.12	2.80	0.06	E	4	88.81	0.22	0.09
Climatic factors
¤ is a function of spring warming rate, relative humidity and temperature	I	8	85.77	5.45	0.02	I	6	92.07	3.49	0.02
¤ is a function of humidity and temperature	H	7	89.32	9.00	0.00	H	5	89.97	1.39	0.05
¤ is a function of spring warming rate and temperature	G	7	83.52	3.20	0.05	G	5	90.00	1.42	0.05
¤ is a function of spring warming rate and relative humidity	F	7	85.40	5.08	0.02	F	5	91.27	2.69	0.03
¤ is a function of temperature	E	6	87.30	6.98	0.01	E	4	88.81	0.22	0.09
¤ is a function of relative humidity	D	6	87.84	7.52	0.01	D	4	89.15	0.56	0.08
¤ is a function of spring warming rate	C	6	83.12	2.80	0.06	C	4	91.07	2.48	0.03
Best model from below	B	5	85.65	5.33	0.02	B	3	89.20	0.62	0.08
Spatial factors
¤ is a function of forest district	B	5	85.65	5.33	0.02	B	3	89.20	0.62	0.08
Best model from below	A	3	88.64	8.32	0.00	–/–a
Seasonal dynamics
¤ is a function of season	A	3	88.64	8.32	0.00	A	3	98.91	10.33	0.00

Models are grouped from bottom (simple models) to top (more complex). We first addressed the fundamental extrinsic factors season and forest district and then included further extrinsic factors (climatic factors, vegetation, roe deer density and rodent density) and intrinsic factors (age, body mass and sex). In each case, the best model from the set below was chosen based on minimum AICc values. P indicates the number of parameters used for fitting each model, Δ is the difference in AICc and w is the AICc weight based on all models (Burnham and Anderson 2002). The models with most support are in italics

aThe effect of season was not significant (p > 0.10) and hence was not included in further models

Presence of Ixodes spp. nymphs on M. glareolus was best explained by model M (Table 3). This model suggested that forest district (Dieburg < Beerfelden, p = 0.39, Lampertheim < Beerfelden, p = 0.46), temperature during the vegetation period (coefficient = 9.56, p = 0.43) and indexed vole density (coefficient = −0.18, p = 0.14) affected the probability that a nymphal Ixodes spp. was present on a bank vole. Yet all model parameters were insignificant (p values > 0.05), suggesting cautious treatment of this model.

Prevalence of nymphal Dermacentor spp. was very low. Since only 2% (2/106) of the captured A. flavicollis was parasitized by Dermacentor spp. nymphs, we did not analyse this host-tick system statistically. The two infested yellow-necked mice infested with a Dermacentor spp. nymph tended to have higher Ixodes spp. larvae burdens (57 ± 45 vs. 18 ± 2) and Dermacentor spp. burdens (39 ± 1 vs. 2 ± 1) than conspecifics without a Dermacentor spp. nymph. Six percent (9/144) of the captured M. glareolus individuals showed a prevalence of Dermacentor spp. nymphs, whereas the maximum per capita count was three nymphs. Individuals parasitized by Dermacentor spp. nymphs fed similar numbers of Ixodes spp. ticks (7 ± 6 vs. 6 ± 1) but on average higher numbers of Dermacentor spp. larvae (40 ± 10 vs. 1 ± 1) compared to conspecifics not feeding a nymph. The selected model N (Table 4) explaining presence of nymphal Dermacentor spp. on M. glareolus should be regarded conservatively. In this model, only one variable, indexed mice density (showing a positive correlation), reached statistical significance (p = 0.01). Other variables in the model (forest district, relative humidity, host age and host body mass) were insignificant (p > 0.05).

Table 4

Support for models explaining the presence (¤) of nymphal ticks (Dermacentor spp.) on bank vole (M. glareolus)

Model	Myodes glareolus
	Model letter	P	AICC	Δi	w i
Individual characteristics
¤ is a function of age, sex and body mass	Q	8	34.30	2.23	0.12
¤ is a function of sex and body mass	P	7	36.02	3.95	0.05
¤ is a function of sex and age	O	7	38.23	6.16	0.02
¤ is a function of age and body mass	N	7	32.07	0.00	0.36
¤ is a function of sex	M	6	36.26	4.19	0.04
¤ is a function of body mass	L	6	34.31	2.24	0.12
¤ is a function of age	K	6	36.02	3.95	0.05
Best model from below	H	5	34.14	2.07	0.13
Rodent host density
¤ is a function of rodent density	J	5	35.40	3.33	0.07
¤ is a function of vole density	I	5	36.34	4.27	0.04
¤ is a function of mice density	H	5	34.14	2.07	0.13
Best model from below	C	4	41.08	9.01	0.00
Roe deer density
¤ is a function of roe deer density	G	5	43.18	11.11	0.00
Best model from below	C	4	41.08	9.01	0.00
Vegetation
¤ is a function of vegetation cover	F	7	42.15	10.08	0.00
Best model from below	C	4	41.08	9.01	0.00
Climatic factors
¤ is a function of relative humidity and temperature	E	5	43.18	11.11	0.00
¤ is a function of temperature	D	4	41.36	9.29	0.00
¤ is a function of relative humidity	C	4	41.08	9.01	0.00
Best model from below	B	3	44.84	12.77	0.00
Spatial factors
¤ is a function of forest district	B	3	44.84	12.77	0.00
Best model from below	–/–a
Seasonal dynamics
¤ is a function of season	A	3	58.91	26.84	0.00

Models are grouped from bottom (simple models) to top (more complex). We first addressed the fundamental extrinsic factors season and forest district and then included further extrinsic factors (climatic factors, vegetation, roe deer density and rodent density) and intrinsic factors (age, body mass and sex). In each case, the best model from the set below was chosen based on minimum AICc values. P indicates the number of parameters used for fitting each model, Δ is the difference in AICc and w is the AICc weight based on all models (Burnham and Anderson 2002). The models with most support are in italics

aThe effect of season was not significant (p > 0.10) and hence was not included in further models

Discussion

After controlling for season and unexplained spatial variation (forest district entered as factor), we found that several extrinsic and intrinsic factors influence the investigated rodent hard-tick systems.

Factors affecting mean Ixodes spp. larvae burdens

It was expected that larval Ixodes spp. burdens were higher in A. flavicollis than in M. glareolus (cf. Matuschka et al. 1991; Humair et al. 1993; Boyard et al. 2008; Paziewska et al. 2010) because M. glareolus may acquire resistance against I. ricinus, the dominant tick species in central Europe (Dizij and Kurtenbach 1995). Highest mean larval burdens were observed in July, whereas we cannot exclude that larval activity might have peaked in June (cf. Randolph 2004). Mean infestation levels of A. flavicollis were also influenced by relative humidity during the vegetation period and by vegetation cover. The effect of relative humidity appears inconclusive: On the one hand, Ixodes ticks require high relative humidity [e.g. they are inactive at relative humidity <70% (Aeschlimann 1972)]; on the other hand, we observed that larval tick burden slightly declined with increasing relative humidity. This finding is in contrast to tick feeding experiments (Randolph and Storey 1999) but in accordance with results from a field study where the numbers of larval Ixodes spp. were also higher under drier conditions (82–89% vs. 91–98% relative humidity) (Boyard et al. 2008). To complicate this issue, vegetation cover (which is usually positively correlated with relative humidity on the forest floor) affected the mean number of larval Ixodes spp. The underlying factors causing these apparently contradictory findings need further experimental clarification. Unambiguously, our findings lend support to the ‘dilution hypothesis’ (Schmidt et al. 1999; Brunner and Ostfeld 2008), i.e. with increasing rodent densities, the mean per capita larval burden of A. flavicollis declines. Further on, adult individuals had higher levels of infestations than sub-adults, which is in line with other studies of tick–rodent systems (Brunner and Ostfeld 2008; Harrison et al. 2010), given the probable correlation between age and body mass (cf. Brunner and Ostfeld 2008).

Mean infestation levels of M. glareolus with Ixodes spp. larvae were affected by similar variables also affecting infestation levels of A. flavicollis. The ‘dilution effect’ was, however, related to indexed vole density (vs. combined rodent density index in A. flavicollis). Further on, adult bank voles and heavier bank voles fed on average more larval Ixodes spp. than younger and lighter conspecifics which largely supports the ‘body size’ hypothesis (Harrison et al. 2010).

In both rodent species, we found no support for the ‘sex-bias’ hypothesis, which is in contrast to similar tick–rodent systems in Europe (Harrison et al. 2010; Boyard et al. 2008) or in the USA (e.g. Brunner and Ostfeld 2008). Further on, it was unexpected that roe deer densities were not correlated with Ixodes spp. larval counts on forest rodents (cf. Gilbert 2010). Potentially, the failure to detect a significant effect of roe deer density was due to the different scales at which roe deer density estimation and rodent trapping were conducted. Moreover, the range of roe deer density indices (2.0–9.9 deer per square kilometre) might not be wide enough to detect a significant relationship between roe deer density and Ixodes larvae density.

Factors affecting mean Dermacentor spp. larvae burdens

Dermacentor spp. ticks were found in all forest districts but in one rodent trapping grid in the Rhine valley, this tick species was very abundant, supporting the hypothesis that Dermacentor spp. is associated with and/or expands along rivers (Bullová et al. 2009; Zygner et al. 2009). The relative abundance of Dermacentor spp. varied considerably at small scale whereas ticks of this genus appear to prefer areas with higher average temperatures; this might be important for potential further expansion of Dermacentor species and associated diseases with respect to climate change scenarios. In contrast to other studies which identified M. glareolus as the main host species (e.g. Randolph et al. 1999; Paziewska et al. 2010), we did not detect an apparent preference for a certain rodent species. Male A. flavicollis, which usually have a relative large home range (Schwarzenberger and Klingel 1994), were disproportionally infested with larval Dermacentor spp.

Dispersion of larval ticks

Larval ticks were highly aggregated on their rodent hosts, and in the case of Ixodes spp. the observed patterns broadly confirmed the ‘20/80’ rule (Woolhouse et al. 1997). In Dermacentor spp., the level of aggregation was even more pronounced, possibly due to the spatial clumping of these ticks. Whereas we found several factors affecting the variation in mean infestation levels, we failed to detect variables affecting the level of aggregation. This was rather disappointing since the level of aggregation is as important as the mean infestation level. Our approach suggests that other undocumented variables might influence the dispersion of larval ticks in rodent populations. Recent findings suggest that individual space use, which is not necessarily correlated with attributes such as sex, age or host density, affects the distribution of ticks among their hosts (Boyer et al. 2010). This would offer a mechanistic explanation for differences in mean tick loads but also for different aggregation levels. Incorporating individual space use of hosts as an explanatory variable is, however, not feasible in a removal study and would necessitate a capture–recapture design.

Factors affecting co-feeding

Overall, we found very few nymphs infesting forest rodents. The typical rodent individual infested with an Ixodes spp. nymph was an A. flavicollis of high body mass, captured in July. Co-feeding of Ixodes spp. ticks in M. glareolus appeared to be a rather erratic event. We found empirical evidence that the spring warming rate was positively correlated with co-feeding presence. A fast increase in spring temperatures relative to January temperatures promotes seasonal synchrony of larval and nymphal activity peaks (Randolph and Sumilo 2007). Rodents that fed at least one nymph were also disproportionally infested with larval ticks, being in accordance with Perkins et al. (2003). In contrast to Perkins et al. (2003), our data do not suggest that “sexually mature males of high body mass” were mainly feeding nymphs. With our data (which only contain few co-feeding aggregations), we only revealed an effect of body mass. Potentially, drawing upon a larger sample size would identify also further intrinsic variables such as host sex. Given a longer time horizon and a larger sample size of tick-infested rodents, it would be interesting to test the effect of time lags and to test explicitly whether high rodent densities in a given year translate into high nymph densities in the following year (Ostfeld et al. 2006; Rosa et al. 2007).

Presence of Dermacentor spp. nymphs was also rather inconsistent. Nymphs of this genus were predominantly found on M. glareolus, being in accordance with Paziewska et al. (2010). Only one variable (indexed mice density) was statistically associated with nymphal presence. Given the fast life cycle of D. reticulatus (and D. marginatus) (Hillyard 1996) and the strong association between larvae and rodents, high rodent densities during spring time might boost nymphal Dermacentor spp. densities in early summer. As an analogue to the Ixodes spp. system, an advanced study drawing upon a larger sample size of rodents infested with Dermacentor spp. nymphs should explicitly test the effect of spring rodent density on the prevalence/abundance of nymphs at a later stage (i.e. during early summer).

Conclusion

Multiple factors appear to influence tick burdens on forest rodent species. We provide evidence for the ‘dilution’ and for the ‘body size’ hypotheses but find little support for the ‘sex-bias’ hypothesis. Co-feeding aggregations which are essential for tick-borne disease transmission (especially tick-borne encephalitis virus) were mainly found in yellow-necked mice of high body mass trapped in areas showing a fast increase in spring temperatures. Whereas Ixodes spp. is the dominant tick genus in woodlands of our study area, Dermacentor spp. is locally very abundant. Its occurrence and its epidemiological role should be monitored closely.