# Species Distribution Modeling of Ixodes ricinus (Linnaeus, 1758) Under Current and Future Climates, with a Special Focus on Latvia and Ukraine

Volodymyr Tytar, Iryna Kozynenko, Mihails Pupins, Arturs Škute, Andris Čeirāns, Jean-Yves Georges, Oksana Nekrasova

► To cite this version:

Volodymyr Tytar, Iryna Kozynenko, Mihails Pupins, Arturs Škute, Andris Čeirāns, et al.. Species Distribution Modeling of Ixodes ricinus (Linnaeus, 1758) Under Current and Future Climates, with a Special Focus on Latvia and Ukraine. Climate, 2024, 12 (11), pp.184-189. http://dx.doi.org/10.3390/cli12110184. hal-04777598

HAL Id: hal-04777598

https://hal.science/hal-04777598v1

Submitted on 12 Nov 2024

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
# Species Distribution Modeling of *Ixodes ricinus* (Linnaeus, 1758) Under Current and Future Climates, with a Special Focus on Latvia and Ukraine

Volodymyr Tytar<fnref n="1" />, Iryna Kozynenko<fnref n="1" />, Mihails Pupins<fnref n="2" />, Arturs Škute<fnref n="2" />, Andris Čeirāns<fnref n="2" />, Jean-Yves Georges<fnref n="3" /> and Oksana Nekrasova<fnref n="1" />,<fnref n="2" />,<fnref n="3" />,*

**Abstract:** This study assesses the impact of climate change on the distribution of *Ixodes ricinus*, which transmits Lyme disease, a growing public health concern. Utilizing ensemble models from the R package 'flexsdm' and climate data from WorldClim, ENVIREM, and CliMond, we project habitat suitability changes for the focus species. The models, validated against Lyme disease incidence rates, predict a 1.5-fold increase in suitable habitats in Latvia, contrasted with a 4.5-fold decrease in suitable habitats within Ukraine over the coming decades. SHAP values are analyzed to determine the most influential climatic features affecting tick distribution, providing insights for future vector control and disease prevention strategies. The optimal bioclimatic environment for *I. ricinus* seems to be an intricate balance of moderate temperatures, high humidity, and sufficient rainfall (bio7, 14, 18, 29). Also, radiation during the wettest quarter (bio24) significantly influences tick distribution in northern countries. This implies an increased presence of ticks in Scandinavian countries, Baltic states, etc. These findings largely coincide with our projections regarding bioclimatic suitability for ticks in Latvia and Ukraine. These shifts reflect broader patterns of vector redistribution driven by global warming, highlighting the urgent need to adapt public health planning to the evolving landscape of vector-borne diseases under climate change.

**Keywords:** *Ixodes ricinus*; Lyme disease; species distribution modeling; SHAP values; climate change; Eastern Europe

## 1. Introduction

Climate change is widely recognized as a key factor influencing biodiversity [1,2]. There is a specific interest in evaluating the potential shifts in habitat suitability for medically or economically important species as a result of climate change. According to reports from the IPCC, climate change is expected to alter the distribution of certain vectors and vector-borne diseases [3,4]. Tick-borne illnesses have emerged as a significant health concern in recent years. Ticks belonging to the suborder Ixodida are found across the globe and are highly sensitive to climate conditions due to their reliance on a complex interplay of climate factors for their development and survival [5]. Ticks serve as carriers for various pathogens, such as protozoa, viruses, bacteria, and nematodes [6,7], which are known to cause several diseases, including Lyme borreliosis. In Europe, Lyme borreliosis is among the most prevalent tick-borne diseases in humans [8], with *Ixodes ricinus* (Linnaeus, 1758) being the primary vector species. The disease is particularly widespread in Central Europe, with the highest infection rates recorded in ixodid ticks [9]. A recent study reported that Central Europe has the highest rate of residents with Lyme disease compared to North America (21% and 9%, respectively) [10].

*Ixodes ricinus* is the most prevalent and widely distributed tick species in Europe [11]. It parasitizes a diverse range of wild vertebrates, as well as incidental hosts like humans, livestock, and pets [12]. In Europe, *I. ricinus* is found across a broad geographical area including Scandinavia, the British Isles, Central Europe, France, Spain, Italy, the Balkans, Eastern Europe, and Ukraine [13] (ECDC 2020). Global average temperatures have increased by 0.7 °C over the past century, with a projected additional increase of 1.1 °C by the end of the 21st century [14]. Climate change can have a dual impact on tick populations and the transmission of tick-borne diseases by influencing vector biology. On the one

<footnote n="1">¹ I. I. Schmalhausen Institute of Zoology, National Academy of Sciences of Ukraine, 01030 Kyiv, Ukraine; vtytar@gmail.com (V.T.); kozinenko@gmail.com (I.K.)</footnote>

<footnote n="2">² Institute of Life Sciences and Technologies, Daugavpils University, LV5400 Daugavpils, Latvia; mihails.pupins@gmail.com (M.P.); arturs.skute@du.lv (A.Š.); cuskisa@gmail.com (A.Č.)</footnote>

<footnote n="3">³ Université de Strasbourg, CNRS, IPHC UMR7178, F-67000 Strasbourg, France; jean-yves.georges@iphc.cnrs.fr</footnote>

<footnote>* Correspondence: oneks22@gmail.com; Tel.: +380-672171970</footnote>
hand, studies have demonstrated that climate change contributes to the expansion of *I. ricinus* in northern regions and/or increases their local abundance [15–17]. Conversely, high temperatures can lead to increased mortality in ticks due to desiccation [6]. Additionally, ticks are less active in seeking hosts during periods of high temperatures, which can result in higher mortality rates by reducing their success in finding hosts [18]. Therefore, global warming may negatively affect the species' habitat suitability in the southern parts of its range of distribution [19] and may potentially lead to geographical contractions of its range in areas where the combination of climate variables for development and survival turns out to become less or entirely not favorable. Under these circumstances, there is a need to construct the current climatic niche of *I. ricinus* and project it into future conditions to assess how climatic change will affect the geographic distribution of the vector. Ultimately, climate-induced changes in vector distribution may affect the epidemiology of vector-borne diseases [20] and are expected to be of profound character [21].

With the global incidence of Lyme disease and other tickborne diseases on the rise, there is a critical need for data-driven tools to assess the magnitude of this problem and provide scientifically grounded guidelines for public health decision makers. Moreover, there is an interest in identifying potential areas of distribution of vector species in order to assess the future infection risks with vector-borne diseases and improve surveillance efforts [22]. In this respect, species distribution modeling (SDM) can be considered as a data-driven numerical tool that combines field observations of species occurrence or abundance with environmental estimates (usually climatic) [23]. Modeling the ecological requirements of species to explore future disease transmission patterns has been found to be challenging [24]; however, it is essential to assess the potential impacts of climate change on vector distribution, which can provide a theoretical basis for the prevention and control of vector-borne diseases [25,26].

In this study, we used an SDM approach to investigate the distribution of the Lyme disease vector *I. ricinus* under the current and forecasted (by 2030 and 2050) climate and to rank habitat suitability areas for the tick at a country scale with a focus on Latvia and Ukraine. Of particular concern is Russia's military action in eastern Ukraine since 2022, which has resulted in the substantial relocation of Ukrainian citizens to mainly western and central regions of the country [27,28]. Therefore, this raises the question of how safe these regions are in terms of disease transmission. This could substantially help health officials in managing and monitoring tick distributions, assessing their overlap and potential contact with humans because it is vital to decrease the risk of zoonotic disease transmission.

## 2. Materials and Methods

### 2.1. Tick input Data

As input, SDMs require georeferenced biodiversity observations. Faced with limited data for Ukraine, we built species distribution models using European and adjacent occurrences. Presence data were retrieved from online public databases, which were accessed using the R package 'rgbif', an interface to the Global Biodiversity Information Facility [29,30] (GBIF.org, 2022), supplemented from the literature [31–34], including Ukrainian and Latvian sources [35–37]. The conducted search for occurrence data yielded a total of 12,095 non-duplicate georeferenced records of *I. ricinus* across Europe and adjacent areas. Using the pre-modeling functions from the 'flexsdm' R package, the number of presence records was reduced to 1827 [38].

### 2.2. Climate Data

SDMs are predominantly influenced by climate, with the variables used in their development typically reflecting climatic factors [39,40]. This is logical as climate plays a key role in determining environmental suitability [41]. Bioclimatic parameters were gathered as raster layers from three climatic databases and utilized separately to construct the projected SDMs and assess their performance.

Nineteen bioclimatic variables indicating general patterns in precipitation and temperature, including extremes and temperature seasonality, were obtained at a resolution of 2.5' from the WorldClim website (https://www.worldclim.org/; accessed on 21 November 2022; [42]; Table S1; Miroc6 ssp126).

For the modeling in this study, a set of 16 climatic and 2 topographic variables (from the ENVIREM dataset, downloaded from http://envirem.github.io; accessed on 26 November 2022; Table S2) at a resolution of 2.5' was utilized. These variables were selected based on a recent reassessment of their biological significance, with many likely to have direct relevance to ecological or physiological processes that influence species distributions [43]. These variables are important to consider in species distribution modeling applications, particularly as some (such as potential evaporation) are directly linked to processes crucial for species ecology. The inclusion of topographic variables is also significant as they can alter the impacts of climate descriptors.
Datasets from CliMond v.1.2 were downloaded from https://www.climond.org/ at a resolution of 10' (accessed on 28 November 2022; Table S1; A1B). These datasets comprised the core set of 19 bioclimatic variables (temperature and precipitation) and an additional set of 16 variables (solar radiation and soil moisture) [44,45].

In previous studies, we used various climate change scenarios to predict future conditions [3,40]. We believe that scenarios utilizing average or minimum climate change values are the most optimal for our purposes. After reviewing the literature, we concluded that scenarios such as A1B (CliMond) and Miroc6 (ssp126-WorldClim) are the most suitable [33,34,39,40,42,44,45]. These scenarios assume a balanced emphasis on all factors and represent a middle path between various possible climate changes, making them useful for our forecasts and analysis.

Climate variables often exhibit high collinearity, and most SDM approaches necessitate selecting one variable from highly correlated ones [46]. To address this, the “removeCollinearity” function in the “virtualspecies” R package was utilized [47]. This function evaluates the correlation among variables in the provided stack of environmental variables, identifying variables that are not collinear and grouping them based on their collinearity levels. Collinearity was assessed across Europe and neighboring regions to avoid shifts in collinearity when projecting training extents [48]. Given that climate variables frequently have skewed distributions or outliers, the Spearman correlation method was applied with a threshold of 0.8 [49].

Removing high collinearity and the selection of one amongst strongly correlated predictors gave the following results, presented in Table 1.

<table><caption>**Table 1.** Results of the selection of uncorrelated predictors (Spearman's r < 0.8).</caption><thead><tr><th>Predictor Dataset</th><th>Selected Subsets of Predictors</th></tr></thead><tbody><tr><td>WorldClim v.2</td><td>Annual mean temperature (bio1), Mean diurnal range (bio2), Temperature seasonality (bio4), Maximum temperature of warmest month (bio5), Annual precipitation (bio12), Precipitation seasonality (bio15), Precipitation of driest quarter (bio17)</td></tr><tr><td>ENVIREM</td><td>Annual PET *, Continentality, Emberger's pluviothermic quotient, mean monthly PET of coldest quarter, mean monthly PET of driest quarter, PET seasonality, mean monthly PET of wettest quarter, terrain roughness index</td></tr><tr><td>CliMond v.1.2</td><td>Mean diurnal temperature range (bio2), Minimum temperature of coldest week (bio6), Annual temperature range (bio7), Mean temperature of driest quarter (bio9), Annual precipitation (bio12), Precipitation of driest week (bio14), Precipitation seasonality (bio15), Precipitation of warmest quarter (bio18), Radiation of wettest quarter (bio24), Radiation of warmest quarter (bio26), Radiation of coldest quarter (bio27), Highest weekly moisture index (bio29)</td></tr></tbody></table>

## 2.3. Modeling Procedure

To prepare the input data, pre-modeling functions from the 'flexsdm' R package were utilized [38]. The calibration area was delineated by creating 500 km buffers around presence points. Occurrence data filtering was implemented to mitigate sample bias by randomly eliminating points that were densely clustered (oversampling) in both environmental and geographical space. The modeling functions within the package facilitated the fitting and evaluation of various modeling approaches, encompassing individual algorithms, tuned models, and ensemble models. Seven machine learning (ML) SDM (species distribution models) methods were applied, including 'Generalized Additive Models', 'Gaussian Process Models', 'Generalized Linear Models', 'Maximum Entropy', 'Artificial Neural Networks', 'Random Forest', and 'Support Vector Machine'. This led to the development of a comprehensive ensemble model, where predictions from individual algorithms were amalgamated to generate a consensus distribution, thereby reducing model uncertainty and enhancing model transferability [50]. To create a global model across Europe using SDMs, we applied the Maxent method based on Maximum Entropy (see Supplementary Materials, Figures 4 and S3, Tables 3 and S1, https://biodiversityinformatics.amnh.org/open_source/maxent/; accessed on 28 November 2022) [51,52].

Initially, the models were assessed using the area under the receiver operating characteristic curve (AUC) and the true skill statistic (TSS)[53,54]. AUC scores range from 0 to 1, with 0 indicating consistently incorrect model predictions and 1 signifying consistently accurate model predictions. AUC values between 0.7 and 0.8 are deemed acceptable, while

<footnote>* PET = potential evapotranspiration.</footnote>
values exceeding 0.8 are considered good to excellent. TSS values range from -1 to +1, with -1 representing systematic inaccuracies and +1 indicating systematic correctness. TSS values below 0.4 are considered poor, those between 0.4 and 0.8 are considered useful, and values above 0.8 are classified as good to excellent [25]. In addition to AUC, the continuous Boyce index from the 'flexsdm' package was employed to address the limitations of the AUC. This index, which only requires presence data, evaluates how closely model predictions align with the observed distribution of presences across prediction gradients [55]. It ranges from -1 to +1, with positive values indicating consistent predictions, values near zero suggesting similarity to a random model, and negative values indicating contradictory predictions [56].

Subsequently, it was assumed that maps depicting habitat suitability for the Lyme disease vector could serve as proxies for the potential distribution of the disease agent itself. The performance of each ensemble model was assessed by correlating the derived habitat suitability values with the reported incidence of Lyme disease (expressed as the number of reported cases per 100,000 people) sourced from the Johns Hopkins Lyme and Tick-Borne Disease Dashboard (https://www.hopkinslymetracker.org; accessed on 23 December 2022). To address spatial dependencies, a modified t-test was employed to calculate the statistical significance of the correlation coefficient (a corrected Pearson's correlation) based on geographically effective degrees of freedom, as implemented in the 'SpatialPack' package [57,58].

Finally, the 'best' ensemble SDM was selected based on the aforementioned performance criteria and the correlation between *I. ricinus* habitat suitability values and reported Lyme disease incidence. The top-performing SDM was categorized into areas of low, medium, and high potential habitat suitability. These thresholds were defined using Jenks natural breaks, a method that optimizes the grouping of numerical variables by minimizing the deviation within each class while maximizing the deviation from the means of other groups. Jenks natural breaks provide a standardized approach to determining class intervals for continuous numerical variables [59].

## 2.4. Feature Importance

A critical challenge in the adoption of ML models is their inherent lack of interpretability, which is often referred to as the “black box” problem. Shapley Additive exPlanations (SHAP) [60–62] represents a significant advancement in addressing this concern by providing a framework for understanding how ML models arrive at their predictions. SHAP offers various advantages over other feature importance methods, including its model-agnostic nature. Leveraging game theory concepts, SHAP provides a robust framework for feature importance attribution in ML models, despite the number of used variables. The SHAP value quantifies the magnitude and direction (positive or negative) of the feature’s influence on the prediction [60,61]. In our case, the relative importance of predictors was measured with mean absolute Shapley values, which can be interpreted as the magnitude of the relative contribution of a predictor (feature) towards a model output [63]. The R package ‘shap-values’ (https://github.com/pablo14/; author Pablo Casas; accessed on 10 December 2023) was used in a modified version to calculate SHAP values for a selected model. A summary plot offers a comprehensive view of the most influential features in a model and ranks features based on their effect on the model’s predictions. In other words, SHAP values can be integrated into global explanations such as variable importance, but using a completely different method from the various feature selection methods widely applied so far [64]. Until only now has the use of SHAP for explaining which environmental factors (including climatic) influence the spatial distribution of species occurrences started to be explored more widely [65–69].

Maps of habitat suitability in the GeoTIFF format were processed and visualized in SAGA GIS v.9.3.0 [70]. Statistical data were analyzed using the PAST v.4.04 software package [71,72] and/or the R environment [73].

## 3. Results

As a result of the modeling, we obtained models with a sufficiently high evaluation (Table 2). SDMs built using various selected subsets of predictors from the WorldClim v.2, ENVIREM, and CliMond v.1.2 datasets performed satisfyingly, with the AUC values considered to mostly correspond to good and excellent, and sound TSS values. In this respect, more weight should be given to the continuous Boyce index, which is considered one of the most appropriate metrics for assessing model predictions applied to presence-only datasets and which is a more reliable metric than the AUC [74-77]. Also, the continuous Boyce index shows a high uniformity of performance of the obtained SDMs and it now reveals a challenge to picking the 'best' model. Importantly, values reported in Table 2 indicate models in which present predictions are highly consistent (>0.9) with the distribution of presences in the datasets.

The 'flexsdm' package supports model evaluation based on a number of the abovementioned performance metrics. Amongst the threshold criteria, maximum training sensitivity plus specificity was considered a reasonable choice [78,79].
<table><caption>**Table 2.** Evaluation metrics for SDMs built using various selected subsets of predictors.</caption><thead><tr><th rowspan="2">Predictor Dataset/Subset</th><th colspan="6">Evaluation Metrics</th></tr><tr><th>Area Under Curve (AUC)</th><th>SD * of AUC</th><th>True Skill Statistic (TSS)</th><th>SD of TSS</th><th>Continuous Boyce Index (BOYCE)</th><th>SD of BOYCE</th></tr></thead><tbody><tr><td>WorldClim v.2</td><td>0.84</td><td>0.01</td><td>0.59</td><td>0.02</td><td>0.91</td><td>0.07</td></tr><tr><td>ENVIREM</td><td>0.84</td><td>0.01</td><td>0.56</td><td>0.02</td><td>0.92</td><td>0.06</td></tr><tr><td>CliMond v.1.2 current</td><td>0.79</td><td>0.06</td><td>0.50</td><td>0.12</td><td>0.91</td><td>0.07</td></tr><tr><td>CliMond v.1.2 A1B scenario for 2030</td><td>0.76</td><td>0.06</td><td>0.46</td><td>0.13</td><td>0.87</td><td>0.18</td></tr><tr><td>CliMond v.1.2 A1B scenario for 2050</td><td>0.74</td><td>0.08</td><td>0.41</td><td>0.12</td><td>0.69</td><td>0.16</td></tr></tbody></table>

Further, the ensemble models for current and projected climate scenarios were clipped to Latvia and Ukraine and classified using the corresponding module in SAGA GIS, in accordance with Jenks natural breaks, into three categories, i.e., “high”, “medium” and “low” habitat suitability (Figures 1–3).

<figure><img src="image_2.png" /><figcaption>**Figure 1.** Jenks natural breaks maps of habitat suitability (HS) for *Ixodes ricinus* under current bioclimatic conditions. (a) Ukraine: red, green, and blue, respectively, for HS of high (0.52–0.71), medium (0.29–0.52), and low (0.08–0.29) value (oblasts: 1—Zakarpatska, 2—Lvivska, 3—Volynska, 4—Rivnenska, 5—Ternopilska, 6—Ivano-Frankivska, 7—Zhytomyrska, 8—Khmelnytska, 9—Chernivetska, 10—Vinnytska, 11—Crimea). (b) Latvia: red, green, and blue, respectively, for HS of high (0.62–0.77), medium (0.46–0.62), and low (0.32–0.46) values (regions: 1—Kurzeme, 2—Pieriga, 3—Vidzeme). Administrative units according to OECD, 2022.</figcaption></figure>

<footnote>* SD = standard deviation.</footnote>
<figure><img src="image_3.png" /><figcaption>**Figure 2.** Jenks natural breaks maps of habitat suitability (HS) for *I. ricinus* as predicted by our models by 2030. (a) Ukraine: red, green, and blue, respectively, denote areas with high (0.61–0.77), medium (0.36–0.61), and low (0.1–0.36) HS. Oblasts: 1—Zakarpatska, 2—Lvivska, 3—Ternopilska, 4—Ivano-Frankivska, 5—Chernivetska. (b) Latvia: red, green, and blue, respectively, denote areas with high (0.62–0.76), medium (0.48–0.62), and low (0.34–0.48) HS. Regions: 1—Kurzeme, 2—Pieriga, 3—Vidzeme.</figcaption></figure>

<figure><img src="image_4.png" /><figcaption>**Figure 3.** Jenks natural breaks maps of habitat suitability (HS) for *I. ricinus* as predicted by our models by 2050. (a) Ukraine: red, green, and navy blue, respectively, denote areas with high (0.63–0.73), medium (0.45–0.63), and low (0.33–0.45) HS. Oblasts: 1—Lvivska, 2—Ivano-Frankivska, 3—Odeska. (b) Latvia: red, green, and navy blue, respectively, denote areas with high (0.69–0.8), medium (0.51–0.69), and low (0.33–0.51) HS. Regions: 1—Kurzeme, 2—Pieriga, 3—Vidzeme, 4—Zemgale.</figcaption></figure>

To gain a comprehensive understanding of tick distribution in Europe, we developed a global European model of tick distribution. Our analysis revealed that factors such as the annual temperature range (bio7) and the precipitation of the driest week (bio14) also influence tick distribution (Table 3, Figures 4, 5, S3, and S4).

<table><caption>**Table 3.** Analysis of variable contributions in the SDM Maxent model for *I. ricinus* in Europe (Climond Dataset, see Table S1).</caption><thead><tr><th>Variable</th><th>Percent Contribution %</th><th>Permutation Importance %</th></tr></thead><tbody><tr><td>bio7 Annual temperature range (Bio05-Bio06) (°C)</td><td>34.6</td><td>21.8</td></tr><tr><td>bio14 Precipitation of driest week (mm)</td><td>15.6</td><td>2.1</td></tr><tr><td>bio24 Radiation of wettest quarter (W m<sup>-2</sup>)</td><td>9.1</td><td>8.6</td></tr><tr><td>bio1 Annual mean temperature (°C)</td><td>7.7</td><td>1.7</td></tr></tbody></table>

For the spread of ticks in Ukraine, the most significant parameters are bio18—precipitation of the warmest quarter—and bio7—annual temperature range. In Latvia, the most important factors are bio29—the highest weekly moisture index—and bio24—the radiation of the wettest quarter (Figures 5, S1, S2, and S4).
<figure><img src="image_5.png" /><figcaption>**Figure 4.** Results of the analysis of SDM Maxent for *I. ricinus* in Europe (current, CliMond). The average training AUC for the replicate runs is 0.849, and the standard deviation is 0.001.</figcaption></figure>

<figure><img src="image_6.png" /><figcaption>**Figure 5.** Absolute summary plots of the datasets for Ukraine (a) and Latvia (b), where the average absolute value of the SHAP</figcaption></figure>

values for each variable is taken in order to obtain a bar chart as a function of the contribution of each variable to the prediction of the CliMond current model. The variables are ordered from most (top) to least (bottom) important. The y-axis represents the variances used in the study, which refer to Bio6—minimum temperature of coldest week; bio7—annual temperature range; bio9—mean temperature of the driest quarter; bio15—precipitation seasonality; bio18—precipitation of the warmest quarter; bio24—radiation of the wettest quarter; bio29—highest weekly moisture index. The x-axis represents the corresponding SHAP value.
## 4. Discussion

The distribution of ticks, especially Ixodes ticks, depends on numerous environmental factors that influence their life cycle, survival (including host seeking), and distribution. Key factors include temperature and humidity, which play a critical role in the physiological processes of ticks. These parasites can exist outside of hosts for extended periods, and their sensitivity to climate changes makes them particularly vulnerable to fluctuations in temperature and humidity levels. To better understand pan-European trends, we conducted a study within Europe. Our analysis revealed that factors such as the annual temperature range (bio7) and the precipitation of the driest month (bio14) also influence tick distribution. Bio14, the precipitation of the driest month, is determined by identifying the driest month for each year in the period, i.e., the month with the least precipitation, and then averaging the precipitation across all of the driest months. It is important to note that the driest month may vary from year to year. The annual temperature range (bio7) measures the temperature range between the most extreme (warmest and coldest) months (Table 3, Figure 4). Sharp temperature fluctuations and insufficient humidity can also affect the viability and distribution of ticks [80-83].

In our study, we analyzed SHAP values to understand feature importance governing current tick distribution in Ukraine and Latvia by calculating the average absolute values for each feature. This permitted us to provide an idea of which features are generally most important in influencing predictions within the area of our countries of interest (Figure 5). Given their extended off-host existence as generalist parasites, ixodid ticks, including the castor bean tick, exhibit a marked sensitivity to both temperature and humidity [84]. This sensitivity arises from the critical influence these environmental factors exert on physiological processes and desiccation rates, both of which are essential for tick survival and life cycle completion [33,88]. In our study, SHAP analysis suggests that rainfalls during the warmest quarter (bio18 in Figure 5a) have a particularly strong influence on predictions for Ukraine. Checking the correlation between bio18 and the ensemble SDM for the current climate demonstrated a strong relationship too, which could be explained by the fact that 'precipitation during the warmest quarter' acts as a trade-off factor for other environmental variables, including temperature and humidity. Overall, the optimal bioclimatic environment for the tick species seems to be an intricate balance of moderate temperatures, high humidity, and sufficient rainfall to maintain that humidity, despite the fact that the correlation between humidity and precipitation is not straightforward [82-83].

In addition to climate, the distribution of *I. ricinus* may also be affected by factors such as land cover, particularly vegetation. Ticks require moisture and humidity to survive and dense vegetation provides these conditions by shading the ground and preventing moisture loss [83]. Additionally, ticks in some stages of the tick lifecycle feed on the blood of various animals that frequent forested and other vegetated areas [12]. Because in our modeling we focused on bioclimatic predictors, strictly speaking, vegetation features were not used; however, we tested the connection between bio18 and NDVI (Normalized Difference Vegetation Index). NDVI is a measure of vegetation health and density. Corresponding rasters for the warmest quarter (June–August) were downloaded from the EDIT geoplatform (http://edit.csic.es/GISdownloads.html; accessed on 4 April 2023). Previous studies have shown that areas with higher NDVI values tend to be more suitable habitats for *I. ricinus* ticks [81]. Indeed, we found fairly strong correlations for summer months (June, July, and August). In our opinion, this once again highlights the importance of the 'precipitation of warmest quarter' factor and/or, associated with it, factors in shaping the tick's ecological niche and distribution, specifically in Ukraine.

In fact, for Latvia we have an analogous situation, where the most important features in influencing predictions within the country are suggested to be bio29 ('highest weekly moisture index') and bio24 ('radiation of wettest quarter') (Figure 5b). Both characterize the moisture regime, with a direct or indirect emphasis on the warm time of the year. The study area experiences its highest weekly moisture index in July. Secondly, the wettest quarter occurs from July to September, with peak rainfall in July.

As already implied, there may be another approach, tentatively tagged "biological", to finding a model that would meet our objective aiming to aid health officials in managing and monitoring tick and tick-borne disease distributions. We assume this can be accomplished by correlating the obtained habitat suitability values of each SDM with the reported incidence of Lyme disease and hypothesize that a closer correlation would point towards a model that fits our needs the most. Existing studies in Europe mainly focus on acarological risk assessment, with very limited investigations exploring Lyme disease occurrence in humans [84]. As an example, incidences of Lyme disease in Ukraine are scarcely reported, so proper statistical analyses cannot be achieved.

Therefore, we used data from Romania since this neighboring country presents a fairly similar climate compared to Ukraine, especially if accounting for warm-season temperatures when tick activity is the highest [85]. Annual incidences of Lyme disease were downloaded for 41 counties and the Bucharest Municipality from the Johns Hopkins Lyme and Tick-Borne Disease Dashboard for a 10-year period (from 2012 to 2021) and averaged. Corresponding averaged habitat suitability values for each administrative entity, represented by a polygon shapefile, were obtained
using the 'grid statistics for polygons' module in SAGA GIS. Firstly, using the R package 'trafo', suitable transformations depending on statistical requirements and the data being analyzed were assessed by checking assumptions of normality, homoscedasticity, and linearity. By log-transforming the annual incidences of Lyme disease using the formula $\ln(x + 1)$, all specified assumptions were shown to be acceptable ($p > 0.05$). Secondly, using log-transformed annual incidences of Lyme disease and the averaged habitat suitability values regarding the considered administrative units in Romania, corrected Pearson's correlation coefficient ($r$) controlling for spatial non-independence were calculated. Results concerning each of the SDMs reveal the following: for data extracted from the SDM based on the WorldClim v.2 subset, $r = 0.3815$ and $p$-value = 0.0584, meaning no confirmation of statistical significance; in the case of the SDM based on the ENVIREM subset, $r = 0.3989$ and the $p$-value is 0.0422; and in the case of the SDM based on the CliMond v.1.2 subset, $r = 0.59150$ and the $p$-value is 0.0015. For the last two subsequent cases, conclusions regarding correlations are of sound statistical significance ($p < 0.05$); however, a stronger correlation coefficient was found for the SDM built on the CliMond v.1.2 subset. Therefore, we can conclude that this particular model is the best fit for our purpose.

In this study, besides the investigation of the present distribution of the tick *I. ricinus*, we aimed to assess how future climatic changes predicted for 2030 and 2050 could affect tick distribution in order to be prepared for proactive management actions involving both the vector and spread of Lyme disease. For this, we employed the CliMond v.1.2 dataset projection of future climate for 2030 and 2050 generated for the A1B scenario for emissions of greenhouse gases and sulfate aerosols (Figures 4, 5, and S4 [39]), which envisages a balanced emphasis on all energy sources. The produced ensemble models using the 'flexsdm' algorithm performed fairly well (Table 2), with acceptable AUC values, consistent TSS, and a high value of the continuous Boyce index. Somewhat poorer was the performance of the 2050 model, but it was within acceptable limits.

In terms of the current bioclimate, such as for Latvia, areas of high habitat suitability for *I. ricinus* occur on 22% of the territory of the country, mainly in the west and, to a lesser extent, in the north (Figure 1b). In Ukraine, areas of high habitat suitability for *I. ricinus* occupy around 20% of the country, with locations found predominantly in the west (Figure 1a). Other smaller patches of habitat suitability were found around Kyiv and in the south of Crimea, where ticks are the most abundant representative of the family [86]. Russia's military action in Ukraine since 2022 has resulted in the internal displacement of over 3.5 million people (The UN Refugee Agency; https://www.unhcr.org/emergencies/ukraine-emergency; accessed on 4 April 2023), with potential health consequences [87–90]. Moreover, while refugee movements are associated with increased risks of infectious disease transmission and are likely to affect zoonotic disease risks, it remains unclear how forced migrations affect disease dynamics [90]. Human susceptibility to disease during forced migration might increase due to exhaustion, malnutrition, and stress arising from displacement, magnified by crowded and substandard living conditions [71,72].

In Ukraine, by 2030, our models predict a reduction in surface area of around two-fold of high habitat suitability for *I. ricinus*, decreasing to <10% of the country area (Figure 2a). These are expected to be located exclusively in the west of Ukraine, whereas patches of habitat suitability around Kyiv and in the south of Crimea are predicted to most likely disappear. By 2050, highly suitable habitats for *I. ricinus* are predicted to continue their contraction (to only 4.5% of the country area), with potential reduction to Lvivska and Ivano-Frankivska oblasts, despite a patch of high habitat suitability that is predicted to appear in the southern region of Odesa (Figure 3a).

On the contrary, in Latvia, areas of high habitat suitability for *I. ricinus* are predicted to expand and occupy up to 27% of the territory of the country by 2030, and up to 33% in 2050, engulfing the entire western region of Kurzeme (Figure 2b). Upon analyzing the models, we can predict that in the north and south of its range, the potential distribution of *I. ricinus* may take different directions: for instance, in Ukraine, in the south of its range, predicted increased fragmentation may lead to a significant decrease in tick populations. On the contrary, in the northern part of Ukraine, the models predict an increase in the potential suitable territory for the spread of the species.

In this context, it is important to consider changes in the distribution of host species, which can significantly affect the spread of ticks such as *I. ricinus*. These ticks feed on the blood of various animals, including rodents, birds, foxes, hares, and larger animals such as deer, moose, humans, and domestic animals. Changes in the distribution of host species, especially in the context of the large-scale migration of people and domestic animals associated with military actions, can be a significant factor influencing the spread of ticks. This is an important topic for future research, especially for those dealing with the medical aspects of tick distribution and related diseases. Furthermore, understanding the ecological niche of *I. ricinus* is crucial for predicting its distribution under current and future climatic conditions. The ecological niche encompasses the range of environmental conditions that support the survival and reproduction of the species. Climate change scenarios, such as those used in our study, help model potential changes in these niches, allowing us to anticipate shifts in tick distribution. Species distribution modeling (SDM) primarily focuses on the theoretical or potential niche, also known as the fundamental niche [40]. This is because SDM typically uses environmental
and species presence data to predict where a species could potentially exist in the absence of constraints such as competition or predation.

Our results regarding the prospects of bioclimatic suitability for the castor bean tick in Ukraine and Latvia are consistent with findings from previous studies. Indeed, previously published SDMs have predicted that Ixodes ticks will increase their distribution and abundance under global warming, with the most generally envisaged scenario being an expansion of the geographical range for ticks, including *I. ricinus*, with potential increases in abundance in some regions [91–93]. Bioclimatic models also predict a shift in the distribution of the tick species, particularly at higher altitudes and latitudes, like Scandinavia, the Baltic states, and Belarus. Conversely, tick populations are expected to decline in areas like the Alps, central and western Italy, and northwestern Poland. This means an increase in tick presence for Scandinavian countries (Sweden, Norway, and Finland), the Baltic states (Estonia, Latvia, and Lithuania), Denmark, and Belarus [94-97].

**5. Conclusions**

Global warming is one of the most important factors in the redistribution of tick populations and the outbreak of tick-borne diseases, and this is supported by recent studies [96–98]. SDMs built for estimating such redistribution using various selected subsets of predictors from the WorldClim v.2, ENVIREM, and CliMond v.1.2 datasets performed satisfyingly. In particular, the continuous Boyce index shows a high uniformity of performance of the obtained SDMs. In another approach correlating the obtained habitat suitability values of each SDM with the reported incidence of Lyme disease, we found a stronger correlation coefficient for the SDM built on the CliMond v.1.2 subset; therefore, we concluded that this particular model is the best fit for our purpose.

SDMs employing the CliMond v.1.2 dataset projection for future climates (2030 and 2050) generated for the A1B scenario for emissions of greenhouse gases and sulfate aerosols performed fairly well; somewhat poorer was the performance of the 2050 model, but this was within acceptable limits. In terms of the current bioclimate, areas of high habitat suitability for *I. ricinus* in Ukraine are predicted to occupy around one-fifth of the country, with locations found predominantly in the west, thus highlighting here the elevated risk of disease transmission. Regarding Latvia, areas of high habitat suitability for the tick are predicted to occur mainly in the west and, to a lesser extent, in the north of the country. Modeling suggests a significant reduction in highly suitable habitats for *I. ricinus* within Ukraine over the coming decades, whereas for Latvia the expansion of such suitable habitats is expected to continue, aggravating future health concerns accordingly.

We report that the relative importance of predictors in our case study was successfully measured with mean absolute Shapley values, underlining the pivotal significance of moisture conditions in the warm months of the year for creating highly suitable habitats for the tick species. Therefore, we believe that the use of multiple types of climatic factors can significantly improve the quality of predictions.

**Supplementary Materials:** The following supporting information can be downloaded at www.mdpi.com/xxx/s1. Table S1: Bioclimatic variables (35) from the WorldClim v.2 (http://www.worldclim.com/version2; accessed on 21 November 2022; [42]) and the CliMond (https://www.climond.org/ at a resolution of 10' (accessed on 28 November 2022; [44,45]), Table S2: Environmental variables (18) from ENVIREM (https://envirem.github.io/; accessed on 26 November 2022; [43]) datasets. Figure S1: Graphs of the most significant bioclimatic variables (CliMond) from different parts of the tick range: a – Latvia; b – Ukraine. Figure S2: Results of niche space exploration (CliMond): 3D ecological envelope: a – Latvia; b – Ukraine. Figure S3: Results of the analysis of SDM Maxent for *I. ricinus* in Europe: (A) current; (B) 2050 (2041–2060) (green-yellow color > 0.5, WorldClim dataset). The average training AUC for the replicate runs is 0.824, and the standard deviation is 0.001. Figure S4: Changes in key bioclimatic variables important for the spread of ticks in Ukraine and Latvia (CliMond, Figures 5 and S1).

**Author Contributions:** Conceptualization, V.T., I.K., M.P., and O.N.; methodology, V.T. and O.N.; software, V.T. and O.N.; validation, V.T. and O.N.; formal analysis, V.T. and O.N.; investigation, V.T., O.N., M.P., and A.Č.; resources, V.T., I.K., O.N., and M.P.; data curation, V.T., I.K., O.N., M.P., and A.Š.; writing—original draft preparation, V.T. and O.N.; writing—review and editing, V.T., O.N., I.K., M.P., A.Č., J.-Y.G., and A.Š.; visualization, V.T.; supervision, V.T., O.N., I.K., M.P., A.Š., and J.-Y.G.; project administration, V.T., O.N., M.P., J.-Y.G., and A.Š. All authors have read and agreed to the published version of the manuscript.

**Funding:** This research was funded through the Emys-R project (https://emysr.cnrs.fr), funded through the 2020–2021 Biodiversa+ and Water JPI joint call for research projects, under the BiodivRestore ERA-NET Cofund (GA N°101003777), with the EU and the funding organizations Agence Nationale de la Recherche (ANR, France, grant ANR-21-BIRE-0005), Bundesministerium für Bildung und Forschung (BMBF, Germany, grant 16LW015), the State Education Development Agency (VIAA, Latvia, grant ES RTD/2022/2), and the National Science Center (NSC, Poland, grant 2021/03/Y/NZ8/00101), and by the project "Ecological and socioeconomic thresholds as a basis for defining adaptive management triggers in Latvian pond aquaculture" (Izp-2021/1-0247). O.N. was supported by Collège de France and Agence Nationale de la Recherche (ANR) through the PAUSE ANR Ukraine program (grant ANR-23-PAUK-0074).
**Data Availability Statement:** Data derived from public domain resources: 1. Species presence records *Ixodes ricinus* (Linnaeus, 1758): GBIF.org (22 May 2024) GBIF Occurrence Download, https://doi.org/10.15468/dl.bf6zjp. 2. Environmental data: The WorldClim v2 website (https://www.worldclim.org/; accessed on 21 November 2022). The ENVIREM dataset was downloaded from http://envirem.github.io (accessed 26 November 2022). The CliMond v1.2 (current, 2030, 2050) datasets were downloaded from https://www.climond.org/ at 10' resolution (accessed on 28 November 2022). 3. Software: Maxent, https://biodiversityinformatics.amnh.org/open_source/maxent/ (accessed on 28 November 2022); the R package 'rgbif', https://github.com/ropensci/rgbif; the R package 'flexsdm', https://sjevelazco.github.io/flexsdm/; https://besjournals.onlinelibrary.wiley.com/doi/full/10.1111/2041-210X.13874; the R package "virtualspecies", https://borisleroy.com/virtualspecies/; https://cran.r-project.org/web/packages/trafo/index.html; the R package 'shap-values', https://github.com/pablo14/; the R-package 'trafo', https://github.com/cran/trafo; the R-package 'SpatialPack', http://spatialpack.mat.utfsm.cl/; the PAST software (v.4.04; accessed on 28 November 2022) package and/or the R environment [72-74], https://www.nhm.uio.no/english/research/resources/past/; the Saga GIS software, https://saga-gis.sourceforge.io/en/index.html (v.9.3.0; accessed on 28 November 2022).

**Acknowledgments:** We are grateful to project 16-00-F02201-000002 for the use of the mobile complex of scientific laboratories for research purposes.

**Conflicts of Interest:** The authors declare no conflicts of interest.

## References

1. Román-Palacios, C.; Wiens, J.J. Recent responses to climate change reveal the drivers of species extinction and survival. *Proc. Natl. Acad. Sci. USA* **2020**, 117, 4211-4217. https://doi.org/10.1073/pnas.1913007117.

2. Weiskopf, S.R.; Rubenstein, M.A.; Crozier, L.G.; Gaichas, S.; Griffis, R.; Halofsky, J.E.; Hyde, K.J.; Morelli, T.L.; Morisette, J.T.; Muñoz, R.C.; et al. Climate change effects on biodiversity, ecosystems, ecosystem services, and natural resource management in the United States. *Sci. Total. Environ.* **2020**, 733, 137782. https://doi.org/10.1016/j.scitotenv.2020.137782.

3. Parry, M.L.; Canziani, O.F.; Palutikof, J.P.; van der Linden, P.J.; Hanson, C.E. *Contribution of Working Group II to the Fourth Assessment Report of the Intergovernmental Panel on Climate Change*; Cambridge University Press: Cambridge, UK.; New York, NY, USA, 2007.

4. Tytar, V.; Nekrasova, O.; Pupins, M.; Skute, A.; Kirjušina, M.; Gravele, E.; Mezaraupe, L.; Marushchak, O.; Čeirāns, A.; Kozynenko, I.; et al. Modeling the Distribution of the Chytrid Fungus *Batrachochytrium dendrobatidis* with Special Reference to Ukraine. *J. Fungi* **2023**, 9, 607. https://doi.org/10.3390/jof9060607.

5. Cumming, G. Comparing climate and vegetation as limiting factors for species ranges of African ticks. *Ecology* **2002**, 83, 255-268.

6. Sonenshine, D.; Roe, R.M. *Biology of Ticks*, 2nd ed.; Oxford University Press: Oxford, UK, 2013; 496p.

7. Čeirāns, A.; Pupins, M.; Kirjusina, M.; Gravele, E.; Mezaraupe, L.; Nekrasova, O.; Tytar, V.; Marushchak, O.; Garkajs, A.; Petrov, I.; et al. Top-down and bottom-up effects and relationships with local environmental factors in the water frog-helminth systems in Latvia. *Sci. Rep.* **2023**, 13, 8621. https://doi.org/10.1038/s41598-023-35780-7.

8. Wormser, G.P.; Dattwyler, R.J.; Shapiro, E.D.; Halperin, J.J.; Steere, A.C.; Klempner, M.S.; Krause, P.J.; Bakken, J.S.; Strle, F.; Stanek, G.; et al. The clinical assessment, treatment, and prevention of Lyme disease, human granulocytic anaplasmosis, and babesiosis: clinical practice guidelines by the Infectious Diseases Society of America. *Clin. Infect. Dis.* **2006**, 43, 1089-1134. https://doi.org/10.1086/508667.

9. Derdáková, M.; Lencakova, D. Association of genetic variability within the Borrelia burgdorferi sensu lato with the ecology, epidemiology of Lyme borreliosis in Europe. *Ann. Agric. Environ. Med.* **2005**, 12, 165.

10. Dong, Y.; Zhou, G.; Cao, W.; Xu, X.; Zhang, Y.; Ji, Z.; Yang, J.; Chen, J.; Liu, M.; Fan, X.; et al. Global seroprevalence and sociodemographic characteristics of Borrelia burgdorferi sensu lato in human populations: A systematic review and meta-analysis. *BMJ Glob. Health* **2022**, 7, e007744.

11. Gray, J.S.; Dautel, H.; Estrada-Peña, A.; Kahl, O.; Lindgren, E. Effects of climate change on ticks and tick-borne diseases in Europe. *Interdiscip. Perspect. Infect. Dis.* **2009**, *2009*, 593232. https://doi.org/10.1155/2009/593232.

12. Estrada-Peña, A.; de la Fuente, J. The ecology of ticks and epidemiology of tick-borne viral diseases. *Antivir. Res.* **2014**, *108*, 104-128. https://doi.org/10.1016/j.antiviral.2014.05.016.

13. ECDC. *Ixodes ricinus*—Factsheet for Experts. 2020. Available online: https://www.ecdc.europa.eu/en/disease-vectors/facts/tick-factsheets/ixodesricinus (accessed on 23 November 2022).

14. Intergovernmental Panel on Climate Change (IPCC). Summary for Policymakers. In: *Global Warming of 1.5°C: IPCC Special Report on Impacts of Global Warming of 1.5°C above Pre-Industrial Levels in Context of Strengthening Response to Climate Change, Sustainable Development, and Efforts to Eradicate Poverty*. Cambridge University Press; 2022:1-24. https://doi.org/10.1017/9781009157940.001.

15. Jore, S.; Vanwambeke, S.O.; Viljugrein, H.; Isaksen, K.; Kristoffersen, A.B.; Woldehiwet, Z.; Johansen, B.; Brun, E.; Brun-Hansen, H.; Westermann, S.; et al. Climate and environmental change drives *Ixodes ricinus* geographical expansion at the northern range margin. *Parasites Vectors* **2014**, *7*, 11.

16. Lindgren, E.; Talleklint, L.; Polfeldt, T. Impact of climatic change on the northern latitude limit and population density of the disease-transmitting European tick *Ixodes ricinus*. *Environ. Health Perspect.* **2000**, *108*, 119. https://doi.org/10.2307/3454509.

17. Tomkins, J.L.; Aungier, J.; Hazel, W.; Gilbert, L. Towards an evolutionary understanding of questing behaviour in the tick *Ixodes ricinus*. *PLoS ONE* **2014**, *9*, e110028. https://doi.org/10.1371/journal.pone.0110028.
18. Randolph, S.E.; Storey, K. Impact of microclimate on immature tick-rodent host interactions (Acari: Ixodidae): Implications for parasite transmission. *J. Med. Entomol.* **1999**, *36*, 741–748. https://doi.org/10.1093/jmedent/36.6.741.

19. Porretta, D.; Mastrantonio, V.; Amendolia, S.; Gaiarsa, S.; Epis, S.; Genchi, C.; Bandi, C.; Otranto, D.; Urbanelli, S. Effects of global changes on the climatic niche of the tick *Ixodes ricinus* inferred by species distribution modelling. *Parasites Vectors* **2013**, *6*, 271. https://doi.org/10.1186/1756-3305-6-271.

20. Gage, K.L.; Burkot, T.R.; Eisen, R.J.; Hayes, E.B. Climate and vector borne disease. *Am. J. Prev. Med.* **2008**, *35*, 436–450.

21. Levy, S. Ticking time bomb? Climate change and *Ixodes scapularis*. *Environ. Health Perspect.* **2014**, *122*, A168. https://doi.org/10.1289/ehp.122-a168.

22. Curriero, F.C.; Wychgram, C.; Rebman, A.W.; Corrigan, A.E.; Kvit, A.; Shields, T.; Aucott, J.N. The Lyme and Tickborne Disease Dashboard: A map-based resource to promote public health awareness and research collaboration. *PLoS ONE* **2021**, *16*, e0260122. https://doi.org/10.1371/journal.pone.0260122.

23. Franklin, J. *Mapping Species Distributions: Spatial Inference and Prediction*; Cambridge University Press: Cambridge, UK, 2009; pp. 1–340.

24. Patz, J.A.; Graczyk, T.K.; Geller, N.; Vittor, A.Y. Effects of environmental change on emerging parasitic diseases. *Int. J. Parasitol.* **2000**, *30*, 1395–1405. https://doi.org/10.1016/s0020-7519(00)00141-7.

25. Zhang, L.; Liu, S.; Sun, P.; Wang, T.; Wang, G.; Zhang, X.; Wang, L. Consensus forecasting of species distributions: the effects of niche model performance and niche properties. *PLoS ONE* **2015**, *10*, e0120056. https://doi.org/10.1371/journal.pone.0120056.

26. Zhang, L.; Ma, D.; Li, C.; Zhou, R.; Wang, J.; Liu, Q. Projecting the potential distribution areas of *Ixodes scapularis* (Acari: Ixodidae) driven by climate change. *Biology* **2022**, *11*, 107. https://doi.org/10.3390/biology11010107.

27. Ukrainian Crisis: Situational Analysis, 5 March 2024. Data Friendly Space. 2024, 49p. Available online: https://reliefweb.int/report/ukraine/ukrainian-crisis-situational-analysis-05-march-2024 (accessed on 5 March 2024).

28. Ukrainian Refugee Crisis (2022–Present). Available online: https://w.wiki/A2Hi (accessed on 9 May 2024).

29. Chamberlain, S.; Oldoni, D.; Barve, V.; Desmet, P.; Geffert, L.; Mcglinn, D.; Ram, K.; rOpenSci; Waller, J. rgbif: Interface to the Global Biodiversity Information Facility API. 2022, R Package Version 3.7.4. Available online: https://CRAN.R-project.org/package=rgbif (accessed on 9 March 2023).

30. GBIF.org. GBIF Occurrence Download. Available online: https://doi.org/10.15468/dl.bf6zjp (accessed on 22 May 2024).

31. Estrada-Peña, A. Distribution, Abundance, and Habitat Preferences of *Ixodes ricinus* (Acari: Ixodidae) in Northern Spain. *J. Med. Entomol.* **2001**, *38*, 361–370. https://doi.org/10.1603/0022-2585-38.3.361.

32. Estrada-Pena, A.; Venzal, J.M.; Acedo, C.S. The tick *Ixodes ricinus*: Distribution and climate preferences in the western Palearctic. *Med. Vet. Entomol.* **2006**, *20*, 189–197. https://doi.org/10.1111/j.1365-2915.2006.00622.x.

33. Alkishe, A.A.; Peterson, A.T.; Samy, A.M. Climate change influences on the potential geographic distribution of the disease vector tick *Ixodes ricinus*. *PLoS ONE* **2017**, *12*, e0189092. https://doi.org/10.1371/journal.pone.0189092.

34. Noll, M.; Wall, R.; Makepeace, B.L.; Newbury, H.; Adaszek, L.; Bødker, R.; Estrada-Peña, A.; Guillot, J.; da Fonseca, I.P.; Probst, J.; et al. Predicting the distribution of *Ixodes ricinus* and Dermacentor reticulatus in Europe: A comparison of climate niche modelling approaches. *Parasites Vectors* **2023**, *16*, 384. https://doi.org/10.1186/s13071-023-05959-y.

35. Yemchuk Ye.M. Fauna Ukrainy. *Ixodid ticks*]; Academy of Sciences of the Ukrainian SSR Press; 1960; Volume 25, 168p. [In Ukrainian].

36. Akimov, I.A.; Nebogatkin, I.V. Distribution of *Ixodes ricinus* (Arachnida, Ixodidae) in Ukraine in the context of tick hazard, and factors favoring its persistence in conditions of fast-going environmental change. *Zoodiversity* **2022**, *56*, 429–434. https://doi.org/10.15407/zoo2022.05.429.

37. Capligina, V.; Seleznova, M.; Akopjana, S.; Freimane, L.; Lazovska, M.; Krumins, R.; Kivrane, A.; Namina, A.; Aleinikova, D.; Kimsis, J.; et al. Large-scale countrywide screening for tick-borne pathogens in field-collected ticks in Latvia during 2017–2019. *Parasites Vectors* **2020**, *13*, 351. https://doi.org/10.1186/s13071-020-04219-7.

38. Velazco, S.J.E.; Rose, M.B.; de Andrade, A.F.A.; Minoli, I.; Franklin, J. flexsdm: An r package for supporting a comprehensive and flexible species distribution modelling workflow. *Methods Ecol. Evol.* **2022**, *13*, 1661–1669. https://doi.org/10.1111/2041-210x.13874.

39. Kriticos, D.J. Regional climate-matching to estimate current and future sources of biosecurity threats. *Biol. Invasions* **2012**, *14*, 1533–1544. https://doi.org/10.1007/s10530-011-0033-8.

40. Nekrasova, O.D.; Tytar, V.M.; Kuybida, V.V. GIS Modeling of Climate Change Vulnerability of Amphibians and Reptiles in Ukraine. In *NAS of Ukraine*; Schmalhausen Institute of Zoology NAS: Kyiv, Ukraine, 2019; 204p, ISBN 978-966-02-8956-7. (In Ukrainian).

41. Schrodt, F.; Santos, M.J.; Bailey, J.J.; Field, R. Challenges and opportunities for biogeography—What can we still learn from von Humboldt? *J. Biogeogr.* **2019**, *46*, 1631–1642. https://doi.org/10.1111/jbi.13616.

42. Fick, S.E.; Hijmans, R.J. WorldClim 2: New 1-km spatial resolution climate surfaces for global land areas. *Int. J. Climatol.* **2017**, *37*, 4302–4315. https://doi.org/10.1002/joc.5086.

43. Title, P.O.; Bemmels, J.B. ENVIREM: An expanded set of bioclimatic and topographic variables increases flexibility and improves performance of ecological niche modeling. *Ecography* **2018**, *41*, 291–307. https://doi.org/10.1111/ecog.02880.

44. Kriticos, D.J.; Webber, B.L.; Leriche, A.; Ota, N.; Macadam, I.; Bathols, J.; Scott, J.K. CliMond: Global high-resolution historical and future scenario climate surfaces for bioclimatic modelling. *Methods Ecol. Evol.* **2012**, *3*, 53–64. https://doi.org/10.1111/j.2041-210x.2011.00134.x.

45. Kriticos, D.J.; Jarošik, V.; Ota, N. Extending the suite of bioclim variables: A proposed registry system and case study using principal components analysis. *Methods Ecol. Evol.* **2014**, *5*, 956–960. https://doi.org/10.1111/2041-210x.12244.
46. Braunisch, V.; Coppes, J.; Arletta, R.; Suchant, R.; Schmid, H.; Bollmann, K. Selecting from correlated climate variables: A major source of uncertainty for predicting species distributions under climate change. *Ecography* **2013**, *36*, 971-983.
https://doi.org/10.1111/j.1600-0587.2013.00138.x.

47. Leroy, B.; Meynard, C.N.; Bellard, C.; Courchamp, F. virtualspecies, an R package to generate virtual species distributions. *Ecography* **2016**, *39*, 599-607. https://doi.org/10.1111/ecog.01388.

48. Feng, X.; Park, D.S.; Liang, Y.; Pandey, R.; Papeş, M. Collinearity in ecological niche modeling: Confusions and challenges. *Ecol. Evol.* **2019**, *9*, 10365-10376. https://doi.org/10.1002/ece3.5555.

49. Shrestha, N. Detecting multicollinearity in regression analysis. *Am. J. Appl. Math. Stat.* **2020**, *8*, 39-42.
https://doi.org/10.12691/ajams-8-2-1.

50. Araujo, M.; New, M. Ensemble forecasting of species distributions. *Trends Ecol. Evol.* **2007**, *22*, 42-47.
https://doi.org/10.1016/j.tree.2006.09.010.

51. Peterson, A.T. Uses and requirements of ecological niche models and related distributional models. *Biodivers. Informatics* **2006**, *3*.

52. Peterson, A.T.; Soberón, J.; Pearson, R.G.; Anderson, R.P.; Martínez-Meyer, E.; Nakamura, M.; Araújo, M.B. *Ecological Niches and Geographic Distributions*; Princeton University Press: Oxford, UK, 2011; 314p.

53. Allouche, O.; Tsoar, A.; Kadmon, R. Assessing the accuracy of species distribution models: Prevalence, kappa and the true skill statistic (TSS). *J. Appl. Ecol.* **2006**, *43*, 1223-1232. https://doi.org/10.1111/j.1365-2664.2006.01214.x.

54. Lobo, J.M.; Jimenez-Valverde, A.; Real, R. AUC: A misleading measure of the performance of predictive distribution models. *Glob. Ecol. Biogeogr.* **2008**, *17*, 145-151.

55. Boyce, M.S.; Vernier, P.R.; Nielsen, S.E.; Schmiegelow, F.K. Evaluating resource selection functions. *Ecol. Model.* **2002**, *157*, 281-300. https://doi.org/10.1016/s0304-3800(02)00200-4.

56. Hirzel, A.H.; Le Lay, G.; Helfer, V.; Randin, C.; Guisan, A. Evaluating the ability of habitat suitability models to predict species presences. *Ecol. Model.* **2006**, *199*, 142-152. https://doi.org/10.1016/j.ecolmodel.2006.05.017.

57. Medina, L.; Kreutzmann, A.-K.; Rojas-Perilla, N.; Castro, P. The R Package trafo for transforming linear regression models. *R J.* **2019**, *11*, 99-123. https://doi.org/10.32614/rj-2019-054.

58. Osorio, F.; Vallejos, R.; Cuevas, F. SpatialPack: Computing the association between two spatial processes. *arXiv* **2016**, arXiv:1611.05289.

59. Jenks, G.F.; Caspall, F.C. Error on choroplethic maps: Definition, measurement, reduction. *Ann. Assoc. Am. Geogr.* **1971**, *61*, 217-244.

60. Lundberg, S.M.; Lee, S.I. A unified approach to interpreting model predictions. In *Proceedings of the 31st International Conference on Machine Learning* 2017, Long Beach, CA, USA, 4-9 December 2017; pp. 3547-3555.

61. Lundberg, S.M.; Nair, B.; Vavilala, M.S.; Horibe, M.; Eisses, M.J.; Adams, T.; Liston, D.E.; Low, D.K.-W.; Newman, S.-F.; Kim, J.; et al. Explainable machine-learning predictions for the prevention of hypoxaemia during surgery. *Nat. Biomed. Eng.* **2018**, *2*, 749-760.
https://doi.org/10.1038/s41551-018-0304-0.

62. Lin, K.; Gao, Y. Model interpretability of financial fraud detection by group SHAP. *Exp. Syst. Appl.* **2022**, *210*, 118354.

63. Bhattacharyay, S.; Milosevic, I.; Wilson, L.; Menon, D.K.; Stevens, R.D.; Steyerberg, E.W.; Nelson, D.W.; Ercole, A.; the CEN-TER-TBI investigators participants. The leap to ordinal: Detailed functional prognosis after traumatic brain injury with a flexible modelling approach. *PLoS ONE* **2022**, *17*, e0270973. https://doi.org/10.1371/journal.pone.0270973.

64. Song, L.; Estes, L. itsdm: Isolation forest-based presence-only species distribution modelling and explanation in r. *Methods Ecol. Evol.* **2023**, *14*, 831-840. https://doi.org/10.1111/2041-210x.14067.

65. Cha, Y.; Shin, J.; Go, B.; Lee, D.-S.; Kim, Y.; Kim, T.; Park, Y.-S. An interpretable machine learning method for supporting ecosystem management: Application to species distribution models of freshwater macroinvertebrates. *J. Environ. Manag.* **2021**, *291*, 112719. https://doi.org/10.1016/j.jenvman.2021.112719.

66. Ghafarian, F.; Wieland, R.; Lütschwager, D.; Nendel, C. Application of extreme gradient boosting and Shapley Additive explanations to predict temperature regimes inside forests from standard open-field meteorological data. *Environ. Model. Softw.* **2022**, *156*, 105466. https://doi.org/10.1016/j.envsoft.2022.105466.

67. Bourhis, Y.; Bell, J.R.; Shortall, C.R.; Kunin, W.E.; Milne, A.E. Explainable neural networks for trait-based multispecies distribution modelling—A case study with butterflies and moths. *Methods Ecol. Evol.* **2023**, *14*, 1531-1542.
https://doi.org/10.1111/2041-210x.14097.

68. Lehnen, S.E.; Lombardi, J.V. Climate envelope modeling for ocelot conservation planning: Peering inside the black box. *Ecosphere* **2023**, *14*, e4477. https://doi.org/10.1002/ecs2.4477.

69. Scavuzzo, C.M.; Scavuzzo, J.M.; Campero, M.N.; Anegagrie, M.; Aramendia, A.A.; Benito, A.; Periago, V. Feature importance: Opening a soil-transmitted helminth machine learning model via SHAP. *Infect. Dis. Model.* **2022**, *7*, 262-276.
https://doi.org/10.1016/j.idm.2022.01.004.

70. Conrad, O.; Bechtel, B.; Bock, M.; Dietrich, H.; Fischer, E.; Gerlitz, L.; Wehberg, J.; Wichmann, V.; Böhner, J. System for Automated Geoscientific Analyses (SAGA) v. 2.1.4. *Geosci. Model Dev.* **2015**, *8*, 1991-2007. https://doi.org/10.5194/gmd-8-1991-2015.

71. Hammer, C.C.; Brainard, J.; Hunter, P.R. Risk factors and risk factor cascades for communicable disease outbreaks in complex humanitarian emergencies: A qualitative systematic review. *BMJ Glob. Health* **2018**, *3*, e000647.

72. Hammer, Ø.; Harper, D.A.T.; Ryan, P.D. PAST: Paleontological Statistics software package for education and data analysis. *Palaeontol. Electron.* **2001**, *4*, 9. Available online: http://palaeo-electronica.org/2001_1/past/issue1_01.htm (accessed on 28 November 2022).
73. R Core Team 2020. *R: A Language and Environment for Statistical Computing*; R Foundation for Statistical Computing: Vienna, Austria, 2020. Available online: https://www.R-project.org/ (accessed on 28 November 2022).

74. Cianfrani, C.; Le Lay, G.; Hirzel, A.H.; Loy, A. Do habitat suitability models reliably predict the recovery areas of threatened species? *J. Appl. Ecol.* **2010**, *47*, 421–430.

75. Manzoor, S.A.; Griffiths, G.; Lukac, M.; Manzoor, S.A.; Griffiths, G.; Lukac, M. Species distribution model transferability and model grain size—Finer may not always be better. *Sci. Rep.* **2018**, *8*, 7168. https://doi.org/10.1038/s41598-018-25437-1.

76. Liu, C.; Berry, P.M.; Dawson, T.P.; Pearson, R.G. Selecting thresholds of occurrence in the prediction of species distributions. *Ecography* **2005**, *28*, 385–393. https://doi.org/10.1111/j.0906-7590.2005.03957.x.

77. Shabani, F.; Kumar, L.; Ahmadi, M. Assessing accuracy methods of species distribution models: AUC, specificity, sensitivity and the true skill statistic. *Glob. J. Hum. -Soc. Sci. B Geogr. Geo-Sci. Environ. Sci. Disaster Manag.* **2018**, *18*, 1–13. Available online: https://globaljournals.org/GJHSS_Volume18/2-Assessing-Accuracy-Methods.pdf (accessed on).

78. Macleod, J. *Ixodes ricinus* in relation to its physical environment: II. The factors governing survival and activity. *Parasitology* **1935**, *27*, 123–144.

79. Alasmari, S.; Wall, R. Metabolic rate and resource depletion in the tick *Ixodes ricinus* in response to temperature. *Exp. Appl. Acarol.* **2021**, *83*, 81–93. https://doi.org/10.1007/s10493-020-00568-1.

80. Gethmann, J.; Hoffmann, B.; Kasbohm, E.; Süss, J.; Habedank, B.; Conraths, F.J.; Beer, M.; Klaus, C. Research paper on abiotic factors and their influence on *Ixodes ricinus* activity—Observations over a two-year period at several tick collection sites in Germany. *Parasitol. Res.* **2020**, *119*, 1455–1466. https://doi.org/10.1007/s00436-020-06666-8.

81. Zając, Z.; Kulisz, J.; Bartosik, K.; Woźniak, A.; Dzierżak, M.; Khan, A. Environmental determinants of the occurrence and activity of *Ixodes ricinus* ticks and the prevalence of tick-borne diseases in eastern Poland. *Sci. Rep.* **2021**, *11*, 15472. https://doi.org/10.1038/s41598-021-95079-3.

82. Wongnak, P.; Bord, S.; Jacquot, M.; Agoulon, A.; Beugnet, F.; Bournez, L.; Cèbe, N.; Chevalier, A.; Cosson, J.-F.; Dambrine, N.; et al. Meteorological and climatic variables predict the phenology of *Ixodes ricinus* nymph activity in France, accounting for habitat heterogeneity. *Sci. Rep.* **2022**, *12*, 7833. https://doi.org/10.1038/s41598-022-11479-z.

83. 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. https://doi.org/10.1186/1756-3305-6-1.

84. Fu, W.; Bonnet, C.; Septfons, A.; Figoni, J.; Durand, J.; Frey-Klett, P.; Rustand, D.; Jaulhac, B.; Métras, R. Spatial and seasonal determinants of Lyme borreliosis incidence in France, 2016 to 2021. *Eurosurveillance* **2023**, *28*, 2200581. https://doi.org/10.2807/1560-7917.ES.2023.28.14.2200581.

85. Qviller, L.; Grøva, L.; Viljugrein, H.; Klingen, I.; Mysterud, A. Temporal pattern of questing tick *Ixodes ricinus* density at differing elevations in the coastal region of western Norway. *Parasites Vectors* **2014**, *7*, 179. https://doi.org/10.1186/1756-3305-7-179.

86. Evstafiev, I.; Simferopol, U.C.S.-E.S. Results of a 30-years-long investigation of small mammals in Crimea. Part 3. Parasites and epizootiology. *Proc. Thériol. Sch.* **2017**, *2017*, 111–135. https://doi.org/10.15407/ptt2017.15.111.

87. Nekrasova, O.; Marushchak, O.; Redinov, K.; Pupins, M.; Čeirāns, A.; Skute, A.; Tytar, V.; Moysiyenko, I.; Theissinger, K.; Georges, J.-Y. Assessing ecocide impacts for developing a conservation strategy in Ukraine. In Proceedings of the World Biodiversity Forum, Davos, Switzerland, 16–21 June 2024; pp. 95–96.

88. Marushchak, O.; Nekrasova, O.; Zinenko, O.; Drohvalenko, M.; Mykytynets, H.; Suriadna, N.; Kotserzhynska, I.; Kotserzhynska, S.; Brusentsova, N.; Kuzmenko, Y.; et al. Herpetofauna at the frontline: So many ways to die. *Responsible Herpetoculture J.* **2024**, *14*, 114–128.

89. Moysiyenko, I.I.; Khodosovtsev, O.Y.; Vasyliuk, O.V.; Parkhomenko, V.V.; Rusin, M.Y.; Viter, S.H.; Kuzemko, A.A.; Drapalyuk, A.M.; Biatov, A.P.; Sadogurska, S.S.; et al. Consequences of the Russian terrorist attack on the Kakhovka HPP for wildlife. In *Traditions of Reserve Management, Modern Problems of Conservation and Post-War Restoration of Protected Areas: Collection of Scientific Papers Based on the Materials of the All-Ukrainian round Table Dedicated to the 160th Anniversary of Friedrich Falz-Fein*; Shapoval, V.V., Ed.; "Hlyboki Balyky", Balyko-Shchuchynka village, Ukraine; Series: "Conservation Biology in Ukraine"; Issue 32; Druk Art: Chernivtsi, Ukraine, 2023; pp. 158–163. (In Ukrainian)

90. Tarnas, M.C.; Desai, A.N.; Lassmann, B.; Abbara, A. Increase in vector-borne disease reporting affecting humans and animals in Syria and neighboring countries after the onset of conflict: A ProMED analysis 2003–2018. *Int. J. Infect. Dis.* **2021**, *102*, 103–109. https://doi.org/10.1016/j.ijid.2020.09.1453.

91. Boeckmann, M.; Joyner, T.A. Old health risks in new places? An ecological niche model for *I. ricinus* tick distribution in Europe under a changing climate. *Health Place* **2014**, *30*, 70–77. https://doi.org/10.1016/j.healthplace.2014.08.004.

92. Li, S.; Gilbert, L.; Harrison, P.A.; Rounsevell, M.D.A. Modelling the seasonality of Lyme disease risk and the potential impacts of a warming climate within the heterogeneous landscapes of Scotland. *J. R. Soc. Interface* **2016**, *13*, 20160140. https://doi.org/10.1098/rsif.2016.0140.

93. Ogden, N.H.; Ben Beard, C.; Ginsberg, H.S.; Tsao, J.I. Possible effects of climate change on Ixodid ticks and the pathogens they transmit: Predictions and observations. *J. Med. Entomol.* **2021**, *58*, 1536–1545. https://doi.org/10.1093/jme/tjaa220.

94. Voyiatzaki, C.; Papailia, S.I.; Venetikou, M.S.; Pouris, J.; Tsoumani, M.E.; Papageorgiou, E.G. Climate changes exacerbate the spread of *Ixodes ricinus* and the occurrence of Lyme Borreliosis and tick-borne encephalitis in Europe-How climate models are used as a risk assessment approach for tick-borne diseases. *Int. J. Environ. Res. Public Health* **2022**, *19*, 6516. https://doi.org/10.3390/ijerph19116516.
95. Bisanzio, D.; Amore, G.; Ragagli, C.; Tomassone, L.; Bertolotti, L.; Mannelli, A. Temporal variations in the usefulness of nor- malized 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. https://doi.org/10.1093/jmedent/45.3.547.

96. Bouchard, C.; Dibernardo, A.; Koffi, J.; Wood, H.; Leighton, P.A.; Lindsay, L.R. Increased risk of tick-borne diseases with climate and environmental changes. *Can. Commun. Dis. Rep.* **2019**, *45*, 83–89. https://doi.org/10.14745/ccdr.v45i04a02.

97. Tsoumani, M.E.; Papailia, S.I.; Papageorgiou, E.G.; Voyiatzaki, C. Climate change impacts on the prevalence of tick-borne diseases in Europe. *Environ. Sci. Proc.* **2023**, *26*, 18. https://doi.org/10.3390/environsciproc2023026018.

**Disclaimer/Publisher's Note:** The statements, opinions and data contained in all publications are solely those of the individual author(s) and contributor(s) and not of MDPI and/or the editor(s). MDPI and/or the editor(s) disclaim responsibility for any injury to people or property resulting from any ideas, methods, instructions or products referred to in the content.