The Crimean‐Congo haemorrhagic fever tick vector Hyalomma marginatum in the south of France: Modelling its distribution and determination of factors influencing its establishment in a newly invaded area

2.1. Study area and time of sampling

We chose to sample ticks in equestrian structures because they are most likely to be infested in areas where H. marginatum is established, as horses seem to show higher infestation rates by H. marginatum adult stages than any other domestic ungulates (Grech‐Angelini et al., 2016; Santos‐Silva et al., 2011). By equestrian structure, we refer to any facility that houses horses; this definition includes stables for horse owners, leisure structures offering occasional horse riding, equestrian centres that provide lessons and usually participate in competitions, horse‐breeding structures and traditional structures rearing horses used in taurine training.

We sampled structures across 14 French continental departments forming a link on the Mediterranean coast between the borders of Spain and Italy, which are two countries where H. marginatum is known to be present. A department corresponds to a NUTS 3 (Nomenclature of Territorial Units for Statistics) administrative division. Some of the selected departments have a Mediterranean environment, while others are characterised by a transitional zone where climate tends to become much cooler and humid and is assumed to be less suitable for tick establishment (Joly et al., 2010).

We first conducted tick collections in 2016, in the department of Bouche‐du‐Rhône (Desjardins et al., 2018); subsequent samplings dedicated specifically to detecting H. marginatum were then carried out in 2017, 2018, 2019 and 2021. The sample from 2021 was only used for model validation as it occurred after the analysis. To ensure data homogeneity, all tick samplings were conducted in spring (from April to July, mainly in May), which is the peak period of activity regarding the H. marginatum adult stage.

2.2. Tick sampling and identification

We aimed to sample ticks on at least eight horses per structure and if possible, from horses using different pastures to have a good representation of infestation in the structure. If the structures had fewer than eight horses, we sampled all available individuals and considered their infestation status as representative of the visited structure if they have access to the whole environment around the structure.

Tick collection consisted of a thorough full body search for any ticks on the animal. When feeding on a host, H. marginatum is usually attached around the genitals and anus, between the thighs, between the udders, and on the abdomen, in rare instances. While exploring its host, prior to attachment at the definitive engorgement site, the tick usually waits behind the pasterns or near the hooves (Apanaskevich, 2004). On horses, other tick species can be found attached on dorsal side of the host, or near the mane or tail. All ticks were removed manually and placed in plastic tubes. Each evening, ticks were sorted by sex and developmental stage, and the species was determined according to its morphological characteristics using relevant species descriptions and identification keys (Apanaskevich & Horak, 2008; Estrada‐Peña et al, 2004; Claudine, 2007). After identification, ticks were frozen in liquid nitrogen. Once returning from the field, ticks were stored at −80°C in the laboratory for later use.

2.3. Equestrian structures selection

We included structures that used as little acaricide treatment as possible and which provided horses with access to vegetated areas (Tirosh‐Levy et al., 2018) to favour horses’ exposure to the environment and therefore to ticks. A structure was included in our sample once we obtained owner consent. Furthermore, before our visit, we specifically asked owners not to apply any parasite treatment and, upon visiting, we verified the absence of such treatment (acaricides, but also anthelminthic ivermectin that may have an effect on ticks). Horses that had recently received treatment (less than 14 days) were excluded from the study. During the visit, GPS coordinates of each horse structure or pasture (when far from the main structure) were recorded. In addition, we filled out a form compiling information about the structure (owner, activity, nature and number of animals), the management of horses (last movements, last treatments, any unusual observations), and the scope of the visit (number of horses examined, number of horses parasitised by ticks, and specifically by H. marginatum). Since horses’ movement within their respective areas should directly impact the number of encounters with ticks and therefore infestation rate and probability of presence, this information was used to determine horse‐movement classes, categorised as followes: (1) Limited (level‐1): horses staying in their paddocks and rarely going on strolls for a day or less; (2) Medium (level‐2): horses trekking for a few days, but staying within the study zone; (3) Important (level‐3): Horses going out of the study zone for competitions, breeding or trade.

We reported in each structure the mean number of ticks found on examined horses. We assumed that the tick detection error on an infested horse was negligible since the tick attachment location on horses is well known for H. marginatum. In addition, most of the examinations were conducted jointly by two observers. We considered that H. marginatum was absent from a structure when no tick was found on a sufficient or representative number of horses, and after ensuring the absence of acaricide treatments; if treatments were conducted within the last 14 days, absence was validated if some specimens from other tick species were found on horses. We considered that the tick was established as a permanent population in a structure when at least two ticks were found on one or two horses, provided that they had not recently moved. For H. marginatum, a single tick found on a horse may be due to an engorged nymph having been introduced by a migratory bird (Capek et al., 2014; Černỳ & Balát, 1957; Jameson et al., 2012), then moulting locally into an adult and finally attaching to a horse, without indicating a local development cycle. This type of single detection has already been reported in Germany, the United Kingdom and Sweden (Chitimia‐Dobler et al., 2019; Grandi et al., 2020; Hansford et al., 2019; Lesiczka et al., 2022). Therefore, when a single H. marginatum was found in the whole structure, it was labelled as a “single presence” and was not considered to be sufficient reliable evidence of H. marginatum establishment (i.e., the existence of a permanent population in the area) to be used in the analysis. As female ticks engorge for 7–12 days and males can remain attached for at least a month, an observed infestation could also result from a recent horse movement from already‐infested areas. Collected information regarding recent horse movements helped to detect such doubtful detection cases, and only tick count data from horses that did not move from another department or country during the month prior to sampling were considered.

2.4. Environmental data collection and treatment

We used a combination of CORINE Land Cover (CORINE Land Cover 2018, Version 2020_20u1) and BD forêt (BD Forêt® version 2.0) to characterise the habitats. CORINE Land Cover (CLC) is a biophysical inventory of land use according to a 44‐item nomenclature. This inventory was produced at a 1/100,000 scale by visual interpretation of satellite images in 39 European countries and was updated in 2018. Pixels are 100‐m wide (area of 10,000 m²). The BD Forêt (BDF) version 2.0 is a reference database for forested areas and semi‐natural environments, which is compiled and made available by the Institut National de l'Information Géographique et Forestière (IGN) (French Institution of Geographic and Forest Information). It describes the forest and natural plant formations by a land cover approach in pixels wider than 70 m (areas larger than 5,000 m²). It is elaborated by visual interpretation of infrared images and is delivered according to a departmental segmentation, with updates ranging from 2014 to 2019. In addition, BDF contains information on a finer scale than CLC. For these reasons, we were much more confident in using BDF despite being limited to forested areas. Therefore, we first used CLC and then transferred BDF values to CLC in overlapping pixels that were classified as forests using the packages “raster” (Hijmans, 2020) and “rgdal” (Bivand et al., 2020). This resulted in a composite raster with a pixel resolution of 100^*100 m. From this composite raster and for each data point, we extracted land cover proportions within a radius of 5 km around the GPS coordinates to have information on the environment around the equestrian structures, which likely reflects the environment where horses can become infested by H. marginatum. Among the several kilometre radii that we tested within the 1 to 10 km range, we chose 5 km as a sensible trade‐off between bias and variance, which produced the best area under the receiver operating characteristic curve (AUC) using the leave‐one‐out method. AUC is a threshold‐independent measure of a model's ability to discriminate between sites where a species is present, versus sites where it is absent (Hanley & McNeil, 1982). To avoid model overfitting, the resulting 76 land‐cover classes from BDF and CLC fusion were reclassified into 7 generic classes: Humid areas such as inland marshes and rice fields, open natural areas such as sclerophyllous vegetation and sparsely vegetated areas, open forests, broad‐leaved forests, coniferous forests, urban or peri‐urban areas and agricultural areas. These classes were chosen as they reflected the main categories of habitats encountered within the study region, with a distinct capacity to retain humidity and allow for sun exposure, as well as a diverse availability of hosts. For each sampling point, the resulting seven variables represent the relative proportion of each habitat in a 5‐km radius. Regression with such compositional explanatory variables, which sum to unity, can be misleading due to their non‐Euclidean geometry (Aitchison 1982), which causes invalid estimation of the regression parameters of the model inference. Therefore, via the package “compositions” (Van den Boogaart et al., 2020), we used an isometric log‐ratio (ILR) transformation (Hron et al., 2012), which maps the original domain into a 6‐dimensional Euclidean space. The resulting variables, called balance coordinates, represent the coefficients of the original data with respect to an orthonormal Euclidean basis on the transformed space.

Regarding climate, we used the SAFRAN (Système d'Analyse Fournissant des Renseignements Adaptés à la Nivologie) analysis module of the French National Meteorological Service, from 2000 to 2018. SAFRAN provides both ground and altitude weather estimates from a numerical interpolation model fitted to weather measurements from the French network of meteorological stations. Weather estimates are presented in a regular 8‐km grid covering France. To characterise the climatic conditions of each sampling point, we collected monthly means of daily maximal and minimal temperatures, precipitation, potential evapotranspiration and relative humidity that are common features for describing climate, averaged across the 18 selected years to remove inter‐annual variation and reflect the climate of the last decades coinciding with the timeframe when the establishment of H. marginatum may have occurred. Principal component analysis (PCA) was performed on climate data to reduce dimensionality and avoid multicollinearity in our regression model. PCA allows for the description of a complex dataset using a small number of uncorrelated variables called principal components (PCs) that retain most of the variability and ecological meanings from the original data. For defining PCs, climatic data were extracted from each grid cell within the extent of the selected prediction area (south of France from the coastline to the latitude of 45°N), using the packages “ade4” (Dray & Dufour, 2007), “sf” (Pebesma, 2018) and “factoextra” (Kassambara & Mundt, 2020). Then, for each tick sampling point, PC coordinates were extracted from the cell where the point was located, using the package “rgeos” (Bivand & Rundel, 2020).

2.5. Statistical analysis

Out of the 169 sampling points used to develop the model (data from 2021 excluded), 113 had no ticks and 13 were single observations, corresponding to about two‐thirds of count observations being zeros. Thus, a zero‐inflated Poisson mixture (ZIP) model (Lambert, 1992) was used to estimate the probability of occurrence and abundance of H. marginatum. The ZIP model assumes that the data are generated by two underlying processes, which results in two distinct layers in the model: The “zero layer” and the “count layer”. The first process is modelled by a Bernoulli model, where the probability π _i of tick absence, resulting in an observed outcome of 0, may depend upon the covariates. The second process models the abundance of ticks, conditional to their presence, as a Poisson regression of mean μ _i. The response variable was the number of H. marginatum collected in each sampling point with the number of horses examined as an offset. The explanatory variables were the linear and quadratic functions of the three principal components of the climatic variables, the ILR balance coordinates of habitat and the level of horse mobility. Quadratic functions of principal components of climatic variables were included as it was hypothesised that the relationship between tick abundance and climate was not linear, with an expected optimum beyond which the population density will decline (Estrada‐Peña et al., 2011; Ogden et al., 2004; Randolph, 1997).

Model selection was based on Akaike's information criteria (AIC) and likelihood ratio tests (LRTs). The first model included all predictors in both the Binomial and Poisson components; we then conducted a backward stepwise elimination process to select the most predictive variables in the final model (Akaike, 1974; Woolf, 1957). At each step, we compared the model with and without each variable and removed the variable for which the LRT test was non‐significant and for which the AIC value was the lowest. The selection stopped when for all remaining variables, the LRT test was significant and the difference between AIC values was lower than 2. The resulting model had the lowest value of AIC and best goodness of fit.

Using the most parsimonious model to interpret the effect of each habitat under ILR transformation on the probability of presence and abundance of H. marginatum, we fitted multiple, equivalent but differently parametrised ZIP models such as the first balance coordinate representing one habitat type over the geometric means of all others as it is easily interpretable.

To assess the potential geographical expansion of H. marginatum in all of southern France, we then drew a raster map of the predicted probability of occurrence (presence/absence) based on the coefficients of the selected model, with a resolution of 100^*100 m. To do this, we sampled all the points within the prediction zone such as the resulting raster drawn from these points have a resolution of 100^*100 m and then we extracted climate and habitat data for each point. The predictive performance of the model was evaluated using three data sources: (i) The training dataset using the “leave one out” method, which consists of successively fitting the model without one observation and testing the prediction of the model at the removed observation; (ii) the collection of testimonies by people made aware of the H. marginatum existence thanks to several reports in the media (only testimonies confirmed with photos and with precise locations were included). As these data were presence‐only, and given that we made a buffer of 5 km around the sampling points, we sub‐divided the study area into 10‐km grids and took the centroid of these grids as a pseudo‐absence point so that each point was spaced 5 km apart and that there was no overlap between their buffer. (iii) The presence–absence data obtained from the last 2021 field survey, which were not included in the training dataset. We used AUC to assess the concordance between presence–absence records and the model predictions. For statistical analyses, we used the R statistical software R version 4.0.2 with the packages MASS (Venables & Ripley, 2002), lmtest (Zeiles & Hothom, 2002) and pscl. (Jackman et al., 2007)

3.1. H. marginatum distribution in the south of France

Out of 2588 ticks collected from horses among the whole study area, the most frequently identified species were Rhipicephalus bursa with 1147 individuals (44%), including 572 males, 557 females and 18 nymphs, followed by H. marginatum with 931 (36%) individuals, including 524 males and 407 females. The other tick species consisted of 244 Ixodes ricinus, including 18 males, 173 females and 53 nymphs (9%), 97 Dermacentor marginatus, including 47 males and 50 females (4%), 69 Dermacentor reticulatus, including 30 males and 39 females (3%), 63 Haemaphysalis punctata, including 24 males, 33 females and 6 nymphs (3%), and 37 Rhipicephalus of the sanguineus sensu lato (s.l.) group, including 13 males and 24 females (1%). Where H. marginatum was present, it was predominant (75%) with 931 individuals among the 1234 collected ticks in the sampling points.

A total of 188 equestrian structures, i.e., sampling points, were visited. The sampling points covered the Mediterranean eco‐climatic zone as well as surrounding climates but their distribution was uneven, with a few sampling sites in some departments at the border of the study area, such as Aveyron, Lozère and Ariège, while the other departments were sampled on at least 6 sampling points. An average of 10 horses per equestrian structure was examined (ranging from 2 to 35). Horses from 82 of these structures stayed in their paddock and rarely went for a stroll (level‐1 movement), animals from 59 structures were kept in pastures or moved regularly for a few days (level‐2) and horses of the 47 remaining structures participate in competition or breeding events (level‐3).

In total, 63 of the 188 visited structures were infested by H. marginatum, including 13 single presences, while the remaining 125 were considered negative. The mean parasitic load by H. marginatum (mean number of H. marginatum per examined horse) ranged from 0 to 11.43, depending on the structure. The tick was reported in 7 of the 14 sampled departments, but its distribution was not homogeneous over this zone. There were “clusters” of presence around Perpignan in Pyrénées‐Orientales, in Corbières hills in the south of Aude, in the hinterland of Montpellier in Hérault, in the north of Gard, the southern part of Ardèche and the extreme south of Drôme, as well as around the Maures mountainous massif in Var (Figure 1). Absences were mainly reported from Alpes‐Maritimes, Drôme and Bouches‐du‐Rhône to the east, and in the western part of the sampling zone. Although single presences were mainly located around presence clusters, a few were reported in areas presenting only absences such as in Bouches‐du‐Rhône, Aude and Tarn.

3.2. Correlations between climatic variables

The first three axes of PCA accounted for 87.7% of the variability of climatic data. They were therefore all included in the model.

The first principal component (PC1) accounted for 50.1% of the variability, with minimal and maximal monthly mean temperature scores strongly and negatively correlated with PC1 (scores > 0.86; Figure 2A). Similarly, monthly mean potential evapotranspiration scores had a moderate‐to‐strong negative correlation with PC1 (scores = 0.55–0.78; PC1, Figure 2A). Precipitation during summer (June to August) was correlated heavily and positively with PC1 (scores > 0.76; PC1, Figure 2A). Relative humidity during spring and summer (May to August) had a moderate and positive correlation with PC1 (scores = 0.43–0.49; PC1, Figure 2A). PC1 thus discriminates climatic conditions that are hot all year round (high temperature and potential evapotranspiration) and dry during summer (low precipitation and relative humidity) from climatic conditions that are cooler all year round and wetter during summer. PC1 can be interpreted as a temperature and summer precipitation index.

The second principal component (PC2) explained 28.7% of the variability, with monthly mean relative humidity negatively and strongly correlated with PC2 (scores > 0.77; PC2, Figure 2B), as well as a similar but moderate relationship with winter precipitation (January to March) (scores 0.61–0.72; PC2, Figure 2B) and a positive moderate relationship with potential evapotranspiration over the year (scores 0.38–0.65; PC2, Figure 2B). PC2 thus discriminates climatic conditions that are relatively dry from climatic conditions that are more wet. PC2 can thus be interpreted as a humidity index.

The third principal component (PC3) explained 8.9% of the variability and was mostly defined by precipitation during autumn (September to November), with scores strongly and positively correlated with PC3 (scores ≥ 0.73; PC3, Figure 2B). PC3 can thus be interpreted as an autumn precipitation index.

3.3. Habitats sampling effort

We measured H. marginatum abundance on structures in seven different habitats. Their distribution in our sample according to various criteria is presented in Table 1.

3.4. Effects of climate and habitat on H. marginatum distribution

The results of the ZIP model are presented in Table 2 and Figure 3. Although the “zero layer” of the model estimates the effect of explanatory variables on the probability of absence, results will be interpreted in terms of probability of presence by switching the sign of the estimates. The terms that remained in the most parsimonious model after the backward selection for the “zero layer” of the model were the three principal components, the squared PC2, and the ILR balance coordinates of habitats. Most of the terms were retained for the “count layer”, except for the squared PC2.

PC1 had a significant negative effect on the probability of presence of H. marginatum but showed a significant positive effect on the parasitic load (Table 2; Figures 3A and 3B). This means that H. marginatum is likely present in areas with climatic conditions that are relatively hot all year round and dry during the summer (Figures 1, 2 and 3 in supplementary data). Regarding the parasitic load, over the range of PC1 values on observed presence, the relationship is convex. With increasing PC1 values, a decrease in the predicted parasitic load is observed until a minimum is reached, followed by an exponential increase that rapidly exceeds the infestation rates observed on horses and thus does not reflect the reality. We noticed that the range of PC1 values on the observed presences are exclusively negative, thus corresponding to hot and dry conditions in summer. Considering the plots drawn from the model when covariates values are fixed to their mean, a 50% probability of presence is reached as soon as the mean annual temperature is 16.3°C and mean monthly summer precipitation is below 32 mm. For a maximal parasitic load of 12 ticks/horse as observed in the field, mean annual temperature remains above 17°C and mean monthly summer precipitation below 32 mm.

PC2 had a negative effect on the parasitic load and a positive effect on the probability of tick presence up to a point where it becomes negative (Table 2; Figure 3D). This means that H. marginatum is likely present under fairly dry conditions (Figure 4 and 5 in supplementary data) until a certain threshold after which conditions seem to become too dry. Conversely, increasing humidity is favourable to tick abundance. Nevertheless, for tick abundance, the range of PC2 values on the observed presences is exclusively positive, corresponding to dry conditions only. Maximal probability of presence is predicted for 71% mean annual relative humidity and mean monthly winter precipitation of 51 mm. An effective maximal parasitic load of 12 ticks/horse is obtained under a relative humidity of 73% and a mean monthly winter precipitation of 57 mm.

PC3 had a positive effect on the probability of tick presence. On tick parasitic load, in the observed presence range of PC3's values, the relationship is slightly convex. (Figures 3E and 3F). However, the relationship is nearly flat suggesting little involvement of PC3 on tick parasitic load. Thus, H. marginatum is likely to be present in areas presenting high precipitation in autumn (Figure 6 in supplementary data). A 50% probability of presence is reached as soon as mean monthly autumn precipitation is above 187 mm.

Open natural habitats had a significant positive effect on the probability of H. marginatum presence whereas urban or peri‐urban and humid habitats had significant negative effects, meaning that this tick is likely to be present in non‐flood‐prone open natural areas where human activities are scarce. Regarding the “count layer” of the model, most habitat parameters had significant effects on the tick parasitic load. In the same manner, as for probability of tick presence, open natural areas were favourable to large tick populations and humid habitats were unsuitable. Open and coniferous forests were considered unsuitable while higher proportions of agricultural and urban or peri‐urban habitats had a positive effect on the number of ticks found on horses. Finally, horse movement only had a significant effect on tick abundance. Horses are more likely to be exposed to a larger amount of ticks when strolling around or trekking for a few days within the sampling zone (level‐2) than horses staying in their pastures and paddocks (level‐1), but also more than horses that are used to traveling far away for competitions, breeding or trade (level‐3).

3.5. Prediction of H. marginatum distribution and model validation

Figure 4 shows the predicted probability of H. marginatum presence in southern France. Its presence was predicted only in coastal areas along the Mediterranean Sea with a small incursion towards the north along the Rhône valley. Within this zone, the highest probabilities were localised on the western margin from Pyrénées‐Orientales to Ardèche while areas to the east of the Rhône River, except Var, were considered less suitable for the establishment of H. marginatum. A certain heterogeneity was also observed within the suitable zone, with areas such as Camargue, Etang de Berre in Bouches‐du‐Rhône, the plains of Béziers, Perpignan in Hérault and Pyrénée‐Orientales, respectively, which were predicted as unsuitable by the model and sometimes confirmed by observed absences. Outside this zone, the western part of the country is predicted to be unsuitable once reaching the plains of Aude, Tarn and Haute‐Garonne, as well as the mountainous areas of Aveyron, Lozère, northern Gard and Ardèche. In the eastern part, the limit is less clear in Drôme, Vaucluse, Bouches‐du‐Rhône and Alpes‐de‐Haute‐Provence. While all observed absences of H. marginatum were well predicted by the model in the western part of the study area, some minor discrepancies between predictions and observations were observed, with suitability overestimation at the eastern margin of Camargue and the hinterland of Var, as well as suitability underestimation in the northern part of Ardèche. This resulted in moderate efficiency of the model using the “leave one out” validation method (AUC of 0.658 in ROC analysis). However, results improved using the testimonies and the occurrences of the 2021 survey for model validation (AUC of 0.761 and 0.833, respectively).