Introduction
Human intervention in ecosystems, including globalization, urbanization, rapid transportation, land use and climate change, has led to a global increase in the range, distribution and transmission rate of numerous pathogens and vector organisms that serve as reservoirs for these pathogens (Kovats et al., Reference Kovats, Campbell-Lendrum, McMichael, Woodward and Cox2001; Andersen and Davis, Reference Andersen and Davis2017; Carvalho et al., Reference Carvalho, Rangel and Vale2017; Semenza and Suk, Reference Semenza and Suk2018; Aguilar-Domínguez et al., Reference Aguilar-Domínguez, Moo-Llanes, Sánchez-Montes, Becker, Feria-Arroyo, de León and Romero-Salas2021; Leder et al., Reference Leder, Openshaw, Allotey, Ansariadi, Barker, Burge, Clasen, Chown, Duffy, Faber, Fleming, Forbes, French, Greening, Henry, Higginson, Johnston, Lappan, Lin, Luby, McCarthy, O'Toole, Ramirez-Lovering, Reidpath, Simpson, Sinharoy, Sweeney, Taruc, Tela, Turagabeci, Wardani, Wong and Brown2021). While the most severe impacts of climate change are expected to affect ectotherm populations in the tropics, population growth rates of insects and other arthropods are also anticipated to increase in mid-to-high latitude regions due to the influence of warmer temperatures on growth rates (Deutsch et al., Reference Deutsch, Tewksbury, Huey, Sheldon, Ghalambor, Haak and Martin2008; Bonebrake and Deutsch, Reference Bonebrake and Deutsch2012; Rocklöv and Dubrow, Reference Rocklöv and Dubrow2020). In addition to the effects of temperature on the biology of disease-transmitting vector organisms, increased precipitation and extreme precipitation events in some regions may lead to an increase in vector abundance, whereas longer dry periods might increase tick mortality (Githeko et al., Reference Githeko, Lindsay, Confalonieri and Patz2000; Rocklöv and Dubrow, Reference Rocklöv and Dubrow2020).
Ticks are the vectors of numerous critical pathogens and pose a significant threat to animal and public health (Jongejan and Uilenberg, Reference Jongejan and Uilenberg2004). The species Hyalomma marginatum Koch, 1844 presents a substantial health threat as the main vector of Crimean–Congo haemorrhagic fever (CCHF) and also the vector of babesiosis and Rickettsia in the Mediterranean region, Africa and Asia (Ergönül, Reference Ergönül2006; Vatansever et al., Reference Vatansever, Estrada-Pena and Ergonul2007; Gale et al., Reference Gale, Estrada-Peña, Martinez, Ulrich, Wilson, Capelli, Phipps, De La Torre, Muñoz, Dottori, Mioulet and Fooks2010; Ros-García et al., Reference Ros-García, M'ghirbi, Bouattour and Hurtado2011; Bonnet et al., Reference Bonnet, Vourc'h, Raffetin, Falchi, Figoni, Fite, Hoch, Moutailler and Quillery2022; Sultankulova et al., Reference Sultankulova, Shynybekova, Kozhabergenov, Mukhami, Chervyakova, Burashev, Zakarya, Nakhanov, Barakbayev and Orynbayev2022). Effective and sustainable control of ticks, like many arthropod pests and vectors, can be achieved with a solid understanding of the species' biology, ecology and distribution (Bonnet et al., Reference Bonnet, Vourc'h, Raffetin, Falchi, Figoni, Fite, Hoch, Moutailler and Quillery2022). Knowledge of H. marginatum's current distribution is crucial for projecting areas at risk for CCHF expansion. The Mediterranean basin holds special importance as the main distributional area and reservoir for H. marginatum. The Mediterranean is also among the areas projected to be most negatively affected by climate change (Ulbrich et al., Reference Ulbrich, May, Li, Lionello, Pinto and Somot2006; Bardsley and Edwards-Jones, Reference Bardsley and Edwards-Jones2007; Newbold et al., Reference Newbold, Oppenheimer, Etard and Williams2020). Studies have indicated an increase in CCHF transmission in the Mediterranean basin (Maltezou and Papa, Reference Maltezou and Papa2010) and a tendency for the northward expansion of CCHF's range in this region (Estrada-Peña and Venzal, Reference Estrada-Peña and Venzal2007; Williams et al., Reference Williams, Cross, Crump, Drost and Thomas2015; Andersen and Davis, Reference Andersen and Davis2017; Fernández-Ruiz and Estrada-Peña, Reference Fernández-Ruiz and Estrada-Peña2021). Recent reports have shown that the vector of CCHF, H. marginatum Koch, 1844 might have expanded its geographic range and shifted northwards into previously unoccupied areas (Vial et al., Reference Vial, Stachurski, Leblond, Huber, Vourc'h, René-Martellet, Desjardins, Balança, Grosbois, Pradier, Gély, Appelgren and Estrada-Peña2016; Bah et al., Reference Bah, Grosbois, Stachurski, Muñoz, Duhayon, Rakotoarivony, Appelgren, Calloix, Noguera, Mouillaud, Andary, Lancelot, Huber, Garros, Leblond and Vial2022). Furthermore, the expansion of suitable areas has been documented in central Europe and Balkans in addition to Mediterranean countries (Fernández-Ruiz and Estrada-Peña, Reference Fernández-Ruiz and Estrada-Peña2021). Recent detections have highlighted that northern latitudes might be becoming more suitable for H. marginatum activity and survival as a result of climate change; however, most of these are records of individual tick detections and yet there is not sufficient evidence of establishment of northern populations (Duscher et al., Reference Duscher, Hodžić, Hufnagl, Wille-Piazzai, Schötta, Markowicz, Estrada-Peña, Stanek and Allerberger2018; Chitimia-Dobler et al., Reference Chitimia-Dobler, Schaper, Rieß, Bitterwolf, Frangoulidis, Bestehorn, Springer, Oehme, Drehmann, Lindau, Mackenstedt, Strube and Dobler2019; Grandi et al., Reference Grandi, Chitimia-Dobler, Choklikitumnuey, Strube, Springer, Albihn, Jaenson and Omazic2020; McGinley et al., Reference McGinley, Hansford, Cull, Gillingham, Carter, Chamberlain, Hernandez-Triana, Phipps and Medlock2021).
Thanks to advances in machine-learning algorithms in ecology, such as MaxEnt (Phillips and Schapire, Reference Phillips and Schapire2004; Elith et al., Reference Elith, Phillips, Hastie, Dudík, Chee and Yates2011; Crisci et al., Reference Crisci, Ghattas and Perera2012; Merow et al., Reference Merow, Smith and Silander2013) and the advent of high-resolution climatic datasets such as Chelsa (Karger et al., Reference Karger, Conrad, Böhner, Kawohl, Kreft, Soria-Auza, Zimmermann, Linder and Kessler2017) and WorldClim (Fick and Hijmans, Reference Fick and Hijmans2017), ecological niche modeling ENM has become a fundamental method for investigating possible current and future distributions of disease-transmitting vectors over the past 2 decades (Chalghaf et al., Reference Chalghaf, Chemkhi, Mayala, Harrabi, Benie, Michael and Ben Salah2018; Raghavan et al., Reference Raghavan, Heath, Lawrence, Ganta, Peterson and Pomroy2020; Aguilar-Domínguez et al., Reference Aguilar-Domínguez, Moo-Llanes, Sánchez-Montes, Becker, Feria-Arroyo, de León and Romero-Salas2021; Alkishe et al., Reference Alkishe, Raghavan and Peterson2021; Moo-Llanes et al., Reference Moo-Llanes, de Oca-Aguilar, Romero-Salas and Sánchez-Montes2021). Correlative ecological niche models primarily rely on combining georeferenced species records with predictive environmental variables (e.g. climatic, topographic, soil) to build a coefficient matrix representing the organism's multidimensional niche (Warren and Seifert, Reference Warren and Seifert2011; Peterson and Soberón, Reference Peterson and Soberón2012; Sillero and Barbosa, Reference Sillero and Barbosa2021). ENMs have been widely used to predict the potential future distributions of vectors under various global circulation scenarios of climate, in addition to their current potential niche (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; Wu et al., Reference Wu, Li, Liu and Bao2022).
The main objective of this study was to predict climatically suitable areas for H. marginatum under both present-day conditions and future climate scenarios. To achieve this aim, bioclimatic parameters for present and near-future scenarios (2011–2040 and 2041–2070) obtained from the Chelsa database and H. marginatum locations from the literature and previous field surveys by the authors were used to build ecological niche models with MaxEnt to show suitable areas for present and future projections. The findings of this study will help predict potential new suitable areas in the Mediterranean basin and Europe that H. marginatum may colonize in the future and enhance surveillance efforts in areas identified as high risk.
Materials and methods
Occurrence records
Local data from Turkey were derived from 2 sources: (1) from published studies by first author (O. H.) until 2021, which used both morphological and molecular methods to identify tick samples (Hekimoglu and Ozer, Reference Hekimoglu and Ozer2017; Hekimoglu et al., Reference Hekimoglu, Sahin, Ergan and Ozer2021) and additional field collected samples from 2018 (n = 13 locations) and 2021 (n = 7 locations) which have been identified using 16s rDNA. Morphological identification was conducted using taxonomic keys (Apanaskevich and Horak, Reference Apanaskevich and Horak2008; Estrada-Peña et al., Reference Estrada-Peña, Mihalca and Petney2018). Since H. marginatum populations show morphological variations (Apanaskevich and Horak, Reference Apanaskevich and Horak2008) and Hyalomma species have been one of the most misidentified taxa using morphological methods (Estrada-Pena and De La Fuente, Reference Estrada-Pena and De La Fuente2016), conducting a molecular step is extremely important. For this purpose, mitochondrial 16S rDNA sequences (Mangold et al., Reference Mangold, Bargues and Mas-Coma1998) were generated for samples of localities where molecular data were not generated before. In total, 66 geographical points from Turkey were included in the analyses.
To represent the distribution area of the species worldwide, we included geographic locations of H. marginatum from different parts of the world from a previous compilation by Estrada-Pena and De La Fuente (Reference Estrada-Pena and De La Fuente2016), composed of literature reviews between 1970 and 2014. Firstly, the geographic information of H. marginatum was recorded on a separate data sheet. Secondly, these raw records were cleaned and reduced by removing localities, (1) where tick collection had been conducted from birds since finding ticks on birds does not mean that ticks can establish populations in these areas (Estrada-Peña et al., Reference Estrada-Peña, Martínez Avilés and Muñoz Reoyo2011; Fernández-Ruiz and Estrada-Peña, Reference Fernández-Ruiz and Estrada-Peña2021); and (2) sampling information from Cyprus was removed as molecular data do not confirm that H. marginatum exists in Cyprus (Hekimoglu and Ozer, Reference Hekimoglu and Ozer2017). After these steps, a total of 565 geographic coordinates were obtained.
Environmental variables
For environmental predictors, Chelsa V2.1 dataset was used (Karger et al., Reference Karger, Conrad, Böhner, Kawohl, Kreft, Soria-Auza, Zimmermann, Linder and Kessler2017, Reference Karger, Lange, Hari, Reyer, Conrad, Zimmermann and Frieler2022) which is available at https://chelsa-climate.org, a relatively new high-resolution (30 arcsec, ~1 km) climate dataset that includes additional important microclimatic variables in addition to the counterparts (Chelsa-Bioclim) of popular WorldClim bioclim variables (Hijmans et al., Reference Hijmans, Cameron, Parra, Jones and Jarvis2005; Fick and Hijmans, Reference Fick and Hijmans2017). These include variables related to microclimate, water content of air and humidity which are important for H. marginatum (Estrada-Peña et al., Reference Estrada-Peña, Martínez Avilés and Muñoz Reoyo2011; Estrada-Peña, Reference Estrada-Peña2023). Furthermore, while the WorldClim dataset used for current distribution predictions uses extrapolations made between 1970 and 2000, the Chelsa dataset uses a dataset created for the period between 1980 and 2010. The additional important variables (Chelsa-Bioclim+) included near surface relative humidity (hurs), vapour pressure deficit (vpd), climate moisture index (cmi), growth degree days above 0°C, net primary productivity (npp) and surface downwelling shortwave radiation (rsds). Growth season temperature (gst) and growth season precipitation (gsp) were also included, which are indicative of H. marginatum, whose activity coincides greatly with the growth season (Estrada-Peña et al., Reference Estrada-Peña, Martínez Avilés and Muñoz Reoyo2011). Bioclim variables 8, 9, 18 and 19 from the Worldclim dataset (Fick and Hijmans, Reference Fick and Hijmans2017) are reported to have significant spatial artefact problems visible as anomalous discontinuities between neighbouring pixels and thus it is generally advised to remove them before carrying out analyses (Escobar et al., Reference Escobar, Lira-Noriega, Medina-Vogel and Peterson2014; 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). In this study, these variables were removed from the Chelsa dataset since it was observed that the same artefacts are also present in Chelsa-Bioclim data. All remaining variable rasters were clipped to the area of interest, which is between latitudes −20° and 60° and longitudes 20° and 60° WGS84. All geographic computations were carried out with QGIS (QGIS Geographic Information System, 2022), GDAL library in Python (Open Source Geospatial Foundation, 2022) and R version 4.2.2 (R Core Team, 2022).
Future projections
For future predictions, environmental variables from different global circulation models (GCM) of ssp370 and ssp585 scenarios for CMIP6 (Eyring et al., Reference Eyring, Bony, Meehl, Senior, Stevens, Stouffer and Taylor2016) were used. ‘Shared socioeconomic pathways’ (SSPs) complement the ‘representative concentration pathways’ (RCPs) by involving socio-economic factors such as population growth, urbanization and technological advances, ssp585 represents the worst-case scenario while ssp370 represents a middle way between the worst and more optimistic scenarios (Hausfather, Reference Hausfather2020). These bad-to-worst case scenarios were chosen for a more prudent risk prediction as anthropogenic climate risks are becoming much more difficult to manage and complicated all around the world (Kemp et al., Reference Kemp, Xu, Depledge, Ebi, Gibbins, Kohler, Rockström, Scheffer, Schellnhuber, Steffen and Lenton2022). All of the available scenarios for Chelsa dataset were included for a better estimation. These include 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) for the Chelsa dataset, and also available at https://chelsa-climate.org.
Ecological niche modelling
Before building the models, occurrence records that fell outside environmental raster pixels were removed. In order to reduce spatial autocorrelation, the locations for H. marginatum were thinned to a 5 km radius using the spThin package for R (Aiello-Lammens et al., Reference Aiello-Lammens, Boria, Radosavljevic, Vilela and Anderson2015). The resulting occurrence data consisted of 470 locations (Fig. 1) which were split into a training set (70%: n = 335), a preliminary test set (cross-validation set) for evaluating the candidate models (25%: n = 111) and another secondary independent test set for evaluation of final models (5%: n = 24). To simulate the accessible M space (Soberon and Peterson, Reference Soberon and Peterson2005; Barve et al., Reference Barve, Barve, Jiménez-Valverde, Lira-Noriega, Maher, Peterson, Soberón and Villalobos2011) for H. marginatum, a buffered minimum convex polygon with a buffer area of 100 km around the occurrence records was created and environmental variable rasters representing the M space were clipped to this buffer area before building models. Thinning and the creation of the buffer zone were carried out with ellipsenm package for R, available at https://github.com/marlonecobos/ellipsenm (Cobos et al., Reference Cobos, Osorio-Olvera, Soberón, Peterson, Barve and Barve2020). Three sets were built by setting correlation thresholds between variables using the vifcor function in the R package usdm (Naimi, Reference Naimi2017). The vifcor function selects the variables with lower variance inflation factors (vif) from correlation pairs to build the sets. The correlation thresholds were 0.9, 0.8 and 0.75, respectively. For the first model dataset sets 2 and 3 were identical, so only set 1 and set 2 were used to build the models. For the Chelsa dataset some of the Bioclim+ variables are not available for future projections (vpd, hurs, rsds and cmi). Thus, 3 different sets were built not including these variables to build a separate second model (model 2) for future projections.
All occurrence points of Hyalomma marginatum after cleaning and thinning.
ENMs were built with the maximum entropy algorithm (Phillips and Schapire, Reference Phillips and Schapire2004) using MaxEnt 3.4.4 (Steven et al., Reference Steven, Miroslav and Schapire2021) implemented in the Kuenm package (Cobos et al., Reference Cobos, Townsend Peterson, Barve and Osorio-Olvera2019) for R. Using the Kuenm package, calibration models were built with all combinations of MaxEnt features L (linear), Q (quadratic), P (product) and H (hinge). All models were repeated with regularization multiplier parameters 0.1, 0.2, 0.3, 0.4, 0.5, 0.6, 0.7, 0.8, 0.9, 1, 2, 3, 4 and 5. The performance of all these models were evaluated primarily with the partial receiver operating characteristics (pROC) test and an omission threshold rate of 5% (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), secondary evaluations were done with Akaike information criterion corrected for small sample size (AICc). The selected eventual calibration models were statistically significant, that have omission rates below 5% and with ΔAICc values below 2 (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).
After model calibrations, final models were created with the selected parameters in the resulting calibration models using all the same occurrences of training and testing data with 10 bootstrap replicates. Then, these final models were additionally evaluated with the independent location data (n = 24) that were not used for building and selecting these models.
Map binarization was done by using the average logistic thresholds of maximum training sensitivity plus specificity of 10 replicates (Liu et al., Reference Liu, White and Newell2013). Consensus maps of the niche models for the future projections were created by taking the average of GCMs. In order to show the effects of model projections under novel future conditions, we transferred models with allowed projections in novel climate conditions (extrapolation with clamping to the GCM scenarios) (Cobos et al., Reference Cobos, Townsend Peterson, Barve and Osorio-Olvera2019). To see the places where conditions are more extreme compared to the calibration area of the models, a mobility-oriented parity metric (MOP) analysis was carried out and extrapolation risk of transfer regions was calculated with nearest 10% reference (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). Because of memory limitations, resolution of calibration area and transfer rasters were changed from 30 s to 2.5 min resolution and MOP analysis were carried out with the mop function included in the kuenm package (Cobos et al., Reference Cobos, Townsend Peterson, Barve and Osorio-Olvera2019). Again, consensus maps for MOP were created taking the average of all GCMs for each ssp scenario and 4 MOP maps were created in total.
Results
Model parameters
Environmental parameters selected based on correlation thresholds and vif for variable sets used in model calibrations are shown in Table 1. In total 420 candidate models were created for the first model (model 1). All of 420 candidate models were statistically significant (pROC test P < 0.05); however, of these only 82 fulfilled the omission rate criteria of 5% (omission rate = 3.6%). Only one of these met the AICc criteria of ⩽2 (ΔAICc = 0.90). This final selected model is a linear quadratic model with a regularization multiplier parameter of 0.2, which uses variables included in set 1 with a mean area under the curve (AUC) ratio of 1.21. The final evaluation with the independent test set resulted in an omission rate of 8.3%. Although this value is larger than the 5% threshold set for evaluations with the cross-validation sets, the independent test set had a sample size of n = 24, and 8.3% omission rate means that the final model classified 22 of 24 test locations correctly which is a pretty good result for a machine learning algorithm. The percent contribution and permutation importance of 12 predictors are shown in Table 2. Out of these, the most important parameter was cmi with a contribution of 33.1%, followed by Bio14 (precipitation amount of the driest month) with a 24.2% contribution and rsds with an 8.4% contribution. Bio6 (minimum air temperature of the coldest month) was the fourth most important predictor (7% contribution).
Environmental predictors used for different sets in model 1 and model 2 with explanation of the predictors in the first column
Percent contribution of environmental predictors to model 1 and model 2
Out of 630 candidate models for the second Chelsa model (model 2) built for future projection, 625 were statistically significant (pROC test P < 0.05). Out of these, 201 met the 5% omission rate criteria (≈ 4.5%), and only 2 of them met the AICc criteria: a PH (product hinge) model with a regularization multiplier of 2 (ΔAICc ≈ 0), and a QPH (quadratic product hinge) model with a regularization multiplier of 2 (ΔAICc ≈ 1.27). Both of these final calibration models used the set 1 variables created for the second model. Mean AUC values were 1.23 and 1.20 for the product hinge and quadratic product hinge models, respectively. The final evaluation with the independent test set showed omission rates 8.3% and 12.5% for the PH and QPH models, respectively, indicating that the PH model is better for predictions. For the second model, there were 11 predictors. The biggest contributors were Bio13 (precipitation amount of wettest month) with 31.3% followed by Bio14 (19.4% contribution) and Gdd0 (growth degree days above 0 °C) with a 13% contribution. Like model 1, Bio6 was the fourth important predictor for model 2 with 8.2% contribution. Threshold values for creating binary maps were 0.409 and 0.483, respectively.
Current and future predictions
Suitable areas for present conditions are shown in Fig. 2. According to the first model, the current climatically suitable region for H. marginatum stretches out from Iberian Peninsula to Anatolia and the Caspian Sea. This region includes most of the Mediterranean basin. This pattern mostly coincides with the reported locations of H. marginatum. The predictions of 2 models are mostly compatible with each other. Compared to the first model, the prediction of the second model for the current distribution shows a wider suitable area in some regions, especially in North Africa and near Hungary. Also, the first model shows wider suitable regions in Trans-Caucasus and Spain compared to the second model. This might be due to the fact that the second model does not include parameters like cmi and rsds which are present in the first model. On the contrary, the second model includes another important parameter, gdd0, which was eliminated from the first model in the correlation and vif selection procedure.
Maps of predicted suitable areas from the ENM results. (A) Red areas show the suitable regions under current conditions according to model 1. (B) Green areas show the suitable regions under current conditions for model 2.
The predicted changes in suitable areas under GCM scenarios compared to present conditions are also shown in Fig. 3. Future projections under ssp370 scenarios for 2011–2040 show a significant widening of the climatically suitable area. The direction of the expansion is northwards in Europe, westwards in Anatolia and southwards in North African regions. Interestingly, the projections for the 2041–2070 period do not show such widening compared to the 2011–2040 period in eastern Europe and Anatolia; for example, the patches of new barely suitable regions in Baltic regions and Russia are lost in the 2041–2070 scenarios. In addition to that, there are significant declines in some regions in both ssp370 and ssp585 scenarios; these areas of suitability decline are especially distinct in central Anatolia and some areas in Balkans and central Europe. While both ssp370 and ssp585 show similar patterns in new areas for the 2041–2070 period, ssp585 scenarios present a wider new region of low suitability in France and eastern Germany. The results of MOP analysis done for between calibration area and future projections are shown in Fig. 4. High extrapolation risk areas are located relatively far from the current distribution and projected shifting areas of H. marginatum and the regions with high or strict extrapolation start from the north and east of Caspian Sea, Middle East (Arabian Peninsula, Israel and Egypt) and Africa (excluding northwest Africa).
Maps of predicted suitable areas for future average of 5 GCM scenarios with differing degrees of loss and gain compared to current conditions for model 2. (A) For ssp370 in the 2011–2040 period. (B) For ssp370 in the 2041–2070 period (C) for ssp585 in the 2011–2040 period (D) for ssp585 in the 2041–2070 period.
Extrapolation risk in 4 GCM projections of H. marginatum with MOP10%. Green to black scale shows increasing risk extrapolation where black areas are regions with strict extrapolation. (A) For ssp370 in the 2011–2040. (B) For ssp370 in the 2041–2070 period. (C) For ssp585 in the 2011–2040 period. (D) For ssp585 in the 2041–2070 period.
Discussion
The outputs from the first model show that the main climatic limiting factors for H. marginatum are cmi and precipitation, followed by surface radiation and minimum temperatures. Relative humidity or moisture is crucial for ticks, as they are highly susceptible to cuticular water loss in their off-host state (Knülle and Rudolph, Reference Knülle and Rudolph1982; Benoit and Denlinger, Reference Benoit and Denlinger2010; Estrada-Peña et al., Reference Estrada-Peña, Ayllón and De La Fuente2012; Requena-García et al., Reference Requena-García, Cabrero-Sañudo, Olmeda-García, González and Valcárcel2017; Leal et al., Reference Leal, Zamora, Fuentes, Thomas and Dearth2020). The importance of water content of air for H. marginatum has also been demonstrated in previous habitat suitability models (Estrada-Peña et al., Reference Estrada-Peña, Martínez Avilés and Muñoz Reoyo2011; Fernández-Ruiz and Estrada-Peña, Reference Fernández-Ruiz and Estrada-Peña2021). Another important parameter is incoming surface radiation, which is a relatively neglected parameter in ENMs for ticks despite the fact that it is a critical parameter for small arthropods (Davis et al., Reference Davis, Farrey and Altizer2005; Battisti et al., Reference Battisti, Marini, Pitacco and Larsson2013; Barrett and O'Donnell, Reference Barrett and O'Donnell2023). Solar radiation might affect many different parameters related to ticks. Firstly, it is the main driver of microclimate temperature, and higher levels of solar radiation might increase tick abundance and questing activity (Kiewra et al., Reference Kiewra, Kryza and Szymanowski2014; Del Fabbro et al., Reference Del Fabbro, Gollino, Zuliani and Nazzi2015). However, due to ultraviolet-B (UVB) radiation, very high levels of solar radiation might also have deleterious effects on small ectotherms like acari (Sakai et al., Reference Sakai, Sudo and Osakabe2012; Sudo and Osakabe, Reference Sudo and Osakabe2015).
The predicted area for the current conditions in model 1 is mostly compatible with previous habitat suitability models created for H. marginatum (Williams et al., Reference Williams, Cross, Crump, Drost and Thomas2015; Fernández-Ruiz and Estrada-Peña, Reference Fernández-Ruiz and Estrada-Peña2021). Additionally, the model output shows potential suitable areas in Trans-Caucasus and Iran, especially in northern Iran near the coast of the Caspian Sea, despite that the dataset of locations used for the models did not include any locations from this region. Although no exact coordinates are presented in the literature, CCHF cases and also H. marginatum ticks carrying the virus have been reported from this area (Shemshad et al., Reference Shemshad, Rafinejad, Kamali, Piazak, Sedaghat, Shemshad, Biglarian, Nourolahi, Beigi and Enayati2012; Sofizadeh et al., Reference Sofizadeh, Telmadarraiy, Rahnama, Gorganli-Davaji and Hosseini-Chegeni2014; Telmadarraiy et al., Reference Telmadarraiy, Chinikar, Vatandoost, Faghihi and Hosseini-Chegeni2015; Sedaghat et al., Reference Sedaghat, Sarani, Chinikar, Telmadarraiy, Moghaddam, Azam, Nowotny, Fooks and Shahhosseini2017).
Model 2 output shows a larger climatically suitable area compared to the more restrictive predictions of model 1. Important parameters like moisture index (cmi) and surface radiation (rsds) were left out for this model, as they are not available for future scenarios. Moisture index is replaced by precipitation of the wettest period (Bio13) as the most important parameter in model 2, which was in sixth place in the previous model. Although net precipitation amount is not as strong as water deficit or moisture index as a predictor, this shows the importance of the water balance for climatic suitability. In addition, surface radiation is replaced by gdd0 (growth degree days above 0 °C), which was eliminated from model 1 due to correlation threshold with other parameters. Fulfilling the required degree days is an important necessity for H. marginatum populations to establish in a new region (Estrada-Peña, Reference Estrada-Peña2023). It has been reported that a yearly accumulated temperature of 3000–4000 °C is necessary for this species (Estrada-Peña et al., Reference Estrada-Peña, Martínez Avilés and Muñoz Reoyo2011); a northern limit roughly coincides with 47°N.
With the increased temperature in recent decades, new reports of H. marginatum have also been increasing in Europe (Hornok and Horváth, Reference Hornok and Horváth2012; Duscher et al., Reference Duscher, Hodžić, Hufnagl, Wille-Piazzai, Schötta, Markowicz, Estrada-Peña, Stanek and Allerberger2018; Bah et al., Reference Bah, Grosbois, Stachurski, Muñoz, Duhayon, Rakotoarivony, Appelgren, Calloix, Noguera, Mouillaud, Andary, Lancelot, Huber, Garros, Leblond and Vial2022; Lesiczka et al., Reference Lesiczka, Daněk, Modrý, Hrazdilová, Votýpka and Zurek2022). Gillingham et al. (Reference Gillingham, Medlock, Macintyre and Phalkey2023), who modelled the probability of survival and establishment of H. marginatum populations depending on cumulative degree days for the United Kingdom, forecasted increased risk in southern UK in the future. Also, the prediction maps for the climatically suitable area for ssp370 and ssp585 scenarios point to a northward expansion in Europe in the future. Depending on the models trained with trends in climate and tick records between 1970 and 2018, Fernández-Ruiz and Estrada-Peña (Reference Fernández-Ruiz and Estrada-Peña2021), assumed that the suitable climatic area in Europe would expand while maintaining the original Mediterranean distribution. Additionally, Williams et al. (Reference Williams, Cross, Crump, Drost and Thomas2015) also predicted a northward expansion under AR5 climate scenarios in the future. Another parameter that should be considered is minimum temperatures, which occupied fourth place in both models 1 and 2. In addition to the above-mentioned effect of temperature (degree days) on development time, overwintering survival is another important factor for H. marginatum, as most adults overwinter in the field (Valcárcel et al., Reference Valcárcel, González, González, Sánchez, Tercero, Elhachimi, Carbonell and Sonia Olmeda2020). An increase in the minimum temperatures in the field might increase winter survival and contribute to the ongoing and possible future expansion in the northern limits of H. marginatum populations (Estrada-Peña et al., Reference Estrada-Peña, Ayllón and De La Fuente2012). Projections for future distributions are compatible with these predictions, showing a very similar pattern of new suitable regions in North Spain, France, the Balkans and western Anatolia. The projections for 2041–2070 point to a continuation of this expansion; however, the projection also points to significant declines in some regions, mostly in central Anatolia and the Balkans. This is most probably due to the decrease in precipitation, as climatic simulations predict that the impact of decreased precipitation will be much higher in the eastern Mediterranean, especially in Turkey (Hoerling et al., Reference Hoerling, Eischeid, Perlwitz, Quan, Zhang and Pegion2012; Turkes et al., Reference Turkes, Turp, An, Ozturk and Kurnaz2020). It was previously suggested that the northern limit for H. marginatum is thought to be mainly determined by temperature, while the southern limit is determined by precipitation and humidity (Gray et al., Reference Gray, Dautel, Estrada-Peña, Kahl and Lindgren2009), and expansion in the north and contraction in the south is an expected outcome in the future.
It also has to be considered that in some regions, the change in suitability might be more complex than this pattern. With projected ongoing climate change, an interchange of climatically suitable and unsuitable areas can be seen in models (Gillingham et al., Reference Gillingham, Medlock, Macintyre and Phalkey2023). A previous model using RCP8.5 scenarios carried out locally for Romania depicted an increase in suitable areas until 2050, followed by a decreasing suitability until 2070 (Domşa et al., Reference Domşa, Sándor and Mihalca2016), which was also detected by future predictions in the current study (Fig. 3). This rather complex pattern in suitability might be due to the 2 main factors (temperature and precipitation) determining suitability. For instance, present-day models show a highly suitable region for H. marginatum in the central and central-north parts of Turkey. Contrasting with the relatively more northward expansion in Europe, future projections indicate that after 2040, this suitable area might shift towards coastal areas where required precipitation would be available compared to the now more arid continental regions.
Assessing the potential distribution range of ticks is important to predict the risk of emergence and re-emergence of tick-borne diseases (Estrada-Peña, Reference Estrada-Peña2008; Zhao et al., Reference Zhao, Wang, Fan, Ji, Liu, Zhang, Li, Zhou, Li, Liang, Liu, Yang and Fang2021). In this sense, the interaction between ticks and climate is extremely important, as the survival and biological functions of ticks depend strongly on the external micro-abiotic environment (Gray et al., Reference Gray, Dautel, Estrada-Peña, Kahl and Lindgren2009). It should also be kept in mind that finding a suitable host is one of the most important drivers, which affects the distribution of tick species in an area. On the contrary, H. marginatum is a 2-host species and uses a wide range of vertebrates at different stages of its life cycle; this complex behaviour of this species is another reason for the difficulties in predicting the future distributions of these vectors (Apanaskevich and Horak, Reference Apanaskevich and Horak2008; Valcárcel et al., Reference Valcárcel, González, González, Sánchez, Tercero, Elhachimi, Carbonell and Sonia Olmeda2020; Bonnet et al., Reference Bonnet, Vourc'h, Raffetin, Falchi, Figoni, Fite, Hoch, Moutailler and Quillery2022). For instance, an area might be classified as suitable for a species by an ENM in future, but the absence or scarce distribution of suitable hosts would prevent the establishment of populations in the area. Conversely, unexpected novel hosts in new environments would provide the necessary link for dispersal in new environments (Bakkes et al., Reference Bakkes, Ropiquet, Chitimia-Dobler, Matloa, Apanaskevich, Horak, Mans and Matthee2021). Additionally, landscape cover, configuration and habitat fragmentation are other potential factors on the distribution of ticks and transmission of tick-borne diseases (Suzán et al., Reference Suzán, Esponda, Carrasco-Hernández, Aguirre and Ostfeld2012; Perez et al., Reference Perez, Bastian, Agoulon, Bouju, Durand, Faille, Lebert, Rantier, Plantard and Butet2016). A previous study showed that in Turkey, the number of CCHF cases was associated with high landscape fragmentation and connectivity (Estrada-Pena et al., Reference Estrada-Peña, Vatansever, Gargili and Ergönul2010). On the contrary, Bah et al. (Reference Bah, Grosbois, Stachurski, Muñoz, Duhayon, Rakotoarivony, Appelgren, Calloix, Noguera, Mouillaud, Andary, Lancelot, Huber, Garros, Leblond and Vial2022) showed that in southern France, sclerophyllous vegetated or sparsely vegetated open natural areas present favourable habitats for H. marginatum rather than more humid or urbanized areas, a similar pattern is also present in Anatolia where the highest distribution is in the central-north while H. marginatum is significantly absent both in humid and densely forested Black Sea coast and hot and dry southern regions. Ecological niche models are generally more deterministic models that build the representative fundamental niche using presently available abiotic parameters and locations, then project this assumed niche space to assumed future abiotic parameter data calculated by predictive global circulation simulations (Escobar, Reference Escobar2020; Sillero et al., Reference Sillero, Arenas-Castro, Enriquez-Urzelai, Vale, Sousa-Guedes, Martínez-Freiría, Real and Barbosa2021). Thus, these projections are for predicting the climatically suitable areas in the predicted future under assumed scenarios. However, these kinds of models are among the most useful tools presently available for estimating the potential risk of vector-borne diseases and deciding the necessary precautions (Medlock and Leach, Reference Medlock and Leach2015, Ortega-Guzmán et al., Reference Ortega-Guzmán, Rojas-Soto, Santiago-Alarcon, Huber-Sannwald and Chapa-Vargas2022). MOP analysis is one of the tools to interpret the potential risk of uncertainty for transfer of models to novel conditions (Owens et al., Reference Owens, Campbell, Dornak, Saupe, Barve, Soberón, Ingenloff, Lira-Noriega, Hensz, Myers and Peterson2013; Cobos et al., Reference Cobos, Townsend Peterson, Barve and Osorio-Olvera2019). MOP analysis showed that the current and predicted future distribution areas did not include areas with high risk of extrapolation. Uncertainty maps for H. marginatum are similar to the Europe and North Africa projections of a previous uncertainty analysis of worldwide projections of the distribution of another tick species (Rhipicephalus sanguineus) to RCP scenarios done by Alkishe et al. (Reference Alkishe, Cobos, Peterson and Samy2020). The main difference is that the present MOP results of projections to ssp scenarios show new and wider areas with high extrapolation in North Africa and Middle East compared to the previous paper by Alkishe et al. (Reference Alkishe, Cobos, Peterson and Samy2020).
With the ongoing human-induced changes in the biotic and abiotic environment we have been observing a significant increase in dispersal and abundance of H. marginatum and also an increase in the transmission of pathogens carried by H. marginatum. Future distribution projections indicate a significant increase in potential risk due to these factors, so it is important to build both statistical and process-based mechanistic models with new data. Also publishing coordinates is strongly suggested, as it is observed that there are many studies which report species presence without any coordinates that can help future coordination between researchers and provide more reliable models. As H. marginatum is the vector of CCHF, change in predicted distributional range might point out potential new risk areas or widening of current endemic areas. This study may be helpful to forecast new risk areas and therefore to expand awareness and to start well-adapted prevention strategies against CCHF in these areas. Similarly, our models show that the climatically suitable region for H. marginatum matches the current distributional area, which can be interpreted as the need to strengthen and to maintain control measures to CCHF in the future.
Data availability
Locations used in this study are available at https://data.mendeley.com/datasets/n8b3kps7nk/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 3 anonymous reviewers for reviewing the paper and for their helpful suggestions.
Author's contributions
O. H. and A. C. K. designed the study; O. H. prepared and cleaned the data; A. C. K. and C. E. carried out the analysis and interpreted the results; A. C. K. and C. E. prepared the tables and maps; A. C. K. prepared the original draft; O. H., A. C. K. and C. E. reviewed and edited the manuscript.
Financial support
This research received no specific grant from any funding agency, commercial or not-for-profit sectors.
Competing interest
None.
Open access
