Modelling habitat suitability for occurrence of human tick-borne encephalitis (TBE) cases in Finland

The numbers of reported human tick-borne encephalitis (TBE) cases in Europe have increased in several endemic regions (including Finland) in recent decades, indicative of an increasing threat to public health. As such, it is important to identify the regions at risk and the most influential factors associated with TBE distributions, particularly in understudied regions. This study aimed to identify the risk areas of TBE transmission in two different datasets based on human TBE disease cases from 2007 to 2011 (n=86) and 2012–2017 (n =244). We also examined which factors best explain the presence of human TBE cases. We used ensemble modelling to determine the relationship of TBE occurrence with environmental, ecological, and anthropogenic factors in Finland. Geospatial data including these variables were acquired from several open data sources and satellite and aerial imagery and, were processed in GIS software. Biomod2, an ensemble platform designed for species distribution modelling, was used to generate ensemble models in R. The proportion of built-up areas, field, forest, and snow-covered land in November, people working in the primary sector, human population density, mean precipitation in April and July, and densities of European hares, white-tailed deer, and raccoon dogs best estimated distribution of human TBE disease cases in the two datasets. Random forest and generalized boosted regression models performed with a very good to excellent predictive power (ROC=0.89–0.96) in both time periods. Based on the predictive maps, high-risk areas for TBE transmission were located in the coastal regions in Southern and Western Finland (including the Åland Islands), several municipalities in Central and Eastern Finland, and coastal municipalities in Southern Lapland. To explore potential changes in TBE distributions in future climate, we used bioclimatic factors with current and future climate forecast data to reveal possible future hotspot areas. Based on the future forecasts, a slightly wider geographical extent of TBE risk was introduced in the Åland Islands and Southern, Western and Northern Finland, even though the risk itself was not increased. Our results are the first steps towards TBE-risk area mapping in current and future climate in Finland.


Introduction
High-latitude regions (> 60°) in the Northern Hemisphere are undergoing rapid changes associated with climate warming. Climate change interacts with the global change through atmospheric circulation and biogeophysical and biogeochemical effects; such changes are particularly affected by rising temperatures at high latitudes (Groisman and Soja, 2007). As a result of these climatic changes, distribution of invasive species is expanding to new regions and environmental conditions have become more suitable for the circulation of several vectorborne viruses (Koch et al., 2017;Hobbs, 2000;Jore et al., 2014;Sutherst et al., 2000). Ensemble modelling and geographical https://doi.org/10.1016/j.ttbdis.2020.101457 Received 2 December 2019; Received in revised form 22 April 2020; Accepted 24 April 2020 information systems (GIS) have been used to understand the connections between vectors, their habitats, and vector-borne diseases (Honig et al., 2019;Stefanoff et al., 2018;Gama et al., 2017;Sun et al., 2017;Deka and Morshed, 2018). The distribution of vector-borne diseases remains understudied in many regions. GIS and ensemble modelling approaches used for identification of influential environmental factors and estimation of vector-borne disease risks can improve knowledge on disease prevention.
Tick-borne encephalitis (TBE) is a zoonotic disease caused by the TBE virus (TBEV). TBEV infection typically induces an influenza-like illness. In one third of cases, the initial illness may be followed by fever, meningitis, or meningoencephalitis. Neurological sequelae, including paresis, may occur. Death occurs in 1-2% of cases (World Health Organization (WHO, 2019). TBE occurs focally in endemic areas across large regions of the temperate and boreal forest regions of Europe and Asia (European Centre for Disease Prevention and Control (ECDC, 2015;Charrel et al., 2004;Woolhouse et al., 2001). The annual number of TBE cases in the European Union and the European Free Trade Association varies between 2000 and 3500 cases; most cases occur between June and September (Beauté et al., 2018; European Centre for Disease Prevention and Control (ECDC, 2015). The Finnish Government has defined TBE as one of the major risks to public health in changing climate and has called for preparedness and risk assessment (Tuomenvirta et al., 2018).
Although ticks are both the vectors and the main reservoir for the virus, vertebrate species can also be infected by TBEV. Small mammals can harbor the virus and transmit it to their offspring. TBEV can be transmitted directly from an infected vertebrate host to its progeny or occasionally transovarially (Alekseev and Chunikhin, 1990). Transmission can also occur by co-feeding, when uninfected ticks feed simultaneously close to infected ticks (Labuda et al., 1993;Randolph et al., 1996). TBEV is transmitted by Ixodes spp. ticks from Western Europe through Russia to the Far East (Lindquist and Vapalahti, 2008). Finland lies in the zone where the two tick species (Ixodes ricinus and I. persulcatus) overlap and they can transmit both the European and Siberian virus subtypes (Süss, 2011;Tonteri et al., 2015;Öhman, 1961;Jääskeläinen et al., 2006Jääskeläinen et al., , 2010Kuivanen et al., 2018). Ixodes ricinus is predominant in Southern Finland while I. persulcatus prevails in Northern Finland (Laaksonen et al., 2017). In Finland and in the neighboring countries Sweden and Russia, the host tick species and the distribution of TBEV have moved to more northern latitudes (Tokarevich et al., 2011;Jaenson et al., 2012;Laaksonen et al., 2017).
The natural cycles of TBEV are sensitive to various environmental and ecological factors, such as climate (Daniel et al., 2015;Semenza and Menne, 2009;Brabec et al., 2017), microclimate (Haider et al., 2017;Randolph et al., 2001), and population density of tick hosts (Heyman et al., 2010;Brugger et al., 2017). Changes in weather conditions influence the distribution of vector ticks and host and reservoir animals and therefore affects pathogen transmission and incidence of human disease cases (Alkishe et al., 2017;Lindgren, 1998a). The likelihood of tick co-feeding (and consequently co-feeding transmission of TBEV), tick host densities, and human-tick encounters are all influenced by climatic associations. Seasonal tick activity starts when the daily temperature rises above 5°C but ceases in temperatures above 25°C (Daniel et al., 2015(Daniel et al., , 2018. Daylight length and relative humidity (RH), particularly long periods of low RH, affect tick activity (Gray et al., 1998;Daniel et al., 2015). In addition, RH also affects tick survival (Gray et al., 1998). The presence of suitable host animals also influences tick abundance. Adult I. ricinus and I. persulcatus normally feed on medium-sized and large-sized animals such as deer and hares, while nymphal ticks mainly feed on small to medium-sized animals such as rodents, birds and hares, but also on large-sized animals (Gray et al., 2016;Brugger et al., 2017). TBE has been found to correlate with the abundance of deer and hares in Sweden (Jaenson et al., 2018), and in Finland with white-tailed deer density (Dub et al., unpublished results).
Habitat suitability refers to condition, when the combination of abiotic environmental variables at the site is included in the environmental conditions that a species needs to grow and maintain viable populations (Hutchinson, 1992). This is the basis for habitat suitability modelling approach which can improve knowledge on the distribution of pathogens, vector-borne diseases, vectors and hosts in the changing climate (Guisan and Zimmermann, 2000;Guisan et al., 2017). In this study, we use the definition of habitat suitability to refer to the ability of a habitat to promote the occurrence of TBE, and to the probability of TBE occurrence. Recently, ensemble modelling has been widely used to predict distribution of vectors (Miller et al., 2018;Chalghaf et al., 2018;Uusitalo et al., 2019) and vector-borne diseases (Deribe et al., 2018;Deka and Morshed, 2018;Gama et al., 2017;Eneanya et al., 2018). This tool also has the advantage of producing more robust decision making in the face of uncertainty in comparison with single-model forecasts (Araújo and New, 2007). Modelling studies of vector-borne diseases and their associations with climate change have been conducted since the 1990s (Rogers and Packer, 1993;Haines et al., 2000). The Southern Finland coastline has been predicted to be suitable for TBEV circulation a few decades ago (Randolph and Rogers, 2001), and several distinct foci have since been found in the region (Jääskeläinen et al., 2016;Smura et al., 2019). However, recent spatial modelling studies on TBE in Finland are missing. The aims of this study are to predict the distribution of TBE in the identified foci (Tonteri et al., 2015) in Finland based on known factors that affect TBE incidence; to discover which environmental, anthropogenic, and ecological factors best explain the probability of TBE disease case occurrence; and to analyze the effect of predicted climate change on TBE distribution based on bioclimatic variables. To reach these aims, we chose to apply biomod2, an ensemble platform for species distribution modelling (Thuiller et al., 2016).

TBE occurrence data and study datasets
Finland (59°50′N, 20°38′E, 70°09′N, 31°30′E) is located in Northern Europe between Sweden and Russia ( Fig. 1.). TBE occurrence data include serologically confirmed human TBE cases by municipality from 2007 to 2017 obtained from the National Infectious Diseases Register (NIDR). In total, 488 TBE cases were reported in 2007-2017. Patients infected abroad and whose location of exposure or date of onset were unknown, were excluded from this study (n = 158) ( Fig. 1.).
Overall, the number of TBE cases has increased annually and the distribution of exposures has spread to new areas from 2007 to 2017 (Figs. 1 and 2). The geographical distribution of TBE cases is mainly focused in coastal and Southern Finland (including the Åland Islands). In the past 6 years, the distribution of TBE cases has also expanded to Central, Eastern, and Northern Finland, excluding northernmost Finland. The annual number of TBE cases has almost quadrupled in 11 years from 2007 (n = 20) to 2017 (n = 72).
TBE datasets were built based on NIDR data. Data were split into two different datasets (2007-2011 and 2012-2017), to identify differences in predictors and TBE risk between the two time periods. The first dataset includes an area of 35,344 km 2 and consists of 24 municipalities with 86 TBE cases and 24 control municipalities without TBE cases (Fig. 3.). The second dataset includes an area of 83,720 km 2 , and consists of 51 (presence) municipalities with 244 TBE cases and 51 control (absence) municipalities (Fig. 3). Absence data were randomly selected from municipalities located in the vicinity of TBE municipalities.

Environmental, ecological and anthropogenic data
The selection of explanatory data were based on those introduced in the existing literature. Environmental data for Finland were obtained directly from satellite imagery, or were derived from the satellite imagery in ESRI ArcGIS (version 10.3.1). Ecological data included tick data and hunting data of selected game animals. Anthropogenic data consisted of data on human population density and people working in primary production. Altogether, explanatory data included 50 predictors (excluding separate variables for each month e.g. temperature) before modelling analysis (Table 1.; Table A.1).

Data analysis and habitat suitability modelling
All the geospatial datasets including environmental and other attributes, were set to the same spatial extent, geographic coordinate system (EUREF FIN TM35FIN), and resolutions (1000 m × 1000 m). Mean values of explanatory data per each municipality were calculated. The compiled dataset consisted of the presence-absence data of TBE and environmental, ecological, and anthropogenic data. Prior to modelling analysis, multicollinearity of the variables was tested using Variance Inflation Factors (VIFs) in R (Besley et al., 1980). The VIFs of the suite of environmental variables were calculated and correlated variables were excluded in a stepwise procedure at a generally accepted threshold value of 5 (Dormann et al., 2013;Ringle et al., 2015). Final explanatory variables used in the modelling consisted of 14 variables in the first dataset and 13 variables in the second dataset (Table 2.).
To model the habitat suitability of TBE in Finland, we used ensemble modelling in the biomod2 platform (version 3.3-7) in RStudio computing software (version 1.2.5033; RStudio Inc.). We fitted our data using the following eight predictive modelling techniques: generalized linear models (GLM), generalized additive models (GAM), classification tree analysis (CTA), artificial neural networks (ANNs), multivariate adaptive regression splines (MARS), generalized boosting models (GBM), random forest (RF), and maximum entropy (MAXENT). Flexible discriminant analysis (FDA) and surface range envelope (SRE) were excluded due to methodological weaknesses (Zhao and Gao, 2015;Elith et al., 2006;Pearson et al., 2006). Models were processed mainly using the default settings of the biomod2 with the following exception: we used the function bam mgcv package in GAM model in R due to its advantages regarding use of random term smoothers, specification of weights, offset, and ability to handle large datasets with good estimates (Wood, 2017).
We performed an iteration of 10 model runs for each of the eight algorithms, (80 model runs total). Model calibration was performed using randomly divided model training (70%) set, and model evaluation was performed using the remaining 30% of the data over the 10 model replicate runs (Guisan and Zimmermann, 2000). The area under the curve of a receiver operating characteristic (ROC) value (Fielding and Bell, 1997) and true skill statistics (TSS) (Peirce, 1884) were produced based on each model to estimate the predictive power of the model. Ensemble model was generated using the best-performing model algorithms (0.7 < ROC > 1.0) (Drew et al., 2010). Sensitivity (the proportion of observed presences) and specificity (the proportion of observed absences) were calculated to quantify the omission errors (Fielding and Bell, 1997). The contribution of variables to the model was calculated by correlating the fitted values of the full models with those from the model in which the values of the predictor have been randomly permuted. The relative importance of variables was averaged over the eight model algorithms and compared with the variable contributions in the models to define the most powerful variables and their relative magnitude. All habitat suitability maps were created using ArcGIS software. Spatial autocorrelation (SAC) of the predictor variables was measured using Moran's Index (Moran, 1950).

Ensemble distribution model performance
The model performance of eight model algorithms with mean ROC and TSS values of 10 model runs in both datasets is presented in Fig. 4. Altogether, six out of the eight models provided reliable estimates for TBE distributions resulting in area under the ROC value of 0.74 at minimum (ROC = 0.74-0.91) in the first dataset. In the second dataset, seven out of the eight model algorithms resulted in an ROC value of 0.76 at minimum, suggesting good to excellent predictive power. RF and GBM models were the best performing model algorithms (ROC = 0.89-0.96) in both datasets.
Ensemble model performance in both datasets was classified as excellent, based on the ensemble median TSS value (0.988) in the first dataset and (0.956) in the second dataset. The ensemble model correctly predicted 98.8% of TBE presences (i.e. sensitivity) and 100% of TBE absences (i.e. specificity) in the first dataset. In the second dataset, sensitivity was 97.5% and specificity was 98.0%. SAC of the predictor variables for both datasets indicated that several covariates were highly spatially autocorrelated (0.5 ≤ Moran's I ≥ 0.9) for short distances (p < 0.05) but not at longer distances. (Figs. A.1. and A.2.).

The relative importance of predictors
The relative contribution of influential variables (%) in both datasets showed that several factors resulted in values close to or greater than the mean importance value; these are thus important predictors (Fig. 5.). Proportion of built-up areas, forest, and people working in the primary sector, mean precipitation in April and July, and white-tailed deer and European hare density were the most important predictors (19-41%) in the first dataset. In the second dataset, the most important predictors (12-27%) were proportion of built-up areas, field, human population density, people working in the primary sector, snow-covered land in November, and raccoon dog density.

TBE risk based on suitability maps
The produced habitat suitability map of the first dataset suggests that the majority of coastal municipalities, municipalities in the Åland Islands, and several municipalities in Central and Southwestern Finland are suitable for TBE transmission (Fig. 6.). The habitat suitability map of TBE in the second dataset, strengthens the suggestion that coastal municipalities widely reflect the geographic distribution of the actual TBE risk and already notified TBE cases (Fig. 1.). Furthermore, the habitat suitability map from 2012 to 2017 period suggests that there is a considerable TBE transmission risk in several municipalities in Central Finland and in Eastern municipalities close to the Russian border. Inland municipalities in Northern Finland are estimated to have a lower risk of TBE transmission in both datasets.

Predicted current and future distribution of TBE based on the bioclimatic data
Our results showed an influence of a suite of environmental, anthropogenic, and ecological factors on probability of TBE disease case occurrence in Finland during two time periods. The habitat suitability maps indicated certain hotspot areas for TBE distribution. However, based on these results we cannot reliably estimate how the distribution of TBE may change in future climate in Finland. Therefore, we decided to use current and future global climate data including 19 bioclimatic variables obtained by WorldClim 2.0 and 1.4 (    (Wan et al., 2015) in April, May, June, July, August and October, normalized difference vegetation index (NDVI) (Didan, 2015) and mean snow cover (Hall and Riggs, 2015) in March, April, November and December were calculated in ArcGIS.    Hijmans, 2017;Hijmans et al., 2005), to explore potential changes in the distribution of TBE in the future within the study area of the second dataset. In this substudy, we ignored all the other variables included in the major study as future prediction data were not available for them. For future climate scenarios, for reasons of feasibility only one global climate model (GCM) called IPSL-CM5 (IPSL Earth system model for the 5th Intergovernmental Panel on Climate Change (IPCC) Assessment (2014)) was used (Dufresne et al., 2013). For greenhouse gas scenarios, we chose representative concentration pathways (RCPs) of both medium (RCP 4.5) and high (RCP 8.5) change for the years 2041-2060 and 2061-2080 to reveal possible climate change influences in two different time periods. Temperature in RCP 4.5 scenario ranges from 2.0°to 4.5°and in RCP 8.5 scenario from 3.5°to 4.5° (van Vuuren et al., 2011;Nazarenko et al., 2015). After multicollinearity of bioclimatic variables was run with VIF, six variables for modelling analyses were included (Table A.2). Modelling analyses were performed by using the same workflow as described in Section 2.3. The performance of eight models was similar within all climate conditions (Table A.3). RF and GBM performed best in current and future climate conditions with medium and high change scenarios and had predictive power from good to excellent (ROC = 0.76-0.94). A composition of influential factors in the second dataset is presented in Fig. 7.A.-E.
The habitat suitability map of TBE in current climate conditions indicates a high TBE risk in the Åland Islands, in Southern and Southeastern Finland, and in Ostrobothnia and Northern municipalities (Fig. 8.). In ensemble forecasts for 2041-2060 and 2061-2080 with medium and high change scenarios, there is a slightly wider geographical extent of TBEV transmission risk in Northern Finland in several inland municipalities compared to current climate conditions. Southern Finland is predicted to be at high risk for TBE in all climatic conditions. No increase in TBE risk is predicted for inland municipalities in Central and Eastern Finland.     R. Uusitalo, et al. Ticks and Tick-borne Diseases 11 (2020) 101457 4. Discussion

Validity of the study
In this study, we used the biomod2 ensemble platform in R to create ensemble models to identify areas with suitable habitat conditions for TBE under present-day and future climate in Finland. Ensemble forecasting yields more accurate estimates over single-model estimates (Araújo and New, 2007) and is commonly used to estimate the potential distributions of species and vector-borne diseases under current and future climate conditions (Deka and Morshed, 2018;Eneanya et al., 2018). Two different methods of accuracy, ROC and TSS were used in order to give more comprehensive view of the model performance. The ROC uses all possible thresholds for classifying the scores into confusion matrices and obtains each matrix' sensitivity and specificity; comparing sensitivity against the corresponding proportion of false positives (Fielding and Bell, 1997). The TSS is independent of prevalence and an intuitive method of performance measurement of SDMs when predictions are expressed as presence-absence maps (Allouche et al., 2006). Even though ROC is widely used in modelling species distributions, TSS might be more realistic and practical method (Shabani et al., 2018). The majority of the models used herein produced moderate to high ROC and TSS values in both datasets indicating reliable estimations. ROC scores range from 0 to 1, where models higher than 0.5 predict better than random draws, and TSS values range from −1 to 1, where models higher than zero indicate a performance better than random (Ruete and Leynaud, 2015). In this study, RF and GBM models had the highest predictive power of all models and MAXENT the lowest (Fig. 4.). GBM and RF have been widely used to predict disease distributions (Eneanya et al., 2018;Elith et al., 2008;Breiman, 2001a) and the use of them avoids over-fitting (Bhatt et al., 2013). In the model performance, eight model algorithms performed generally better for the second dataset than the first dataset. This is probably due to the lower amount of TBE presence-absence data (n = 110) in the first dataset than in the second dataset (n = 295).
In our study, data from human TBE cases were available at the municipality level and produced results at the same resolution. Explanatory data were mainly available at higher resolution (e.g. temperature 1000m × 1000m), and because explanatory data were calculated to correspond the resolution of response variable, some information was lost. Due to the narrow geographical range of TBE cases, it was not possible to predict TBE risk through the whole country. The SAC of some covariates may produce some uncertainty in the study results even though spatial autocorrelation does not necessarily generate bias (Diniz-Filho et al., 2003). Game animal density data were not real density data but hunting data, and this may cause bias in the study results even though hunting data correlates with animal densities (Cattadori et al., 2003;Jore et al., 2014;Jaenson et al., 2018). Excluding the Åland Islands, hunting data were performed in game management area (GMA) level, which is larger than municipality boundaries (the mean area of GMA is approximately 25 000 km 2 ) and thus, are not real density data per municipality. Although tick species presence-absence data were available but due to the weakened model performance, data were excluded from the final modelling analysis. The control municipalities were randomly selected from the neighboring municipalities of TBE municipalities and this may cause overmatching to the study results, particularly concerning models with climate-based factors. The same amount of control municipalities (absence locations) and municipalities with TBE (presence locations) were chosen for the modelling analysis based on recommendations of producing reliable species distribution models (Barbet-Massin et al., 2012). Effective TBE vaccination in risk localities after a few cases may make the risk in the nature appear low or even nonexistent. However, data on vaccination coverage were not available. Vaccination coverage is particularly high in the Åland Islands. Furthermore, there are likely other influential variables related to TBE distributions that were not included in the modelling process, such as small mammal densities and microclimate situations. These factors have been confirmed to affect tick abundance and distribution of TBE but suitable data of these variables in Finland are not currently available (Hubálek and Rudolf, 2012;Randolph and Rogers, 2000).

Influential factors of TBE occurrence
Our study suggests that the distribution of TBE is affected by environmental factors such as mean precipitation in April and July, proportion of snow-covered land in November, forest, and field; ecological factors such as white-tailed deer, European hare, and raccoon dog density; and anthropogenic factors including proportion of built-up areas, human population density, and people working in the primary sector. Consistent with previous research from other global locations, these factors have been important drivers for large-scale distributional patterns of TBE (Randolph et al., 2008;Czupryna et al., 2016;Brugger et al., 2017;Gray et al., 2009;Brabec et al., 2017). Both I. ricinus and I. persulcatus are vulnerable to desiccation and consequently require high relative humidity (> 80%) in their microhabitats to be able to quest and survive (Gray, 1998). The increased bursts of humidity provided by more rainfall help in maintaining adequately humid shelters for ticks on the ground floor and reduce moisture loss during questing, improving tick survival and lengthening questing periods. Precipitation in April has previously been found to correlate with TBE incidence (Czupryna et al., 2016). Dry periods at the beginning of the tick season has been found to lead to tick mortality and reduced late-season populations for Ixodes scapularis and I. ricinus ticks (Berger et al., 2014;Perret et al., 2000). July is typically seen as a "hiatus" month in the activity of I. persulcatus adults (Gray et al., 2016;Korenberg, 2000), and I. ricinus nymphs and adults, when ticks cease questing and withdraw to humid microhabitats, either to escape unfavorable ambient conditions or due to diapause (Cayol et al., 2017;Sormunen et al., in press).
Higher precipitation in July might not only allow for continued tick questing during this period, but also ensure that ticks in shelters survive to later activity periods in August-September (Sormunen et al., 2016). Snow cover acts as an insulating blanket over ground litter and tends to further insulate ticks from the frigid winter air temperatures (Lundkvist et al., 2011;Vollack et al., 2017). The proportion of snow-covered land in November may positively affect tick activity and survival because snow cover protects ticks from exposure to freezing and frequent temperature shifts when the air temperature decreases (Dautel et al., 2016).
Medium-sized and large animals such as deer, hares and raccoon dogs are potential hosts for I. ricinus nymphs and adults, and I. persulcatus adults (Gray et al., 2016;Klemola et al., 2019). The population sizes of white-tailed deer and raccoon dogs have significantly increased in Finland recently (Natural Resources Institute of Finland, 2018) and may possibly be connected with higher I. persulcatus and I. ricinus density. White-tailed deer density and hare density have earlier been confirmed to correlate with I. ricinus abundance and consequently TBE distribution (Brugger et al., 2017;Jaenson et al., 2018). The proportion of field and forest area in municipalities have indirect effects on TBE transmission. Forests are typical habitats for many important tick host animals, such as deer and hares, and higher amounts of forests therefore typically increase host animal and consequently tick abundance. On the other hand, the increasing proportion of field area often means that a more fragmented habitat mosaic is formed, wherein the amount of boundary areas between different habitats increases. These increasing edge effects allow for greater biodiversity and often higher animal densities (Tack et al., 2012;Czupryna et al., 2016;Nadolny and Gaff, 2018). High human population density and a high proportion of builtup areas are associated with the large number of population and naturally increase TBE risk. People who are working in the primary sector spend relatively more time outdoors than other sectors and consequently have a higher risk for getting tick bites (Randolph et al., 2008). In our future forecasts, mean temperature of the warmest month and wettest quarter and temperature seasonality were the most influential bioclimatic factors on TBE distribution. High spring and summer temperatures and mild winter temperatures are drivers of new tick establishment and higher TBE risk at high-latitudes in Northern Europe (Gray et al., 2009;Randolph and Rogers, 2000).

Conclusions
Our results confirm the influence of game animal densities and anthropogenic and environmental factors on distribution of human TBE disease cases in Finland. The proportion of built-up areas, field, forest, and snow-covered land, people working in the primary sector, human population density, mean precipitation in April and July, white-tailed deer density, and European hare and raccoon dog density were found correlated with the occurrence of TBE. In future forecasts, temperature variables were the most influential drivers for higher TBE risk. Based on habitat suitability maps of 2007-2011 and 2012-2017, high-risk areas of TBEV transmission were estimated to be in the Åland Islands, the coastal regions of Southern, Western, and Northern Finland, and several municipalities in Central and Eastern Finland. In future forecasts for 2041-2060 and 2061-2080 climate, a slightly wider geographical extent of TBE risk was introduced in the Åland Islands and Southern, Western and Northern Finland, even though the risk itself was not significantly increased (Fig. 7). Identified risk areas were consistent with previous study results, in which Southern Finland (including the Åland Islands) and Southeastern Finland were estimated to be suitable for TBE transmission in the 2020 forecast, and risk areas were suggested to expand up to Central Finland in the 2080 forecast (Randolph and Rogers, 2000). Higher TBE risk in Northern regions will be reasonable as temperature increase is greater in the northernmost latitudes (Trenberth et al., 2007), which makes the region more favorable for tick activity (Soucy et al., 2018;Medlock et al., 2013;Lindgren and Gustafson, 2001;Gray et al., 2009). Our findings provide insight into the identification of risk areas and influential factors for TBE in Finland and can be applied to other regions located at high latitudes in the Northern Hemisphere. Our goal in future studies is to combine larger and more detailed datasets of human TBE cases from Scandinavia and create predictions across Northern Europe.

Declaration of Competing Interest
The authors have declared that no competing interests exist.