- Research
- Open Access
- Published:

# Estimating *Ixodes ricinus* densities on the landscape scale

*International Journal of Health Geographics*
**volume 14**, Article number: 23 (2015)

## Abstract

### Background

The study describes the estimation of the spatial distribution of questing nymphal tick densities by investigating *Ixodes ricinus* in Southwest Germany as an example. The production of high-resolution maps of questing tick densities is an important key to quantify the risk of tick-borne diseases. Previous *I. ricinus* maps were based on quantitative as well as semi-quantitative categorisations of the tick density observed at study sites with different vegetation types or indices, all compiled on local scales. Here, a quantitative approach on the landscape scale is introduced.

### Methods

During 2 years, 2013 and 2014, host-seeking ticks were collected each month at 25 sampling sites by flagging an area of 100 square meters. All tick stages were identified to species level to select nymphal ticks of *I. ricinus*, which were used to develop and calibrate Poisson regression models. The environmental variables height above sea level, temperature, relative humidity, saturation deficit and land cover classification were used as explanatory variables.

### Results

The number of flagged nymphal tick densities range from zero (mountain site) to more than 1,000 nymphs/100 m^{2}. Calibrating the Poisson regression models with these nymphal densities results in an explained variance of 72 % and a prediction error of 110 nymphs/100 m^{2} in 2013. Generally, nymphal densities (maximum 374 nymphs/100 m^{2}), explained variance (46 %) and prediction error (61 nymphs/100 m^{2}) were lower in 2014. The models were used to compile high-resolution maps with 0.5 km^{2} grid size for the study region of the German federal state Baden-Württemberg. The accuracy of the mapped tick densities was investigated by leave-one-out cross-validation resulting in root-mean-square-errors of 227 nymphs/100 m^{2} for 2013 and 104 nymphs/100 m^{2} for 2014.

### Conclusions

The methodology introduced here may be applied to further tick species or extended to other study regions. Finally, the study is a first step towards the spatial estimation of tick-borne diseases in Central Europe.

## Background

The study describes the estimation of the spatial distribution of questing tick densities using the example of *Ixodes ricinus* in Southwest Germany. *Ixodes ricinus* is the most well-known and studied European tick species transmitting various arthropod-borne diseases [1, 2]. Suitable habitats are found all over in Germany [3]. The knowledge of the spatial distribution of tick densities is a key factor in quantifying the risk of tick-borne diseases such as tick-borne encephalitis, Lyme borreliosis, anaplasmosis, rickettsiosis, and others.

Previous studies on the spatial distribution of ticks mainly comprise the determination of the presence and absence of single or multiple species. Frequently, these tick data were exclusively published for selected sites, but recently more and more models for spatial interpolation or analysis have been developed. These spatial distribution models, also known as ecological niche models, are generally applied on a continental scale. Therefore, various statistical models to describe species habitats have been developed. An overview on models as well as methodological caveats in the modelling with examples for *I. ricinus* was described by Estrada-Peña and co-authors [4]. Several applications focus on the projection of climate niche models to estimate future scenarios. These studies provide maps of suitability indices and presence/absence maps of both the present and the future distribution of *I. ricinus*, covering Europe and Northern Africa [4, 5].

Calculations of the spatial distribution of tick densities, however, are rare. The few papers available comprise the mapping of the nymphal density of *I. pacificus* in California using GIS and satellite derived data [6], the mapping of the density of *I. scapularis* in the eastern United States using geo-statistical interpolation [7] and a zero-inflated negative binominal model [8]. In Europe, the spatial distribution of tick densities was exclusively estimated on local scales. Studies comprise the mapping of *I. ricinus* densities at two sites in the Czech Republic [9] and in a German nature reserve [10]. The two latter, however, are simply based on semi-quantitative categories of tick densities classified according to the number of ticks collected at different study sites, which were characterised by typical vegetation cover. No statistical models were developed for these density maps [9, 10]. In a further local scale analysis [11], generalized estimating equations were applied for mapping of *I. ricinus* densities in an Italian nature reserve (model domain: 5 km^{2}). Although other statistical models for *I. ricinus* were developed, these were not used to compile density maps [12]. In contrast, the study presented here, is a contribution to the quantitative spatial density estimation of ixodid ticks on the landscape scale (35,750 km^{2}). Here, the landscape scale is defined for model domains of 100–50,000 km^{2} (Fig. 1), where the environmental variables land-use, topography and climate are affecting the distribution of ticks. Thus, biotic interactions expressed for example by vector-to-host ratios and other local phenomena were not considered.

A total of 25 sampling sites in Southwest Germany, in the federal state Baden-Württemberg, were selected according to their vegetation cover and tick habitat suitability. *Ixodes ricinus* was flagged over an area of 100 m^{2} each month from March to October for each site. These field data were used to develop a generalised linear model (Poisson regression) describing the total number of nymphal ticks at the different sampling sites, which were collected during 1 year by monthly flagging an area of 100 m^{2}. Air temperature, air humidity, height above sea level and land cover classification were used as explanatory variables. These were preselected according to their biological relevance. Air temperature was selected, since thermal conditions determine the rate of tick development from stage to stage and activity patterns, which in turn influence survival success and population density [10, 14, 15]. This, and potentially other synergetic processes, leads to decreasing tick densities with increasing altitude [14]. Ticks and other terrestrial vertebrates must maintain their body water balance to survive. In a subsaturated atmosphere, unavoidable body water losses in ticks occur, e.g., by discharging faeces and urine, by gas exchange, and via the cuticle. Saturation deficit, the difference between the saturation vapour pressure at 100 % relative humidity (at a given temperature) and the given vapour pressure, is to a large extent a measure of the drying power of the atmosphere and as such largely responsible for the rate of body water loss of a given tick in a given situation [16]. Unfed ticks and some engorged ticks are capable of compensating suffered body water losses by active water vapour uptake when the ambient relative humidity surpasses a certain threshold, the so-called critical equilibrium humidity [17]. As a consequence, both saturation deficit and relative humidity must be taken in consideration when investigating the influence of humidity on ticks (Olaf Kahl 2015, personal communication). The habitat type influences not only the microclimatic conditions, but also the host spectrum and abundance. *Ixodes ricinus* prefers woodland habitats such as broadleaved, mixed or coniferous forests [10, 15].

As a result, predictive models were applied to extrapolate the tick densities to the entire region of Baden-Württemberg. A spatial resolution of 30 arc seconds, corresponding to about 0.5 km^{2}, was selected to depict maximum features of the tick density taking into account simultaneously the resolution of the climate data. Although the tick density maps were exclusively developed for *I. ricinus*, the method may be applied for other tick species as well.

## Methods

### Study area

The model domain covers the entire region of the German federal state Baden-Württemberg with an area of 35,750 km^{2}. It is located in the southwest of Germany (Fig. 2a) and has a wide ecological variability with low lands and mountains up to 1,500 m height above sea level (Fig. 2b). The historical natural coverage was woodland with European beech (*Fagus sylvatica*) as the main tree species. Today the natural landscape has been widely replaced by areas of anthropogenic utilization. The main tree species is the economically important European spruce (*Picea abies*) with 38 % followed by European beech with 21 % of the total forest area. Following the well-known Köppen-Geiger climate classification [18], the climate conditions in Baden-Württemberg as well as in the entire region of Germany were characterized by Cfb climate (C = warm temperate climate, f = fully humid, b = warm summers). Even under climate change conditions the Cfb climate in Baden-Württemberg will be preserved [19]. Total annual precipitation varies between 600 and 2,000 mm.

### Tick sampling and in situ measurements

From March to October (period 2013–2014) host-seeking ticks were collected at 25 sites (Table 1) by dragging a 1 m^{2} cotton flag over 100 m^{2} of woodland scrub and litter. The cloth was checked visually every 10 m^{2} to remove adult and nymphal ticks, which were collected in separate vials for every drag. The standardised sampling took place every month on dry days (no rain) at comparable daytime. All tick stages were identified to species level [20–22], but only nymphal ticks of *I. ricinus* were used in this study. Questing tick density was calculated as the total number of nymphs monthly collected per 100 m^{2} during 1 year, following Schulz and co-authors [15]. At each sampling site a weather station was operated throughout the year. From these continuous measurements, mean annual values of air temperature and relative humidity were calculated. The saturation deficit was estimated by applying the Magnus formula [23], which provides a functional relationship between saturation vapour pressure, temperature and relative humidity. Main vegetation cover was determined according to the *coordination of information on the environment* (CORINE) land classification [24]. Finally, geographical coordinates and the heights above sea level were documented (Table 1).

### Gridded environmental data

The gridded data are variables available on a high-resolution raster covering the entire study area. They were used to predict the spatial distribution of the tick density. Exclusively variables with biological significance for the life cycle of ticks were used [15, 25–28]. Figure 2 depicts these explanatory variables comprising the height above sea level, the CORINE land cover classification as well as the air temperature, relative humidity and saturation deficit.

Gridded temperature fields for Germany were provided by the climate data centre of the German weather service [29]. This dataset covers the period from 1881 up to present with a temporal resolution of 1 day and a spatial resolution of 1 km. The mean annual temperature for the years 2013 (Fig. 2d) and 2014 (not shown) were calculated and applied. As humidity measures are not included in this dataset, they were taken from the HYRAS dataset [30]. The HYRAS dataset comprises gridded data with a temporal resolution of 1 day and a spatial resolution of 5 km for the period 1951–2006. As this dataset does not cover the period 2013–2014, climate means (period 1976–2005) of the relative humidity (Fig. 2e) and the saturation deficit (Fig. 2f) were used. Additionally, the vegetation cover was classified using the CORINE land cover dataset [24]. It comprises a total of 44 land cover classes derived from satellite imagery. The recent version of the CORINE data was issued 2013 providing a spatial resolution of 3 arc seconds (<100 m). CORINE level I is divided into five classes: *artificial surfaces*, *agricultural areas*, *forest and semi natural areas*, *wetlands* and *water bodies*. For this study, agricultural areas were denoted as A. Further, from the forest and semi natural areas only the CORINE level III subclasses broad-leafed forest (B), coniferous forest (C), and mixed forest (M) were used (Fig. 2c). No tick densities were estimated for urban areas and water bodies. Thus, the n = 25 sampling sites were divided into 4 land cover categories. Ideal for model calibration the land cover classes were uniformly distributed over the sampling sites, which is reasonably met for 7 sampling sites characterised by class A, 7 by class B, 6 by class C, and 5 by class M (Table 2). Note that the largest and enclosed woodland area is located in the mountainous region of the Black Forest (Fig. 2b). All other regions are characterised by a fragmented mixture of agricultural areas, forests and residential areas.

The spatial resolution of the density map was selected to be 30 arc seconds, which corresponds to grid lengths of about 0.9 km (latitude) and 0.6 km (longitude). Climate variables described above have slightly lower resolutions and were disaggregated using the bi-linear interpolation of the R package *raster* [31, 32]. Categorical CORINE data, however, were provided with the higher resolution of 3 arc seconds. Therefore, they were aggregated to the 30 arc sec grid by selecting the most frequent land class within each grid box.

### Predictive statistical model

To predict the nymphal densities (N) in units of nymphs/100 m^{2}, Poisson regression models were applied. The environmental variables height above sea level (H), temperature (T), relative humidity (RH), saturation deficit (SD) and the 4 land cover classes (LC) were selected as explanatory variables.

While the height and the climate variables were integrated in the model as numerical values, the land cover is considered as a categorical variable (classes A, B, C and M). For all other land cover classes the tick density is known to be very low or even zero. They comprise urban areas and water bodies for which no tick density was estimated. McFadden’s pseudo R
^{2}_{p}
(the coefficient of determination R^{2}, adapted for additional categorical variables such as LC) and the root mean square error (RMSE) were selected as goodness-of-fit measures. The final model was evaluated by leave-one-out cross-validation (LOOCV).

All analyses and graphical presentations were conducted using the open-source statistical computing environment R [31]. The gridded data were processed with the R-package raster [32].

## Results

The results of the field study were summarised in Table 2. For each of the 25 sampling locations described in Table 1, the total number of *I. ricinus* nymphs flagged monthly during 1 year is given together with the air temperatures for 2013 and 2014 and the long-term mean of the relative humidity, both extracted from gridded data provided by the German weather service (Deutscher Wetterdienst, DWD). Additionally, the calculated saturation deficit and land cover classifications are listed. Table 2 shows a wide range of values for the flagged nymphs. In 2013 for example, nymphal densities range from N = 0 at the mountain site FB (Feldberg) or N = 3 at AH (Allerheiligen) to N = 1,066 at BT (Botnang). Thus, for an optimal representation nymphal densities were plotted on a logarithmic scale. Solving the Poisson regression model to predict nymphal densities for 2013 and 2014 using the gridded explanatory variables results in the scatterplots depicted in Fig. 3. The Poisson regression model performs slightly better with the explorative variables of the year 2013. Looking at McFadden’s pseudo R
^{2}_{p}
, a total of 71.7 % of the variation in the observed nymphal density was explained by the modelled nymphal density in 2013. In contrast, only 46.1 % of the explained variance was estimated for 2014. The regression coefficients β for both models are summarised in Table 3. While for 2013 all of the explanatory variables contribute highly significantly to the performance of the model, for 2014 only three variables depict a significant contribution. Further, a collinearity of some variables was estimated. Although this collinearity does not influence the explained variances, it leads to regression parameters, which are not generally applicable (see the discussion in the next section).

Poisson regression models with the parameters depicted in Table 2 were used to construct nymphal density maps of Baden-Württemberg for 2013 (Fig. 4) and 2014 (Additional file 1). Each map depicts the spatial distribution of the density of nymphs and should be interpreted as the amount of *I. ricinus* nymphs that may be collected by monthly flagging an area of 100 m^{2}. Green grids represent regions with low nymphal densities (N = 0–50 nymphs/100 m^{2}), while red gradations indicate regions with high nymphal densities (N = 50–1,000 nymphs/100 m^{2}). Low densities of *I. ricinus* are restricted to higher altitudes in the Black Forest characterised by slightly fragmented coniferous forest habitats (Fig. 2c, class C). Moderate densities were estimated for all other hilly countryside with heights around 300–800 m characterised by shorter growing periods and/or forested areas with high proportions of coniferous forest. Very high densities were estimated for the warmest parts at altitudes below 400 m. This includes, in particular, the regions along the river Rhine and Neckar with surrounding areas, as well as the ambience of Lake Constance. The estimates for all other parts with warm to moderate climatic conditions tend to high nymphal densities. Urban areas and water bodies were excluded from the analysis (Fig. 4, yellow and blue areas). Comparing the maps for 2013 and 2014 depicts a similar distribution of the nymphal density, which is generally lower in 2014 (Fig. 5). These lower densities are caused by extraordinary high temperatures in 2014. The difference of the mean temperatures listed in Table 2 is 1.6 °C, which is of the same order as the temperature increase predicted by climatologists for the next 100 years. It seems that higher temperatures at lower altitudes are responsible for the observed decrease of the *I. ricinus* nymphal densities, while an increase of the generally lower temperatures at higher altitudes, e.g. in the black forest, cause an increase of the nymphal densities (Fig. 5). The frequency distribution of the difference between the nymphal densities mapped for 2013 and 2014 are depicted in Fig. 6. According to this frequency distribution about 50 % of the study region depicts only minor changes in the nymphal density (below ± 25 nymphs/100 m^{2}). For the other 50 % of the area of Baden-Württemberg considerably lower nymphal densities were estimated for 2014.

To verify the Poisson regression model leave-one-out cross-validation (LOOCV) was applied [33]. It is a special case of the well-known p-out cross-validation with p = 1 as appropriate for the sampling size of n = 25 (flagging sites). LOOCV calculates all possible combinations of models calibrated with 24 sampling sites, which were validated with the remaining 1 site not used for calibration. Typically, LOOCV results were expressed by verification measures such as the prediction residual error sum of squares (PRESS) or the root mean square error (RMSE). For 2013, a value of RMSE = 227 nymphs/100 m^{2} was calculated providing a more realistic error measure than RMSE = 110 nymphs/100 m^{2} calculated with the whole set of sampling sites (Fig. 3). Due to generally lower observed nymphal densities, LOOCV results in lower errors of RMSE = 104 nymphs/100 m^{2} in 2014.

## Discussion

So far only few landscape scale tick density maps, such as those of the nymphal density of *I. pacificus* in California [6], were published. The *I. ricinus* maps presented here (Fig. 4, Additional file 1) contribute to the development of tick density maps compiled with statistical models. The nymphal densities are to be interpreted as the number of nymphs obtained by monthly flagging of an area of 100 m^{2}. As a consequence of the flagging method only active questing ticks are sampled. Further, the proportion of questing ticks compared to the overall population depends on weather conditions [15]. Therefore, the questing tick densities are naturally lower than the actual densities [34, 35] and variable to a certain extend. This has to be taken into account by interpreting the results of the modelled nymphal densities. However, investigating only forests and forest-like habitats has a positive effect on the comparability of the results, since forest habitats provide comparable microclimatic conditions in contrast to other habitats such as meadows or clearances. Since forests offer very good survival conditions for *I. ricinus* ticks in Baden-Württemberg, modelled nymphal densities have to be interpreted as the upper end of the population density scale. Therefore, lesser population densities are to be expected in other habitats.

Uncertainties appeared in the modelling process, especially in the selection of explanatory variables. As suggested by experts [4], only biologically meaningful variables were considered instead of the frequently applied approach of reducing significant variables from a large set of possible explanatory variables. As a side effect, collinearities (e.g. between altitude and temperature) were estimated, which should be avoided in generally applicable models [4]. Due to the low number of explanatory variables used in the model, this is, however, not possible without reducing the model’s performance. Thus, the biological significance was higher valuated than statistical features and no generally applicable regression parameters were given. Instead, the regression parameters were separately fitted for each year and must be fitted again if further years will be investigated. In so doing, the model is demonstrably usable for mapping, i.e. spatial interpolation, but could not be used to predict future scenarios.

By the choice of the generalized linear model (GLM) it was possible to include the categorical variable land cover. As a result, the performance of the statistical model increased and the tick maps could be calculated with the very high spatial resolution of 0.5 km^{2}. The mean nymphal tick densities of the four land cover (LC) classes are depicted in Fig. 7. For the climatological unremarkable year 2013, a doubling rule may be derived from it. By the transition from one LC class to another, in the order of the classes C-A-M-B, the nymphal density of *I. ricinus* doubles. As expected, the statistically estimated density of *I. ricinus* is minimal in coniferous forests and maximal in broad-leafed forests. For the extraordinary hot year 2014, however, the nymphal densities of the LC classes M and B were not significantly different from class A. Only the nymphal densities in LC class C were significantly lower than in class A as depicted by the p-values of the regression coefficients in Table 2. In succession this might be responsible for the lower explained variance of the Poison regression model applied to the data of 2014. In this context it should be noted that A is the default LC in the Poisson model (as described above). Further, the definition of class A is notable, which contains arable land, permanent crops, pastures and heterogeneous agricultural areas. However, for a final statement concerning nymphal densities and LC longer time series of observations are needed. Although it is well known that *I. ricinus* is present even in city parks [36], no tick densities were estimated for urban areas. Tick habitats in urban areas were treated as sub-scale and not resolved by the maps presented here.

The model performance was expressed by explained variances of R
^{2}_{p}
= 71.7 % for the model fitted with the data of 2013 and R
^{2}_{p}
= 46.1 % for the model of 2014 (Fig. 3), indicating a high reliability of the *I. ricinus* density maps (Figs. 4, 5, Additional file 1). The results should be evaluated with regard to the uncertainties in observed tick densities as well as the low number of explanatory variables used for modelling. As a verification measure independent from observations used to calibrate the models, LOOCV errors of RMSE = 227 nymphs/100 m^{2} for 2013 and RMSE = 104 nymphs/100 m^{2} for 2014 were estimated. These errors are of the order of the mean nymphal density or 20 % of the maximum nymphal density, thus twice of the error estimated during the calibration process (RMSE = 110 nymphs/100 m^{2} for 2013 and RMSE = 61 nymphs/100 m^{2} for 2014). For high nymphal densities these errors are comparable to the observational error, i.e. the accuracy of flagging. It should be noted that it is also possible to fit the model with in situ measurements (Additional file 2), although gridded variables must be used to compile geographical maps.

Since the original research project was designed to estimate the density and distribution of ticks on the one hand, and to estimate the pathogen infestation on the other hand, the previous site selection had a focus on habitats and regions where ticks are likely to be found. Therefore, only a few sites with low and zero tick numbers were selected. This leads to uncertainties in observed and modelled nymphal densities, respectively. The inclusion of additional field data from land cover classes and altitudes with low tick densities would improve future model predictions. Varying abundance of several hosts, i.e. deer as the dominant host species for female ticks and small mammals as hosts for immature stages, are also likely to influence tick densities [34]. Future model improvements are expected by incorporating these data.

## Conclusions

Two years of *I. ricinus* observations from 25 sampling sites in the federal state Baden-Württemberg were used to develop a Poisson regression model. This predictive model was applied to compile very high-resolution maps depicting the spatial distribution of nymphal tick densities for the study area of 35,750 km^{2}. The accuracy of the mapped tick densities was described by explained variances for the model calibration and prediction errors for the estimation of un-sampled sites. For the latter LOOCV was applied, indicating the reliability of the density maps. Although the methodology introduced here was shown to be suitable for compiling tick density maps, it may not be applied to other study regions without deriving new regression coefficients. For the same reason an application to predict climate change scenarios may not be recommended. Although field data were collected monthly (see the seasonal cycles in Additional file 3), only annually accumulated tick densities were considered to reduce the influence of observational errors and overdispersion [37], i.e. high nymphal densities due to local aggregation not resolved at the landscape scale (sub-grid scale effect).

Finally, the study is a first step towards a risk estimation for tick-borne diseases in Central Europe. Following the theory of the basic reproduction number for tick-borne diseases [38], this risk depends—beside other parameters—on the vector density or more exactly on the vector to host ratio. Mapping the risk of tick-borne diseases is therefore a challenge for the future. For other arthropod-borne diseases, however, vector density and disease risk maps are well established as described for Bluetongue [39] or Rift Valley Fever [40].

## References

Petney TN, Pfäffle MP, Skuballa JD. An annotated checklist of the ticks (Acari: Ixodida) of Germany. System Appl Acarol. 2012;17:115–70.

Estrada-Peña A, Venzal JM, Sánchez Acedo C. The tick

*Ixodes ricinus*: distribution and climate preferences in the western Palaearctic. Med Vet Entomol. 2006;20:189–97.Rubel F, Brugger K, Monazahian M, Habedank B, Dautel H, Leverenz S, et al. The first German map of georeferenced ixodid tick locations. Parasit Vectors. 2014;7:477.

Estrada-Peña A, Estrada-Sánchez A, Estrada-Sánchez D. Methodological caveats in the environmental modelling and projections of climate niche for ticks, with examples for

*Ixodes ricinus*(Ixodidae). Vet Parasitol. 2015;208:14–25.Porretta D, Mastrantonio V, Amendolia S, Gaiarsa S, Epis S, Genchi C, et al. Effects of global changes on the climatic niche of the tick

*Ixodes ricinus*inferred by species distribution modelling. Parasit Vectors. 2013;6:217.Eisen RJ, Eisen L, Lane RS. Predicting density of

*Ixodes pacificus*nymphs in dense woodlands in Mendocino county, California, based on geographic information systems and remote sensing versus field-derived data. Am J Trop Med Hyg. 2006;74:632–40.Diuk-Wasser MA, Gatewood AG, Cortinas R, Yaramych-Hamer S, Tsao J, Kitron U, et al. Spatiotemporal patterns of host-seeking

*Ixodes scapularis*nymphs (Acari: Ixodidae) in the United States. J Med Entomol. 2006;43:166–76.Diuk-Wasser MA, Vourch G, Cislo P, Hoen AG, Melton F, Hamer SA, et al. Field and climate-based model for predicting the density of host-seeking nymphal

*Ixodes scapularis*, an important vector of tick-borne disease agents in the eastern United States. Global Ecol Biogeogr. 2010;19:504–14.Daniel M, Zitek K, Danielová V, Kriz B, Valter J, Kott I. Risk assessment and prediction of

*Ixodes ricinus*tick questing activity and human tick-borne encephalitis infection in space and time in the Czech Republic. Int J Med Microbiol. 2006;296(S1):41–7.Schwarz A, Maier WA, Kistemann T, Kampen H. Analysis of the distribution of the tick

*Ixodes ricinus L.*(Acari:Ixodidae) in a nature reserve of western Germany using Geographic Information Systems. Int J Hyg Environ Health. 2009;212:87–96.Bisanzio D, Amore G, Ragagli C, Tomassone L, Bertolotti L, Mannelli A. Temporal variations in the usefulness of Normalized Difference Vegetation Index as a predictor for

*Ixodes ricinus*(Acari: Ixodidae) in a*Borrelia lusitaniae*focus in Tuscany, central Italy. J Med Entomol. 2008;45:547–55.Ruiz-Fons F, Fernández-de-Mera IG, Acevedo P, Gortázar C, de la Fuente J. Factors driving the abundance of

*Ixodes ricinus*ticks and the prevalence of zoonotic*I. ricinus*-borne pathogens in natural foci. Appl Environ Microbiol. 2012;78:2669–76.Pearson RG, Dawson TP. Predicting the impacts of climate change on the distribution of species: are bioclimate envelope models useful? Global Ecol Biogeogr. 2003;12:361–71.

Materna J, Daniel M, Metelka L, Harčarik J. The vertical distribution, density and the development of the tick

*Ixodes ricinus*in mountain areas influenced by climate changes (The Krkonoše Mts., Czech Republic). Int J Med Microbiol. 2008;298(S1):25–37.Schulz M, Mahling M, Pfister K. Abundance and seasonal activity of questing

*Ixodes ricinus*ticks in their natural habitats in southern Germany in 2011. J Vector Ecol. 2014;39:56–65.Perret JL, Guigoz E, Rais O, Gern L. Influence of saturation deficit and temperature on

*Ixodes ricinus*tick questing activity in a Lyme borreliosis-endemic area (Switzerland). Parasitol Res. 2000;86:554–7.Kahl O, Knülle W. Water vapour uptake from subsaturated atmospheres by engorged immature ixodid ticks. Exp Appl Acarol. 1988;4:73–83.

Kottek M, Grieser J, Beck C, Rudolf B, Rubel F. World Map of the Köppen-Geiger climate classification updated. Meteorol Z. 2006;15:259–63.

Rubel F, Kottek M. Observed and projected climate shifts 1901-2100 depicted by world maps of the Köppen-Geiger climate classification. Meteorol Z. 2010;19:135–41.

Arthur DR. British ticks. London: Butterworths; 1963.

Pérez-Eid C. Les tiques. Identification, biologie, importance médicale et vétérinaire. Monographies de microbiologie collection dirigée par Jean-Paul Larpent, Lavoisier. 2007.

Hillyard PD. Ticks of north-west Europe. Shrewsbury: Field Studies Council; 1996.

Murray FW. On the computation of saturation vapour pressure. J Appl Meteorol. 1967;6:203–4.

European Environment Agency: Corine Land Cover 2006 raster data 2013. Version 17, available at http://www.eea.europa.eu/data-and-maps/data/corine-land-cover-2006-raster-3.

Burri C, Cadenas MF, 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 Zoon Dis. 2007;7:50–8.Berger KA, Ginsberg HS, Dugas KD, Hamel LH, Mather TN. Adverse moisture events predict seasonal abundance of Lyme disease vector ticks

*(Ixodes scapularis)*. Parasit Vectors. 2014;7:181.Qviller L, Grøva L, Viljugrein H, Klingen I, Mysterud A. Temporal pattern of questing tick

*Ixodes ricinus*density at differing elevations in the coastal region of western Norway. Parasit Vectors. 2014;7:179.Gray JS, Dautel H, Estrada-Peña A, Kahl O, Lindgren E. Effects of climate change on ticks and tick-borne diseases in Europe. Interdiscipl Persp Inf Dis. 2009 (Article ID 593232).

Climate Data Centre of the German Weather Service. Gridded fields of air temperature (daily mean). 2014.

Frick C, Steiner H, Mazurkiewicz A, Riediger U, Rauthe M, Reich T, et al. Central European high-resolution gridded daily data sets (HYRAS): mean temperature and relative humidity. Meteorol Z. 2014;23:15–32.

R Core Team. R: A language and environment for statistical computing, Vienna, Austria. R Foundation for Statistical Computing 2014, available at http://www.R-project.org/.

Hijmans RJ. raster: Geographic analysis and modeling. R package version 2.3-12, 2014. Available at http://CRAN.R-project.org/package=raster.

Hawkins DM, Basak SC, Mills D. Assessing model fit by cross-validation. J Chem Inf Model. 2003;43:579–86.

Dobson ADM. History and complexity in tick-host dynamics: discrepancies between ‘real’ and ‘visible’ tick populations. Parasit Vectors. 2014;7:231.

Dobson ADM. Ticks in the wrong boxes: assessing error in blanket-drag studies due to occasional sampling. Parasit Vectors. 2013;6:344.

Tappe J, Strube C. Anaplasma phagocytophilum and

*Rickettsia*spp. infections in hard ticks (*Ixodes ricinus*) in the city of Hanover (Germany): Revisited. Ticks and Tick-borne Dis. 2013;4:432–8.Petney TN, van Ark H, Spickett AM. On sampling tick populations: the problem of overdispersion. Onderstepoort J Vet Res. 1990;57:123–7.

Randolph SE, Chemini C, Furlanello C, Genchi C, Hails RS, Hudson PJ, et al. The ecology of tick-borne infections in wildlife reservoirs. In: Hudson PJ et al. editors. The Ecology of Wildlife Diseases, Oxford University Press. 2001.

Brugger K, Rubel F. Bluetongue disease risk assessment based on observed and projected

*Culicoides obsoletus*spp. vector densities. PLoS One. 2013;8(4):e60330.Barker CM, Niu T, Reisen WK, Hartley DM. Data-driven modeling to assess receptivity for Rift Valley fever virus. PLoS Negl Trop Dis. 2013;7(11):e2515.

## Author’s contributions

The project was designed by TP and SN, the field data were collected and determined by DB, MP, NL, PS, RG and RO, the statistical analysis and modelling was designed by FR, the model was implemented and tested by DB, KB, KL, JR, MW and FR, and the paper written by DB, KB, SN, TP, MP, NL, KL and FR. All authors read and approved the final manuscript.

### Acknowledgements

The study was financed by the Ministry of the Environment, Climate Protection and the Energy Sector with their support programme BWPLUS under project numbers BWZ 11001, 11005, 11006 and 11007, and the Graduate School for Climate and Environment (GRACE). We are grateful to Florian Hogewind for technical assistance and to Sebastian Schmidtlein for institutional support.

### Compliance with ethical guidelines

**Competing interests** The authors declare that they have no competing interests.

## Author information

### Authors and Affiliations

### Corresponding author

## Additional files

### Additional file 1:

**Figure S1.**
*Ixodes ricinus* nymphal ticks per 100 m^{2} for 2014. Map of the total number of nymphal ticks monthly flagged during 2014 and interpolated to the entire region of Baden-Württemberg, Germany. Sampling locations are marked by a circle showing both the observed (left half) and the modelled (right half) tick density.

### Additional file 2.

Poisson regression model with in-situ measurements. **Table S1**. Environmental variables observed at the 25 sampling sites. Sites specific environmental variables from in-situ measurements comprising temperature T in °C, relative humidity RH in %, saturation deficit SD in hPa as well as land cover classes A (agricultural land), B (broad-leaved forest), C (coniferous forest) and M (mixed forest) for 2013 and 2014. **Table S2**. Summary of regression models for 2013 and 2014 using in-situ variables. For each explanatory variable the regression coefficient b, the standard error SE, the z-value (test statistics) and the p-value (significance) are given. Note that land cover classifications A, B, C and M are categorical variables set to 0 (false) or 1 (true), from which class A was selected as default (b = 0). **Figure S2**. Observed vs. modelled *Ixodes ricinus* nymphs per 100 m^{2}. Comparison of observed vs. modelled nymphal densities using in-situ measured explanatory variables for 2013 (left) and 2014 (right). The model performance is expressed by explained pseudo variances R
^{2}_{p}
and root mean square errors (RMSE).

### Additional file 3:

**Figure S3.** Time series of monthly Ixodes ricinus nymphal ticks per 100 m^{2}. Selected sites for land classes C, A, M and B.

## Rights and permissions

**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.

## About this article

### Cite this article

Boehnke, D., Brugger, K., Pfäffle, M. *et al.* Estimating *Ixodes ricinus* densities on the landscape scale.
*Int J Health Geogr* **14**, 23 (2015). https://doi.org/10.1186/s12942-015-0015-7

Received:

Accepted:

Published:

DOI: https://doi.org/10.1186/s12942-015-0015-7

### Keywords

- Ixodid ticks
- Generalized linear model
- Population density
- Climate
- Land cover classification