Ecological niche and potential geographic distributions of Dermacentor marginatus and Dermacentor reticulatus (Acari: Ixodidae) under current and future climate conditions

Dermacentor marginatus is a vector disease of both humans and animals and transmits the causative agents of Q fever (Coxiella burnetii) and the spotted fever group (Rickettsia raoultii and R. slovaca), as well as of Crimean–Congo hemorrhagic fever. Dermacentor reticulatus can transmit various pathogens such as Francisella tularensis, Babesia spp., tick encephalitis virus, Coxiella burnetii, Omsk hemorrhagic fever virus, and Rickettsia spp. and can cause serious skin lesions. Herein, ecological niche modeling (ENM) is used to characterize the niches of these two ticks and describe their potential distributional patterns under both current and future climate conditions, as a means of highlighting geographic distributional shifts that may be of public health importance. We assessed distributional implications of five general circulation models (GCMs), under two shared socio-economic pathways (SSP245 and SSP585) for the period 2041–2060. Predictions for D. marginatus showed broad suitable areas across western, central, and southern Europe, with potential for expansion in northern and eastern Europe. Dermacentor reticulatus has suitable areas across western, central, and northern Europe. Under future scenarios, new expansions were observed in parts of northern and eastern Europe and highland areas in central Europe. Despite broad overlap between the niches of the ticks, D. marginatus has a broader niche, which allows it to show greater stability in the face of the changing climate conditions. Areas of potential geographic distributional expansion for these species should be monitored for actual distributional shifts, which may have implications for public health in those regions.

1 Introduction

Ticks are the second most dangerous disease vectors (after mosquitoes) for humans, but for animals, they are the first (Ahmed et al., 2007). They can cause various health problems for humans and considerable economic losses by transmitting different pathogens to livestock, including bacteria (Rickettsia spp.) (Parola and Raoult, 2001), viruses (tick-borne encephalitis virus) (Dumpis et al., 1999), and protozoa (Babesia divergens, B. microti, and Theileria parva) (Homer et al., 2000). More than a billion cattle (80% of the global population) are at risk of tick-borne diseases, and the global annual loss estimate is around USD 7 billion (Estrada-Peña and Salman, 2013).

The ixodid ticks Dermacentor marginatus (ornate sheep tick) and D. reticulatus (ornate cow tick) are two Eurasian tick species of veterinary and public health importance (Rubel et al., 2016). The two species have different geographic ranges except for some areas of overlap in central Europe (Buczek et al., 2015). The two species have similar life cycles that last 1–2 years; D. marginatus adults are active during spring and fall and can be active during winter in southern areas; D. reticulatus adults are active from late August through May, and only winter weather interrupts their questing activities (Dautel et al., 2006; Földvári et al., 2016).

Dermacentor marginatus is a vector of the causative agents of Q fever (Coxiella burnetii), spotted fevers (Rickettsia raoultii and R. slovaca) (Walter et al., 2016), and Crimean–Congo hemorrhagic fever (virus family Bunyaviridae) (Nuttall et al., 1994). Dermacentor marginatus uses a considerable variety of domestic animals (sheep, goats, horses, cattle, dogs) and wild animals (deer, wolf, hare, hedgehog, wild boar) as their main hosts (Masala et al., 2012). This tick species can be found in lowland, forest, alpine steppe, and semi-desert areas; its geographic distribution extends from Europe (Spain, Italy, Switzerland, France, Germany, Poland) to north Africa (Morocco, Tunisia) and western Asia (Iran, Kazakhstan) and includes much of central Asia (Walter et al., 2016; Rubel et al., 2015).

Dermacentor reticulatus transmits the causative agents of various diseases such as tularemia (Francisella tularensis), babesiosis (Babesia canis, B. caballi), tick-borne encephalitis virus, Q fever (Coxiella burnetii), Omsk hemorrhagic fever virus, and Rickettsia spp., and the bite of this species can cause extensive skin lesions in animals (Abdullah et al., 2016; Buczek et al., 2002; De Marco et al., 2017; Földvári et al., 2016; Jongejan et al., 2015; Wójcik-Fatla et al., 2011). This tick species uses more humid habitats in urban and suburban areas such as fallow land and parks (Rubel et al., 2016). Dermacentor reticulatus infests a broad range of mammal hosts (dog, cattle, deer) and bite humans only incidentally (Pfäffle et al., 2015). This tick species is well documented in central Europe and is distributed mostly in temperate and cooler regions (Rubel et al., 2016). Physiological experiments indicate that adults are cold-tolerant and can survive at temperatures of −10 ^∘C for up to 150 d (Medlock et al., 2017).

Europe and the Mediterranean region are expecting to face considerable change in rainfall and temperature; the latter is forecasted to increase by 1.5–2.5 ^∘C in coming decades (Gray et al., 2009; Eea, 2019). Those abiotic factors (temperature and precipitation) are crucial for ticks' survival on the ground (knowing that 99 % of tick life is off-host) (Sonenshine and Roe, 2013). Given the global warming impact on shaping the geographic distribution of D. marginatus and D. reticulatus, our aim here is to characterize their ecological niches and potential geographic distributions under climate change scenarios – understanding likely distributional changes in these species will be useful in guiding public health interventions as tick-borne disease risks evolve. We used ecological niche modeling tools to understand the species' niches in current environmental space and project them onto geography under two shared socio-economic pathways (SSP245 and SSP585), using five general circulation models (GCMs) in the period 2041–2060. This analysis has allowed us to explore these species' niches in environmental dimensions, to delineate their current potential geographic distributions, and to predict possible range expansions and reductions in the future.

2 Methods

2.1 Data preparation

We obtained occurrence data from various sources for D. marginatus and D. reticulatus: Global Biodiversity Information Facility (GBIF; http://www.gbif.org, last access: 15 January 2022; 190 and 514 (old data) and 2778 (new updated data)), VectorMap (http://vectormap.si.edu/, last access: 30 November 2020; 342 and 45), and the scientific literature (571 and 520) (sources summarized in Fig. S1 in the Supplement). We filtered the data to avoid duplicates, errors, localities with missing coordinates, and other georeferencing errors (Cobos et al., 2018). Based on the number of data points and environmental heterogeneity in the area and to avoid problems with spatial autocorrelation, we rarefied the data spatially based on a 50 km distance filter, considering the spatial resolution of the climate data (see below) of 10^′ or ∼17 km using the spThin R package (Aiello-Lammens et al., 2015) – this spatial resolution was chosen in light of the rather imprecise spatial resolution of the input occurrence data available to us. After these data cleaning steps, we had 236 occurrence points for D. marginatus and 198 for D. reticulatus. We split the cleaned occurrence data randomly into two sets of equal size for model calibration and evaluation steps involved in model calibration and used the entire set for producing final models.

To avoid negative impacts of inappropriate background choices for model outcomes (Anderson and Raza, 2010), we used explicit hypothesis of accessible areas for calibrating models (M; Barve et al., 2011). That is, we assumed a symmetric, long-term dispersal away from known populations and assumed that dispersal would likely be mediated by movements of host animals. In this regard, M areas were delimited based on 200 km buffer areas around the occurrence records for each species (Fig. 1).

Figure 1 Occurrence points and calibration areas for Dermacentor marginatus (red dots, red buffer), and D. reticulatus (blue dots, blue buffer).

We used bioclimatic variables for current and future (2041–2060) scenarios from WorldClim version 2.1, at 10^′ spatial resolution (Hijmans et al., 2005; available at http://www.worldclim.org, last access: 11 February 2022); variables 8, 9, 18, and 19 were removed because they are known to present abrupt discontinuities in areas not characterized by geographic breaks, which have been considered spatial artifacts (Escobar et al., 2014; Bede-Fazekas and Somodi, 2020). Variables were masked to fit the calibration areas (equivalent to the accessible area, M), and a principal component analysis (PCA) was used to reduce dimensionality and multicollinearity. PCAs were performed separately for temperature and precipitation variables to allow us to assess which factors contribute most to our models. We followed Cobos et al. (2019) in creating distinct sets of predictors (Table 1) and testing them in concert with other parameter settings to find the best models and sets of variables during model calibration (see below).

Table 1 Variables included in the five sets of predictors tested during model calibration. X indicates using a principal component (PC) variable.

To represent future climatic conditions, we used bioclimatic variables from five general circulation models (GCMs) in the context of two shared socio-economic pathways (SSP245 and SSP585). For future climate data layers, the GCMs used were (1) the high-resolution version of the Beijing Climate Center Climate System Model (BCC-CSM2-HR), (2) a Centre National de Recherches Météorologiques model (CNRM-CM6-1), (3) the Model for Interdisciplinary Research on Climate (MIROC6), (4) an Institut Pierre-Simon Laplace model (IPSL-CM6A-LR), and (5) the Canadian Earth System Model version 5 (CanESM5). GCM choice was based on availability of data for both SSP scenarios and on frequency of use in other such research applications.

2.2 Ecological niche modeling and model transfers

The combination of five sets of environmental variables, 15 feature classes (all combinations of linear, L; quadratic, Q; product, P; hinge, H), and 15 regularization multiplier values (0.1 to 1 at intervals of 0.3 and 2 to 6 at intervals of 1) resulted in 1125 candidate models per species. We evaluated candidate models based first on statistical significance (partial ROC (receiver operating characteristic), pROC, p≤0.05; Peterson et al., 2008) and then based on predictive performance (omission rates < 5 %; Anderson et al., 2003). Finally, we applied a criterion of minimum complexity, the Akaike Information Criterion corrected for small sample sizes (AICc; Warren and Seifert, 2011); differences between particular AICc values and minimum values are referred to as ΔAICc, and candidate models with ΔAICc < 2 (Warren and Seifert, 2011) were selected to produce final models.

Final models were created using the complete set of occurrences and the parameterizations selected during model calibration. We performed 10 bootstrap replicates and transferred the models globally to current and future scenarios (representations of results were restricted to an area of direct influence that includes Europe, north Africa, and western Asia). We calculated medians of final predictions for each calibration area in which final models were produced to summarize model results. We binarized models using a threshold of the allowable omission error rate (E) of 5 %, assuming this value as a percentage of data that may have included errors that misrepresented environments used by the species. We calculated differences in suitable areas between current and the two future scenario SSPs (SSP245 and SSP585). To represent changes in suitable areas, we used the agreement of changes (stable, gain, and loss) among the five GCMs per SSP scenario. Specifically, for each SSP scenario, we compared all projections to future conditions based on distinct GCMs against the current projection and quantified the agreement of gain and loss of suitable areas, as well as the stability of suitable and unsuitable conditions (see Campbell et al., 2015). All modeling was done in R 4.2.1 (R Core Team, 2020) using Maxent 3.4.1 (Phillips et al., 2017) with the kuenm package (Cobos et al., 2019).

We assessed strict extrapolation risk using the mobility-oriented parity metric (MOP, considering the nearest 5 % of the reference point cloud of environmental space in the calibration areas compared to environmental conditions in projection areas; Owens et al., 2013). We estimated variance contributions from distinct sources (replicates, parameter settings, GCMs, and SSPs) in our model projections (Peterson et al., 2018). Variability is calculated as the variance in all layers that correspond to distinct groups considering the sources of variation mentioned above. MOP and model variability were represented geographically following Owens et al. (2013) and Cobos et al. (2019), respectively. Note that D. reticulatus does not present all sources of variation because only one set of parameters was selected for final models.

2.3 Ecological niche characterizations in environmental space

We characterized the ecological niches of the two species using three distinct approaches: (1) a niche similarity test (Schoener's D index) (Broennimann et al., 2012) and (2) a niche overlap test based on ellipsoids (e.g., Banks et al., 2021; Nuñez-Penichet et al., 2021), both to compare the two species' niches, and (3) a visualization of ellipsoidal niches and changes in suitability in a three-dimensional environmental space. Niche similarity tests were performed with 1000 replicates, using the test proposed by Broennimann et al. (2012) implemented in the ecospat R package (Di Cola et al., 2017). This test describes how similar the niches are considering the environmental values of occurrence records and areas that are relevant for the analyses for both species (i.e., our calibration areas). Environmental variables used for these analyses were the first PCA of temperature variables and the first PCA for precipitation variables.

We measured the overlap between ellipsoidal envelopes for the two species in two environmental spaces (precipitation and temperature) separately. We used the first three PCAs of temperature and precipitation in two distinct analyses to obtain the environmental information of the occurrence points for each species and fit three-dimensional minimum volume ellipsoids (Van Aelst and Rousseeuw, 2009; Osorio-Olvera et al., 2020). To measure overlap using these ellipsoids, we used two different approaches. First, we measured full overlap, which considered the whole set of environments that overlap between the two ellipsoids (note that some environments may not be observed in the M areas). To this end, a cloud of points uniformly distributed in the multi-dimensional environmental space that covers the two ellipsoids is used to detect how those envelopes overlap (Qiao et al., 2016). Additionally, we measured background overlap, which considers only the subset of the full overlap corresponding to environments manifested in the M areas of the two species of interest (Nuñez-Penichet et al., 2021). The overlap is measured using the Jaccard index (Mammola, 2019), where the value of overlap is the ratio between the number of points within the intersection of two ellipsoids (E1 ∩ E2) and the total number of points contained by two ellipsoids (E1 ⋃ E2): J = E1 ∩ E2$/$E1 ⋃ E2. To measure full overlap between the ellipsoid envelopes, we used 10⁶ uniformly distributed points. To perform analyses of overlap using the environmental combinations accessible to the species, we used the values of the PCs corresponding to the M of each species. To test statistical significance of our observed values of overlap, we compared the observed results to a null distribution of 1000 overlap values derived from comparing ellipsoid envelopes created from points randomly sampled from each species' M (N is the total number of occurrences for each species). The null hypothesis in this test is that the observed degree of niche overlap reflects M-area environmental similarity, versus an alternative hypothesis of niche differentiation, so it is key to consider the environments available to each species within its M. If the observed niche overlap value falls inside the upper 95 % of the null distribution, the null hypothesis cannot be rejected. As the full overlap analysis uses points that are uniformly distributed, significance tests (which rely on available conditions) were performed only for analyses of background overlap. Analyses of niche overlap were performed in R using the package ellipsenm (Cobos et al., 2020).

To understand the dynamics of changes in suitable areas in future scenarios, we used representations in three dimensions for precipitation and temperature variables. These plots were constructed to show the species' ecological niches as in ellipsoid envelope models and the environmental values in future conditions. Future conditions were classified as never suitable, always suitable (stable), gain, and loss to represent changes in the same way they were shown in geographic projections (Campbell et al., 2015). Plots for each GCM and SSP combination were produced instead of searching for an agreement of changes as in geographic projections because detecting such agreements is complicated in environmental space. Niche overlap analyses and three-dimensional representations of niches and future conditions of suitability were performed in R using the packages ellipsenm and rgl (Adler et al., 2019).

Figure 2 Left panels: potential geographic distribution of Dermacentor marginatus based on current conditions (in blue and gray) and future conditions (blue denotes no longer suitable; red denotes newly suitable). Right panels: agreement in strict extrapolation areas among the five general circulation models. Model variability: median of variance coming from GCMs, replicates, and SSP in future projections.

3 Results

3.1 Model calibration results

From 1125 candidate models for each species, D. marginatus and D. reticulatus, 996 and 1033 were significantly better than random expectations, respectively (pROC test, p≤0.05). Of these models, 17 and 81 met the omission rate (OR) criteria (i.e., OR ≤ 0.05), and only 1 and 2 models were selected as the best models for D. marginatus and D. reticulatus, respectively, based on AICc. For D. marginatus, models performed better with the variables in Set 4 (temperature PC1 and PC2 and precipitation PC1, PC2, and PC3), whereas for D. reticulatus variables in Set 3 (temperature PC1, PC2, and PC3 and precipitation PC1) were selected (Tables 1 and 2).

Table 2 Parameter settings selected in the process of model calibration, to produce final models for Dermacentor marginatus and D. reticulatus. Reg. denotes regularization; AUC denotes area under the curve; ROC denotes receiver operating characteristic; AICc denotes Akaike information criterion corrected for small sample sizes.

3.2 Changes in potential geographic distribution

Current predictions for D. marginatus showed broad suitable areas across western, central, and southern Europe, including southern Great Britain, Germany, France, the Netherlands, Belgium, Italy, Greece, and Turkey. Some areas in northern Europe and north Africa were also predicted as suitable. However, these areas were small, including parts of southern Sweden, Morocco, Algeria, and Tunisia (Fig. 2).

The model transfer to future conditions showed stable suitable areas (i.e., suitable both now and in the future) across central Europe and in restricted areas of north Africa. Range reduction (loss) was observed in southern parts of the current suitable areas; gains were detected especially in northern Europe (Sweden, Denmark, Finland), eastern Europe (Russia), and northwestern Europe (Fig. 2). In general, more losses and fewer gains were detected in the SSP585 scenario compared to SSP245. MOP results showed high agreement in strict extrapolative areas among future scenarios, with most such areas concentrated in Asia (Fig. 2). Model variability results showed no variation coming from parameters, and low variation came from replicates and GCMs SSPs (Fig. 2), but there was high contribution to variation from SSP in different parts of Europe, Turkey, and northern Africa (Fig. 2).

Current predictions for D. reticulatus showed different suitable areas compared to those for D. marginatus: large parts of southern Europe were considered not suitable for this species. Suitable areas extended across western, central, eastern, and parts of southern Europe, including Great Britain, Germany, France, the Netherlands, Belgium, Italy, northern Spain, Turkey, and eastern Europe including Russia. Models also identified some areas of northern Europe, in southern Sweden, Finland, and Norway (Fig. 3).

Figure 3 Left panels: potential geographic distribution of Dermacentor reticulatus based on current (in blue and gray) and future (blue denotes no longer suitable; red denotes newly suitable) conditions. Right panels: agreement in strict extrapolation areas among the five general circulation models. Model variability: median of variance coming from GCMs, replicates, SSP, and parameters in future projections.

Under future conditions, predictions for the two SSP scenarios showed similar patterns of range stability, expansion, and contraction (Fig. 3). Stable suitable areas were distributed across western, central, and eastern Europe. Range reduction was anticipated across large areas of southern Russia, Ukraine, Moldova, Romania, Slovakia, Italy, France, Spain, and Portugal. Gains of suitable area were in northern Europe, mainly in Sweden, Denmark, Finland, Russia, and Great Britain, and in highland areas across central Europe (Fig. 3). Strict extrapolative areas were found in both southern and northern Europe, as well as in parts of Asia and north Africa, in both scenarios (Fig. 3). Model variability came mainly from GCMs and SSPs and was concentrated in central and northern Europe; variation deriving from GCMs and SSPs was low (Fig. 3).

Global projections showed higher potential of D. marginatus to invade new areas outside its native range, compared to D. reticulatus. Suitable areas for D. marginatus were mainly in North America, South America, sub-Saharan countries, and Asia; only a few were in South America and Australia (Fig. S1). For D. reticulatus, in contrast, few suitable areas were detected in North America and China (Fig. S2). However, areas such as sub-Saharan countries, South America, and northern Asia were detected as strict extrapolation areas by MOP analysis for the species D. marginatus, whereas strict extrapolation areas were detected in areas that are not predicted as suitable for D. reticulatus (Figs. S1 and S2).

3.3 Characteristics of ecological niches

The niche similarity test showed a value of Schoener's D of 0.268 and a p value of 0.389, so we could not reject the null hypothesis of niche similarity (Fig. S7). The volume of the ellipsoidal ecological niche for D. marginatus was larger than that for D. reticulatus (41.5 against 7.4 units using precipitation PCs, 101.3 versus 33.0 units using temperature PCs). Considering precipitation, the full overlap between two ellipsoidal niches was 0.22 and the overlap using available environments was 0.76 (Fig. 4). As regards temperature variables, full overlap resulted in 0.21 and the overlap using available conditions resulted in 0.39 (Fig. 4). The null hypothesis of niche overlap (using the background) was rejected for both precipitation and temperature dimensions (Fig. 4). The change in conditions seen in environmental space showed that temperature had a larger impact in terms of changes in suitability compared to precipitation, for both species (Figs. 5 and S3–S6). Changes in suitability in environmental space also showed a unidirectional pattern of changes.

Figure 4 (a, b) Species niche overlap based on ellipsoid envelope models (background overlap). Blue and red ellipsoids represent ecological niches of D. marginatus and D. reticulatus, respectively. Warmer colors indicate conditions that are closer to the centroid of D. marginatus. (c, d) Comparison of observed values of niche overlap against the null distribution of 1000 overlap values of random ellipsoids. Solid blue lines show observed values. Dashed blue lines show the 5 % confidence interval.

Figure 5 Changes in environmental conditions from the species niches' perspective, shown for one GCM (CanESM5). Blue dots indicate sites losing suitability. Red dots indicate sites gaining suitability. Gray dots indicate stable suitable sites.

4 Discussion

This study perforce focused on data points available to us, which were concentrated in Europe and parts of western Asia and north Africa (e.g., D. marginatus). Recently, updated GBIF data served to include more occurrence points for D. reticulatus, including parts of eastern Europe and Russia, thanks to iNaturalist becoming a new data contributor to GBIF. In this regard, this study is the first in terms of characterizing ecological niches and modeling the potential geographic distributions of D. marginatus and D. reticulatus under current and future climate conditions and provides a comparison between older data and a broad set of occurrence points for D. reticulatus. In this way, our results highlight the importance of data sharing and having a complete a set of data as possible with which to assess the potential geographic distribution of a species.

Characterizations in environmental and geographic spaces showed that D. marginatus has a broader niche than D. reticulatus. Although niche similarity could not be rejected using the test from Broennimann et al. (2012), the rejection of similarity using the overlap analysis in terms of precipitation and temperature dimensions and the broader niche of D. marginatus could indicate the higher stability of suitable areas for this species under future conditions (Figs. 2–4). We included detailed assessments of model uncertainty to detect areas with strict extrapolation and variation from different sources in model projections. Those two manifestations of model uncertainty are crucial to identifying areas that may lead to misinterpretations of the potential geographic distribution of those tick vectors (Alkishe et al., 2020). For instance, in future projections, areas of strict extrapolation are places with non-analogous conditions to those of the current period; therefore, whether suitability was detected or not in such areas depends entirely on the way the model extrapolates. Variability in our results plays yet another critical role, allowing us to detect areas where predictions can be more safely trusted (i.e., variability is low).

The study by Rubel et al. (2016) discussed the geographic distributions of Dermacentor species, assessed the ranges of the two species treated herein based on georeferenced sampling sites, and found that the limited distribution ranges for D. marginatus and D. reticulatus are 33–51^∘ N in latitude and 41–57^∘ N in latitude, respectively. A recent study by Okely et al. (2020) showed similar results for D. marginatus under current climate conditions, where models showed potential geographic distributional areas across Europe, north Africa, and Asia. For D. reticulatus, our initial results showed a more restricted potential geographic distribution across central and northern Europe before incorporating newer data from iNaturalist that represented more eastern areas more completely (Figs. 2 and 6d). However, these results changed to include broader areas in eastern Europe and Russia when we used new occurrences (Fig. 6b and c).

Figure 6 Comparison among three different models for Dermacentor reticulatus. (a) Old occurrence points in green dots and new occurrence points in blue. (b) Model transfer globally under current and future conditions using WorldClim version 2.1 and iNaturalist data from GBIF. (c) Model transfer globally under current and future conditions using WorldClim version 1.4 and iNaturalist data from GBIF. (d) Model transfer globally under current and future conditions using WorldClim version 1.4 and old data from GBIF (without iNaturalist data).

Our results added more detail about suitable areas under current conditions and allowed us to anticipate potential changes in suitable areas in the future by projecting models globally. The trend of gains seen for these tick species concurs with previous findings, obtained with similar methods, for other species in Europe (Dermacentor marginatus, Haemaphysalis punctata, Hyalomma marginatum, Rhipicephalus annulatus, R. bursa, and Ixodes ricinus) (Williams et al., 2015; Alkishe et al., 2017) and other parts of the world (James et al., 2015; Raghavan et al., 2019). Our global model showed fewer suitable areas in North America, compared to the species' native distributional area in Europe and Asia for D. marginatus, whereas D. reticulatus showed more restricted suitable areas (Figs. S1 and S2).

Major changes in suitable areas were associated with changes in temperature variables, especially for D. reticulatus, which coincides with the ecological niche differentiation of the two species in these dimensions. Our analyses and representations allowed us to detect that increases in temperatures will likely lead to significant losses and gains of suitable areas (Figs. 5 and S3–S6). These results support the findings of other studies, in which changes in the distribution of species have been mainly associated with changes in temperature (Buczek et al., 2015; Földvári et al., 2016; Walter et al., 2016). Although temperature conditions will change more or less evenly across the region, D. reticulatus is anticipated to see more losses of suitable areas (due to more areas becoming warmer) than D. marginatus in view of its narrower ecological niche that restricts it to higher latitudes and elevations. This idea has been previously suggested by Thuiller et al. (2005) and can be appreciated clearly when exploring the two species' niches and distributions in environmental space.

Our models help in reconstructing potential geographic-scale distributions of species based on environmental factors but do so without considering important factors of microenvironments for the ticks. Biological rhythms for ticks corresponding to daily and seasonal cycles (Belozerov, 1982) and the complexities related to distinct life stages, for instance, are among the most important factors that we do not consider. Temperature and relative humidity have strong influences on tick survival and development rates (Mehlhorn, 2012). Given direct impacts of climate change, particularly in the form of milder winters, new distributional patterns for some Ixodes ticks and tick-borne diseases have been observed, including new latitudinal and elevational records in the Northern Hemisphere in different temperate and boreal regions (Caminade et al., 2019).

Ticks depend on hosts, and physiological features of these hosts can affect their behavior and activity (e.g., feeding and attachment period) (Carroll, 2003; McCoy et al., 2013). In other words, tick population dynamics also depend on the impacts of environmental changes on their host, with potentially complex interactions (Perret et al., 2004; Wikel, 2018). Animal hosts also play an important role by facilitating tick dispersal into new regions (Randolph, 2004; Medlock et al., 2013). In lab experiments, competition between D. marginatus and D. reticulatus could also determine survival, where female D. marginatus show more intensive blood feeding than female D. reticulatus (Buczek et al., 2015). Although we did not include those biotic factors in our models, given the diffuse, one-to-many nature of tick–host interactions in these two species, these aspects may be highly relevant to understanding future scenarios of distribution and interactions of the two ticks studied here (Estrada-Peña et al., 2016). It is worth mentioning that our predictions do not represent estimated densities of species in areas but rather just an approximation of the suitable areas. We believe that these predictions will support public health organizations in anticipating geographic expansions and reductions and help them to take the appropriate steps to investigate and control the spread of those ticks in new areas.