Introduction
As anthropogenic climate change becomes increasingly apparent through rising mean annual temperatures and the significant increase in both the frequency and magnitude of extreme weather events, investigations of the relationships between vector populations and climate have garnered increased prominence (Kilpatrick and Randolph, Reference Kilpatrick and Randolph2012; Dantas-Torres, Reference Dantas-Torres2015; Semenza and Suk, Reference Semenza and Suk2018; Rocklöv and Dubrow, Reference Rocklöv and Dubrow2020; Gilbert, Reference Gilbert2021). Consequently, studying how vectors interact with their climatic environment, with a focus on mechanistic and correlative models, has gained additional importance for human health and welfare (Estrada-Peña, Reference Estrada-Peña2008; Li et al., Reference Li, Heyman, Cochez, Simons and Vanwambeke2012; Zhao et al., Reference Zhao, Wang, Fan, Ji, Liu, Zhang, Li, Zhou, Li, Liang, Liu, Yang and Fang2021).
Ixodes ricinus is one of the most important tick species in Europe because of its role as a vector for Lyme disease and tick-borne encephalitis (Jaenson and Lindgren, Reference Jaenson and Lindgren2011; Ostfeld and Brunner, Reference Ostfeld and Brunner2015; Estrada-Peña et al., Reference Estrada-Peña, Alexander and Wint2016). The abiotic factors that affect the abundance and distribution of I. ricinus have long been studied. Specifically, humidity and vegetation litter layers have been noted as the most important requirements for the survival and host-seeking activity of the species (Knülle and Rudolph, Reference Knülle, Rudolph, Obenchain and Galun1982; Gray et al., Reference Gray, Kahl, Lane, Levin and Tsao2016; Van Gestel et al., Reference Van Gestel, Matthysen, Heylen and Verheyen2022). Geographical variation in temperature is another significant factor that has been suggested to influence the timing of host-seeking activity and molting to the next life stage (Jaenson and Lindgren, Reference Jaenson and Lindgren2011; Gilbert et al., Reference Gilbert, Aungier and Tomkins2014). As a result of global warming, changes in precipitation patterns and increasing temperatures could increase the mortality rate of I. ricinus, particularly due to the potential increase in droughts and cold periods (Dautel et al., Reference Dautel, Kämmer, Kahl, Braks, Van Wieren, Takken and Sprong2016). Conversely, the occurrence of milder winters could create new distribution areas with altitudinal shifts, such as mountainous regions and higher latitudes in Europe (Materna et al., Reference Materna, Daniel, Metelka and Harčarik2008; Jaenson and Lindgren, Reference Jaenson and Lindgren2011; Martello et al., Reference Martello, Mannelli, Ragagli, Ambrogi, Selmi, Ceballos and Tomassone2014; Hvidsten et al., Reference Hvidsten, Frafjord, Gray, Henningsson, Jenkins, Kristiansen, Lager, Rognerud, Slåtsve, Stordal, Stuen and Wilhelmsson2020).
Ecological niche modelling (ENM) is a fundamental correlative method for investigating the distribution of disease transmission vectors, including ticks (Raghavan et al., Reference Raghavan, Heath, Lawrence, Ganta, Peterson and Pomroy2020; Alkishe et al., Reference Alkishe, Raghavan and Peterson2021; Moo-Llanes et al., Reference Moo-Llanes, de Oca-Aguilar, Romero-Salas and Sánchez-Montes2021). ENMs combine the locations of species with environmental variables to construct a representative multidimensional niche of the species (Warren and Seifert, Reference Warren and Seifert2011; Peterson and Soberón, Reference Peterson and Soberón2012). In addition to predicting the current potential niches of the species, ENMs are utilized to project possible future distributions under global circulation model (GCM) scenarios (Carvalho et al., Reference Carvalho, Rangel and Vale2017; Aguilar-Domínguez et al., Reference Aguilar-Domínguez, Moo-Llanes, Sánchez-Montes, Becker, Feria-Arroyo, de León and Romero-Salas2021; Alkishe and Peterson, Reference Alkishe and Peterson2022). Although most distribution models use macroclimatic parameters gathered from standardized weather stations, the conditions of the microclimate can differ greatly from those of the macroclimate (Maclean et al., Reference Maclean, Duffy, Haesen, Govaert, De Frenne, Vanneste, Lenoir, Lembrechts, Rhodes and Van Meerbeek2021; Marcin et al., Reference Marcin, Raschmanová, Miklisová and Kováč2021); recently, there have been significant advances in including microclimatic parameters in distribution modelling (Lembrechts et al., Reference Lembrechts, Nijs and Lenoir2019; Stark and Fridley, Reference Stark and Fridley2022). The microclimate might be more decisive than the macroclimate for small arthropods such as ticks because they are much more exposed to conditions close to the ground and under litter, and conditions in the microclimate could greatly affect the survival and activity of ticks (Randolph and Storey, Reference Randolph and Storey1999; Lauterbach et al., Reference Lauterbach, Wells, O'Hara, Kalko and Renner2013; Boehnke et al., Reference Boehnke, Gebhardt, Petney and Norra2017; Volk et al., Reference Volk, Lubelczyk, Johnston, Levesque and Gardner2022).
Several modelling studies have suggested possible changes in the distribution areas of I. ricinus in the future, which could increase the risk of dissemination of related diseases. The common conclusion drawn from these studies is that the range has expanded to a large part of Europe, North Africa and the Middle East (Porretta et al., Reference Porretta, Mastrantonio, Amendolia, Gaiarsa, Epis, Genchi, Bandi, Otranto and Urbanelli2013; Alkishe et al., Reference Alkishe, Peterson and Samy2017; Cunze et al., Reference Cunze, Glock, Kochmann and Klimpel2022). However, a significant portion of the studies included specimens from North Africa, which were later identified as Ixodes inopinatus (Estrada-Peña et al., Reference Estrada-Peña, Nava and Petney2014). Therefore, reevaluating previous occurrence records is extremely necessary. Additionally, occurrence records from countries such as Turkey, where I. ricinus was observed but generally was neglected in these analyses, will provide important information about the current and future distributions of the species. Additionally, most previous studies except Cunze et al. (Reference Cunze, Glock, Kochmann and Klimpel2022) on I. ricinus have used CMIP5 and earlier GCMs to predict possible future climatic suitabilities, and updated projections that use the latest GCMs will provide better predictions.
The primary goal of this study was to predict climatically favourable areas for I. ricinus under both current and future climate conditions via the Maxent approach by incorporating microclimatic and macroclimatic variables. Tick presence data were gathered from the current literature, and a Macroclimate suitability model was built using the CHELSA bioclimatic dataset, whereas a Mixed (microclimate and macroclimate) suitability model was built by combining Soiltemp and WorldClim datasets. The second aim was to estimate the possible future distributions under projected GCMs using suitability models. To this end, SSP3-7.0 and SSP5-8.5 GCMs for the 2011–2040, 2041–2070 and 2071–2100 periods, respectively, were used to carry out future projections with the Macroclimate model. The results of this study will be informative for evaluating the risk of tick-borne diseases and predicting future threats posed by emerging diseases transmitted by I. ricinus.
Materials and methods
Species distribution data
To infer the current and future distributions of I. ricinus, various resources have been utilized to gather occurrence data: VectorMap (www.vectormap.org), the Global Biodiversity Information Facility (www.gbif.org) and scientific literature (Boehnke et al., Reference Boehnke, Brugger, Pfäffle, Sebastian, Norra, Petney, Oehme, Littwin, Lebl, Raith, Walter, Gebhardt and Rubel2015; Estrada-Pena and De La Fuente, Reference Estrada-Pena and De La Fuente2016; Krawczyk et al., Reference Krawczyk, Röttjers, Fonville, Takumi, Takken, Faust and Sprong2022). Occurrence points from Turkey were primarily derived from Hekimoğlu (Reference Hekimoğlu2022), who utilized both morphological and molecular methods for sample identification. Additionally, field sampling data from 2023 (n = 2, İstanbul) were included in the final local dataset (n = 19).
The raw records were cleaned and refined by excluding localities (1) where tick collection was conducted from birds, as the presence of ticks on birds does not necessarily indicate that ticks establish populations in those areas; (2) samples from North Africa, as these have been redetermined as I. inopinatus (Rollins et al., Reference Rollins, Margos, Brachmann, Krebs, Mouchet, Dingemanse, Laatamna, Reghaissia, Fingerle, Metzler, Becker and Chitimia-Dobler2023), whereas I. inopinatus could be distributed in southern Spain and Portugal, the absence of clear evidence distinguishing them from I. ricinus led us to retain the geographical points of these countries in the analyses. (3) Locations from the Mediterranean and Aegean regions of Turkey were excluded because of a lack of confirmation in the scientific literature. Following these steps, a preliminary set of 4668 geographic location points was obtained.
Environmental variables
Microclimatic predictors from the Soiltemp (available at https://zenodo.org/records/7134169) dataset were used to build the first model (Lembrechts et al., Reference Lembrechts, Aalto, Ashcroft, De Frenne, Kopecký, Lenoir, Luoto, Maclean, Roupsard, Fuentes-Lillo, García, Pellissier, Pitteloud, Alatalo, Smith, Björk, Muffler, Ratier Backes, Cesarz, Gottschall, Okello, Urban, Plichta, Svátek, Phartyal, Wipf, Eisenhauer, Pușcaș, Turtureanu, Varlagin, Dimarco, Jump, Randall, Dorrepaal, Larson, Walz, Vitale, Svoboda, Higgens, Halbritter, Curasi, Klupar, Koontz, Pearse, Simpson, Stemkovski, Graae, Sørensen, Høye, Fernández Calzado, Lorite, Carbognani, Tomaselli, Forte, Petraglia, Haesen, Somers, Van Meerbeek, Björkman, Hylander, Merinero, Gharun, Buchmann, Dolezal, Matula, Thomas, Bailey, Ghosn, Kazakis, de Pablo, Kemppinen, Niittynen, Rew, Seipel, Larson, Speed, Ardö, Cannone, Guglielmin, Malfasi, Bader, Canessa, Stanisci, Kreyling, Schmeddes, Teuber, Aschero, Čiliak, Máliš, De Smedt, Govaert, Meeussen, Vangansbeke, Gigauri, Lamprecht, Pauli, Steinbauer, Winkler, Ueyama, Nuñez, Ursu, Haider, Wedegärtner, Smiljanic, Trouillier, Wilmking, Altman, Brůna, Hederová, Macek, Man, Wild, Vittoz, Pärtel, Barančok, Kanka, Kollár, Palaj, Barros, Mazzolari, Bauters, Boeckx, Alonso, Zong, Di Cecco, Sitková, Tielbörger, van den Brink, Weigel, Homeier, Dahlberg, Medinets, Medinets, De Boeck, Portillo-Estrada, Verryckt, Milbau, Daskalova, Thomas, Myers-Smith, Blonder, Stephan, Descombes, Zellweger, Frei, Heinesch, Andrews, Dick, Siebicke, Rocha, Senior, Rixen, Jimenez, Boike, Pauchard, Scholten, Scheffers, Klinges, Basham, Zhang, Zhang, Géron, Fazlioglu, Candan, Sallo Bravo, Filip Hrbacek, Laska, Cremonese, Haase, Moyano, Rossi and Nijs2020) with 30 arcsec (~1 km) resolution. The Soiltemp dataset includes bioclimatic temperature variables similar to those of the WorldClim and CHELSA datasets. However, unlike these datasets, which use interpolated data from standardized weather stations, the environmental data in the Soiltemp dataset are interpolated from both observational and experimental data gathered from dataloggers placed in microclimates (in soil) worldwide. Another advantage of Soiltemp is that the most intense data gathering was carried out in Europe, the main distributional area of I. ricinus. As the Soiltemp data included only temperature-related parameters (the first 11 bioclimatic variables), they were combined with precipitation parameters (bioclimatic vars 11–17), vapour pressure and surface radiation parameters from the WorldClim dataset (Fick and Hijmans, Reference Fick and Hijmans2017) available at https://www.worldclim.org to construct the input dataset. The WorldClim parameters were chosen for combination because CHELSA and Soiltemp have some compatibility issues in Maxent models. For the Macroclimate model and projections, data were gathered from the CHELSA dataset, which is a ready-to-use model output of estimates of ERA-Interim climatic reanalysis downscaled with the CHELSA algorithm and is available at https://CHELSAclimate.org (Karger et al., Reference Karger, Lange, Hari, Reyer, Conrad, Zimmermann and Frieler2022) with a 30 arcsec, (~1 km) resolution. Climatic variables for the 1981–2010 period were used to build the base Macroclimate model for the present. In addition to featuring measurements for a more recent period, the CHELSA dataset includes additional parameters that are biologically more relevant to Ixodes ticks, such as growth degree days above 5°C (Gdd5) (Jaenson and Lindgren, Reference Jaenson and Lindgren2011; Gray et al., Reference Gray, Kahl and Zintl2021), net primary productivity (Npp), which is related to vegetation cover important for ticks, growth season temperature (Gst) and growth season precipitation (Gsp). Additionally, since bioclimatic factors 8, 9, 18 and 19 were found to include some inconsistencies related to spatial artefacts first found in the WorldClim dataset and present in the CHELSA dataset (Escobar, Reference Escobar2020; Aguilar-Domínguez et al., Reference Aguilar-Domínguez, Moo-Llanes, Sánchez-Montes, Becker, Feria-Arroyo, de León and Romero-Salas2021; Alkishe and Peterson, Reference Alkishe and Peterson2022), these variables were not used in datasets (Soiltemp, WorldClim and CHELSA). As the Mixed model dataset also included some macroclimatic variables (precipitation variables), the term ‘Mixed’ was chosen for differentiation from the other dataset, which included only macroclimatic variables. For the two datasets, variables were clipped to the area of interest between latitudes −20° and 65° and longitudes 20° and 65° WGS84. Geographic computations were performed with QGIS (QGIS Geographic Information System, 2022), the GDAL library in Python (Open Source Geospatial Foundation, 2022), and R version 4.2.2 (R Core Team, 2022).
Future projections
As future projections are not available for microclimatic variables, GCMs available for the CHELSA dataset for three periods were used to construct future projections with the Macroclimate model selected for the CHELSA dataset: present-near future (2011–2040), near future (2041–2070) and distant future (2071–2100). The scenarios included SSP3-7.0 and SSP5-8.5 projections of UKESM1–0-LL (Sellar et al., Reference Sellar, Jones, Mulcahy, Tang, Yool, Wiltshire, O'Connor, Stringer, Hill, Palmieri, Woodward, de Mora, Kuhlbrodt, Rumbold, Kelley, Ellis, Johnson, Walton, Abraham, Andrews, Andrews, Archibald, Berthou, Burke, Blockley, Carslaw, Dalvi, Edwards, Folberth, Gedney, Griffiths, Harper, Hendry, Hewitt, Johnson, Jones, Jones, Keeble, Liddicoat, Morgenstern, Parker, Predoi, Robertson, Siahaan, Smith, Swaminathan, Woodhouse, Zeng and Zerroukat2019), MRI-ESM2-0 (Oshima et al., Reference Oshima, Yukimoto, Deushi, Koshiro, Kawai, Tanaka and Yoshida2020), MPI-ESM1-2 (Mauritsen et al., Reference Mauritsen, Bader, Becker, Behrens, Bittner, Brokopf, Brovkin, Claussen, Crueger, Esch, Fast, Fiedler, Fläschner, Gayler, Giorgetta, Goll, Haak, Hagemann, Hedemann, Hohenegger, Ilyina, Jahns, Jimenéz-de-la-Cuesta, Jungclaus, Kleinen, Kloster, Kracher, Kinne, Kleberg, Lasslop, Kornblueh, Marotzke, Matei, Meraner, Mikolajewicz, Modali, Möbis, Müller, Nabel, Nam, Notz, Nyawira, Paulsen, Peters, Pincus, Pohlmann, Pongratz, Popp, Raddatz, Rast, Redler, Reick, Rohrschneider, Schemann, Schmidt, Schnur, Schulzweida, Six, Stein, Stemmler, Stevens, von Storch, Tian, Voigt, Vrese, Wieners, Wilkenskjeld, Winkler and Roeckner2019), IPSL-CM6-LR (Boucher et al., Reference Boucher, Servonnat, Albright, Aumont, Balkanski, Bastrikov, Bekki, Bonnet, Bony, Bopp, Braconnot, Brockmann, Cadule, Caubel, Cheruy, Codron, Cozic, Cugnet, D'Andrea, Davini, de Lavergne, Denvil, Deshayes, Devilliers, Ducharne, Dufresne, Dupont, Éthé, Fairhead, Falletti, Flavoni, Foujols, Gardoll, Gastineau, Ghattas, Grandpeix, Guenet, Guez, Guilyardi, Guimberteau, Hauglustaine, Hourdin, Idelkadi, Joussaume, Kageyama, Khodri, Krinner, Lebas, Levavasseur, Lévy, Li, Lott, Lurton, Luyssaert, Madec, Madeleine, Maignan, Marchand, Marti, Mellul, Meurdesoif, Mignot, Musat, Ottlé, Peylin, Planton, Polcher, Rio, Rochetin, Rousset, Sepulchre, Sima, Swingedouw, Thiéblemont, Traore, Vancoppenolle, Vial, Vialard, Viovy and Vuichard2020) and GFDL-ESM4 (Dunne et al., Reference Dunne, Horowitz, Adcroft, Ginoux, Held, John, Krasting, Malyshev, Naik, Paulot, Shevliakova, Stock, Zadeh, Balaji, Blanton, Dunne, Dupuis, Durachta, Dussin, Gauthier, Griffies, Guo, Hallberg, Harrison, He, Hurlin, McHugh, Menzel, Milly, Nikonov, Paynter, Ploshay, Radhakrishnan, Rand, Reichl, Robinson, Schwarzkopf, Sentman, Underwood, Vahlenkamp, Winton, Wittenberg, Wyman, Zeng and Zhao2020). Shared socioeconomic pathway (SSP) scenarios improve representative concentration pathways (RCPs) by including important socioeconomic parameters, such as urbanization, population growth and climate change (Hausfather, Reference Hausfather2020). These datasets are also available at https://CHELSA-climate.org.
Ecological niche models
Locations falling outside environmental rasters were removed before building the models. To reduce spatial sampling bias and to remove any unnoticed duplicates, the locations of I. ricinus were thinned to 10, 20, 30, 40 and 50 km, and five different location sets were obtained via the spThin package for R (Aiello-Lammens et al., Reference Aiello-Lammens, Boria, Radosavljevic, Vilela and Anderson2015). The resulting datasets were split into a training set (70%), a cross-validation set for evaluating the candidate models (25%), and an independent test set for evaluating the final models (5%), which is crucial for revealing the real predictive power of a model. To better predict the distribution of a species, an accessible M space is needed (Soberon and Peterson, Reference Soberon and Peterson2005; Barve et al., Reference Barve, Barve, Jiménez-Valverde, Lira-Noriega, Maher, Peterson, Soberón and Villalobos2011). To construct the M space, a minimum convex polygon of a 100 km buffer area around the location points was created for the models using the ellipsenm package for R, which is available at https://github.com/marlonecobos/ellipsenm (Cobos et al., Reference Cobos, Townsend Peterson, Barve and Osorio-Olvera2019), and environmental rasters were cut to this buffer area. For the Mixed model and the Macroclimate model, three sets were created for each of them by setting 0.9, 0.8 and 0.7 intervariable correlation thresholds using the vifcor function in the usdm package for R (Naimi, Reference Naimi2017), which eliminates variables from correlated pairs depending on the variance inflation factor (VIF), and by consulting the evaluations of the jackknife results of the preliminary Maxent models. An additional dataset was prepared for the Macroclimate model, depending on the previous models and jackknife results, for a total of three variable sets for the Mixed model and four sets for the Macroclimate model (Table 1).
Table 1. The environmental predictors used in different sets in the microclimate model and macroclimate model are explained in the first column
Data Sources.
a Soiltemp.
b CHELSA.
c WorldClim.
ENMs were built according to the maximum entropy algorithm (Phillips et al., Reference Phillips, Anderson and Schapire2006) and with the Kuenm R package (Cobos et al., Reference Cobos, Townsend Peterson, Barve and Osorio-Olvera2019), which includes Maxent version 3.4.4 (Steven et al., Reference Steven, Miroslav and Schapire2021). The first calibration models were built with combinations of linear, quadratic, product and hinge features (L, Q, P, H) with regularization parameters of 0.1, 0.2, 0.3, 0.4, 0.5, 0.6, 0.7, 0.8, 0.9, 1, 2, 3, 4 and 5. To evaluate the models, a minimum threshold of 0.7 was set for the area under the curve (AUC), a maximum threshold of 0.05 was set for the omission rate, and a partial receiver operating characteristic (pROC) test was used (Peterson and Soberón, Reference Peterson and Soberón2012; Aguilar-Domínguez et al., Reference Aguilar-Domínguez, Moo-Llanes, Sánchez-Montes, Becker, Feria-Arroyo, de León and Romero-Salas2021). Further postanalysis model selection was carried out by the Akaike information criterion corrected for small sample sizes (AICc), and models with ΔAICc values less than 2 were selected as the final models (Warren and Seifert, Reference Warren and Seifert2011; Nuñez-Penichet et al., Reference Nuñez-Penichet, Osorio-Olvera, Gonzalez, Cobos, Jiménez, DeRaad, Alkishe, Contreras-Díaz, Nava-Bolaños and Utsumi2021). These final models were used to construct the final projections with 10 bootstrap replicates, and they were also tested with the independent testing locations that were set aside before model training and cross-validation.
Raster plots, which show the suitability of the distributions, were binarized using the 10th percentile training presence logistic threshold, which was taken from the Maxent results (Liu et al., Reference Liu, Berry, Dawson and Pearson2005; Kramer-Schadt et al., Reference Kramer-Schadt, Niedballa, Pilgrim, Schröder, Lindenborn, Reinfelder, Stillfried, Heckmann, Scharf, Augeri, Cheyne, Hearn, Ross, Macdonald, Mathai, Eaton, Marshall, Semiadi, Rustam, Bernard, Alfred, Samejima, Duckworth, Breitenmoser-Wuersten, Belant, Hofer and Wilting2013). For the future projections, 6 consensus maps were calculated by taking the median of the GCM scenarios for SSP3-7.0 and SSP5-8.5 for the periods 2011–2040, 2041–2070 and 2071–2100, respectively. To minimize possible deviations due to the novelty of the conditions compared with the original distribution area, the models were transferred to future conditions without any projection under extreme novel conditions (no extrapolation in the GCM scenarios) (Cobos et al., Reference Cobos, Townsend Peterson, Barve and Osorio-Olvera2019). For the Mixed and Macroclimate models, the range (max val-min val) of the suitability maps of 10 bootstrap replications was used to create uncertainty maps. To map the uncertainty for future projections, we calculated the range of suitabilities between the different GCMs. Additionally, to determine the places where conditions are more extreme than in the calibration area of the models, a mobility-oriented parity (MOP) analysis was carried out, and the extrapolation risk of transfer regions was calculated with the nearest 10% reference via the MOP package for R (Owens et al., Reference Owens, Campbell, Dornak, Saupe, Barve, Soberón, Ingenloff, Lira-Noriega, Hensz, Myers and Peterson2013; Alkishe et al., Reference Alkishe, Cobos, Peterson and Samy2020; Flores-López et al., Reference Flores-López, Moo-Llanes, Romero-Figueroa, Guevara-Carrizales, López-Ordoñez, Casas-Martínez and Samy2022). Three types of MOP maps were created, (1) representing how many variables in the region of projection are non-analogous to those in the calibration area and where non-analogous conditions were found towards (2) high values and (3) low values of each variable were created for each GCM. All the maps were created with QGIS software and R software.
Results
Model parameters
As models that used locations thinned to a greater distance than 20 km had AUC values less than 0.7, models created with 20 km thinned locations were selected and had a total of n = 1001 locations (Fig. 1). Among these locations, ~70% (n = 713) were used to train the models, ~25% (n = 238) were used for cross-validation, and ~5% (n = 50) were set aside for the final evaluation tests. For the Mixed model, out of all 630 models, 627 were statistically significant, while only 24 passed the 0.05 omission rate criterion, and only one model from these fulfilled the AICc. The final model was a linear product (LP) model with a regularization multiplier of 0.5; this model used the second set, with intercorrelations less than 0.8 and had an omission rate of 0.046, a pROC value of ~0 and an AUC ratio of 1.16. The final evaluations performed on the independent test set showed an omission rate of 0.06, which is marginally above 0.05; on the other hand, this means that the Mixed model correctly classified approximately 47 points out of 50 points from the independent test locations. The contributions of the 6 environmental variables to the final Mixed model are shown in Table 2. The most important environmental variable was Bio4, the temperature seasonality near the soil, with a contribution of 40.6%, followed by vapour pressure, which contributed 36%, and Bio14 (precipitation amount in the wettest month) and Bio15 (precipitation seasonality), both of which contributed 9.6%. Bio6, the near-soil minimum temperature of the coldest month, had a 3.3% contribution.
Figure 1. All occurrence points of Ixodes ricinus after cleaning and thinning.
Table 2. Percent contribution of the environmental predictors to the models
Data Sources.
a Soiltemp.
b CHELSA.
c WorldClim.
For the Macroclimate dataset, the AUC also fell below 0.7 beyond a 20 km thinning length in the models generated using the CHELSA 1981–2010 dataset; thus, the 20 km thinned locations were also used in the models. Out of the 840 candidate models, 757 were statistically significant, 49 passed the 0.05 omission rate criterion, and only one fulfilled both the omission rate and AICc conditions. The final Macroclimate model was a linear quadratic product hinge (LQPH) model combining all features with a regularization multiplier of 2 that used the 4th variable set (Table 2), which had a 0.7 upper intercorrelation threshold. The model had an omission rate of 0.048 and a pROC of ~0, with an AUC of 1.14. Final evaluations with the independent test revealed that the Macroclimate model had better results in predicting the independent test results, with a final evaluation omission rate of 0.041, which is even better than the cross-validation omission rate of the same model. Although this result is better than that of the final evaluations of the Mixed model, the values are actually very close, as the 0.041 omission rate means that it correctly assigned 48 out of the 50 independent test locations correctly. According to the Macroclimate model, the largest contributor was Bio7 (annual range of air temperature), which contributed 52% of the total contribution. This was followed by Gsp (growth season precipitation), which contributed 26.5%, whereas Bio3 (isothermality) contributed 5.6%.
Current and future predictions
For the current prediction maps, the binary thresholds were 0.27 and 0.28, respectively, with very similar 10th percentile training presence logistic threshold values for both the Mixed and Macroclimate models. The predicted areas of suitability and uncertainty are shown in Fig. 2. Both the Mixed and Macroclimate models predicted a wide distributional area that covers much of Europe except for gaps in mountainous areas and parts of southern Europe, especially Spain and Greece. The climatically suitable regions mostly coincided with the current distribution of I. ricinus. On the other hand, in western Asia, suitable regions are located on the coast of the Black Sea. The Mixed model revealed a wider suitable region in the Caucasus, whereas the Macroclimate model for 1981–2010 predicted a smaller suitable region in the southern Caspian Sea. Both models also predicted a small suitable region on the coast of North Africa, with the Mixed model showing a slightly larger area in North Africa. On the other hand, the level of uncertainty is much greater in southern Europe and North Africa.
Figure 2. Maps of predicted suitable areas for I. ricinus and uncertainty from the ENM results. (A) Red areas show the suitable regions under current conditions according to the Macroclimate model. (B) Uncertainty values of the Macroclimate model. (C) Green areas show the suitable regions under current conditions according to the Mixed model. (D) Uncertainty values of the Mixed model.
The consensuses of the different GCM future projections are shown in Figs 3 and 4, and each individual future projection is included in Supplementary File 1. The predictions of SSP3-7.0 for the 2011–2040 period reveal significant expansion in eastern and northern Europe, whereas there are some declines in southern Europe, the Balkans and the Black Sea region. Interestingly, despite being a more extreme scenario, the SSP5-8.5 projection for 2011–2040 showed a smaller area of expansion, while the pattern of decline was similar to that of the SSP3-7.0 scenario. The level of uncertainty increases from western Europe to the peripheral regions where most of the expansion and decline occurs. Future projections for 2041–2070 suggest the continuation of increasing suitability in northeastern Europe, with the SSP3-7.0 scenario showing a spread that is slightly greater than that of SSP5-8.5. Conversely, the decreases in the southern region are greater in the SSP5-8.5 projection for 2041–2070. Distant future projections for 2071–2100 show that this trend of decreased suitability in south and increase in the north is expected to continue in both scenarios, with the SSP5-8.5 projection depicting a very extreme decline in the south, where most of the Mediterranean and Balkans will not be suitable for I. ricinus.
Figure 3. Side-by-side predicted suitable areas for I. ricinus and uncertainty values for the median of 5 GCM scenarios for SSP3-7.0 with differing degrees of loss and gain compared with the current conditions (1981–2010) for the Macroclimate model.
Figure 4. Side-by-side predicted suitable areas for I. ricinus and uncertainty values for the median of 5 GCM scenarios for SSP5-8.5 with differing degrees of loss and gain compared with the current conditions (1981–2010) for the Macroclimate model.
Uncertainty and MOP results
For the present conditions, the uncertainty of the mixed model is relatively low within the main distributional area of I. ricinus, whereas the highest levels of uncertainty are observed at the margins (Fig. 2), particularly in the eastern Mediterranean, southern Caucasus (Georgia), North African Coast, southwestern Spain and western Scandinavia. Uncertainty attains the highest values in the eastern Mediterranean and North Africa, outside the distribution of I. ricinus. Therefore, the high suitability results in these areas are not very reliable. For the Macroclimate model, although the maximum level of uncertainty is greater than that of the Mixed model, the overall level of uncertainty is lower, as shown in Fig. 2. The general pattern of uncertainty is similar and increases towards the margins. Specifically, east of the Mediterranean, North Africa, west of Spain and western Scandinavia present higher levels of uncertainty. Additionally, the uncertainty in future projections is greater than that in present projections. However, the pattern of uncertainty remains similar to that at present, with western Europe (including France, the United Kingdom and parts of Germany) experiencing the lowest uncertainty levels, whereas higher uncertainty is evident around the Mediterranean basin, the Balkans, western Spain and western Scandinavia. The general level of uncertainty is very high in the 2071–2100 projections, particularly in the SSP5-8.5 projections; additionally, there are more areas with extreme uncertainty, suggesting that distant future projections (especially the worst-case scenarios) are not very reliable for tick distribution models (see Figs 3 and 4).
The MOP results closely mirror other uncertainties and are included in Supplementary File 2. In the 2011–2040 and 2041–2070 SSP3-7.0 and SSP5-8.5 future projections, areas with significant extreme extrapolation begin at the south of Anatolia to Iran and the Middle East, North Africa and east of the Caspian Sea. Almost all regions with significant extrapolations fall outside the current and predicted future distribution range of I. ricinus, except for small areas in southern Anatolia and North Africa. For the 2071–2100 projections, in accordance with the uncertainty maps, the areas with extreme extrapolation are also dispersed to the coast of Mediterranean regions, including Spain, Italy, Greece and Turkey, which is more prominent in the SSP5-8.5 projection.
Discussion
Ixodes ricinus has high adaptability to European environments and a tendency to expand its geographic range towards northern altitudes, which is mostly explained by climate change (Materna et al., Reference Materna, Daniel, Metelka and Harčarik2008; Jaenson and Lindgren, Reference Jaenson and Lindgren2011; Martello et al., Reference Martello, Mannelli, Ragagli, Ambrogi, Selmi, Ceballos and Tomassone2014; Garcia-Vozmediano et al., Reference Garcia-Vozmediano, Krawczyk, Sprong, Rossi, Ramassa and Tomassone2020; Hvidsten et al., Reference Hvidsten, Frafjord, Gray, Henningsson, Jenkins, Kristiansen, Lager, Rognerud, Slåtsve, Stordal, Stuen and Wilhelmsson2020). This raises the question of whether this expansion is likely to give rise to new, distant foci for Lyme and TBE diseases. Environmental parameters play a significant role in influencing the distribution of I. ricinus. Factors such as climatic conditions, habitat characteristics, and landscape heterogeneity have been identified as key determinants affecting the abundance and prevalence of this tick species (Krasnov et al., Reference Krasnov, Stanko and Morand2007; Ruiz-Fons et al., Reference Ruiz-Fons, Fernández-de-Mera, Acevedo, Gortázar and De La Fuente2012; Hauser et al., Reference Hauser, Rais, Morán Cadenas, Gonseth, Bouzelboudjen and Gern2018).
The Mixed model identified temperature seasonality as the primary environmental predictor, whereas the Macroclimate model indicated that the annual temperature range was the most significant contributor. These findings are consistent with those of previous studies. The annual temperature range, annual mean temperature and annual precipitation were identified as the most influential factors by Alkishe et al. (Reference Alkishe, Peterson and Samy2017). Similarly, temperature seasonality is considered a relevant factor according to Cunze et al. (Reference Cunze, Glock, Kochmann and Klimpel2022). The second most influential factor is vapour pressure and growth season precipitation according to the Mixed and Macroclimate models, respectively. The off-host survival of ticks depends strongly on water availability since desiccation is one of the most prominent causes of tick mortality and decreasing questing activity (Perret et al., Reference Perret, Rais and Gern2004; Tagliapietra et al., Reference Tagliapietra, Rosà, Arnoldi, Cagnacci, Capelli, Montarsi, Hauffe and Rizzoli2011; Nolzen et al., Reference Nolzen, Brugger, Reichold, Brock, Lange and Thulke2022; Van Gestel et al., Reference Van Gestel, Matthysen, Heylen and Verheyen2022). Factors such as higher soil moisture content and increased cloud cover have been linked to increased questing activity in I. ricinus (Medlock et al., Reference Medlock, Pietzsch, Rice, Jones, Kerrod, Avenell, Los, Ratcliffe, Leach and Butt2008; Lauterbach et al., Reference Lauterbach, Wells, O'Hara, Kalko and Renner2013). Additionally, the distribution of I. ricinus has expanded in the past three decades due to more favourable biotic and abiotic conditions, which can be influenced by changes in precipitation and an increase in mean winter temperatures (Gray and Ogden, Reference Gray and Ogden2021). Other modelling studies have underscored the importance of precipitation as a significant factor in the distribution of this species in both current and future projections (Alkishe et al., Reference Alkishe, Peterson and Samy2017; Cunze et al., Reference Cunze, Glock, Kochmann and Klimpel2022). Temperature has been recognized in various studies as a crucial determinant of the completion of the life cycle and duration of the host-seeking activity of I. ricinus at northern latitudes (Lindgren et al., Reference Lindgren, Tälleklint and Polfeldt2000; Schwarz et al., Reference Schwarz, Hönig, Vavrušková, Grubhoffer, Balczun, Albring and Schaub2012; Hauser et al., Reference Hauser, Rais, Morán Cadenas, Gonseth, Bouzelboudjen and Gern2018). However, Gst (growth season temperature) made a very small contribution to the final Macroclimate models, and Gdd5 (degree days above 5°C) was not included in the final selected models. Another discrepancy with these previous studies is that the minimum temperature of the coldest month (Bio6) contributed less than other variables in the Mixed model; moreover, in the final selected models for Macroclimate, Bio6 was not included because of correlation thresholds. However, when Bio6 was included in the side models, Bio6 was again behind other parameters, with a contribution smaller than 3%. This might be due to the difference in the CHELSA and Soiltemp datasets used in creating the models. Additionally, the relationship between minimum temperature and tick survival might be more complicated because of the buffering effect of microclimatic factors such as snow cover and leaf litter on tick survival and activity (Van Gestel et al., Reference Van Gestel, Matthysen, Heylen and Verheyen2022). For example, in situ measurements have shown that the humidity and temperature in the understorey of forests where I. ricinus resides are indeed much more mild and stable than those measured from standard weather stations (Boehnke et al., Reference Boehnke, Brugger, Pfäffle, Sebastian, Norra, Petney, Oehme, Littwin, Lebl, Raith, Walter, Gebhardt and Rubel2015). Additionally, the increased winter mortality effect of removing leaf litter and snow cover on Ixodes scapularis was documented in experimental field plots (Volk et al., Reference Volk, Lubelczyk, Johnston, Levesque and Gardner2022).
It might be expected that microclimatic parameters are more accurate for predicting the distribution of ticks because of the proven strong connection of tick survival with microclimatic abiotic factors in the environment; however, the omission rates of the distribution models are very similar to those of the Macroclimate model, which has a very slight edge. This could be attributed to several factors. One explanation might be that the two microclimatic parameters that were included in the final mixed model may not be good predictors; however, on the contrary, Bio4 was the greatest contributor to the model, which supports the predictive power of macroclimatic parameters. Notably, the Mixed model also had very good predictive results. A more plausible explanation might be that the Macroclimate model included additional important parameters that improved the power of the model. Gsp was the second most important contributor to the Macroclimate model, which is not surprising given that most tick species are active in the growing season. In addition, the second most important parameter for the Mixed model was water vapour pressure from WorldClim, a macroclimatic dataset. The Soiltemp dataset includes only temperature-related variables, and precipitation is the same for macro- and microclimate. Thus, the inclusion of additional parameters related to water content and humidity (crucial factors for ticks) at the microclimate level would improve these models since relative humidity often varies with microclimate, at least as much as temperature (Van Gestel et al., Reference Van Gestel, Matthysen, Heylen and Verheyen2022).
The models indicate that the region climatically suitable for I. ricinus covers a large portion of Europe except southern Spain and Greece and the coastal regions of the Black Sea in West Asia. This outcome is consistent with the results of previous models of I. ricinus. On the other hand, model predictions revealed a more constrained distribution in the southern regions of Europe, especially in southern and western Anatolia; this is an expected result, as we know that there is presently no record of I. ricinus in this region (Hekimoğlu, Reference Hekimoğlu2022). The model also revealed a small area of suitability on the North African coast; however, the suitability of this region is not reliable due to the high degree of uncertainty. Although some of the previous studies included I. ricinus from North Africa, there is a high probability that some of these records involved I. inopinatus, as recent genomic analysis clearly demonstrated the presence of I. inopinatus in North Africa (Younsi et al., Reference Younsi, Fares, Cherni, Dachraoui, Barhoumi, Najjar and Zhioua2020; Rollins et al., Reference Rollins, Margos, Brachmann, Krebs, Mouchet, Dingemanse, Laatamna, Reghaissia, Fingerle, Metzler, Becker and Chitimia-Dobler2023). Nevertheless, the absence of current and future distributions of I. ricinus in the Mediterranean region cannot be ascribed solely to I. inopinatus. Although the distribution of I. inopinatus in different European countries, such as Germany, Spain and Portugal, has been suggested (Estrada-Peña et al., Reference Estrada-Peña, Nava and Petney2014; Chitimia-Dobler et al., Reference Chitimia-Dobler, Rieß, Kahl, Wölfel, Dobler, Nava and Estrada-Peña2018; Hauck et al., Reference Hauck, Springer, Pachnicke, Schunack, Fingerle and Strube2019), genomic data have shown that German samples are in fact I. ricinus (Rollins et al., Reference Rollins, Margos, Brachmann, Krebs, Mouchet, Dingemanse, Laatamna, Reghaissia, Fingerle, Metzler, Becker and Chitimia-Dobler2023). In this case, the taxonomic status of Mediterranean populations should be reevaluated, as this could influence the projections of I. ricinus in these areas. Considering the taxonomic reevaluation mentioned above, many recent studies excluded North African populations from their dataset (Cunze et al., Reference Cunze, Glock, Kochmann and Klimpel2022; Noll et al., Reference Noll, Wall, Makepeace, Newbury, Adaszek, Bødker, Estrada-Peña, Guillot, Da Fonseca, Probst, Overgaauw, Strube, Zakham, Zanet and Rose Vineer2023). The findings of these studies largely align with our current and future predictions. A large part of Europe constitutes the distribution area of the species, with the most suitable areas located in Western Europe. While Mediterranean countries such as France and Italy include suitable areas for the species, Portugal and the northern part of Spain also appear suitable in the western Mediterranean. Along the Black Sea coast, the most suitable areas include coastal regions encompassing the northern part of Turkey.
Future projections indicate a spread in northern and eastern Europe, confirming other projections with different bioclimatic datasets (Cunze et al., Reference Cunze, Glock, Kochmann and Klimpel2022; Noll et al., Reference Noll, Wall, Makepeace, Newbury, Adaszek, Bødker, Estrada-Peña, Guillot, Da Fonseca, Probst, Overgaauw, Strube, Zakham, Zanet and Rose Vineer2023). A consensus is also reached with areas that will not be suitable in the future, particularly for the predicted habitat loss in Spain, Greece and the Balkans (Cunze et al., Reference Cunze, Glock, Kochmann and Klimpel2022). Another outcome is that in the present–near-future scenarios (2011–2040), the midway scenario (SSP3-7.0) estimates a slightly wider new distributional area for I. ricinus in northeastern Europe than does SSP5-8.5, a worse-case scenario. This outcome is compatible with the previous future projections by Porretta et al. (Reference Porretta, Mastrantonio, Amendolia, Gaiarsa, Epis, Genchi, Bandi, Otranto and Urbanelli2013) and also Alkishe et al. (Reference Alkishe, Peterson and Samy2017), where midway scenarios also predicted a wider new region of suitability in the near future. While the spread in the 2041–2070 scenarios is similar between the SSP3-7.0 and SSP5-8.5 scenarios, the decline in the southern distributional areas, including the Balkans and Mediterranean regions, is much more prominent in SSP5-8.5 than in the other regions, which is most likely due to the predictions of great declines in precipitation in this region. Distant future scenarios (2071–2100) continue this increase in suitability towards the north, whereas the decline in southern regions becomes much more excessive (especially in the SSP5-8.5 scenario). On the other hand, the level of uncertainty and the degree of extrapolation are also extreme in these regions for this period, and thus, it would be inaccurate to reach any assumption for this period. Several climatic studies in recent decades since the early 2000s have deemed the Mediterranean basin a climate change hotspot, with projections showing increased temperature and aridity, heightened vulnerability to drought and high temperatures and a greater frequency of heat waves (Ulbrich et al., Reference Ulbrich, May, Li, Lionello, Pinto, Somot, Lionello, Malanotte-Rizzoli and Boscolo2006; Giorgi and Lionello, Reference Giorgi and Lionello2008; Naumann et al., Reference Naumann, Alfieri, Wyser, Mentaschi, Betts, Carrao, Spinoni, Vogt and Feyen2018). Furthermore, the impacts of recent climate change on increased drought frequency and magnitude have already been documented in Mediterranean-type climates (Hoerling et al., Reference Hoerling, Eischeid, Perlwitz, Quan, Zhang and Pegion2012; Feng et al., Reference Feng, Thompson, Woods and Porporato2019). Current projections indicate that the most suitable distribution areas for the species in Turkey are the coasts of the Black Sea and the Thrace region. However, while smaller suitable areas are also indicated along the Mediterranean coast of the country, the uncertainty in these predicted areas is again quite high. Future projections suggest a loss of suitability on the southern coast and, additionally, decline in certain localities, especially in the inland areas of Thrace in the north. This pattern indicates that expected climate changes in the Mediterranean Basin in the future will result in a decrease in the distribution areas of this species in Turkey, similar to the projections for Spain.
Ecological niche models are valuable tools for assessing the possible presence and future distribution of parasite vectors; however, several caveats exist, and special care is needed when building distribution models. Major challenges are the potential bias and limitations in the outputs of niche models, especially when projecting to future scenarios (Peterson et al., Reference Peterson, Cobos and Jiménez-García2018); therefore, providing the necessary uncertainty values associated with these models is crucial. Additionally, ectoparasites, including ticks, partially depend on their hosts for survival and transportation, posing other risks in modelling. Host species and their abundance are among the most important factors influencing the distribution of ticks to new areas and establishing populations. Moreover, the large number of migrating birds increases the probability of the geographic spread of I. ricinus and related diseases, emphasizing the importance of considering avian hosts in disease epidemiology (Waldenström et al., Reference Waldenström, Lundkvist, Falk, Garpmo, Bergström, Lindegren, Sjöstedt, Mejlon, Fransson, Haemig and Olsen2007; Ciebiera et al., Reference Ciebiera, Jerzak, Nowak-Chmura and Bocheński2019). Although hosts play a crucial role in the transmission of ticks, several limitations are associated with including hosts in modelling studies. For example, I. ricinus is a three-host species that feeds on different hosts during its life cycle (Hofmeester et al., Reference Hofmeester, Coipan, Van Wieren, Prins, Takken and Sprong2016), which can introduce complexity when multiple hosts are used in modelling studies. Recent papers have presented contrasting hypotheses regarding tick–host associations. Some suggest that ticks select hosts on the basis of the environment, whereas others propose that ticks select environments and feed on any available host within those environments (Nava and Guglielmone, Reference Nava and Guglielmone2013; Zhang et al., Reference Zhang, Zhang and Liu2019; Ginsberg et al., Reference Ginsberg, Hickling, Pang, Tsao, Fitzgerald, Ross, Rulison and Burke2022; Estrada-Peña et al., Reference Estrada-Peña, Guglielmone and Nava2023). Moreover, while it has been suggested that hosts can influence tick distribution at a smaller geographical scale, the distribution of ticks over a wide spatial scale is determined primarily by direct climatic effects rather than by host presence (Cumming, Reference Cumming2002). Considering all these factors, understanding the interactions between hosts, the environment and ticks is essential for developing effective strategies to manage tick-borne diseases and control their spread. On the other hand, a very recent study by Fabri et al. (Reference Fabri, Hofmeester, Ecke, Sprong, Timmermans, Heesterbeek and Cromsigt2024) revealed that I. ricinus depends much more on abiotic factors than on the composition of its hosts.
Conclusion
This study aimed to update previous projections by reviewing coordinates, utilizing new datasets and additionally implementing microclimatic parameters to construct assistive models that complement other distribution models that rely on macroclimatic parameters. According to the results, the Mixed model, combining micro and macro parameters achieved similar levels of accuracy with the Macroclimate model. However, the integration of microclimates in distribution models is still in its starting phase, and the inclusion of additional microclimatic parameters in available datasets, updating the presently available variables with new data and also widening these models to other species will be beneficial. For species distribution models, we still depend on macroclimatic datasets since predictions of very fine-grained parameters such as the microclimate in the future are not yet available, and macroclimate data still provide good predictions at broader scales. Therefore, it is essential to continuously replicate, validate and update previous models to better predict the possible future distributions of disease-transmitting ticks.
Supplementary material
The supplementary material for this article can be found at https://doi.org/10.1017/S003118202400132X
Data availability statement
The coordinates used in this study are available at Mendeley Data. https://data.mendeley.com/datasets/98m45vcfp2/1
Acknowledgements
The numerical calculations reported in this paper were partially performed at TUBITAK ULAKBIM, High Performance and Grid Computing Center (TRUBA resources). We also thank three anonymous reviewers for reviewing the manuscript and for their helpful suggestions.
Author contributions
O. H. and A. C. K. designed the study; O. H. prepared and cleaned the data; A. C. K. carried out the analysis and interpreted the results; A. C. K. prepared the tables and maps; A. C. K. and O.H. prepared the original draft and reviewed and edited the manuscript.
Financial support
This research received no specific grant from any funding agency or from the commercial or not-for-profit sectors.
Competing interest
None.
Ethical standards
Not applicable.
Open access
