Next Article in Journal
A Phenotype Classification of Internet Use Disorder in a Large-Scale High-School Study
Next Article in Special Issue
Using Landscape Analysis to Test Hypotheses about Drivers of Tick Abundance and Infection Prevalence with Borrelia burgdorferi
Previous Article in Journal
Effects of Rural Medical Insurance on Chronically Ill Patients’ Choice of the Same Hospital Again in Rural Northern China
Previous Article in Special Issue
A Predictive Model Has Identified Tick-Borne Encephalitis High-Risk Areas in Regions Where No Cases Were Reported Previously, Poland, 1999–2012
 
 
ijerph-logo
Article Menu

Article Menu

Font Type:
Arial Georgia Verdana
Font Size:
Aa Aa Aa
Line Spacing:
Column Width:
Background:
Article

Effect of Climate and Land Use on the Spatio-Temporal Variability of Tick-Borne Bacteria in Europe

by
Roberto Rosà
1,
Veronica Andreo
1,2,
Valentina Tagliapietra
1,*,
Ivana Baráková
1,3,
Daniele Arnoldi
1,
Heidi Christine Hauffe
1,
Mattia Manica
1,
Fausta Rosso
1,
Lucia Blaňarová
4,
Martin Bona
5,
Marketa Derdáková
3,
Zuzana Hamšíková
3,
Maria Kazimírová
3,
Jasna Kraljik
3,
Elena Kocianová
6,
Lenka Mahríková
3,
Lenka Minichová
6,
Ladislav Mošanský
4,
Mirko Slovák
3,
Michal Stanko
4,
Eva Špitalská
6,
Els Ducheyne
7,
Markus Neteler
8,
Zdenek Hubálek
9,
Ivo Rudolf
9,
Kristyna Venclikova
9,10,
Cornelia Silaghi
11,12,13,
Evelyn Overzier
11,
Robert Farkas
14,
Gábor Földvári
14,
Sándor Hornok
14,
Nóra Takács
14 and
Annapaola Rizzoli
1
add Show full author list remove Hide full author list
1
Department of Biodiversity and Molecular Ecology, Research and Innovation Centre, Fondazione Edmund Mach, 38010 San Michele all’Adige, Italy
2
Department of Earth Observation Science, Faculty of Geo-Information Science and Earth Observation (ITC), University of Twente, 7500 AE Enschede, The Netherlands
3
Institute of Zoology, Slovak Academy of Sciences, 84506 Bratislava, Slovakia
4
Parasitological Institute, Slovak Academy of Sciences, 04001 Košice, Slovakia
5
Department of Anatomy, Pavol Jozef Šafárik University, 04001 Košice, Slovakia
6
Institute of Virology, Biomedical Research Center, Slovak Academy of Sciences, 84505 Bratislava, Slovakia
7
Avia-GIS, Risschotlei 33, 2980 Zoersel, Belgium
8
Mundialis GmbH & Co. KG, 53111 Bonn, Germany
9
Institute of Vertebrate Biology, v.v.i., Academy of Sciences of the Czech Republic, 60365 Brno, Czech Republic
10
Institute of Macromolecular Chemistry CAS, 16206 Prague 6, Czech Republic
11
Comparative Tropical Medicine and Parasitology, Ludwig-Maximilians-Universität, 80802 Munich, Germany
12
Institute of Parasitology, National Centre for Vector Entomology, Vetsuisse-Faculty, University of Zurich, 8057 Zürich, Switzerland
13
Institute of Infectology, Friedrich-Loeffler-Institut, 17493 Greifswald, Germany
14
Department of Parasitology and Zoology, University of Veterinary Medicine, 1078 Budapest, Hungary
*
Author to whom correspondence should be addressed.
Int. J. Environ. Res. Public Health 2018, 15(4), 732; https://doi.org/10.3390/ijerph15040732
Submission received: 27 February 2018 / Revised: 29 March 2018 / Accepted: 10 April 2018 / Published: 12 April 2018

Abstract

:
The incidence of tick-borne diseases caused by Borrelia burgdorferi sensu lato, Anaplasma phagocytophilum and Rickettsia spp. has been rising in Europe in recent decades. Early pre-assessment of acarological hazard still represents a complex challenge. The aim of this study was to model Ixodes ricinus questing nymph density and its infection rate with B. burgdorferi s.l., A. phagocytophilum and Rickettsia spp. in five European countries (Italy, Germany, Czech Republic, Slovakia, Hungary) in various land cover types differing in use and anthropisation (agricultural, urban and natural) with climatic and environmental factors (Normalized Difference Vegetation Index (NDVI), Normalized Difference Water Index (NDWI), Land Surface Temperature (LST) and precipitation). We show that the relative abundance of questing nymphs was significantly associated with climatic conditions, such as higher values of NDVI recorded in the sampling period, while no differences were observed among land use categories. However, the density of infected nymphs (DIN) also depended on the pathogen considered and land use. These results contribute to a better understanding of the variation in acarological hazard for Ixodes ricinus transmitted pathogens in Central Europe and provide the basis for more focused ecological studies aimed at assessing the effect of land use in different sites on tick–host pathogens interaction.

1. Introduction

Human alteration of natural ecosystems is now evident on more of 75% of the Earth’s ice-free land mass as a result of urbanization, agricultural and other land uses, with less than a quarter remaining as intact habitats [1]. Urbanization, in particular, has increased worldwide in recent decades and more than half of the world’s population now lives in urban areas, with the expectation that 66% will live in urban areas by 2050 [2,3].
Human land use and exploitation of natural ecosystems represent one of the major drivers of zoonotic disease emergence by disrupting disease dynamics and cross-species transmission in multi-host, multi-pathogen systems (‘perturbation hypothesis’) and/or by increasing exposure of hosts to novel pathogens (‘novel pathogen pool hypothesis’) [4]. However, understanding the mechanisms by which land use change leads to disease emergence is still rudimentary.
Vector borne diseases, and in particular those transmitted by ticks, appear to be particularly sensitive to land use changes [5] and, therefore, they represent a valuable model to quantify the indirect impact of anthropization and land use changes on human and animal health.
The widespread occurrence of the castor bean tick, Ixodes ricinus in Europe, and the pathogens it transmits, together with the rise of clinical cases of tick-borne diseases (TBDs) in humans and livestock, have made TBDs one of the major One Health issues in recent years [6]. I. ricinus has a broad ecological plasticity and capacity to exploit anthropic landscapes, and its distribution has increased in the last three decades as a result of more favorable biotic and abiotic conditions [5]. Concomitantly, the annual incidences of bacterial diseases such as Lyme borreliosis and rickettsiosis have also increased steadily [7,8]. Several authors have suggested that the prediction of ‘acarological hazard’ in different habitat conditions could be useful to public health authorities for planning and implementing targeted educational and preventive actions for TBDs [9,10]. Hazard is defined in epidemiological risk assessment procedures as “the set of circumstances that could lead to harm” [11], i.e., in this case, the occurrence within natural enzootic cycles of certain tick-borne pathogens (TPBs) with known pathogenic potential. In turn, “risk is the actual exposure of susceptible hosts to TBPs” [11] that takes into account the interaction between infected ticks and humans [11]. The acarological hazard depends on the co-occurrence of TBPs, a competent and active vector (e.g., Ixodes ricinus) and a competent reservoir host enabling the transmission in endemic cycles of the pathogens from infected to uninfected ticks [12].
Predicting acarological hazard is challenging and, despite significant progress in estimating tick occurrence and abundance, the approximation of infection rates in ticks is more complex, since this is related to a series of parameters that are difficult to measure in the field: tick abundance and tick seasonal activity, as well as the contact rate of infected ticks with reservoir hosts (where reservoir capacity varies according to vertebrate species and their immune status) [13].
The seasonal abundance of questing I. ricinus nymphs (considered the most relevant epidemiological developmental stage) and their questing activity pattern are dependent on habitat structure, microclimatic conditions, and the availability of tick-feeding hosts [14,15]. The prediction of spatial and temporal distribution of I. ricinus based on ecological and climatic factors has progressed from studies capturing the short term phenology of the ticks based on ground climate data with simple models [16] to more complex ones based on remote sensing data using correlative [16,17] or modified matrix [18] approaches. Remote sensing (RS) imagery has also proven to be very useful in predicting changes in habitat and climatic conditions at different temporal and spatial scales, and it has been widely used to map the distribution of several disease vectors [19], including I. ricinus [20,21,22].
Here, we attempt to assess the spatio-temporal variability of B. burgdorferi s.l., A. phagocytophilum and spotted fever group (SFG) Rickettsia spp. infection in questing I. ricinus nymphs. Specifically, we tested the associations between the abundance of host-seeking I. ricinus nymphs and the density of infected nymphs (DIN) with a combination of RS parameters such as Land Surface Temperature (LST), Normalized Difference Vegetation Index (NDVI), Normalized Difference Water Index (NDWI), and climatic data obtained from interpolations of meteorological station data, such as accumulated precipitation. To fit our models, we used multi-year field data on the abundance of questing I. ricinus and prevalence of infection with TBPs from five EU countries. Our aim was to provide a first insight on the variables affecting TBDs infection hazard on a wide geographical scale in sites with different land use.

2. Materials and Methods

2.1. Questing Tick Data

Questing ticks were collected monthly during their peak activity in spring (April, May, June) from 2011 to 2013 from 19 sampling sites distributed across Italy, Germany, Slovakia, Czech Republic and Hungary within the framework of the EU FP7 project EDENext (Figure 1).
In each country, urban, agricultural and natural sites were sampled (Table 1), defined according to CORINE land cover classification layers with a spatial resolution of 100 m [23].
Specifically, urban sites belonged to category 1.4.1 i.e., green urban areas (these are areas with vegetation within an urban context, ornamental and recreational character and public accessibility), agricultural sites belonged to category 2.3.1 i.e., pastures (here are included areas with permanent grassland), natural sites belonged to category 3.1.3 i.e., mixed forests (in these forests, shrub and bush understorey is present and neither conifer or deciduous species predominate).
Questing ticks were sampled using the standard dragging method. A 1 m2 white blanket attached to a rod was pulled over the vegetation along three established transects of 100 m2 (i.e., 0.01 ha) per sampling site, covering an area of 300 m2. The blanket was checked for the presence of ticks on both sides every 5 m. Ticks collected were put in vials, separately per transect, and brought to the laboratory for identification and analyses. Although the dragging method preferentially captures I. ricinus nymphs, underestimating the true tick abundance in the area and introducing a bias with regard to instar composition (Mejlon et al., 1997), the goal in this study was to estimate temporal and spatial variation of questing nymph abundance, so larval and adult stages were not included in the analyses. All captured ticks were identified to species and life stage using a stereo-microscope and reference keys [25,26,27,28]. Nymphs were individually placed in vials with 70% ethanol and stored at −20 °C until DNA extraction.
Ixodes ricinus nymph counts were pooled over the sampling periods, and questing nymph density (nymphs/hectare) for each year and site was computed. To assess the relative risk to public health, the estimated DIN per hectare for each pathogen was calculated by multiplying the infection prevalence for each pathogen (%) by the density of questing nymphs for each sampling site.

2.2. Molecular Analyses

Each partner laboratory screening of TBPs was performed by PCR-based methods with previously published primers: in Italy, Borrelia spp., A. phagocytophilum and Rickettsia spp. were detected according to Rijpkema et al. [29], Massung et al. [30], and Reye et al. [31], respectively, with minor modifications (see [32]). In Germany and Slovakia, screening for Borrelia spp. was carried out according to Derdáková et al. [33]. In Germany, Slovakia, Czech Republic and Hungary, detection of A. phagocytophilum and Rickettsia spp. was carried out according to Courtney et al. [34] and Regnery et al. [35], respectively; and described, for Germany in Overzier et al. [36,37] for Slovakia in Blaňarová et al. [38], Svitálková et al. [39], Špitalská et al. [40] and Minichová et al. [41]; for the Czech Republic in Venclíková et al. [42]; and for Hungary in Hornok et al. [43]. The comparability of results by each partner laboratory was guaranteed by an inter-laboratory quality ring test of molecular PCR methods for the detection of I. ricinus transmitted pathogens. Analyses for species/genospecies identification of each pathogen were not carried out by all partner laboratories.

2.3. Climatic and Environmental Data

Environmental predictors were selected based on published evidence of their importance to tick populations [22,44,45,46]. For each sampling site, we obtained climatic and environmental data from RS and interpolated climatic datasets. All environmental data were processed in GRASS GIS 7 [47], and extracted from the spatial database at the locations corresponding to the sampling sites. Sampling transect lengths are below the pixel size for most spatial data, so data averaging over a wider area was not considered appropriate.
Land surface temperature (LST) data were obtained from EuroLST dataset [48]. These data were collated from the Moderate Resolution Imaging Spectroradiometer (MODIS) products MOD11A1 and MYD11A1. The original MODIS LST products were reconstructed (i.e., gap-filled to remove void pixels due to clouds) and downscaled from 1000 to 250 m resolution, a higher resolution Digital Elevation Model (DEM) [48]. For the analyses, daily reconstructed LST data were used to derive monthly and seasonal mean temperature (3-month moving windows) for the current (t) and previous year (t − 1). This 3-month aggregation has already proven to perform better as a predictor when compared with shorter time windows [49]. We also obtained autumnal cooling [50] as the slope of a linear regression between daily LST and DOY (day of year) over three months (August, September and October) [50]. Autumnal cooling is by definition only estimated for the year before each sampling period, while the other variables are estimated both for the sampling year and the year before. Spring warming, on the other hand, was estimated as the slope of a linear regression between daily LST and DOY (day of year) over February, March and April of the current and previous year.
Precipitation was measured as total rainfall and number of days with rainfall per month, and in a 3-month moving window as described above, from the gridded ECA&D dataset (European Climate Assessment & Dataset, Version 13.1 [51]) at an approximately 25 km pixel resolution [52].
We obtained the NDVI from MOD13Q1, and the NDWI from MOD09A1. Both are MODIS composite products (16 days and eight days, respectively) and have a spatial resolution of 500 m. The gaps in NDVI and NDWI time series were filled and outliers removed using a harmonic analysis of each series [53]. These variables were used as proxies for vegetation coverage (NDVI) and for available environmental water (NDWI), which includes surface water as well as vegetation water content. We also aggregated NDVI and NDWI data per month and seasonally (3-month moving window from February to August) for the current year (year of sampling) and the year before samplings using the arithmetic average. January NDVI and NDWI were not included in the analyses due to the presence of snow cover, which can dramatically alter the reliability of satellite acquisition of these parameters [50]. Spring greening and spring wetting were estimated as the slope of a linear regression between NDVI or NDWI, respectively, and DOY over three months (February, March and April) for the sampling and previous year.

2.4. Statistical Analyses

We investigated the relationship between nymph density and DIN with a categorical variable representing the land use category (agricultural, natural, urban) and several environmental and climatic variables (i.e., LST, precipitation, NDVI and NDWI). All statistical analyses were performed using R version 3.4.3 [54].
The code used for statistical analysis is available at [55].
Negative Binomial Generalized Linear Mixed Models (GLMM) (through the glmmTMB package, [56]) were used to investigate the association between nymph density and land use category accounting for environmental and climatic variables. Negative Binomial distribution was chosen to account for overdispersion in the count data. DIN were modeled using Linear Mixed Models (LMM) after log-transforming the DIN (through the lme4 package; [57]). In both analyses, country and year were considered as random effects. All climatic and environmental variables were standardized (i.e., by subtracting their mean and dividing by their standard deviation) before including them in the models. The use of standardized variables was recommended as some explanatory variables, such as total precipitation, had a much larger range of variation compared to the rest of the variables considered.
For all climatic and environmental variables, a preliminary analysis was carried out to ascertain in which period they proved to be the best predictor of both total nymphs counted and DIN. Precisely, for each environmental variable, univariate models with a single covariate (environmental variable) in a specific temporal window were computed in turn. We considered fifteen different temporal windows that were: (i) each of the six months from January to June during the year of collection (as no tick collections were carried out after June); (ii) four wider windows obtained by aggregating a 3-month period from January to June of the same year (i.e., January–February–March, February–March–April, March–April–May, April–May–June); (iii) each quarter of the previous year; and (iv) the whole previous year. Mean values were considered for each variable except for precipitation where the total over each temporal window was computed.
For each variable (LST, precipitation, NDVI and NDWI), the temporal window producing the lowest AIC (Akaike Information Criterion) was selected for inclusion in a subsequent full model. Terms that were not significant for any of the temporal windows were not included in the full model. The Variance Inflation Factor (VIF) was used to test for collinearity among all explanatory variables in the full model, removing variables with VIF > 4 [58]. Following exclusion of collinear and non-significant variables, we developed full models for both nymph density and DIN including the remaining environmental and climatic variables, each measured over its optimum temporal window as selected through preliminary analyses, along with the habitat type. Starting from the initial full model, we carried out a model selection procedure (based on AIC) to find the best model for nymph density and DIN for three different pathogens. Residuals of the best models (lowest AIC) were used to check for model assumptions. For all models, we computed the variance explained (conditional R2) following [59].

3. Results

3.1. Questing Nymph Abundance

In total, 17,832 ticks were collected of which 13,291 were nymphs. The highest value for questing nymph density per country and per year (April–May–June period) was observed in Slovakia (average value = 518.7 nymph/year/300 m2 (3 transects of 100 m2 per site); 95% confidence interval = 352.1–764), while the lowest was recorded in the Czech Republic (79; 57.1–109.3) (Figure 2). The remaining countries showed similar values (Italy: 266.4, 140–506.9; Germany: average = 278.7, 190.8–407.2; Hungary: 222.2, 210.6–234.4). The highest value for questing nymph density per land use per year (April–May–June period) was observed within natural sites (370.9, 256.4–536.5), compared to the urban sites (256.3, 153.5–428) and agricultural ones (208.9, 144.2–302.6) (Figure 2). However, as shown by the following analysis, these differences were not statistically significant.
The results of the Negative Binomial GLMM indicate that questing nymph density was significantly associated with environmental and climatic variables, while no differences were observed among land use categories (Table 2). The best model explaining questing nymph density included a positive effect of average NDVI computed in the period April–May–June (Figure 3), the same period in which sampling was carried out, and a non-significant negative effect of the accumulated precipitation in the last quarter of the previous year (Table 2). The model for questing nymph explained the 32% of the variance (R2 = 0.32).

3.2. Infection Prevalence and DIN

The prevalence of infection with A. phagocytophilum and Rickettsia spp. was assessed for all countries, with values 2.51% (se = 0.32%) and 7.2% (se = 0.65%), respectively. B. burgdorferi s.l. (screened in Italy, Germany and Slovakia) was the pathogen with the highest prevalence of infection, with an overall prevalence of 19.32% (se = 1.42%).
As described within the methods, to assess the risk to public health we carried out the statistical analysis for the DIN, obtained by multiplying the infection prevalence for each pathogen by the density of questing nymphs for each sampling site. LMM models indicate that DIN varied significantly among land use categories for A. phagocytophilum and B. burgdorferi while no differences were observed in DIN for Rickettsia spp. (Table 3, Figure 4). The model for DIN with A. phagocytophilum explained the 59% of the variance (R2 = 0.59), while R2 for models for DIN with B. burgdorferi s.l. and Rickettsia spp. were 0.68 and 0.66 respectively.
Specifically, the natural habitat had the highest DIN values for B. burgdorferi s.l., the urban habitat type showed the highest DIN values for A. phagocytophilum, while, for Rickettsia spp., no differences were observed among habitat types (Figure 4).
Concerning the effect of environmental variables on DIN, a negative effect of the accumulated precipitation during the first three months (Jan–Feb–Mar) of the previous year was observed on DIN for B. burgdorferi s.l. (Table 3). On the other hand, the best model for A. phagocytophilum and Rickettsia spp. included a positive effect of NDVI in the same year. Precisely, NDVI recorded in January was positively related to DIN for A. phagocytophilum and NDVI recorded in March was positively related to DIN for Rickettsia spp. (Table 3).

4. Discussion

Land use change is considered among the most important drivers of TBD emergence. Therefore, in a given area, shifting levels of anthropisation, type of socio-economical activity, and changes in the presence of domestic and wildlife species can also provoke the conditions promoting variation in the acarological hazard for TBPs. Predicting tick population seasonal dynamics together with their pathogen infection rate is useful for an early pre-assessment of the acarological hazard, but the parameters are challenging to estimate. In this study, we used environmental and climatic variables that are known to be related to I. ricinus phenology and can be obtained on a wide geographical scale from RS imagery and different types of land use, to explain the spatio-temporal variability in questing I. ricinus nymphs spring density as well as their pathogen infection rates among different European countries.
NDVI is a controversial predictor of tick abundance: some authors consider its relationship with vegetation water content arguable [46], and do not consider NDVI a suitable predictor of tick abundance [60]. Others confirm the correlation of NDVI with relative humidity in the vegetation layer [61], and therefore, with I. ricinus survival rate and tick abundance and questing behavior [15,22,36,37,62,63,64]. In our study, NDVI recorded in late spring of the sampling year was a very good predictor of questing tick density, supporting the latter hypothesis.
Moreover, it has been suggested that RS of vegetation moisture (e.g., NDWI) may outperform NDVI in modelling the incidence of Lyme borreliosis [46]. In our analyses, however, we did not find any significant association between questing nymph density or DIN with NDWI, with the exception of a marginally significant positive association of March NDWI with Rickettsia spp. DIN in univariate models (results not shown). NDWI is based on canopy reflectance values, which could explain why this parameter fails to predict conditions on the forest floor that are more closely correlated with the occurrence of ticks [62].
NDVI was also positively correlated with DIN of A. phagocytophilum and Rickettsia spp. when measured a few months before or partially overlapping with the tick sampling period, making it a suitable early predictor of annual acarological hazard within a reasonable time frame for these pathogens.
With regards to DIN, we found that the accumulated rainfall in the first trimester of the previous year had a negative effect on the DIN of B. burgdorferi s.l., while no correlation was detected for the other pathogens studied here. The negative effect of rainfall detected for DIN of B. burgdorferi s.l. was found also in the nymphs abundance model. Therefore, one can speculate that the driving effect of rainfall on DIN is probably driven by the effect of rainfall on nymph abundance. This association may be related to the effect of rainfall on tick oviposition rates [65], number of eggs per clutch and/or hatching success rather than on questing activity [66], since heavy rain silts up egg masses [67]. Therefore, some climatic factors during peak months of egg deposition could be critical for subsequent distribution of ticks. The effect of land use category on the DIN was different according to the pathogen considered. DIN of B. burgdorferi s.l. was higher in natural habitats, probably because of the occurrence and abundance of several competent reservoirs for these bacteria, such as rodents and wild birds. Some studies suggest that free-living ungulates, which are not competent reservoir hosts for B. burgdorferi s.l., might act as dilutors, i.e., they reduce the prevalence of this pathogen in ticks [68]. However, it is controversial whether the dilution effect occurs regardless of host species ratio or only at unrealistically high densities of deer [69]. Anyway, in order to fully understand the effect of vertebrate host species composition in driving infection prevalence, detailed data on pathogen species/strain identification would be needed.
Regarding A. phagocytophilum, the highest DIN was detected in urban sites with a particularly high value in Slovakia. In this site, the high roe deer density of ca. 300/1000 ha [70] is believed to support the maintenance of high density of questing nymphs. The overall highest DIN for A. phagocytophilum in urban areas may also be attributed to the common occurrence of hedgehogs: in a previous study, it was observed that 76.1% of Erinaceus roumanicus were infected with this pathogen in a city park of Budapest [71]. The role of urban parks, in particular when their size is relevant and/or their connectivity with natural surrounding areas is possible, is to favor the interaction among vectors, wildlife (including those with large size), domestic animals and humans. For this reason, they are becoming emerging hotspots for vector borne pathogens transmission in Europe [72].
Due to the efficient transovarial transmission of the SFG Rickettsia spp., I. ricinus is thought to be both their main vector and reservoir [73]. DIN with Rickettsia spp. seems to be proportional to density of host-seeking nymphs in the studied countries, with significantly higher values in Germany and Slovakia, while no differences were found between habitat types. The role of vertebrate hosts in Rickettsia transmission remains to be elucidated.

5. Conclusions

In conclusion, our study provides evidence of variation of the acarological hazard for I. ricinus transmitted pathogens in relation to habitat type and climatic condition. These latter, expressed in NDVI, measured in late spring, could be used as a good predictor of questing nymphs density. On the other hand, DIN shows a complex relationship with climatic and land use variables according to the pathogens considered, demonstrating the need to perform more detailed ecological studies aimed at assessing the effect of land use in different habitat types on TBP emergence.

Acknowledgments

The study was funded by the EU grant FP7-261504and is catalogued by the EDENext Steering Committee as EDENext 474 (http://www.edenext.eu); the contents of this publication are the sole responsibility of the authors and do not necessarily reflect the views of the European Commission. As this is a joint paper involving research in five countries, the effort was enormous and many people helped in field and lab work. We would like to acknowledge them all, in particular: Lenka Betášová, Hana Blažejová, Andrea Gomez Chamorro, Margherita Collini, Petra Jedličková, Gábor Majoros, Jan Mendel, Markus Metz, Petra Straková, Sándor Szekeres, Claudia Thiel, and Tim Tiedemann.

Author Contributions

R.R., V.T., I.B., V.A., M.D., E.D., R.F., H.C.H., Z.H., M.N., C.S., M.K., A.R. drafted the paper. V.T., I.B., M.D., D.A., F.S., C.S., M.K., I.R., K.V., S.H., G.F., N.T., Z.H., L.Ma., L.Mo., M.St., E.K., M.Sl., E.Š., L.Mi., L.B., J.K., M.B., E.O. carried out field work and molecular analysis. R.R., V.A. carried out the statistical analysis. M.N., E.D., V.A. collected and analyzed remote sensing data. All authors critically reviewed the manuscript.

Conflicts of Interest

The authors declare no conflict of interest.

References

  1. Ellis, E.C.; Ramankutty, N. Putting People in the Map: Anthropogenic Biomes of the World. Front. Ecol. Environ. 2008, 6, 439–447. [Google Scholar] [CrossRef]
  2. Bradley, C.A.; Altizer, S. Urbanization and the Ecology of Wildlife Diseases. Trends Ecol. Evol. 2007, 22, 95–102. [Google Scholar] [CrossRef] [PubMed]
  3. United Nations, Department of Economic and Social Affairs, Population Division (2014). World Urbanization Prospects: The 2014 Revision, Highlights (ST/ESA/SER.A/352); United Nations: New York, NY, USA, 2014. [Google Scholar]
  4. Allen, T.; Murray, K.A.; Zambrana-Torrelio, C.; Morse, S.S.; Rondinini, C.; Di Marco, M.; Breit, N.; Olival, K.J.; Daszak, P. Global Hotspots and Correlates of Emerging Zoonotic Diseases. Nat. Commun. 2017, 8, 1124. [Google Scholar] [CrossRef] [PubMed]
  5. Medlock, J.M.; Hansford, K.M.; Bormane, A.; Derdakova, M.; Estrada-Peña, A.; George, J.-C.; Golovljova, I.; Jaenson, T.G.T.; Jensen, J.-K.; Jensen, P.M.; et al. Driving Forces for Changes in Geographical Distribution of Ixodes ricinus Ticks in Europe. Parasites Vectors 2013, 6, 1. [Google Scholar] [CrossRef] [PubMed] [Green Version]
  6. Heyman, P.; Cochez, C.; Hofhuis, A.; van der Giessen, J.; Sprong, H.; Porter, S.R.; Losson, B.; Saegerman, C.; Donoso-Mantke, O.; Niedrig, M.; et al. A Clear and Present Danger: Tick-Borne Diseases in Europe. Expert Rev. Anti-Infect. Ther. 2010, 8, 33–50. [Google Scholar] [CrossRef] [PubMed]
  7. Rizzoli, A.; Hauffe, H.; Carpi, G.; Vourc, H.G.; Neteler, M.; Rosa, R. Lyme Borreliosis in Europe. Euro Surveill. 2011, 16. [Google Scholar] [CrossRef]
  8. Brouqui, P.; Parola, P.; Fournier, P.E.; Raoult, D. Spotted Fever Rickettsioses in Southern and Eastern Europe. FEMS Immunol. Med. Microbiol. 2007, 49, 2–12. [Google Scholar] [CrossRef] [PubMed]
  9. Takken, W.; van Vliet, A.J.H.; Verhulst, N.O.; Jacobs, F.H.H.; Gassner, F.; Hartemink, N.; Mulder, S.; Sprong, H. Acarological Risk of Borrelia burgdorferi Sensu Lato Infections across Space and Time in The Netherlands. Vector Borne Zoonotic Dis. 2017, 17, 99–107. [Google Scholar] [CrossRef] [PubMed]
  10. Mannelli, A.; Boggiatto, G.; Grego, E.; Cinco, M.; Murgia, R.; Stefanelli, S.; De Meneghi, D.; Rosati, S. Acarological Risk of Exposure to Agents of Tick-Borne Zoonoses in the First Recognized Italian Focus of Lyme Borreliosis. Epidemiol. Infect. 2003, 131, 1139–1147. [Google Scholar] [CrossRef] [PubMed]
  11. Randolph, S.E.; EDEN-TBD Sub-Project Team. Human Activities Predominate in Determining Changing Incidence of Tick-Borne Encephalitis in Europe. Euro. Surveill. 2010, 15, 24–31. [Google Scholar] [CrossRef] [PubMed]
  12. Cagnacci, F.; Bolzoni, L.; Rosà, R.; Carpi, G.; Hauffe, H.C.; Valent, M.; Tagliapietra, V.; Kazimirova, M.; Koci, J.; Stanko, M.; et al. Effects of Deer Density on Tick Infestation of Rodents and the Hazard of Tick-Borne Encephalitis. I: Empirical Assessment. Int. J. Parasitol. 2012, 42, 365–372. [Google Scholar] [CrossRef] [PubMed]
  13. Ostfeld, R.S.; Levi, T.; Jolles, A.E.; Martin, L.B.; Hosseini, P.R.; Keesing, F. Life History and Demographic Drivers of Reservoir Competence for Three Tick-Borne Zoonotic Pathogens. PLoS ONE 2014, 9, e107387. [Google Scholar] [CrossRef] [PubMed]
  14. Ruiz-Fons, F.; Fernández-de-Mera, I.G.; 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–2676. [Google Scholar] [CrossRef] [PubMed]
  15. Estrada-Peña, A.; Estrada-Sánchez, A.; de la Fuente, J. A Global Set of Fourier-Transformed Remotely Sensed Covariates for the Description of Abiotic Niche in Epidemiological Studies of Tick Vector Species. Parasites Vectors 2014, 7, 302. [Google Scholar] [CrossRef] [PubMed]
  16. Daniel, M.; Vráblík, T.; Valter, J.; Kríz, B.; Danielová, V. The TICKPRO Computer Program for Predicting Ixodes ricinus Host-Seeking Activity and the Warning System Published on Websites. Cent. Eur. J. Public Health 2010, 18, 230–236. [Google Scholar] [PubMed]
  17. Estrada-Peña, A. Geostatistics and Remote Sensing Using NOAA-AVHRR Satellite Imagery as Predictive Tools in Tick Distribution and Habitat Suitability Estimations for Boophilus microplus (Acari: Ixodidae) in South America. National Oceanographic and Atmosphere Administration-Advanced Very High Resolution Radiometer. Vet. Parasitol. 1999, 81, 73–82. [Google Scholar] [PubMed]
  18. Dobson, A.D.M.; Finnie, T.J.R.; Randolph, S.E. A Modified Matrix Model to Describe the Seasonal Population Ecology of the European Tick Ixodes ricinus. J. Appl. Ecol. 2011, 48, 1017–1028. [Google Scholar] [CrossRef]
  19. Palniyandi, M. The Role of Remote Sensing and GIS for Spatial Prediction of Vector-Borne Diseases Transmission: A Systematic Review. J. Vector Borne Dis. 2012, 49, 197–204. [Google Scholar]
  20. Estrada-Peña, A.; Estrada-Sánchez, D. Deconstructing Ixodes ricinus: A Partial Matrix Model Allowing Mapping of Tick Development, Mortality and Activity Rates. Med. Vet. Entomol. 2014, 28, 35–49. [Google Scholar] [CrossRef] [PubMed]
  21. 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. [Google Scholar] [CrossRef] [PubMed]
  22. Estrada-Peña, A.; Alexander, N.; Wint, G.R.W. Perspectives on Modelling the Distribution of Ticks for Large Areas: So Far so Good? Parasites Vectors 2016, 9, 179. [Google Scholar] [CrossRef] [PubMed]
  23. CLC Illustrated Nomenclature Guidelines. Available online: https://land.copernicus.eu/user-corner/technical-library/corine-land-cover-nomenclature-guidelines/html/ (accessed on 27 February 2018).
  24. Global Multi-Resolution Terrain Elevation Data 2010 (GMTED2010). Available online: https://lta.cr.usgs.gov/GMTED2010 (accessed on 27 February 2018).
  25. Filippova, N.A. Ixodid Ticks of the Subfamily Ixodinae; Fauna USSR-Acarina; Nauka Leningrad: Moscow, Russia, 1977; Volume IV. [Google Scholar]
  26. Siuda, K. Kleszcze Polski (Acari: Ixodida) Systematyka i Rozmieszczenie; Polskie Towarzystwo Parazytologiczne: Warszawa, Poland, 1993. [Google Scholar]
  27. Hillyard, P.D. Ticks of North-West Europe; The Dorset Press: Dorchester, UK, 1996. [Google Scholar]
  28. Estrada-Pena, A.; Bouattour, A.; Camicas, J.-L.; Walker, A.R. Ticks of Domestic Animals in the Mediterranean Region. A Guide to Identification of Species; University of Zaragoza: Zaragoza, Spain, 2004. [Google Scholar]
  29. Rijpkema, S.G.; Molkenboer, M.J.; Schouls, L.M.; Jongejan, F.; Schellekens, J.F. Simultaneous Detection and Genotyping of Three Genomic Groups of Borrelia burgdorferi Sensu Lato in Dutch Ixodes ricinus Ticks by Characterization of the Amplified Intergenic Spacer Region between 5S and 23S rRNA Genes. J. Clin. Microbiol. 1995, 33, 3091–3095. [Google Scholar] [PubMed]
  30. Massung, R.F.; Slater, K.; Owens, J.H.; Nicholson, W.L.; Mather, T.N.; Solberg, V.B.; Olson, J.G. Nested PCR Assay for Detection of Granulocytic Ehrlichiae. J. Clin. Microbiol. 1998, 36, 1090–1095. [Google Scholar] [PubMed]
  31. Reye, A.L.; Stegniy, V.; Mishaeva, N.P.; Velhin, S.; Hübschen, J.M.; Ignatyev, G.; Muller, C.P. Prevalence of Tick-Borne Pathogens in Ixodes ricinus and Dermacentor reticulatus Ticks from Different Geographical Locations in Belarus. PLoS ONE 2013, 8, e54476. [Google Scholar] [CrossRef] [PubMed]
  32. Baráková, I.; Derdáková, M.; Selyemová, D.; Chvostáč, M.; Špitalská, E.; Rosso, F.; Collini, M.; Rosà, R.; Tagliapietra, V.; Girardi, M.; et al. Tick-Borne Pathogens and Their Reservoir Hosts in Northern Italy. Ticks Tick-Borne Dis. 2017, 9, 164–170. [Google Scholar] [CrossRef] [PubMed]
  33. Derdakova, M.; Beati, L.; Pet’ko, B.; Stanko, M.; Fish, D. Genetic Variability within Borrelia burgdorferi Sensu Lato Genospecies Established by PCR-Single-Strand Conformation Polymorphism Analysis of the rrfA-rrlB Intergenic Spacer in Ixodes ricinus Ticks from the Czech Republic. Appl. Environ. Microbiol. 2003, 69, 509–516. [Google Scholar] [CrossRef] [PubMed]
  34. Courtney, J.W.; Kostelnik, L.M.; Zeidner, N.S.; Massung, R.F. Multiplex Real-Time PCR for Detection of Anaplasma phagocytophilum and Borrelia burgdorferi. J. Clin. Microbiol. 2004, 42, 3164–3168. [Google Scholar] [CrossRef] [PubMed]
  35. Regnery, R.L.; Spruill, C.L.; Plikaytis, B.D. Genotypic Identification of Rickettsiae and Estimation of Intraspecies Sequence Divergence for Portions of Two Rickettsial Genes. J. Bacteriol. 1991, 173, 1576–1589. [Google Scholar] [CrossRef] [PubMed]
  36. Overzier, E.; Pfister, K.; Thiel, C.; Herb, I.; Mahling, M.; Silaghi, C. Diversity of Babesia and Rickettsia Species in Questing Ixodes ricinus: A Longitudinal Study in Urban, Pasture, and Natural Habitats. Vector Borne Zoonotic Dis. 2013, 13, 559–564. [Google Scholar] [CrossRef] [PubMed]
  37. Overzier, E.; Pfister, K.; Thiel, C.; Herb, I.; Mahling, M.; Silaghi, C. Anaplasma phagocytophilum in Questing Ixodes ricinus Ticks: Comparison of Prevalences and Partial 16S rRNA Gene Variants in Urban, Pasture, and Natural Habitats. Appl. Environ. Microbiol. 2013, 79, 1730–1734. [Google Scholar] [CrossRef] [PubMed]
  38. Blaňarová, L.; Stanko, M.; Carpi, G.; Miklisová, D.; Víchová, B.; Mošanský, L.; Bona, M.; Derdáková, M. Distinct Anaplasma phagocytophilum Genotypes Associated with Ixodes Trianguliceps Ticks and Rodents in Central Europe. Ticks Tick-Borne Dis. 2014, 5, 928–938. [Google Scholar] [CrossRef] [PubMed]
  39. Svitálková, Z.; Haruštiaková, D.; Mahríková, L.; Berthová, L.; Slovák, M.; Kocianová, E.; Kazimírová, M. Anaplasma phagocytophilum Prevalence in Ticks and Rodents in an Urban and Natural Habitat in South-Western Slovakia. Parasites Vectors 2015, 8, 276. [Google Scholar] [CrossRef] [PubMed]
  40. Špitalská, E.; Stanko, M.; Mošanský, L.; Kraljik, J.; Miklisová, D.; Mahríková, L.; Bona, M.; Kazimírová, M. Seasonal Analysis of Rickettsia Species in Ticks in an Agricultural Site of Slovakia. Exp. Appl. Acarol. 2016, 68, 315–324. [Google Scholar] [CrossRef] [PubMed]
  41. Minichová, L.; Hamšíková, Z.; Mahríková, L.; Slovák, M.; Kocianová, E.; Kazimírová, M.; Škultéty, Ľ.; Štefanidesová, K.; Špitalská, E. Molecular Evidence of Rickettsia spp. in Ixodid Ticks and Rodents in Suburban, Natural and Rural Habitats in Slovakia. Parasites Vectors 2017, 10, 158. [Google Scholar] [CrossRef] [PubMed]
  42. Venclíková, K.; Mendel, J.; Betášová, L.; Blažejová, H.; Jedličková, P.; Straková, P.; Hubálek, Z.; Rudolf, I. Neglected Tick-Borne Pathogens in the Czech Republic, 2011–2014. Ticks Tick-Borne Dis. 2016, 7, 107–112. [Google Scholar] [CrossRef] [PubMed]
  43. Hornok, S.; Meli, M.L.; Gönczi, E.; Halász, E.; Takács, N.; Farkas, R.; Hofmann-Lehmann, R. Occurrence of Ticks and Prevalence of Anaplasma phagocytophilum and Borrelia burgdorferi s.l. in Three Types of Urban Biotopes: Forests, Parks and Cemeteries. Ticks Tick-Borne Dis. 2014, 5, 785–789. [Google Scholar] [CrossRef] [PubMed]
  44. Schulze, T.L.; Jordan, R.A. Meteorologically Mediated Diurnal Questing of Ixodes scapularis and Amblyomma americanum (Acari: Ixodidae) Nymphs. J. Med. Entomol. 2003, 40, 395–402. [Google Scholar] [CrossRef] [PubMed]
  45. Perret, J.-L.; Rais, O.; Gern, L. Influence of Climate on the Proportion of Ixodes ricinus Nymphs and Adults Questing in a Tick Population. J. Med. Entomol. 2004, 41, 361–365. [Google Scholar] [CrossRef] [PubMed]
  46. Barrios, J.M.; Verstraeten, W.W.; Maes, P.; Clement, J.; Aerts, J.M.; Farifteh, J.; Lagrou, K.; Van Ranst, M.; Coppin, P. Remotely Sensed Vegetation Moisture as Explanatory Variable of Lyme Borreliosis Incidence. Int. J. Appl. Earth Obs. Geoinf. 2012, 18, 1–12. [Google Scholar] [CrossRef]
  47. Neteler, M.; Hamish Bowman, M.; Landa, M.; Metz, M. GRASS GIS: A Multi-Purpose Open Source GIS. Environ. Model. Softw. 2012, 31, 124–130. [Google Scholar] [CrossRef]
  48. Metz, M.; Rocchini, D.; Neteler, M. Surface Temperatures at the Continental Scale: Tracking Changes with Remote Sensing at Unprecedented Detail. Remote Sens. 2014, 6, 3822–3840. [Google Scholar] [CrossRef]
  49. Rosà, R.; Marini, G.; Bolzoni, L.; Neteler, M.; Metz, M.; Delucchi, L.; Chadwick, E.A.; Balbo, L.; Mosca, A.; Giacobini, M.; et al. Early Warning of West Nile Virus Mosquito Vector: Climate and Land Use Models Successfully Explain Phenology and Abundance of Culex Pipiens Mosquitoes in North-Western Italy. Parasites Vectors 2014, 7, 269. [Google Scholar] [CrossRef] [PubMed]
  50. Randolph, S.E.; Green, R.M.; Peacey, M.F.; Rogers, D.J. Seasonal Synchrony: The Key to Tick-Borne Encephalitis Foci Identified by Satellite Data. Parasitology 2000, 121 (Pt 1), 15–23. [Google Scholar] [CrossRef] [PubMed]
  51. The European Climate Assessment & Dataset. Available online: https://www.knmi.nl/kennis-en-datacentrum/project/the-european-climate-assessment-dataset (accessed on 27 February 2018).
  52. Haylock, M.R.; Hofstra, N.; Klein, A.M.; Klok, E.J.; Jones, P.D.; New, M. A European Daily High-Resolution Gridded Data Set of Surface Temperature and Precipitation for 1950–2006. J. Geophys. Res. 2008, 113. [Google Scholar] [CrossRef]
  53. Roerink, G.J.; Menenti, M.; Verhoef, W. Reconstructing Cloudfree NDVI Composites Using Fourier Analysis of Time Series. Int. J. Remote Sens. 2000, 21, 1911–1917. [Google Scholar] [CrossRef]
  54. R: The R Project for Statistical Computing. Available online: https://www.R-project.org/ (accessed on 12 October 2017).
  55. IJERPH_Rosa_2018. Available online: www.geodati.fmach.it/IJERPH_Rosa_2018 (accessed on 27 February 2018).
  56. Brooks, M.E.; Kristensen, K.; van Benthem, K.J.; Magnusson, A.; Berg, C.W.; Nielsen, A.; Skaug, H.J.; Maechler, M.; Bolker, B.M. Modeling Zero-Inflated Count Data with glmmTMB. bioRxiv 2017. [Google Scholar] [CrossRef]
  57. Bates, D.; Mächler, M.; Bolker, B.; Walker, S. Fitting Linear Mixed-Effects Models Usinglme4. J. Stat. Softw. 2015, 67. [Google Scholar] [CrossRef]
  58. Pan, Y.; Jackson, R.T. Ethnic Difference in the Relationship between Acute Inflammation and Serum Ferritin in US Adult Males. Epidemiol. Infect. 2008, 136, 421–431. [Google Scholar] [CrossRef] [PubMed]
  59. Nakagawa, S.; Schielzeth, H. A general and simple method for obtaining R2 from Generalized Linear Mixed-effects Models. Methods Ecol. Evol. 2013, 4, 133–142. [Google Scholar] [CrossRef]
  60. Trout Fryxell, R.T.; Moore, J.E.; Collins, M.D.; Kwon, Y.; Jean-Philippe, S.R.; Schaeffer, S.M.; Odoi, A.; Kennedy, M.; Houston, A.E. Habitat and Vegetation Variables Are Not Enough When Predicting Tick Populations in the Southeastern United States. PLoS ONE 2015, 10, e0144092. [Google Scholar] [CrossRef] [PubMed]
  61. Benedetti, R.; Rossini, P. On the Use of NDVI Profiles as a Tool for Agricultural Statistics: The Case Study of Wheat Yield Estimate and Forecast in Emilia Romagna. Remote Sens. Environ. 1993, 45, 311–326. [Google Scholar] [CrossRef]
  62. Hubálek, Z.; Halouzka, J.; Juricová, Z. Host-Seeking Activity of Ixodid Ticks in Relation to Weather Variables. J. Vector Ecol. 2003, 28, 159–165. [Google Scholar] [PubMed]
  63. Alonso-Carné, J.; García-Martín, A.; Estrada-Peña, A. Assessing the Statistical Relationships among Water-Derived Climate Variables, Rainfall, and Remotely Sensed Features of Vegetation: Implications for Evaluating the Habitat of Ticks. Exp. Appl. Acarol. 2015, 65, 107–124. [Google Scholar] [CrossRef] [PubMed]
  64. 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–555. [Google Scholar] [CrossRef] [PubMed]
  65. Olwoch, J.M.; Rautenbach, C.D.W.; Erasmus, B.F.N.; Engelbrecht, F.A.; Van Jaarsveld, A.S. Simulating Tick Distributions over Sub-Saharan Africa: The Use of Observed and Simulated Climate Surfaces. J. Biogeogr. 2003, 30, 1221–1232. [Google Scholar] [CrossRef]
  66. Adejinmi, J.O. Effect of Water Flooding on the Oviposition Capacity of Engorged Adult Females and Hatchability of Eggs of Dog Ticks: Rhipicephalus Sanguineus and Haemaphysalis Leachi Leachi. J. Parasitol. Res. 2011, 2011, 824162. [Google Scholar] [CrossRef] [PubMed]
  67. Sutherst, R.W. An Experimental Investigation into the Effects of Flooding on the Ixodid Tick Boophilus microplus (Canestrini). Oecologia 1971, 6, 208–222. [Google Scholar] [CrossRef] [PubMed]
  68. Rosef, O.; Paulauskas, A.; Radzijevskaja, J. Prevalence of Borrelia burgdorferi Sensu Lato and Anaplasma phagocytophilum in Questing Ixodes ricinus Ticks in Relation to the Density of Wild Cervids. Acta Vet. Scand. 2009, 51, 47. [Google Scholar] [CrossRef] [PubMed]
  69. Rosà, R.; Pugliese, A.; Norman, R.; Hudson, P.J. Thresholds for Disease Persistence in Models for Tick-Borne Infections Including Non-Viraemic Transmission, Extended Feeding and Tick Aggregation. J. Theor. Biol. 2003, 224, 359–376. [Google Scholar] [CrossRef]
  70. Kazimírová, M.; Hamšíková, Z.; Kocianová, E.; Marini, G.; Mojšová, M.; Mahríková, L.; Berthová, L.; Slovák, M.; Rosá, R. Relative Density of Host-Seeking Ticks in Different Habitat Types of South-Western Slovakia. Exp. Appl. Acarol. 2016, 69, 205–224. [Google Scholar] [CrossRef] [PubMed]
  71. Földvári, G.; Jahfari, S.; Rigó, K.; Jablonszky, M.; Szekeres, S.; Majoros, G.; Tóth, M.; Molnár, V.; Coipan, E.C.; Sprong, H. Candidatus Neoehrlichia Mikurensis and Anaplasma phagocytophilumin Urban Hedgehogs. Emerg. Infect. Dis. 2014, 20, 496–498. [Google Scholar] [CrossRef] [PubMed]
  72. Rizzoli, A.P.; Silaghi, C.; Obiegala, A.; Rudolf, I.; Hubálek, Z.; Földvári, G.; Plantard, O.; Vayssier-Taussat, M.; Bonnet, S.; Špitalská, E.; et al. Ixodes ricinus and its transmitted pathogens in urban and peri-urban areas in Europe: New hazards and relevance for public health. Front. Public Health 2014, 2, 1–26. [Google Scholar] [CrossRef] [PubMed]
  73. Raoult, D.; Roux, V. Rickettsioses as Paradigms of New or Emerging Infectious Diseases. Clin. Microbiol. Rev. 1997, 10, 694–719. [Google Scholar] [CrossRef] [PubMed]
Figure 1. Map of the 19 ticks sampling sites in Italy, Germany, Czech Republic, Slovakia and Hungary (see Table 1).
Figure 1. Map of the 19 ticks sampling sites in Italy, Germany, Czech Republic, Slovakia and Hungary (see Table 1).
Ijerph 15 00732 g001
Figure 2. Boxplot of observed questing I. ricinus nymphs collected over the year (April–May–June period) in different countries and habitat types; x-axis = country; y-axis = number of collected nymphs.
Figure 2. Boxplot of observed questing I. ricinus nymphs collected over the year (April–May–June period) in different countries and habitat types; x-axis = country; y-axis = number of collected nymphs.
Ijerph 15 00732 g002
Figure 3. Best models for questing nymph density (Negative Binomial Generalized Linear Mixed Model); values on the x-axis represent the mean normalized difference vegetation index (NDVI) over three months (April, May, June); values on the y-axis represent the number of collected nymphs. The solid black line represents the fitted values (highlighting the relationship between NDVI and the typical country-year) computed by considering the accumulated precipitation in the 4th quarter at its mean value. Dashed lines are the 95% confidence intervals for the fitted values. Coloured lines represent the association between the number of collected nymphs and NDVI within each country. Points are observed values.
Figure 3. Best models for questing nymph density (Negative Binomial Generalized Linear Mixed Model); values on the x-axis represent the mean normalized difference vegetation index (NDVI) over three months (April, May, June); values on the y-axis represent the number of collected nymphs. The solid black line represents the fitted values (highlighting the relationship between NDVI and the typical country-year) computed by considering the accumulated precipitation in the 4th quarter at its mean value. Dashed lines are the 95% confidence intervals for the fitted values. Coloured lines represent the association between the number of collected nymphs and NDVI within each country. Points are observed values.
Ijerph 15 00732 g003
Figure 4. Best Linear Mixed Models for density of infected nymphs (DIN) for A. phagocytophilum (left panel), B. burgdorferi s.l. (central panel), Rickettsia spp. (right panel). Filled circles represent the fitted values, and empty circles are the observed values. Solid lines represent the 95% confidence intervals for the fitted values.
Figure 4. Best Linear Mixed Models for density of infected nymphs (DIN) for A. phagocytophilum (left panel), B. burgdorferi s.l. (central panel), Rickettsia spp. (right panel). Filled circles represent the fitted values, and empty circles are the observed values. Solid lines represent the 95% confidence intervals for the fitted values.
Ijerph 15 00732 g004
Table 1. Description of study sites (see also Figure 1). Elevations were taken from Global Multi-Resolution Terrain Elevation Data 2010 (mn30_grd layer [24]).
Table 1. Description of study sites (see also Figure 1). Elevations were taken from Global Multi-Resolution Terrain Elevation Data 2010 (mn30_grd layer [24]).
CountrySite NumberSampling SiteLand Use CategoryAltitude (m a.s.l.)LatitudeLongitude
Italy1LamarNatural78446.12872611.058944
Italy2CavedineAgricultural71745.98540210.963142
Italy3PietramurataNatural46846.01325810.927981
Italy4TrentoUrban28546.03518711.139236
Germany5TussenhausenNatural64048.11827910.589147
Germany6KerschlachAgricultural72447.91714211.212342
Germany7Englischer GartenUrban51448.15048111.590053
Germany8Berg StarnbergUrban65948.11011710.575944
Germany9Nymphenburger SchlossparkUrban52248.16081411.492586
Germany10Dörnbergpark RegensburgUrban34549.01547812.085803
Czech Republic11PohanskoNatural16248.72713316.902319
Czech Republic12ValticeUrban21548.73491116.753142
Czech Republic13SuchovAgricultural42648.89744217.581928
Slovakia14BratislavaUrban18448.16666717.066667
Slovakia15FúgeľkaNatural38648.36666717.300000
Slovakia16RozhanovceAgricultural28048.75000021.366667
Hungary17PilisszentkeresztNatural46847.70083318.884722
Hungary18CsabrendekAgricultural15947.05388917.323333
Hungary19BudapestUrban10547.55027819.052778
Table 2. Best model for questing nymph density (Negative Binomial Generalized Linear Mixed Model). The columns report the estimated coefficients for explanatory variables, their standard errors, z-values (estimate to standard error ratio) and p-value for z-statistic. Independent variables have been standardized. NDVI = normalized difference vegetation index.
Table 2. Best model for questing nymph density (Negative Binomial Generalized Linear Mixed Model). The columns report the estimated coefficients for explanatory variables, their standard errors, z-values (estimate to standard error ratio) and p-value for z-statistic. Independent variables have been standardized. NDVI = normalized difference vegetation index.
Explanatory VariableEstimateStd. Errorz-ValuePr(>|z|)
Intercept5.4570.26021.026<0.001 ***
NDVI (Apr–May–Jun)0.2640.1242.1310.033 *
Accumulated precipitation (Oct–Nov–Dec, previous year)−0.2650.149−1.7790.075
Signif. codes: *** < 0.001; * < 0.05.
Table 3. Best parsimonious models for density of infected nymphs (DIN) (Linear Mixed Model) with A. phagocytophilum, B. burgdorferi s.l. and Rickettsia spp. The columns report the estimated coefficients for explanatory variables, their standard errors, t-values (estimate to standard error ratio) and p-value for the t-statistic. Reference level is Agricultural for habitat type.
Table 3. Best parsimonious models for density of infected nymphs (DIN) (Linear Mixed Model) with A. phagocytophilum, B. burgdorferi s.l. and Rickettsia spp. The columns report the estimated coefficients for explanatory variables, their standard errors, t-values (estimate to standard error ratio) and p-value for the t-statistic. Reference level is Agricultural for habitat type.
ModelExplanatory VariableEstimateStd. Errort ValuePr(>|t|)
DIN for Borrelia Intercept6.2450.14842.216<0.001 ***
Habitat type Natural0.7790.2091.7240.002 **
Habitat type Urban−0.2150.245−0.8780.395
Accumulated Precipitation (Jan-Feb-Mar, previous year)−0.3100.095−3.2680.006 **
DIN for AnaplasmaIntercept2.4780.6154.0260.006 **
Habitat type Natural0.9800.5401.8140.080
Habitat type Urban1.8820.5913.1850.003 **
NDVI (January)0.7400.2612.8310.008 **
DIN for RickettsiaIntercept5.2140.4054.267<0.001 ***
NDVI (March)0.4880.15329.6030.003 **
Signif. codes: *** <0.001; ** <0.01.

Share and Cite

MDPI and ACS Style

Rosà, R.; Andreo, V.; Tagliapietra, V.; Baráková, I.; Arnoldi, D.; Hauffe, H.C.; Manica, M.; Rosso, F.; Blaňarová, L.; Bona, M.; et al. Effect of Climate and Land Use on the Spatio-Temporal Variability of Tick-Borne Bacteria in Europe. Int. J. Environ. Res. Public Health 2018, 15, 732. https://doi.org/10.3390/ijerph15040732

AMA Style

Rosà R, Andreo V, Tagliapietra V, Baráková I, Arnoldi D, Hauffe HC, Manica M, Rosso F, Blaňarová L, Bona M, et al. Effect of Climate and Land Use on the Spatio-Temporal Variability of Tick-Borne Bacteria in Europe. International Journal of Environmental Research and Public Health. 2018; 15(4):732. https://doi.org/10.3390/ijerph15040732

Chicago/Turabian Style

Rosà, Roberto, Veronica Andreo, Valentina Tagliapietra, Ivana Baráková, Daniele Arnoldi, Heidi Christine Hauffe, Mattia Manica, Fausta Rosso, Lucia Blaňarová, Martin Bona, and et al. 2018. "Effect of Climate and Land Use on the Spatio-Temporal Variability of Tick-Borne Bacteria in Europe" International Journal of Environmental Research and Public Health 15, no. 4: 732. https://doi.org/10.3390/ijerph15040732

APA Style

Rosà, R., Andreo, V., Tagliapietra, V., Baráková, I., Arnoldi, D., Hauffe, H. C., Manica, M., Rosso, F., Blaňarová, L., Bona, M., Derdáková, M., Hamšíková, Z., Kazimírová, M., Kraljik, J., Kocianová, E., Mahríková, L., Minichová, L., Mošanský, L., Slovák, M., ... Rizzoli, A. (2018). Effect of Climate and Land Use on the Spatio-Temporal Variability of Tick-Borne Bacteria in Europe. International Journal of Environmental Research and Public Health, 15(4), 732. https://doi.org/10.3390/ijerph15040732

Note that from the first issue of 2016, this journal uses article numbers instead of page numbers. See further details here.

Article Metrics

Back to TopTop