Grading Habitats for Ticks by Mapping a Suitability Index based on Remotely Sensed Data and Meta® population dataset in Aosta Valley, NW Italy

Abstract

Ticks represent a reservoir of zoonotic pathogens, and their numbers are increasing largely in wildlife. This work is aimed at producing maps of suitable habitats for ticks in Aosta Valley, Italy based on multitemporal EO data and veterinary datasets (tick species and distribution in wild hosts). EO data were processed in Google Earth Engine considering the following inputs: A) Growing Degree Ticks (GDT), B) NDVI from MOD09GA, C) NDVI entropy, D) distance from water bodies, E) topography, F) rainfalls from CHIRPS as monthly composites along the 2020, 2021 and 2022 years. Ticks were collected from hunted, injured, and found-dead wild animals ( Sus scrofa, Capreolus capreolus, Rupicapra rupicapra, Cervus elaphus); they were labeled at species level using taxonomic keys. Between September 2020 and December 2022, a total of 90 ticks were collected from 89 wild animals. Ixodes ricinus was the most prevalent tick species, followed by Dermacentor marginatus and Dermacentor spp. Molecular analyses demonstrated the presence of Anaplasma spp., B. burgdorferi sensu lato and Rickettsia spp. pathogens in infected ticks. To assess human population potential exposure to tick Meta^® population dataset was used. In conclusion this study shows the potentialities of Remote sensing improving the technological transfer to the veterinarian sector.

Introduction

Ticks are a known reservoir of pathogens that can be transmitted to humans and animals, especially wildlife (Estrada-Peña et al., 2008; Estrada-Peña & Salman, 2013). Many pathogens carried by ticks are zoonotic (Giangaspero et al., 2012; Wall & Shearer, 2001), meaning they can be transmitted from animals to humans. Some common zoonotic diseases transmitted are: Lyme disease (Pascucci et al., 2015; Stanek et al., 2012; Steere et al., 1983), Rocky Mountain spotted fever (RMSF) (Dantas-Torres, 2007), Ehrlichiosis (Ismail et al., 2010), Anaplasmosis (Kocan et al., 2000; Ristic, 1981), Babesiosis (Homer et al., 2000) and Tick-borne encephalitis (TBE) (Dumpis et al., 1999; Lindquist & Vapalahti, 2008). Specifically, Lyme disease is caused by the bacterium Borrelia burgdorferi, which is transmitted primarily by the black-legged tick (also known as the deer tick) in the northeastern and north-central United States (Marques et al., 2021), as well as parts of Europe and Asia (Steere et al., 1983). Symptoms can include fever, headache, fatigue, and a characteristic bullseye rash. Rickettsioses are caused by species of Rickettsia, a genus comprising Gram-negative bacteria, characterized by their strictly intracellular location and their association with arthropods (Raoult & Roux, 1997). Symptoms include fever, headache, rash, and in severe cases, organ failure. Then, Ehrlichiosis is caused by several different species of bacteria belonging to the genus Ehrlichia (Ismail et al., 2010). Different species of ticks can transmit these bacteria, including the lone star tick and the black-legged tick (Dagostin et al., 2023; Ebani et al., 2015; Viani et al., 2023). Symptoms can range from mild flu-like symptoms to severe illness involving multiple organ systems. Furthermore, Anaplasmosis is caused by the bacterium Anaplasma phagocytophilum, which is transmitted primarily by the black-legged tick (Arsevska et al., 2024; Kocan et al., 2000). Symptoms include fever, headache, muscle aches, and fatigue. Then, Babesiosis, which is caused by parasites of the genus Babesia, which are transmitted by the black-legged tick. Symptoms can vary widely, ranging from mild flu-like symptoms to severe illness with complications such as anemia and organ failure (Ristic & Kreier, 1984). Finally, TBE is caused by the tick-borne encephalitis virus, which is transmitted by several species of ticks in Europe and Asia. Infection can lead to inflammation of the brain (encephalitis), with symptoms including fever, headache, confusion, and in severe cases, paralysis and death (Estrada-Peña & Salman, 2013; Kalluri et al., 2007; Dantas-Torres et al., 2012; Millán et al., 2016; Schreiber et al., 2014). Pathologies affecting domesticated and companion animals have also received considerable attention, although there has been less focus on those affecting wildlife (Balestrieri et al., 2006; Botti et al., 2013; Viani et al., 2023). However, in certain regions of the world, there has been an increase in research in this area (Arsevska et al., 2024; Ippoliti et al., 2019; Short et al., 2017; Zhao et al., 2021), partly due to the establishment of control centers for game animals intended for trade and consumption, and partly as a reflection of environmental conditions, in line with a One Health approach (Floris et al., 2024; Orusa et al., 2020). In today's climate change and anthropogenic frameworks (Caminade, 2022; Caminade et al., 2019; Napp et al., 2018), the surveillance of diseases and the identification of conditions conducive to potential vectors, along with their distribution and spread, have become paramount for preventive measures and control efforts (Caminade et al., 2017). This underscores the importance of implementing management plans (Orusa et al., 2023; Orusa & Mondino, 2019). Under this premise, the advancements and technological transfers provided by geomatics to the veterinary sector enables the bolstering of monitoring initiatives and facilitate the development of predictive modeling and risk analysis (Porrello et al., 2019; Viani et al., 2023). While applications utilizing satellite data to aid in epidemiology and veterinary medicine are not yet widespread, there exist substantial opportunities for experimentation and application to create valuable tools in this field (Carella et al., 2022; Orusa et al., 2023). In particular, the development of Earth Observation (EO) missions and GIS techniques and more recently Artificial Intelligence (AI) have contributed to the enhancement of epidemiological data analysis capabilities creating new opportunities (Candeloro et al., 2020; Johnson et al., 2018; Orusa et al., 2024a; Orusa & Borgogno Mondino, 2021). In the specific case of vector-borne diseases, the utilization of geomatics and remote sensing enables the establishment of a concrete One Health approach (Anyamba et al., 2001; Lambin et al., 2010; Tran et al., 2013, Tran et al., 2016). Through the analysis of environmental components, a better understanding of the dynamics of interaction between pathogen/vector and host is achievable (Anyamba et al., 2009, Anyamba et al., 2019; Caldwell et al., 2021; Chretien et al., 2015; Ebi et al., 2020), elucidating how variations in space and time can affect the relationship within the disease triangle (Chamaillé et al., 2010; Orusa et al., 2020). In this regard, suitability maps based on locally or district-level geo-referenced ground data facilitate the creation of management plans, preventative actions, risk analysis, or modeling, thereby promoting day to day public health activities. The primary challenge lies not in accessing these data, which is often freely available in the form of open data, but in comprehending and appropriately processing it (Borgogno-Mondino et al., 2022; Farbo et al., 2022; Meloni et al., 2022; De Petris et al., 2021; Sarvia et al., 2021). Conversely, the availability of georeferenced health data collected, for instance, through appropriate apps operated by experts such as veterinarians and foresters, enables the production of moment-specific maps and/or the modeling of future trends (Caminade, 2022; Latini et al., 2020; Latini et al., 2021; Thépaut et al., 2018). Such analyses may involve both straightforward and more intricate models, leveraging neural networks and AI, as observed in other sectors (Conte et al., 2015; Espunyes et al., 2021; Orusa et al., 2024a). Since, the full potential of Geomatics and Remote Sensing capabilities in the ordinary workflows of veterinarians remains largely untapped, in this work an empirical approach is proposed aimed at recognizing and characterizing habitats of ticks through a suitability index based on EO data. The work results from an interdisciplinary collaboration among veterinarians, geomatics, foresters, and biologists and is primarily intended for understanding tick-wildlife host relationships starting from a characterization of tick habitats. The basic idea is the exploitation of multi-temporal remote sensing data and Meta® population (Tiecke et al., 2017; Verhulst et al., 2021) datasets in the framework of a One Health approach. Specifically, the following goals can be listed: i) identification of tick species and related pathogens; ii) detection of tick distribution in wild hosts; iii) mapping and characterization of a suitability index grading habitats for ticks; (iv) intersection of the suitability maps with human population to derive information useful to support risk assessments.

Materials and Methods

Study Area

This study, although scalable worldwide, was developed in the Aosta Valley Region (North Western Italy). This territory hosts the highest peaks in the entire Alpine chain and Europe as a whole. Bordered by France to the west and Switzerland to the north, despite being the smallest region in Italy, its geomorphology renders it one of the most complex from a landscape perspective, situated in the Western Alps (Orusa et al., 2022a; Orusa et al., 2022b) (Figure 1).

Figure. 1. Study area, Aosta Valley Autonomous Region NW Italy. RGB Sentinel-2 composite imagery 2022.

Veterinary workflow and protocols

Ticks were collected from wildlife that were hunted, injured or found dead belonging to the following species: wild boar (Sus scrofa L.), roe deer (Capreolus capreolus L.), alpine chamois (Rupicapra rupicapra L.), red deer (Cervus elaphus L.) within the study area. Analyzed ticks were sourced from animals evenly dispersed across the whole Aosta Valley territory. The number of host animals was consistent, albeit not equally distributed at a spatial level due to variation across different hunting districts. The available dataset covers approximately three years: from September 2020 to December 2022. In this period, 90 ticks were collected and analyzed from 89 wild animals parasitized by the aforementioned species. Characterization and identification of tick species was conducted starting from taxonomic keys according to the approach proposed by (Estrada-Peña, et al. 2018). The presence of tick-borne pathogens was investigated using previously published PCR protocols. In particular, tests were performed for Anaplasma spp., Borrelia burgdorferi s.l., Coxiella burnetii, Rickettsia spp. pathogens. Specifically, these end-point PCR protocols were optimized in a final volume of 25 μl, using Platinum™ Taq DNA Polymerase (Invitrogen, Thermo Fisher Scientific Inc.) and 5 μl of each DNA extract protocols employed to detect the pathogens. The references of PCR protocols used were: C. burnetii (WOAH Manual for Terrestrial Animals 2018, Cap. 3.1.16 par. B.1.3 - Target region: IS1111; WOAH, 2024; Jebara, 2007), Anaplasma spp. (Stuen, Nevland & Moum 2003 – Target gene: 16S rRNA), B. burgdorferi s.l. (Skotarczak, Wodecka & Cichocka 2002– Target gene: FLA) and Rickettsia spp. (Regnery, Spruill & Plikaytis 1991 – Target gene: citrate synthase). Finally, for TBE virus a Real Time PCR was carried out using the previously published method (Schwaiger & Cassinotti, 2003) in a total volume of 20 μl using QuantiTect Multiplex PCR (QIAgen, Germany) and 4 μl of RNA extract.

All positive samples at screening PCRs, were sent to the National Reference Centre for Anaplasma, Babesia, Rickettsia and Theileria (C.R.A.Ba.R.T, Experimental Zooprophylactic Institute of Sicily https://www.izssicilia.it/chi-siamo/centri-di-referenza/c-r-a-ba-r-t last accessed February 6, 2024) for sequencing.

EO data and geospatial workflow

Maps of the suitability index regarding habitats for ticks were generated using various geospatial layers, many of which derived from satellite remote sensing data.

Datasets

The Land Surface Temperature (LST) and the Normalized Difference Vegetation Index (NDVI) Daily collection from MODIS/061/MOD11A1 and MODIS/061/MOD09GA were obtained and processed by Google Earth Engine (GEE) (Gorelick et al., 2017). NDVI data has a native spatial resolution of 500 m: differently, LST data has a resolution of 1 km. Collections were obtained to cover the 3-year period when tick data were available. Data were pre-processed to remove (mask out) “bad pixels” due to clouds, shadows and defective pixels layers provided as Bitmasks (more information available at https://lpdaac.usgs.gov/documents/118/MOD11_User_Guide_V6.pdf and https://lpdaac.usgs.gov/documents/306/MOD09_User_Guide_V6.pdf , last accessed on February 6, 2024). Rainfall data were similarly obtained from the UCSB-CHG/CHIRPS/DAILY through GEE. They are provided as raster layers with a native spatial resolution of 5 km, over-resampled up to 1 km (Funk et al., 2015). Specifically, the Climate Hazards Group InfraRed Precipitation with Station data (CHIRPS) is a 30+ year quasi-global rainfall dataset. CHIRPS incorporate 0.05° resolution satellite imagery with in-situ station data to create gridded rainfall time series for trend analysis and seasonal drought monitoring with daily data. Additionally, a river/water bodies layer was obtained from the Aosta Valley Land cover map (Orusa et al., 2022a; Orusa et al., 2022b) available through the Aosta Valley SCT geoportal (https://geoportale.regione.vda.it/, last accessed on February 6, 2024). This map layer contains water channels and rivers. The distance-to-waterbodies was then calculated through the proximity grid function available in SAGA GIS tool version 8.5.0 (Conrad et al., 2015). Additional terrain features such as slope, aspect, TPI (Topographic Position Index), TRI (Terrain Ruggedness Index), and median yearly light hours duration, were computed using the Aosta Valley DTM (Digital Terrain Model) from SCT with a ground sample distance of 2 m, also within SAGA GIS tool version 8.5.0. The DSM (Digital Surface Model) was used to compute the Sky View factor necessary to compute and map the median yearly light hours. Finally, for population exposure and risk assessment, the HDX Meta® Population dataset (Orusa et al., 2023; Tiecke et al., 2017) was used, having a 30 m grid step. It was obtained for the entire period of analysis from https://dataforgood.facebook.com/ (last accessed on February 6, 2024).

Data processing

LST and NDVI dates were stacked at single year level. Native stacks were then processed at pixel level, removing bad pixels according to the above-mentioned masks. Local temporal profiles were then daily regularized (2nd order polynomial) and filtered by Savitzky-Golay filter (window size = 6, polynomial order=2) (Chen et al., 2004; Press & Teukolsky, 1990; Schafer, 2011).

Since ticks have a suitable active life thermal range well described in scientific literature (Estrada-Peña & Salman, 2013), starting from LST a new measure was defined for pixels that satisfy the condition reported in equation 1. Specifically, during the year, only the images that satisfied equation 1 were considered. Then, only the images in which the condition of equation 1 was true were averaged. In other words, the average LST for each year is calculated considering only the images that satisfy the temperature range for ticks within the stack profile. Therefore, only temperature suitable for ticks (LST_ticks) were considered in the computation of the mean composite (hereinafter called GDT - Growing Degree Ticks):

Then, regarding NDVI profile, only NDVI stacked images that temporally match and overlap the LST condition, were considered. In particular, a threshold of NDVI>0.3 was applied within the profile as reported in equation 2. It is worth to note that NDVI>0.3 represents areas where vegetation is present and therefore in which ticks may be present.

Therefore, NDVI<0.3 were masked out and only images that satisfy NDVI>0.3 (the vegetated ones) and matched with the LST condition (GDT) were considered in the computation of the mean NDVI composite (thNDVI). As previously reported, LST and NDVI data were processed and downloaded from GEE through the Javascript front end. The LST condition and NDVI thresholding were performed within R Studio 2023.12.1. For precipitation data from CHIRPS, a Javascript code was written within GEE to extract cumulative rainfall for each year of analysis. Finally, the entropy (which refers to a measure of uncertainty or randomness within a system and quantifies the disorder or unpredictability present in the spatial distribution of pixel values, indicating the level of complexity within an image) has been computed in order to assess patches in the landscape and in particular vegetation that may affect tick life as reported by (Burrows et al., 2022; Da Re et al., 2019) as well as hosts (Carella et al., 2022). This application calculates Haralick advanced (Haralick et al., 1973; Zhuang et al., 1987), and higher-order texture features for each pixel in the selected channel of the input image, specifically the thNDVI channel. thNDVI entropy (HthNDVI) was identified as a suitable image texture parameter capable of discerning changes in intensity within both natural (such as forests, grasslands, and agriculture) and anthropogenic systems. As demonstrated by (De Marinis et al., 2021), entropy serves as a valuable tool for delineating agricultural production systems as well as health mapping (Carella et al., 2022). In fact, abrupt changes (exceeding 0.05 in HthNDVI) have an impact on suitable tick landscape, and at the same time influencing wildlife movements and the potential for their interactions. Lower values of HthNDVI indicate orderliness in the landscape, while higher values denote a higher degree of disorder and entropy within the system, a perfect condition for tick and wildlife as also reported in other veterinary diseases (Orusa et al., 2020). Starting from thNDVI entropy was computed as reported in equation 3 using the open source Orfeo toolbox vers 8.0.0 (Grizonnet et al., 2017; Inglada & Christophe, 2009; McInerney & Kempeneers, 2015; Gascoin et al., 2019):

Modeling and mapping techniques

As suitability models are heavily influenced by input variables, their selection plays a pivotal role in defining an accurate model. Thus, the choice of certain variables was informed by the evidence presented in the extensive scientific literature and available on the subject matter further supported by empirical evidence (Anyamba et al., 2009; Ebi et al., 2020). While conducting a thorough analysis of variables using principal components and rigorous statistical techniques is preferable before constructing a suitability model, in this study, we relied on existing literature. We introduced ground temperature based on the tick's vital range conditions as a novel and less explored variable. Additionally, factors such as vegetation threshold and entropy were considered, as previous studies have shown their impact on host movement.

Maps of potentially suitable habitats for ticks were generated using the above-mentioned geospatial data. Input datasets are reported in Table I as follow:

Table. I. Input data suggested to derive maps of potentially suitable habitats for ticks per each investigated year upon Ebi et al 2020.

Since variables reported in Table I have different units, in order to generate the yearly map of the suitable habitats for ticks the, TfMSSmall function (eq. 4, Khan & Mohiuddin, 2018) was used to transform and normalize input maps:

where:

σ: standard deviation;

μ: mean;

ϑ: a σ multiplier from 1-10 depending on the weight;

τ: a μ multiplier from 1-10 depending on the weight;

S_n(x,y): input variable.

A Multiple Regression Analysis (MRA) was followed to define Weights (hereinafter call ω). Due to the lack of precise georeferencing for all ground tick survey data points, those localized at the municipality or sub-municipality level were converted into point data by assigning centroids using QGIS version 3.28 (QGIS Development Team, 2018). Consequently, a spatial ground tick dataset with point-level georeferencing, hereinafter referred to as GTP, was generated. We acknowledge that this method may not be the optimal approach for assigning a location due to a higher degree of uncertainty. However, the non-georeferenced points were less than 7% and we have determined the centroid by exclusively considering vegetated areas in the Aosta Valley land cover and relying on information provided by foresters and hunters, which we have assumed to be accurate. MRA was performed involving all the variables reported in Table II and GTP. Equation 5 below reports the computation of ω.

ω_i: input weight;

R²_(x,y): determination coefficient.

Weights obtained for each input variable (Table I) after performing MRA were reported in Table II.

Table. II. Weights obtained based on coefficient of determination by performing equation 5 after MRA.

Finally, the Suitability Index Maps for ticks (SIM) were generated using the function outlined in Equation 6, implemented in ArcGIS Pro version 2.8.0 (Corbin, 2018; Law & Collins, 2019), within the Suitability Modeler toolbox for each year investigated as follows:

where:

S_(x,y): input variable, with i from 1 to n, depending on the number of variables.;

ω: weight of each variable

Subsequently, SIM maps were subjected in SAGA GIS vers.8.5.0 to an ISODATA unsupervised classification clustering (Memarsadeghi et al., 2007) in order to group the suitability of ticks into a maximum permissible number of 12 classes and a minimum of 5 so as to carry out the risk analysis and better represent it, taking into consideration the distribution of the population over the entire territory of the autonomous region of Valle d'Aosta. The algorithm ISODATA stands for Iterative Self-Organizing Data Analysis Techniques. This algorithm which allows the number of clusters to be automatically adjusted during the iteration process by merging similar clusters and splitting clusters with large standard deviations. In particular, the tool is based on Christos Iosifidis' Isodata implementation (http://users.ntua.gr/chiossif/Free_As_Freedom_Software/isodata.c, last accessed on February 16 2024).

Specifically an intersection with HDX Meta® Population dataset was performed to assess the risk exposure per each municipality within the study area through zonal statistics. Although the ISODATA algorithm was created to discriminate the optimal number of classes and since in the initial conditions it was assumed a priori and forcibly to limit the maximum number of classes to 12, we wanted to test the statistical separability of the classes so as to verify the goodness of the initial a priori assumption. Therefore, the separability of each class obtained from ISODATA was checked by computing the Jeffries–Matusita (J-M) distances (Sen et al., 2019) to verify the adequate divisibility of the optimal summary classes already obtained from the unsupervised classification (Wang et al., 2018) to define the exposure to risk for the population. The complete workflow has been summarized in Figure 2. Finally, to map the risk of tick exposure across the municipalities of the Aosta Valley, an ISODATA clustering analysis was conducted, as depicted, dividing the regions into 10 classes ranging from minimum (1) to maximum (10). It is worth noting that the ranges from 1 to 10 express the population vulnerability. Notably, the Meta® dataset was employed for this purpose.

Figure. 2. Flowchart reporting the main steps of the proposed methodology for obtaining proper maps of suitable habitats for ticks.

Results

During the three years of analysis (2020, 2021, and 2022), a total of 90 ticks were detected from 89 wild animals. The taxonomic analysis performed on the sampled ticks revealed their respective taxa. According to detected species, the correspondent percentage frequency was computed (Table III).

Table. III. Tick species detection and percentage frequency over the ground survey collection.

Ixodes ricinus (93.4%) proved to be the most prevalent tick species, followed by Dermacentor marginatus (5.5%) and Dermacentor spp. (1.1%). It is worth noting that the tick bodies could sustain extraction without sustaining damage, enabling a proper recognition of species according to a taxonomic classification. Over the entire GTP within the study area, Ixodes ricinus was the most frequent tick species found on roe deer (Capreolus capreolus L.), alpine chamois (Rupicapra rupicapra L.) and red deer (Cervus elaphus L.), followed by Dermacentor spp., found only on roe deer. Dermacentor marginatus was the dominant parasitic species on wild boar (Sus scrofa L.). Table IV below reports the absolute frequency of each tick species by host over the whole period investigated. Furthermore, pie chart plots have been created from these data (Figure SM 1).

Table. IV. Absolute frequency of each tick species depending on the host.

In general, it can be said that molecular analyses demonstrated the presence of at least one of the analyzed pathogens in all tick samples. In order to better understand the distribution of pathogens on each tick and host species, pie charts in Figures SM 1 and SM 2 were created starting from the positivity count data reported in Tables IV and V.

Table. V. Positivity to pathogens tested carried by tick species.

The PCRs test for ticks screening resulted in the different pathogen positivity as reported in table SM I. In the table it is specified the result of molecular identification of each single pathogen carried out by the National Reference Centre for Anaplasma, Babesia, Rickettsia and Theileria.

Adopting the methodology outlined in the preceding section and visually depicted in schematic format in Figure 2, SIM were generated for the three years of analysis, as reported in Figure 3. It is worth noting that grids SIM were created with 1000 m Ground Sampling Distance (GSD) which is the lowest native spatial resolution of the suitability inputs. In fact, the Cumulative Rainfall (CR) was resampled at 1000 m. The SIM maps for the year 2020, 2021 and 2022 are reported below:

Figure. 3. Habitat Suitability Map for Ticks in 2022 (A), 2021 (B), 2020 (C). EPSG 23032 (UTM WGS84 32 N).

While the tick-prone areas predominantly encompass vegetated terrain rather than urban landscapes, urban areas were not masked within the clustering. This decision was made to account for the population's potential exposure when transitioning from inhabited to neighboring wooded and vegetated regions. In the Table VI, it has been reported the risk and standard deviation values obtained for each cluster class whose statistical separability was tested using the J-M test.

Table. VI. Risk level class per each cluster through the years investigated and the standard deviation per each class.

Since the Meta® data are not freely available for all years, only the 2020 map is reported (mapped data are available on request) and the tabular data is shown in the Table VII following Figure 4.

Figure. 4. Habitat suitability maps for ticks (SIM) involving the following years: (A) 2022, (B) 2021, (C) 2020. EPSG 23032 (UTM WGS84 32 N). Representation scale 1 : 350.000.

Table VII below presents the Aosta Valley population's risk of exposure to ticks over the years based on the SIM. Additionally, in the supplementary material section, the distribution of population within each risk class across the municipalities of Aosta Valley is provided (Tables SM II- SM IV).

Table. VII. Population potentially exposed expressed in numbers of people through the years analyzed.

The data reported in Table VII have also been graphed in a histogram (see Figure 5)

Figure. 5. Population potentially exposed to ticks and relative class risk probability.

Finally, on the basis of the SIM for each of the three years, the possible shift and expansion in altitude of the areas favorable to ticks was analyzed. Therefore, the SIM maps were firstly thresholded according to the SIM value>70 considering only the optimal conditions for tick survival and reproduction. This threshold was chosen on the basis of the discovery data and empirical tests (please look at https://www.epicentro.iss.it/en/ticks/overview last accessed April 10, 2024). It should also be noted that the classes with the greatest exposure fall within this threshold. Then the maps were vectorized and zonal statistics were conducted based on the DTM of the Aosta Valley, calculating the deltas between t₀ (2020) and t_n (2022) as indicated in the Table VIII and Figure 6 below.

Table. VIII. Zonal statistics involving altitude on SIM.

Figure. 6. The most favorable areas identified from the Habitat Suitability Maps for Tick Distribution over the years 2020, 2021, and 2022.

Discussion

Tick-borne diseases pose a growing threat to animal and human health worldwide. To effectively address this challenge, we need to understand where ticks thrive and predict future habitat suitability changes. This is where geomatics and remote sensing step in, offering powerful tools to map tick landscapes and guide a One Health approach that considers the interconnectedness of animal, human and environmental health. Integrated surveillance of wildlife communities, ticks and their pathogens and their possible interactions are topics to be studied in deep, in order to prevent necessary biodiversity restoration from increasing tick-borne disease risks (Kohler et al., 2010).

In this study, we suggest a possible application contributing to the domain of One Health. Zoonotic diseases, those transmissible between animals and humans, pose a growing threat, especially in regions like the Aosta Valley where wildlife frequently interact with domestic animals and people, both through recreation and hunting activities. Climate change might further amplifiy this risk. This study's data have shown that areas suitable for ticks, even at high altitudes, have shifted, increasing the risk for both wildlife and human populations. The study identified three tick species, Ixodes ricinus (93.4%), Dermacentor marginatus (5.5%) and Dermacentor spp. (1.1%), with Ixodes ricinus carrying the highest burden of zoonotic pathogens like Anaplasma phagocitophylumspp., B. burgdorferi sensu lato and different species of Rickettsia spp. A specificity of Rickettsia species for a one-tick host genus (or species) has been observed (Perlman et al., 2006). In this study, we found one R. slovaca in D. marginatus tick: this is in accordance with recent studies (Silva-Pinto et al., 2014; Socolovschi et al., 2009) that consider this pathogen an emerging zoonotic species carried by its main vector Dermacentor tick genus, particularly D. marginatus, which is the most frequent vector for the human transmission of this pathogen. Indeed, a broad range of Rickettsiae has been detected in Ixodes ricinus, the most widespread hard tick species in European countries. This tick mainly appears to be a competent vector for R. helvetica and R. monacensis (Portillo et al., 2015) as we found in our study: 8 ticks positive for R. helvetica and 4 ticks positive for R. monacensis. Interestingly, roe deer and wild boar, the most hunted animals, were also the most frequent wildlife tick hosts: wild animals could have a large impact on tick epidemiology and represent a useful tool as a sentinel species to characterize and monitor tick populations and related risks for tick-borne pathogen spread.

In 2020 and 2021, the distribution of areas most conducive to ticks remained relatively stable. However, in 2022, an expansion of favorable habitats occurred following a heatwave and severe drought, mitigated in Valle d'Aosta by the implementation of irrigation systems (Figure 6). Compared to previous years, there was an increase of over 600 meters in maximum altitude and approximately 400 meters in average altitude. While the optimal altitude limit previously ranged around 900-1000 meters, except for high-altitude environments with specific topographical features, in the last year, this limit has shifted to 1300-1400 meters, with favorable habitats even at higher altitudes, as indicated by maximum values in the statistics. Although 2022 was anomalous, it aligns with projected future climate trends and findings also in other mountains areas (Martello et al., 2014). If water, particularly vegetation moisture, is not a limiting factor, ticks appear to thrive under optimal conditions at altitudes previously considered less hospitable.

In terms of the opportunities, limitations and future prospects of using geomatics and remote sensing in this context, we start with opportunities: the approach presented can be readily scaled from national to global scales. These technologies can identify vulnerable areas, enabling focused territorial management with the development and application of One Health policies. The suggested workflow allows for daily, monthly, or yearly updates, supporting rapid response to observed changing patterns and offering the potential for a population alert and tick monitoring service. Creation of a health databases with georeferenced data from the ground from field surveys would be required. To achieve this goal a collaboration would be necessary by creating common GIS databases and EO Data Health services (similar to Copernicus Health), where diverse specialists can work together for broader impact. Finally, the creation of such a technological infrastructure will empower modeling and forecasting. In fact, Multitemporal data allows for robust modeling, revealing long-term trends and predicting future suitability changes.

Concerning limits and constraints, the spatial and temporal resolution of EO data used in this study is the principal point of the issue. Low-resolution data can skew results especially in geomorphologically complex areas like mountains. The lowest native spatial resolution input deeply affects the final suitability map (SIM). Future high-resolution missions like the space program IRIDE ASI ESA hold great promise in terms of boosting the health technological transfer offered by Earth Observation. Furthermore, the Meta® population data are not entirely free and open-source, therefore limitations in availability and temporal coverage has to be considered using this dataset. Additionally, extensive ground truth data is a crucial key point for validation, highlighting the need for app-based data collection solutions. In this study, we have not properly checked the consistency of our results using a robust validation dataset. It is worth noting that rigorous methods are essential to avoid misleading results; in fact, changes in input variables can significantly impact outcomes. This also applies to wildlife data; incomplete georeferenced wildlife data may affect the SIM. As previously explained in this study, in some cases centroids have been used, and this is a constraint factor due to the lack of correct spatial information.

Furthermore, it has to be considered that georeferencing ground data of wildlife hunted or found dead can provide valuable insights into species distribution, habitat use, and mortality patterns. However, there are several limitations and potential biases associated with this approach. Regarding biases in surveying it has to be considered: a) capture effort; the frequency and intensity of hunting or monitoring efforts can vary spatially and temporally, leading to uneven sampling across the study area. Areas with higher human activity or accessibility may be overrepresented, while remote or inaccessible regions may be underrepresented; b) Time of Year: wildlife mortality rates and hunting activities can fluctuate seasonally due to factors such as breeding, migration, or changes in environmental conditions. Biases may arise if sampling efforts are not evenly distributed throughout the year, resulting in skewed representations of species distribution and mortality patterns; c) Target Species: hunters may selectively target certain species based on factors such as abundance, trophy value, or cultural preferences. This may introduce bias into the dataset, affecting the accuracy of species distribution models or mortality estimates; d) Detection Bias: dead wildlife may not be uniformly detected or reported, leading to underestimation of mortality rates or biased spatial distributions. Factors such as carcass concealment, scavenger activity, or observer detection abilities can influence the likelihood of detecting carcasses; e) Georeferencing ground data relies on accurate GPS coordinates to precisely locate wildlife observations. However, GPS errors, such as signal interference or device limitations, can introduce positional inaccuracies, affecting the reliability of spatial analyses and mapping efforts; f) Habitat Heterogeneity: wildlife mortality events may occur in diverse habitats with complex topography, vegetation cover, or land use patterns. Spatial inaccuracies in georeferenced data may obscure fine-scale habitat associations or hinder the identification of environmental drivers influencing mortality patterns; g) data Compatibility: ground data collected by hunters or wildlife monitoring programs may vary in format, quality, or metadata standards, posing challenges for data integration and harmonization. Inconsistent data formats or missing information may hinder cross-study comparisons or collaborative research efforts.

In summary, while georeferencing ground data of hunted or found dead wildlife offers valuable opportunities for ecological research and conservation, it is crucial to acknowledge and address potential biases, limitations, and ethical considerations associated with surveying methodologies and data interpretation. Implementing standardized protocols, integrating multiple data sources, and conducting rigorous validation procedures can enhance the reliability and robustness of georeferenced wildlife datasets for scientific analysis and management decision-making.

Concerning future perspectives, the aforementioned issues may be easily overcome by developing One Health and wildlife apps for comprehensive data collection to geodatabases. Future missions and higher resolution data may certainly improve SIM. The IRIDE data offers great possibilities combining higher spatial and temporal resolutions and will improve accuracy and relevance of such studies and applications (Orusa et al., 2023). Another interesting potential component to this could be real-time monitoring with fixed-point analysis, similar to long-term ecological research stations (LTER), combined with cloud-based data processing, which can enable near-real-time health monitoring. Moreover, to facilitate a true One Health approach, an aggregated free-of-charge Dynamic population tracking dataset would be very useful. Particularly, dynamic data on human and animal movements (e.g., from Google or better public institutions) holds potential for developing tailored services. At the same time, matching farms with animal numbers can aid veterinarians in targeted interventions. Last but not least, the improvement of machine learning techniques and AI fed with consistent georeferenced tick data, wildlife, domesticated animal and human populations and environmental data offers enormous potential to developing disease forecast models. This will be especially important in alpine regions more sensitive to the influences of climate change. Consequently, in light of what has been illustrated, Geomatics and Remote sensing offer invaluable tools for mapping tick suitability and supporting a One Health approach (Caimotto et al., 2020; Orusa et al., 2023; Orusa et al., 2024b; Viani et al., 2024). Addressing data limitations and leveraging future technological advancements will unlock even greater potential in mitigating tick-borne disease risk. By working collaboratively, these technologies can map not just ticks, but a healthier future for all.

Conclusion

This study has developed a scalable approach to mapping tick suitability area (SIM) and possible monitoring services. At the same time, it highlights a worrying trend: exposure to high-risk tick zones in the Aosta Valley has increased, likely due to rising temperatures. Additionally, "tick friendly" areas have spread to higher altitudes. In fact, the mean Δ_2022-2020has increased by 400 m. Furthermore, the risk of exposure of the population to ticks has increased in the last year due to the more favorable conditions and expansion of the more suitable areas. The majority of the Aosta Valley population lives on the bottom of the valley, in areas normally at higher risk but with the occurrence at high altitude of conditions more favorable to parasites, the classes at greater risk (SIM >70%) have also involved many lateral valleys and the related population. The successful integration of remote sensing and GIS tools, such as Meta® Population, demonstrates the potential for a technological transfer to the veterinary sector, aligning with the One Health approach. Moving forward, this research underscores the need for: Enhanced surveillance with an ongoing monitoring of tick distribution and pathogen prevalence, using remote sensing and GIS tools, to inform targeted interventions. A key point in further application will be the development of an app to collect point geo-referenced samples. Public awareness through the education of the public about tick-borne diseases, high-risk areas, and preventive measures, especially for those living at lower and middle altitudes. Then, veterinary collaboration with experts in disciplines such as geomatics, foresters, planners, doctors and policy makers to develop effective tick control strategies for animals, particularly the most hunted species. Finally, One Health implementation, fostering collaboration between human and animal health sectors, alongside environmental agencies, to create a comprehensive One Health approach for managing tick-borne disease risks like in the Aosta Valley. By embracing these actions, we can leverage the power of technology and collaboration to protect the health of both humans and animals in the face of a changing climate and evolving disease landscapes.

Acknowledgements

A special thanks to the National Centre for Anaplasma, Babesia, Rickettsia and Theileria (C.R.A.Ba.R.T.) of the Experimental Zooprophilactic Institute of Sicily for the support during the work.

Data Availability Statement

The GIS and Remote Sensing datasets and layer for this study may requested to enrico.borgogno@unito.it for the veterinarians to riccardo.orusa@izsto.it

Author Contributions

Conceptualization, A.V., T.O.; Data curation, A.V., T.O. and E.B.M.; Formal analysis, A.V., T.O. M.L.M, and S.R.; Investigation, A.V., T.O.; Methodology, A.V., T.O., R.O., M.L.M., E.B.M.; Project administration, R.O. and E.B.M.; Resources, A.V. and R.O.; Software, T.O. and E.B.M.; Supervision, S.R., M.L.M., E.B.M., R.O.; Validation, A.V., T.O., E.B.M., R.O.; Visualization, A.V., T.O.; Writing—original draft, A.V., T.O.; Writing—review & editing, AV., T.O., C.N., S.R., R.O. and E.B.M. All authors have read and agreed to the published version of the manuscript.

Competing interest statement

The authors declare no conflict of interest.