Model of risk of exposure to Lyme borreliosis and tick-borne encephalitis virus- Infected ticks in the border area of the Czech Republic (South Bohemia) and Germany (Lower Bavaria and Upper Palatinate)

International Journal of Environmental Research and Public Health Article Model of Risk of Exposure to Lyme Borreliosis and Tick-Borne Encephalitis Virus-Infected Ticks in the Border Area of the Czech Republic (South Bohemia) and Germany (Lower Bavaria and Upper Palatinate) Václav Hönig 1,2,3,* , Pavel Švec 4, Lukáš Marek 5, Tomáš Mrkviˇcka 6, ZubrikováDana 7,8 , Maria (Vögerl) Wittmann 7, Ondˇrej Masaˇr 4, Daniela Szturcová4, Daniel R ˚užek 1,2,3 , Kurt Pfister 7and Libor Grubhoffer 1,2 1Institute of Parasitology, Biology Centre, Academy of Sciences of Czech Republic, Branisovska 31, 370 05 Ceske Budejovice, Czech Republic; [email protected] (D.R.); libor[email protected] (L.G.) 2Faculty of Science, University of South Bohemia, Branisovska 31, 370 05 Ceske Budejovice, Czech Republic 3Veterinary Research Institute, Hudcova 296, 621 00 Brno, Czech Republic 4Department of Geoinformatics, VSB—Technical University of Ostrava, 17. listopadu 15, 708 00 Ostrava, Czech Republic; [email protected] (P.S.); [email protected] (O.M.); daniela.sztur[email protected] (D.S.) 5GeoHealth Laboratory, Geospatial Research Institute, University of Canterbury, Private Bag 4800, Christchurch 8140, New Zealand; [email protected] 6Faculty of Economics, University of South Bohemia, Studentska 13, 370 05 Ceske Budejovice, Czech Republic; [email protected] 7Institute of Comparative Tropical Medicine and Parasitology, Ludwig-Maximilians-Universität München, Leopoldstr. 5, 80802 Munich, Germany; [email protected] (Z.D.); [email protected] (M.(V.)W.); [email protected] (K.P.) 8Institute of Parasitology, Slovak Academy of Sciences, Hlinkova 3, 040 01 Kosice, Slovakia *Correspondence: [email protected] Received: 22 February 2019; Accepted: 26 March 2019; Published: 2 April 2019   Abstract: In Europe, Lyme borreliosis (LB) and tick-borne encephalitis (TBE) are the two vector-borne diseases with the largest impact on human health. Based on data on the density of host-seeking Ixodes ricinus ticks and pathogen prevalence and using a variety of environmental data, we have created an acarological risk model for a region where both diseases are endemic (Czech Republic—South Bohemia and Germany—Lower Bavaria, Upper Palatinate). The data on tick density were acquired by flagging 50 sampling sites three times in a single season. Prevalence of the causative agents of LB and TBE was determined. Data on environmental variables (e.g., altitude, vegetation cover, NDVI, land surface temperature) were obtained from various sources and processed using geographical information systems. Generalized linear models were used to estimate tick density, probability of tick infection, and density of infected ticks for the whole area. A significantly higher incidence of human TBE cases was recorded in South Bohemia compared to Bavarian regions, which correlated with a lower tick density in Bavaria. However, the differences in pathogen prevalence rates were not significant. The model outputs were made available to the public in the form of risk maps, indicating the distribution of tick-borne disease risk in space. Keywords: tick; Lyme borreliosis; tick-borne encephalitis; Ixodes ricinus; risk modeling; geographical information systems Int. J. Environ. Res. Public Health 2019,16, 1173; doi:10.3390/ijerph16071173 www.mdpi.com/journal/ijerph Int. J. Environ. Res. Public Health 2019,16, 1173 2 of 17 1. Introduction Vector-borne diseases represent a serious health risk to human and animals in many parts of the world. In Europe, pathogens transmitted by Ixodes ricinus (I. ricinus) ticks, namely agents of Lyme borreliosis (LB) and tick-borne encephalitis (TBE), are the most widespread vector-borne diseases of human [ 1 , 2 ]. To be able to take appropriate measures to limit the level of exposure of human populations to vector-borne pathogens, detailed knowledge of the distribution of the risk of infection in space and time is needed. Lyme borreliosis is caused by spirochetes of the Borrelia burgdorferi sensu lato complex (Spirochaetales,Spirochaetaceae) and TBE by the tick-borne encephalitis virus (Flavivirus,Flaviviridae). Both the diseases are zoonoses maintained in geographically more or less specified areas in nature, called “natural foci”. Humans are accidental, predominantly dead-end hosts, which get infected when entering the natural focus. In these foci, causative agents of vector-borne diseases are maintained by circulation between arthropod vectors and transmission competent hosts. Therefore, the epidemiology of these infections is influenced not only by the intrinsic features of the pathogen, but also by the behavior and ecology of its vectors and hosts. Moreover, all three components of vector-borne pathogen natural cycles are influenced by numerous environmental factors of abiotic (temperature, air humidity, access of sunlight, altitude etc.) [ 3 – 6 ] or biotic (natural host availability, transmission competence of the host, vector-host contact rate, host immunity etc.) [ 7 – 9 ] nature. Because the environment influences all three components of the natural circulation, environmental variables may be used for the prediction of risk of infection for humans entering certain areas [10,11]. The existing models for tick-borne pathogens created for ecological or epidemiological purposes are usually based on the estimation of host-seeking tick density and/or pathogen prevalence in the field, collection of different environmental data, and use of various statistical approaches (multiple linear regression, classification and regression trees, generalized linear models etc.) in order to find a mathematical expression of the relationship between the environmental variables and characteristics of tick and pathogen distribution [6,12–17]. The main aim of this study was to use a set of tick abundance and pathogen prevalence data to prepare summarized, easy to access information on the risk of infection by tick-borne diseases for the general public. Based on the results of an extensive field sampling and laboratory analysis of the tick populations [ 18 , 19 ], we have developed a model for the prediction of the density of host-seeking I. ricinus ticks, LB spirochete and TBEV prevalence, and density of LB spirochete/TBEV-infected ticks using commonly available environmental data. Using tools of geographic information systems (GIS), we have integrated the results of the research, including the risk models, in map outputs which have been made available to the public through an interactive internet-based map portal. 2. Materials and Methods 2.1. Study Sites The data on the distribution of I. ricinus tick density, and TBEV and B. burgdorferi s.l. prevalence in ticks were collected in the Czech-German borderland in 2010–2011. Records of retrospective LB and TBE human disease case incidence were obtained for the same area for the period 2001–2009. The study area consisted of the region of South Bohemia in the Czech Republic and the regions of Lower Bavaria and Upper Palatinate of the German Free State of Bavaria (Figure 1). Geographically, climatically, and biologically, the area is considerably heterogeneous ranging from Danube Valley at the German-Austrian border to the highest sections of Šumava (Bohemian Forest) and Bavarian Forest mountain range (280–1450 m a.s.l.). The total area under survey comprised 30,077 km 2 (10,056 km 2 in the Czech Republic and 20,021 km 2 in Germany). The regions of South Bohemia and Bavaria are recognized as areas of endemic occurrence of LB and TBE [20,21]. Int. J. Environ. Res. Public Health 2019,16, 1173 3 of 17 Figure 1. Surveyed area of South Bohemia (Czech Republic), Lower Bavaria and Upper Palatinate (Bavaria, Germany); Lower Bavaria—Niederbayern, Upper Palatinate (Oberpfalz), Bayern—Bavaria; localization of tick sampling sites indicated by a cross sign. The acarological risk of tick-borne disease infection is determined by the prevalence of the pathogen in ticks and abundance of the tick population. Previous studies have shown that both the parameters may be predicted based on their known relationship to certain environmental factors [ 6 , 10 , 22 , 23 ]. Therefore, in the field, in a network of 50 study sites spread throughout the study area, we have collected data on the density of host-seeking I. ricinus,B. burgdorferi s.l., and TBEV prevalence. In parallel, a large dataset of environmental data was established and subsequently used to generate predictions of host-seeking tick density and pathogen prevalence for every pixel of the study area. 2.2. Field Study The distribution of ticks and tick-borne pathogens was estimated in a network of 50 study sites (30 in South Bohemia, 20 in Bavaria—Lower Bavaria and Upper Palatinate) (Figure 1), as described in [ 18 ]. The density of host-seeking I. ricinus ticks was estimated by “flagging” each sampling site (600 m 2 per site and sampling event) three times per year (May, June/July, September/October), resulting in 150 separate entries (for details, please see [ 18 ]). Only nymphal and adult I. ricinus ticks were sampled, as the density of larval ticks is difficult to estimate due to their highly aggregated distribution, they rarely infest humans [24], and they generally have a low pathogen prevalence [25]. Collected ticks were transported to the laboratory and tested for the presence of DNA of B. burgdorferi and TBEV RNA using conventional PCR for borrelia [ 26 , 27 ], reverse transcription followed by conventional PCR for TBEV (samples from South Bohemia) [ 28 ], and reverse transcription followed by real-time PCR (samples from Bavaria) [ 29 ] (for more detailed information see [ 18 ] and [ 19 ]). The data for both countries were analyzed as a single dataset. Altogether, data from 15,150 and 28,862 ticks were available considering B. burgdorferi s.l. and TBEV infection status, respectively. 2.3. Environmental Data Apart from the data on tick and pathogen distribution, we have used various sources to collect epidemiological (numbers of disease cases), demographic (numbers of inhabitants), and environmental data (altitude, slope, exposition, vegetation cover, temperature etc.). The types of data and their respective sources are summarized in Table 1. Satellite images acquired by the MODIS instrument were used as a source of land surface temperature (LST) and vegetation cover information—normalized difference vegetation index (NDVI). From the NASA Land Processes Distributed Active Archive Center, Int. J. Environ. Res. Public Health 2019,16, 1173 4 of 17 the following products were downloaded for the whole area under survey: Surface Temperature & Emissivity 8-Day L3 Global, MOD11A2, and Vegetation Indices 16-Day L3 Global MOD13Q1. Individual measurements were extracted for the sampling localities and specified dates. Table 1. Environmental and other accessory data used for model construction and testing *. Type of Data Data Description Source Czech Rep. Source Bavaria Epidemiological number LB cases disease cases per municipality National Institute of Public Health, Prague RKI, Berlin; local public health authorities number TBE cases Demographic number of inhabitants number of inhabitants per municipality Czech Statistical Office Federal Statistical Office Germany Physical-geographical altitude digital elevation model ArcCR500, ArcDATA, Prague Vektor 500, ATKIS® Basis DLM, Bayerische Vermessungsverwaltung slope exposition Climatic land surface temperature (LST) MODIS Land Surface Temperature, 1000 m raster NASA, LP DAAC NASA, LP DAAC Vegetation cover NDVI MODIS NDVI, 250 m raster NASA, LP DAAC NASA, LP DAAC forest type CORINE landcover 2006, 100 m raster EEA EEA * LB—Lyme borreliosis; TBE—tick-borne encephalitis; NDVI—normalized difference vegetation index; RKI—Robert Koch Institute; NASA, LP DAAC—National Aeronautics and Space Administration, Land Processes Distributed Active Archive Center; EEA—European Environment Agency. Using GIS, the data were extracted and optimized to enter the analysis. With regard to differences in the resolution of the source data, we have worked with an optimized resolution of the whole dataset of a 250 m pixel size. Minimum, maximum, mean seasonal (May–September), and mean annual LST and NDVI values were calculated for individual sampling sites. In the case of physical geographical characteristics and forest type, a buffer area around the collection site was taken into account (pixel of collection site and eight surrounding pixels). The characteristics were expressed as the proportion of the area of a given quality (e.g., covered by deciduous forests). The density of host-seeking ticks was expressed as the number of ticks per 100 m 2 , and prevalence as portion of samples positive for LB spirochetes or TBEV of total samples tested. Epidemiological data were obtained (National Institute of Public Health, Czech Republic; Robert Koch Institute and local public health authorities, Germany) as numbers of human disease cases registered between 2001–2009 in each of the 626 municipalities in South Bohemia and 536 in the selected regions of Bavaria. Data on TBE and LB were available for South Bohemia, where both diseases are reported mandatory. In the case of Lower Bavaria and Upper Palatinate, only data on TBE could be gathered. Census of population and housing data were used as the source of numbers of inhabitants. Disease incidence was expressed as the number of cases per 100,000 inhabitants. 2.4. Model Construction Independent models were created for the prediction of the density of host-seeking ticks (nymphal and adult I. ricinus), and probability of infection of a tick by LB spirochetes and TBEV. Subsequently, the tick density model was combined with either the LB spirochetes or TBEV infection probability model, giving a prediction of the density of infected ticks as a measure of risk of infection for an unvaccinated person entering a particular geographical area. First, the correlation of all potential explanatory variables was examined by the Spearman ranking test, and cross-correlations among the variables were inspected. Significantly cross-correlated explanatory variables were grouped together. Subsequently, the set of potential explanatory variables was tested for correlation with the dependent variables. From each group of cross-correlated explanatory variables, one representative with the highest correlation coefficient with the dependent variable was chosen for further analysis. Generalized linear models were used for host-seeking tick density prediction. With regard to the character of the data, a log-linear model with a Poisson distribution (Poisson regression model) Int. J. Environ. Res. Public Health 2019,16, 1173 5 of 17 was applied. Individual explanatory variables were evaluated based on a Likely-hood type 3 test. The quality of the model was evaluated using McFadden’s pseudo R2. In the case of the prediction of B. burgdorferi s.l. and TBEV infection in ticks, a logit model was applied. The relationship of a bivariate dependent variable (tick infected/tick not infected) with potential explanatory variables of the same dataset as in the case of host-seeking tick density was examined in the model. First, all the candidate explanatory variables were included in the model and subsequently, the variables that were not statistically significantly important were excluded, as is done using a backwards stepwise approach. Individual variables were excluded based on the comparison of Wald statistics and (cross-) correlation coefficient values. BIC and AIC criteria were controlled, indicating the performance of the whole model. The performance of the two resulting infection probability models was evaluated using ROC (Receiver Operating Characteristic) analysis. An ROC curve evaluates the relationship between the sensitivity (true positive fraction) and specificity (false positive fraction) of the model [ 30 ] and it is a suitable measure for the binary classification. The common way of interpretation of the ROC is quantifying the fraction of an area under the ROC curve (AUC). Generally speaking, the higher the AUC, the better the ability of the model to classify. A meaningful model should have an AUC higher than 0.5, which indicates that it is better than a random predictor [31]. 2.5. Statistical and GIS analyses Statistical analyses were performed using STATISTICA ver. 9 (StatSoft CR, Prague, Czech Republic) and R [ 32 ]. A non-parametric Mann-Whitney test was used for a comparison of disease incidence rates and host-seeking tick densities, and a chi square test for a comparison of the prevalence rate (LB) or minimum infection rate (TBEV). Differences with p< 0.05 were considered statistically significant. All GIS analyses were performed and map outputs generated using ArcGIS 10.X (ESRI, Redlands, USA). The project outputs were integrated in a web-based map portal using ArcGIS for Server (ESRI, Redlands, USA). 3. Results 3.1. Disease Incidence In the period of 2001–2009, altogether, 1453 TBE cases were registered in South Bohemia compared to 417 cases reported in Lower Bavaria and Upper Palatinate. When expressed as incidence, 231 TBE cases per/100,000 inhabitants were reported from the South Bohemian region compared to 18 cases from Bavaria. The differences were found to be statistically significant (Mann-Whitney UTest, p< 0.0001). 3.2. Density of Host-Seeking I. Ricinus Ticks and Pathogen Prevalence All the analyses were based on aggregated data from South Bohemia and Bavaria, resulting in 50 study sites, 150 individual sampling events, 15,150 I. ricinus ticks analyzed for the presence of B. burgdorferi s.l., and 28,862 ticks tested for TBEV. The mean density of host-seeking ticks reached 30.1 nymphs/100 m 2 (range 1.3–110.2), and 2.0 adult ticks/100 m 2 (0–10.7). The mean density of all host-seeking ticks and nymphal ticks in South Bohemian sampling sites reached significantly higher values (37.4 ticks/100 m 2 and 34.9 nymphs/100 m 2 ) than the abundance of host-seeking ticks in Bavarian sites (24.5 ticks/100 m 2 and 22.8 nymphs/100 m 2 ) (Mann-Whitney U test, p< 0.001). In the case of adult ticks, the difference was not statistically significant (2.3 adult ticks/100 m 2 in South Bohemia and 1.7 adults/100 m2in Bavaria). The total prevalence of LB spirochetes reached 12.3% (1865 positive samples/15,150 tested), and the minimum infection rate (MIR) (samples tested in pools) of TBEV reached 0.30% (87/28,862). Slightly higher prevalence rates of B. burgdorferi were observed in samples from Bavaria, whereas higher MIRs of TBEV were found in South Bohemian samples. The only statistically significant differences ( χ2 test; p< 0.01) were found between the Borrelia prevalence in adult ticks, which was Int. J. Environ. Res. Public Health 2019,16, 1173 6 of 17 significantly higher in samples from Bavarian localities (Tables 2and 3). Ticks positive for DNA of B. burgdorferi s.l. were found at each of the 50 study sites. At least one TBEV positive tick was sampled at 21 of 30 (70%) localities in South Bohemia and at 7 of 20 (35%) study sites located in Bavaria. Table 2. Prevalence of Borrelia burgdorferi s.l. in field-collected Ixodes ricinus ticks. Tick Stage Nymphs Adults Total Positive/Tested Prevalence Positive/Tested Prevalence Positive/Tested Prevalence S. Bohemia 1239/9809 12.63% 117/1373 8.52% * 1356/11,182 12.13% Bavaria 458/3599 12.73% 51/369 13.82% * 509/3968 12.83% Total 1697/13,408 12.66% 168/1742 9.64% 1865/15,150 12.31% Statistically significant differences (p< 0.05) are indicated by *. Table 3. Minimum infection rate of TBEV in field-collected Ixodes ricinus ticks. Tick Developmental Stage Nymphs Adults Total Positive/Tested MIR Positive/Tested MIR Positive/Tested MIR S. Bohemia 54/18,829 0.29% 10/1228 0.81% 64/20,057 0.32% Bavaria 20/8203 0.24% 3/602 0.50% 23/8805 0.26% Total 74/27,032 0.27% 13/1830 0.71% 87/28,862 0.30% MIR—minimum infection rate. Statistically significant differences (p< 0.05) are indicated by *. 3.3. Host-Seeking Tick Density Model Based on the Spearman ranking test, a set of potential explanatory variables was identified (Table 4) to be included in the initial model. Altitude as a well-known factor influencing tick density was included in the model, despite its poor correlation with the dependent variable in the test. Furthermore, significant differences were found among the different sampling events spring—summer—autumn (ANOVA; p< 0.01) and thus the sampling event was added to the model as a categorical predictor. Table 4. Correlation of potential explanatory variables with the density of host-seeking Ixodes ricinus ticks *. Variable Spearman R p-Level LST date of collection 0.136 0.096 LST seasonal mean 0.323 0.001 LST seasonal minimum −0.174 0.033 LST seasonal maximum 0.231 0.004 Altitude 0.089 0.277 Deciduous forest proportion −0.073 0.373 Coniferous forest proportion 0.016 0.843 Mixed forest proportion 0.207 0.011 Other veg. cover proportion −0.174 0.034 Proportion slope over 5% −0.028 0.732 Proportion northwards 0.053 0.523 Proportion southwards −0.267 0.001 Average altitude −0.008 0.924 NDVI Date of collection −0.055 0.501 NDVI May −0.047 0.569 NDVI June −0.090 0.276 NDVI July −0.226 0.005 Int. J. Environ. Res. Public Health 2019,16, 1173 7 of 17 Table 4. Cont. Variable Spearman R p-Level Proportion southwards −0.267 0.001 Average altitude −0.008 0.924 NDVI Date of collection −0.055 0.501 NDVI May −0.047 0.569 NDVI June −0.090 0.276 NDVI July −0.226 0.005 NDVI August −0.257 0.002 NDVI September −0.235 0.004 NDVI seasonal mean −0.244 0.003 NDVI July-September mean −0.286 0.001 NDVI seasonal maximum −0.126 0.124 NDVI seasonal minimum −0.310 0.001 * LST—land surface temperature, NDVI—normalized difference vegetation index recorded in the indicated month of the year of tick sampling. Seasonal minimum, maximum, and mean values were calculated for the period May–September. Proportions of area of a particular character and average altitude were calculated from a pixel of sampling site and eight surrounding pixels, pixel size—250 m. Variables included in the initial model are indicated in bold. A log-linear model with a Poisson distribution was used for host-seeking tick density prediction. According to the residual analysis, an outlier was identified and excluded from further analysis. All included variables were found to be significant in the model (Likely-hood type 3 test, p< 0.001). The parameters were estimated, resulting in the following equation (Table 5): ln (host-seeking tick density) = 5.973866 + season + 0.019914 LST Seasonal Mean − 0.000770 Altitude − 0.492528 Portion Southwards −1.046671 NDVI August Table 5. Parameter estimates—host-seeking tick density model *. Parameter Coefficient Standard Error Wald Statistics Lower Conf. Limit (95%) Upper Conf. Limit (95%) p Intercept 5.9739 0.1178 50.7313 5.7429 6.2045 <0.001 LST Seasonal Mean 0.0199 0.0039 5.0860 0.0122 0.0276 <0.001 Altitude −0.0008 0.0001 −9.5683 −0.0009 −0.0006 <0.001 Proportion Southwards −0.4925 0.0208 −23.6768 −0.5333 −0.4518 <0.001 NDVI August −1.0467 0.0601 −17.4058 −1.1641 −0.9284 <0.001 Period—Spring 0.3935 0.0144 27.2851 0.3652 0.4218 <0.001 Period—Summer 0.0968 0.0154 6.2837 0.0666 0.1269 <0.001 * LST—land surface temperature, NDVI August—normalized difference vegetation index recorded in August. Seasonal mean LST values were calculated for the time period of May–September. Proportions of area of a particular character (area of south–facing slopes) were calculated from the total area of the pixel of sampling site and eight surrounding pixels, pixel size—250 m. For the only categorical variable, “season”, the following parameters were estimated: spring = 0.393470 , summer = 0.096751, autumn = 0. Thus, three independent models were created for the different parts of the tick season. The McFadden’s pseudo R 2 of the presented model reached the value of 0.203, indicating a model with very good fit. The comparison of observed and predicted tick counts on individual locations together with its confidence interval is depicted in Figure 2. Int. J. Environ. Res. Public Health 2019,16, 1173 8 of 17 Figure 2. Comparison of observed (determined in field) and predicted (estimated by the model of host-seeking tick density) I. ricinus tick counts on individual sampling sites. Grey lines show the observed values, red crosses indicate the value predicted by the model, and red lines represent the confidence interval. 3.4. B. burgdorferi and TBEV Probability of Infection Model The probability of infection of a tick by B. burgdorferi s.l. or TBEV was predicted using a logit model. The parameter estimates and goodness of fit data for the probability of infection of the tick by B. burgdorferi and TBEV are presented in Tables 6and 7, respectively. Season of collection, the only categorical explanatory variable, was not found to be significant for either Borrelia or for the TBEV risk model. The probability of infection of a tick by spirochetes of Lyme borreliosis was expressed by the following equation (Φstanding for the normal distribution function): Probability of B. burgdorferi infection = Φ (0.546025 − 1.064178 NDVI May − 0.153217 Portion Northwards + 0.000483 Altitude −0.051443 LST Seasonal Mean) Table 6. Parameter estimates for the model of probability of B. burgdorferi infection of ticks *. Parameter Coefficient Standard Error Wald Statistics Lower Conf. Limit 95% Upper Conf. Limit 95% p Intercept 0.5460 0.2657 2.0553 0.0243 1.0674 0.040 NDVI May −1.0642 0.167 −6.3716 −1.3926 −0.7328 <0.001 Proportion northwards −0.1532 0.0422 −3.6290 −0.2360 −0.0713 <0.001 Altitude 0.0005 0.0002 3.0032 0.0002 0.0008 0.003 LST seasonal mean −0.0514 0.0084 −6.159 −0.0679 −0.0351 <0.001 * LST—land surface temperature, NDVI May—normalized difference vegetation index recorded in May. Seasonal mean LST value was calculated for May-September. Proportions of area of a particular character (area of north-facing slopes) were calculated from the total area of the pixel of sampling site and eight surrounding pixels area, pixel size—250 m. The probability of infection of a tick by TBEV was expressed by the following equation ( Φ standing for the normal distribution function): Probability of TBEV infection = Φ ( − 2.275259–0.038585 LST seasonal minimum + 0.218823 Mixed forest portion −2.526153 NDVI May + 2.323306 NDVI seasonal maximum). Int. J. Environ. Res. Public Health 2019,16, 1173 9 of 17 Table 7. Parameter estimates for the model of probability of TBEV infection in ticks *. Parameter Coefficient Standard Error Wald Statistics Lower Conf. Limit 95% Upper Conf. Limit 95% p Intercept −2.2753 0.4913 −4.6309 −3.2658 −1.3322 <0.001 LST seasonal minimum −0.0386 0.0151 −2.5633 −0.0682 −0.0091 0.010 Proportion mixed forest 0.2188 0.0858 2.5497 0.0497 0.3884 0.011 NDVI May −2.5262 0.5601 −4.5103 −3.6222 −1.4000 <0.001 NDVI seasonal maximum 2.3233 0.6853 3.3903 0.9201 3.6734 0.001 * LST—land surface temperature, NDVI May—normalized difference vegetation index recorded in May. LST seasonal minimum and maximum values were calculated for May-September. Proportions of area of a particular character (area covered by mixed forests) were calculated from the total area of the pixel of sampling site and eight surrounding pixels, pixel size—250 m. Validation of the classification ability of the previous two models was evaluated using ROC analysis. ROC curves are depicted in Figure 3, together with their AUC. AUC of the model for B. burgdorferi was 0.58, while the model for TBEV infection reached the level of 0.64. Figure 3. Receiver operating characteristics (ROC) of classification ability of B. burgdorferi infection probability model (a) and TBEV infection probability model (b); AUC—area under ROC curve. By joining the model of host-seeking tick density prediction with either the model of probability of B. burgdorferi or TBEV infection, a risk model predicting the density of infected ticks was created. The mathematical models were integrated into the GIS system and risk maps were generated (Figure 4). Comparison of the classification (prediction of general density of ticks without the influence of season) of the South Bohemian and Bavarian regions is presented in Table 8. Furthermore, all the epidemiological data and risk models were made available for the general public via an online map portal (http://gisak.vsb.cz/klistata/index_en.html). Int. J. Environ. Res. Public Health 2019,16, 1173 16 of 17 40. Zeman, P.; Januska, J. Epizootiologic background of dissimilar distribution of human cases of Lyme borreliosis and tick-borne encephalitis in a joint endemic area. Comp. Immunol. Microbiol. Infect. Dis. 1999,22, 247–260. [CrossRef] 41. Danielova, V. Natural foci of tick-borne encephalitis and prerequisites for their existence. Int. J. Med. Microbiol. 2002,291 (Suppl. 33), 183–186. [CrossRef] 42. Norman, R.; Bowers, R.G.; Begon, M.; Hudson, P.J. Persistence of tick-borne virus in the presence of multiple host species: Tick reservoirs and parasite mediated competition. J. Theor. Biol. 1999 ,200, 111–118. [CrossRef] 43. Hudson, P.J.; Rizzoli, A.; Rosà, R.; Chemini, C.; Jones, L.D.; Gould, E.A. Tick-borne encephalitis virus in northern Italy: Molecular analysis, relationships with density and seasonal dynamics of Ixodes ricinus. Med. Vet. Entomol. 2001,15, 304–313. [CrossRef] [PubMed] 44. Ostfeld, R.S.; Glass, G.E.; Keesing, F. Spatial epidemiology: An emerging (or re-emerging) discipline. Trends Ecol. Evol. (Amst.) 2005,20, 328–336. [CrossRef] [PubMed] 45. Stefanoff, P.; Rubikowska, B.; Bratkowski, J.; Ustrnul, Z.; Vanwambeke, S.O.; Rosinska, M. A predictive model has identified tick-borne encephalitis high-risk areas in regions where no cases were reported previously, Poland, 1999–2012. Int. J. Environ. Res. Public Health 2018,15, 677. [CrossRef] 46. Suss, J. Epidemiology and ecology of TBE relevant to the production of effective vaccines. Vaccine 2003 ,21 (Suppl. 1), S19–S35. [CrossRef] 47. Imhoff, M.; Hagedorn, P.; Schulze, Y.; Hellenbrand, W.; Pfeffer, M.; Niedrig, M. Review: Sentinels of tick-borne encephalitis risk. Ticks Tick-Borne Dis. 2015,6, 592–600. [CrossRef] 48. Kitron, U. Risk maps: Transmission and burden of vector-borne diseases. Parasitol. Today (Regul. Ed.) 2000 , 16, 324–325. [CrossRef] 49. Alonso-Carne, J.; Garcia-Martin, A.; Estrada-Pena, A. Modelling the phenological relationships of questing immature Ixodes ricinus (Ixodidae) using temperature and NDVI data. Zoonoses Public Health 2015 ,63, 40–52. [CrossRef] 50. 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. [CrossRef] [PubMed] 51. Estrada-Pena, A. Distribution, abundance, and habitat preferences of Ixodes ricinus (Acari: Ixodidae) in northern Spain. J. Med. Entomol. 2001,38, 361–370. [CrossRef] [PubMed] 52. Rizzoli, A.; Hauffe, H.C.; Tagliapietra, V.; Neteler, M.; Rosà, R. Forest structure and roe deer abundance predict tick-borne encephalitis risk in Italy. PLoS ONE 2009,4, e4336. [CrossRef] 53. Stanek, G.; Reiter, M. The expanding Lyme Borrelia complex–clinical significance of genomic species? Clin. Microbiol. Infect. 2011,17, 487–493. [CrossRef] [PubMed] 54. Rudenko, N.; Golovchenko, M.; Mokrácek, A.; Piskunová, N.; Ruzek, D.; Mallatová, N.; Grubhoffer, L. Detection of Borrelia bissettii in cardiac valve tissue of a patient with endocarditis and aortic valve stenosis in the Czech Republic. J. Clin. Microbiol. 2008,46, 3540–3543. [CrossRef] [PubMed] 55. McFadden, D. Quantitative Methods for Analyzing Travel Behaviour of Individuals: Some Recent Developments. In Behavioural Travel Modelling; Hensher, D., Stopher, P., Eds.; Croom Helm: London, UK, 1978; pp. 279–318. 56. Dobler, G.; Hufert, F.; Pfeffer, M.; Essbauer, S. Tick-Borne Encephalitis: From Microfocus to Human Disease. In Progress in Parasitology; Mehlhorn, H., Ed.; Parasitology Research Monographs; Springer: Berlin/Heidelberg, Germany, 2011; pp. 323–331. ISBN 978-3-642-21395-3. 57. Fabbro, S.D.; Gollino, S.; Zuliani, M.; Nazzi, F. Investigating the relationship between environmental factors and tick abundance in a small, highly heterogeneous region. J. Vector Ecol. 2015 ,40, 107–116. [CrossRef] [PubMed] 58. Ferquel, E.; Garnier, M.; Marie, J.; Bernède-Bauduin, C.; Baranton, G.; Pérez-Eid, C.; Postic, D. Prevalence of Borrelia burgdorferi Sensu Lato and Anaplasmataceae Members in Ixodes ricinus Ticks in Alsace, a Focus of Lyme Borreliosis Endemicity in France. Appl. Environ. Microbiol. 2006,72, 3074–3078. [CrossRef] 59. Morán Cadenas, F.; Rais, O.; Jouda, F.; Douet, V.; Humair, P.-F.; Moret, J.; Gern, L. Phenology of Ixodes ricinus and infection with Borrelia burgdorferi sensu lato along a northand south-facing altitudinal gradient on Chaumont Mountain, Switzerland. J. Med. Entomol. 2007,44, 683–693. [CrossRef] [PubMed] Int. J. Environ. Res. Public Health 2019,16, 1173 17 of 17 60. Perez, D.; Kneubühler, Y.; Rais, O.; Gern, L. Seasonality of Ixodes ricinus ticks on vegetation and on rodents and Borrelia burgdorferi sensu lato genospecies diversity in two Lyme borreliosis–endemic areas in Switzerland. Vector Borne Zoonotic Dis. 2012,12, 633–644. [CrossRef] 61. Perez, G.; Bastian, S.; Agoulon, A.; Bouju, A.; Durand, A.; Faille, F.; Lebert, I.; Rantier, Y.; Plantard, O.; Butet, A. Effect of landscape features on the relationship between Ixodes ricinus ticks and their small mammal hosts. Parasit. Vector 2016,9, 20. [CrossRef] [PubMed] 62. Paul, R.E.L.; Cote, M.; Le Naour, E.; Bonnet, S.I. Environmental factors influencing tick densities over seven years in a French suburban forest. Parasit. Vectors 2016,9, 309. [CrossRef] [PubMed] 63. Süss, J.; Schrader, C.; Abel, U.; Voigt, W.P.; Schosser, R. Annual and seasonal variation of tick-borne encephalitis virus (TBEV) prevalence in ticks in selected hot spot areas in Germany using a nRT-PCR: Results from 1997 and 1998. Zent. Bakteriol. 1999,289, 564–578. [CrossRef] © 2019 by the authors. Licensee MDPI, Basel, Switzerland. This article is an open access article distributed under the terms and conditions of the Creative Commons Attribution (CC BY) license (http://creativecommons.org/licenses/by/4.0/).