Geomatics and EO Data to Support Wildlife Diseases Assessment at Landscape Level: A Pilot Experience to Map Infectious Keratoconjunctivitis in Chamois and Phenological Trends in Aosta Valley (NW Italy)

: Geomatics and satellite remote sensing o ﬀ er useful analysis tools for several technical-scientiﬁc ﬁelds. This work, with reference to a regional case of study, investigates remote sensing potentialities for describing relationships between environment and diseases a ﬀ ecting wildlife at landscape level in the light of climate change e ﬀ ects onto vegetation. Speciﬁcally, the infectious keratoconjunctivitis (IKC) of chamois ( Rupicapra rupicapra L.) in Aosta Valley (NW Italy) was investigated at the regional level. IKC ( Mycoplasma conjunctivae ) is a contagious disease for domestic and wild ruminants ( Caprinae and Ovinae ). Two types of analysis were performed: one aimed at exploring by remotely sensed data phenological metrics (PMs) and evapotranspiration (ET) trends of vegetation in the area; one investigating the correlation between PMs and ET, versus IKC prevalence. The analysis was based on TERRA MODIS image time series ranging from 2000 to 2019. Ground data about IKC were available for a shorter time range: 2009–2019. Consequently, PMs and ET trend investigations were focused on the whole times range (2000–2019); conversely, correlation analysis was achieved with reference to the reduced 2009–2019 period. The whole study was based on freely available data from public archives. MODIS products, namely MOD13Q1 v.6 and MOD16A2, were used to derive PM and ET trends, respectively. Digital Terrain Model describe local topography; CORINE Land Cover map was adopted to describe land use classes. PMs and ET (as derivable from EO data) proved to signiﬁcantly changed their values in the last 20 years, with a continuous progressive As far as correlation analysis was concerned, ET and some PMs End of Season (EOS) and Length of Season (LOS) condition IKC to results, the proposed methodology can be retained as an e ﬀ ective tool for supporting public health and eco-pathological

. The possible role of remote sensing in the patho-system as represented by the so-called disease triangle.
A good example is represented by atmospheric pollution that was recognized to increase sensitivity to pulmonary diseases, as the last pandemic event, coronavirus (SARS-CoV-2), has suggested [15].
GIS studies about endo-and ecto-parasitoses of veterinary interest, with particular reference to zoonoses agents, represent today the greatest contribution to veterinary and faunistic sectors [16]. For many years, the World Organization for Animal Health located in Paris (OIE) and the World Health Organization (WHO) in Geneva have been underlining the importance of geomatics and remote sensing applications [17] in the One Health perspective.
This work, with reference to a regional case study, investigates remote sensing potentialities for describing relationships between environment and diseases affecting wildlife at landscape level. Moreover, it is intended to describe the effects of climate change onto the vegetation component, with special concern about pastures. The study area corresponds to the entire Aosta Valley Region located in the Italian Western Alps. In particular, a new analysis approach is presented to operate at landscape level to analyze if and how environmental factors could condition the occurrence of infectious keratoconjunctivitis (IKC, Mycoplasma conjunctivae) in chamois. IKC is a contagious disease for domestic and wild ruminants (Caprinae and Ovinae) [18]. In chamois, the disease can be serious [19] and, as in other wild ruminants, blindness can occur [20], with consequent death of the animal from trauma (e.g., fall from cliffs or starvation) [21]. The period of mountain pasture is risky for the potential contact between domestic and wild infected animals; over the years, several outbreaks have been reported in wild ungulates in the Alps [22] and this is the reason that makes monitoring/surveillance plans still active. IKC caused by Mycoplasma conjunctivae is a complex disease of domestic and wild Caprinae, with great variations in the clinic-pathological and epidemiological picture. In wildlife, IKC is sometimes associated with high mortality [23,24]. It has been suggested that the pathogenesis of IKC is influenced by host predispositions, virulence of M. conjunctivae strains, secondary infections, and environmental factors [25]. Sex and age imbalance in affected populations were observed in severe outbreaks [26], indicating that age and social behavior, including sexual segregation, may be important risk factors. Differently, differences in virulence between different strains do not seem to play a major role; mycoplasmal load is obviously associated to the presence and severity of signs. However, the driver of mycoplasmal multiplication in the host is unknown. Environmental factors might have a role, regarding both the expression of the disease in individual cases and the onset of an outbreak in a population [24]. The underlying hypothesis of this work is that remote sensing could support comprehension of the role of environmental patterns in conditioning IKC patho-system, and related pathologies, as for other diseases. Altitude, air quality, and UV light have been discussed as possible predisposing factors for IKC in wild ungulates along with overcrowding [27]. Multiple outbreaks of IKC in Alpine ibex and Alpine chamois populations This work, with reference to a regional case study, investigates remote sensing potentialities for describing relationships between environment and diseases affecting wildlife at landscape level. Moreover, it is intended to describe the effects of climate change onto the vegetation component, with special concern about pastures. The study area corresponds to the entire Aosta Valley Region located in the Italian Western Alps. In particular, a new analysis approach is presented to operate at landscape level to analyze if and how environmental factors could condition the occurrence of infectious keratoconjunctivitis (IKC, Mycoplasma conjunctivae) in chamois. IKC is a contagious disease for domestic and wild ruminants (Caprinae and Ovinae) [18]. In chamois, the disease can be serious [19] and, as in other wild ruminants, blindness can occur [20], with consequent death of the animal from trauma (e.g., fall from cliffs or starvation) [21]. The period of mountain pasture is risky for the potential contact between domestic and wild infected animals; over the years, several outbreaks have been reported in wild ungulates in the Alps [22] and this is the reason that makes monitoring/surveillance plans still active. IKC caused by Mycoplasma conjunctivae is a complex disease of domestic and wild Caprinae, with great variations in the clinic-pathological and epidemiological picture. In wildlife, IKC is sometimes associated with high mortality [23,24]. It has been suggested that the pathogenesis of IKC is influenced by host predispositions, virulence of M. conjunctivae strains, secondary infections, and environmental factors [25]. Sex and age imbalance in affected populations were observed in severe outbreaks [26], indicating that age and social behavior, including sexual segregation, may be important risk factors. Differently, differences in virulence between different strains do not seem to play a major role; mycoplasmal load is obviously associated to the presence and severity of signs. However, the driver of mycoplasmal multiplication in the host is unknown. Environmental factors might have a role, regarding both the expression of the disease in individual cases and the onset of an outbreak in a population [24]. The underlying hypothesis of this work is that remote sensing could support comprehension of the role of environmental patterns in conditioning IKC patho-system, and related pathologies, as for other diseases. Altitude, air quality, and UV light have been discussed as possible predisposing factors for IKC in wild ungulates along with overcrowding [27]. Multiple outbreaks of IKC in Alpine ibex and Alpine chamois populations have been described in literature [28]. Different outbreaks of infectious keratoconjunctivitis (IKC) affecting alpine chamois and ibex in the western and central Swiss Alps and Aosta Valley were recorded in the period 2001-2019 [29]. conjunctival swabs by means of a nested PCR in 27 of the 28 chamois tested. The outbreaks occurred in an area covering 1590 km 2 . Deep valleys acted as a barrier to the spread of the disease. Many of the affected animals were juveniles, and more females than males died of IKC. The disease was more common during the summer and autumn. In some outbreaks, mortality can reach 30 percent, as, for example, in chamois in Italy, France, and Switzerland, and hundreds of chamois may die. Major outbreaks were recorded in 2001-2003 [30] and 2016-2018 [20]. With these premises, in this work, two types of analysis were performed: one aimed at exploring, by remotely sensed data, phenological metrics (PMs) and evapotranspiration (ET) trends of vegetation; one investigating correlation between PMs and ET versus IKC prevalence. PMs/ET analysis was based on TERRA MODIS image time series ranging from 2000 to 2019. Ground data about IKC were available for a shorter time range: 2009-2019. Consequently, PMs and ET trends investigation were done for the whole times range (2000-2019); conversely, correlation analysis was achieved with reference to the 2009-2019 period.

Study Area
The study area corresponds to the entire Aosta Valley, an administrative region located in the Northern West Alps of Italy, close to the border with France and Switzerland. It is an alpine region that hosts the highest peaks in Europe ( Figure 2). It sizes about 3263 km 2 and has a population of about 126,000 people [31]; it is the smallest, least populous, and least densely populated region of Italy. Despite of this, it is one of the most abundant areas in terms of fauna and flora, and biodiversity in general, hosting many protected areas such like Gran Paradiso National Park, Mont Avic Regional Park, and Mont Mars reserve. For this reason, it is an open-cell laboratory especially for environmental and biological sciences. Consequently, wildlife diseases and zoonosis can be proficiently studied in this area. In Figure 3, boundaries of regional protected areas are reported.
have been described in literature [28]. Different outbreaks of infectious keratoconjunctivitis (IKC) affecting alpine chamois and ibex in the western and central Swiss Alps and Aosta Valley were recorded in the period 2001-2019 [29]. Between the years 2001 and 2003, in Switzerland, Mycoplasma conjunctivae was identified from conjunctival swabs by means of a nested PCR in 27 of the 28 chamois tested. The outbreaks occurred in an area covering 1590 km 2 . Deep valleys acted as a barrier to the spread of the disease. Many of the affected animals were juveniles, and more females than males died of IKC. The disease was more common during the summer and autumn. In some outbreaks, mortality can reach 30 percent, as, for example, in chamois in Italy, France, and Switzerland, and hundreds of chamois may die. Major outbreaks were recorded in 2001-2003 [30] and 2016-2018 [20]. With these premises, in this work, two types of analysis were performed: one aimed at exploring, by remotely sensed data, phenological metrics (PMs) and evapotranspiration (ET) trends of vegetation; one investigating correlation between PMs and ET versus IKC prevalence. PMs/ET analysis was based on TERRA MODIS image time series ranging from 2000 to 2019. Ground data about IKC were available for a shorter time range: 2009-2019. Consequently, PMs and ET trends investigation were done for the whole times range (2000-2019); conversely, correlation analysis was achieved with reference to the 2009-2019 period.

Study Area
The study area corresponds to the entire Aosta Valley, an administrative region located in the Northern West Alps of Italy, close to the border with France and Switzerland. It is an alpine region that hosts the highest peaks in Europe ( Figure 2). It sizes about 3263 km 2 and has a population of about 126,000 people [31]; it is the smallest, least populous, and least densely populated region of Italy. Despite of this, it is one of the most abundant areas in terms of fauna and flora, and biodiversity in general, hosting many protected areas such like Gran Paradiso National Park, Mont Avic Regional Park, and Mont Mars reserve. For this reason, it is an open-cell laboratory especially for environmental and biological sciences. Consequently, wildlife diseases and zoonosis can be proficiently studied in this area. In Figure 3, boundaries of regional protected areas are reported.    Table 1 reports altitude ranges and their spatial distribution, respectively. Aosta Valley climate is strongly influenced by topography. High mountains that surround the region prevent the access of humid air masses of Mediterranean or Atlantic origin. This makes the local climate characterized by a high degree of aridity in the central area (center of the main bottom valley), with rainfalls lower than 500 mm·y −1 ; border areas, in particular south-eastern ones and north-western valleys, differently, show an average rainfall above 1400 mm·y −1 [32].
In winter, precipitations are mainly snowy. In summer, convective rainfalls are quite common, determining frequent thunderstorms; spring and autumn are characterized by stratified rainfalls, possibly lasting for several days, with a consequent increase of flood risk.
In force of its highly variegated climatic situation, Aosta Valley land cover is heterogenous. In the present study, particular attention was paid to low natural vegetation where IKC diffusion is highly possible due to the interaction and competition for grasses among animals, in particular wildlife ruminants and breed animals [33][34][35][36]. In Aosta Valley, where breeding activities and ungulates wildlife are very abundant, animal competitive interaction generally occurs in mesotrophic  Table 1 reports altitude ranges and their spatial distribution, respectively. Aosta Valley climate is strongly influenced by topography. High mountains that surround the region prevent the access of humid air masses of Mediterranean or Atlantic origin. This makes the local climate characterized by a high degree of aridity in the central area (center of the main bottom valley), with rainfalls lower than 500 mm·y −1 ; border areas, in particular south-eastern ones and north-western valleys, differently, show an average rainfall above 1400 mm·y −1 [32].
In winter, precipitations are mainly snowy. In summer, convective rainfalls are quite common, determining frequent thunderstorms; spring and autumn are characterized by stratified rainfalls, possibly lasting for several days, with a consequent increase of flood risk.
In force of its highly variegated climatic situation, Aosta Valley land cover is heterogenous. In the present study, particular attention was paid to low natural vegetation where IKC diffusion is highly possible due to the interaction and competition for grasses among animals, in particular wildlife Remote Sens. 2020, 12, 3542 6 of 22 ruminants and breed animals [33][34][35][36]. In Aosta Valley, where breeding activities and ungulates wildlife are very abundant, animal competitive interaction generally occurs in mesotrophic and eutrophic pastures or areas with good pastoral values [37][38][39]. Under this hypothesis, supported only by empirical evidence, an attempt was made to analyze whether a "favorable" type of land cover could play as hotspot for diseases spreading. From an ecological point of view, a greater competitive interaction among animals for the same resource was expected in these areas. Under this scenario, environmental analysis could drive to a better comprehension of the epidemiological relationships linking disease and hosts.

Veterinary Ground Samples
The most of chamois samples were collected at the hunting wildlife control centers of RFD (Regional Forestry Districts) managed by Corpo Forestale della Valle d'Aosta (Forest Guards). Samples were obtained from chamois that were shot by hunters and analyzed by veterinary officers of CeRMAS (National Reference Center for Wildlife Diseases). Map coordinates of locations where animals were collected were not available. Consequently, data from all the analyses were achieved at the regional level.
As far as chamois analysis is concerned, swabs were collected under the third eyelid from both eyes, transported and stored refrigerated until analyzed. Ocular clinical signs of conjunctivitis and keratitis (ocular damage, inflammation or discharged) or signs of blindness as indicated by abnormal behavior such as stumbling, circling, uncertain gait or inability to climb were recorded at sampling. Age, sex, and body condition were also recorded. Some samples came from chamois that were found dead at the ground showing signs of conjunctivitis or keratitis. The sampling period ranges between 2009 and 2019, when several severe outbreaks of chamois IKC occurred in the Alps [24,40].
At the laboratory, eye swabs were placed into sterile tubes with 0.45 mL of lysis buffer (1 M Tris-HCl, pH 8.5, 0.5% Tween 20, 0.24 mg/mL proteinase K) and mixed for 30 s. Cells were lysed for 60 min at 60 • C and, successively, heated up to 95 • C for 15 min in order to inactivate proteinase K. Obtained lysates were tested with a specific real time PCR to assess the presence of M. conjunctivae; according to Vilei et al. 2007 [41] TaqMan real time PCR reactions were performed by using 2.5 µL of test sample, 900 nM of lppS forward primer (5 -CAGCTGGTGTAGCACTTTTTGC-3 ) and lppS reverse primer (5 -TTAACACCTATGCTCTCGTCTTTGA-3 ), 300 nM of lppS probe (5 -TGCTTCGACTACCAAATATGATGGTGATCCTCT-3 with 6FAM reporter dye and TAMRA quencher affixed on the 5 and 3 ends, respectively), and TaqMan Universal PCR Master Mix in a 25 µL volume. An exogenous Internal Positive Control was introduced for all reactions to check for the presence of eventual PCR inhibitors. PCR reactions were run by StepOne Plus instrument (Thermo Fisher) using the following cycling parameters: one step at 50 • C for 2 min and at 95 • C for 10 min, 40 cycles of denaturation at 95 • C for 15 s, and extension at 60 • C for 1 min were performed. Real-time fluorescence measurements were taken for each sample by using the StepOne™ Software v2.3 (Thermo Fisher Scientific: https://www.thermofisher.com/order/catalog/product/4376357#/4376357) and the PCR cycle number at which the fluorescent signal crossed the cycle threshold (set manually) was recorded as CT value. The fluorescence emission baseline was set manually two cycles before the cycle with a significant fluorescence signal. The specificity of the TaqMan assay was evaluated by testing genomic DNA of mycoplasmas other than M. conjunctivae and of other ocular pathogens [41,42].
Prevalence data were finally computed by Equation (1) with reference to the entire regional territory for all the monitored years (2009-2019).
where Pr = disease prevalence (%), C = number of positive disease cases detected by PCR and optical analysis from samples, P = number of examined chamois.

EO and Geographical Digital Data
MOD13Q1 v.6 [43] product from the NASA TERRA Moderate Resolution Imaging Spectroradiometer (MODIS) mission was used to map the Normalized Difference Vegetation Index (NDVI) over the area in the period 18 February 2000-31 December 2019. Four hundred-fifty two MOD13Q1 images were obtained for free by Google Earth Engine [44]. They were stacked into a NDVI time series (NTS) and filtered according to the Pixel Reliability layer, supplied together with the product, to map reliability of each scene pixel. Stacking and filtering (Savitzky-Golay, [45,46]) were operated by a self-developed IDL 8.1 routine [47]. Vegetation Indices from MOD13Q1 v.6 product are a composite one, having a time step of 16 days and a spatial resolution of 250 m. The MOD13Q1 v.6 product provides a Vegetation Index (VI) value at a per pixel basis. It is referred to as a continuity index to the existing one from the National Oceanic and Atmospheric Administration-Advanced Very High Resolution Radiometer (NOAA-AVHRR) data. Composition algorithm selects the best available pixel value from all the acquisitions within the considered 16 days period. Selection criteria are the following: no clouds, low viewing angle, highest NDVI value adopting the classic formula as follow: Equation (2).
where ρNIR and ρRED are the at-the-ground reflectance of the near-infrared (MODIS band 2) and red (MODIS band 1) bands, respectively. Evapotranspiration (ET) maps were retrieved from the MOD16A2 v.6 collection [48] by Google Earth Engine for the period 1 January 2000-31 December 2019 for a total of 868 images.
MOD16A2 v.6 Evapotranspiration/Latent Heat Flux product is an 8-day composite product with a geometric resolution of 500 m. The algorithm adopted for the MOD16A2 data product collection refers to the Penman-Monteith equation [49], which includes inputs of daily meteorological reanalysis data such as albedo, land cover, vegetation property dynamics, and Land Surface Temperature (LST) that in some studies were adopted to evaluate ET [50] or heat fluxes [51]. Pixel values for the Net Evapotranspiration (ET 0 ) is the sum of all 8 days within the composite period expressed in kg m −2 8 d −1 . MOD16A2 v.6 ET layers were stacked into an ET time series (ETS); no filtering was applied. Stacking was operated by a self-developed IDL 8.1 routine.
The Shuttle Radar Topography Mission (SRTM) digital elevation model was used for this work [52,53]. SRTM v3 product (SRTM Plus) is provided by NASA JPL with a grid size of 1 arc-second (approximately 30 m). Native SRTM was resampled at 30 m (SRTM30) to make pixel squared.
All data were converted from native geographical reference systems into the ED50 UTM 32N one. MOD13Q1 v.6 and MOD16A2 v.6 layers were geometrically oversampled up to 30 m to refine area zonation.

Land Cover Data
To properly describe land cover at the regional level, the Corine Land Cover 2018 dataset (hereinafter called CLC2018) was used; CORINE (Coordination of Information on the Environment) Land Cover inventory was initiated in 1985 to standardize data collection on land in Europe to support environmental policy development. The project is coordinated by the European Environment Agency (EEA) in the frame of the EU Copernicus program and implemented by national teams. The number of participating countries has increased over time currently including 33 member countries and six cooperating countries with a total area of over 5.8 million km 2 . CLC2018 specifically, is one of the available datasets produced within the general CORINE frame and refers about land cover/land use status in 2018. The reference year of the first CLC inventory was 1990 and the first update was achieved in 2000. The current update cycle is 6 years. Satellite imagery provides the geometrical and thematic basis of maps; it is integrated with in-situ data as essential ancillary information.  Table 2. LC classes were considered with reference to altitude ranges, in order to give a more comprehensive description on the ongoing processes affecting vegetation and its relationship with IKC prevalence.
As indicated in Table 2, this study only considered the following CORINE land cover types: pastures, natural grassland, moors, and heathlands. These classes were selected since majorly favoring the interaction between breed animals and ungulates. In Aosta Valley, where breeding activities and ungulates wildlife are very abundant, this competitive interaction generally occurs in mesotrophic and eutrophic pastures or areas with good pastoral values. Under this hypothesis, supported only by empirical evidence, an attempt was made to analyze whether a favorable type of land cover could play as a hotspot for diseases spreading. From an ecological point of view, a greater competitive interaction among animals for the same resource was expected in these areas. Under this scenario, environmental analysis could drive to a better comprehension of the epidemiological relationships linking disease and hosts.
CLC2018 and SRTM30 were used as descriptors of environmental conditions possibly conditioning occurrences of IKC. Consequently, IKC was preliminarily tested against favorable land cover type extent with respect to 3 different altitude classes (Table 3).

Methodology
To investigate possible relationships of IKC prevalence with vegetation-related factors, two analyses were performed. One (hereinafter called Analysis 1) was aimed at testing if any significant climatic trend could be recognized affecting phenological metrics (PMs, see forward on) and ET (as measured by NTS and ETS, respectively), in the period 2000-2019.
A second investigation (hereinafter called Analysis 2) was, conversely, addressed at verifying if any significant correlation could be recognized between IKC prevalence and some of PMs that showed significant trends in the previous analysis.
On the basis of the IKC disease dataset, prevalence was calculated by Equation (1). IKC yearly prevalence computation relied on the entire chamois population inventory as reported by the local faunistic centers: Corpo Forestale della Valle d'Aosta and hunter committees. Veterinary analyses were performed by IZS PLV SC Aosta and CeRMAS. Faunistic season was assumed to start in September and terminate at the end of the next year August.
Since interaction between breed animals and ungulates is more likely in good quality grass areas, rangelands could represent hotspots for the spreading of the disease.
CLC2018 classes 231 = pastures, 321 = natural grasslands, 322 = moors and heathland were a-priori assumed as favorable land cover classes to test IKC prevalence against to (see Figure 4). They were aggregated into a single macro-class, hereinafter called "FAV" (see Figure 5).

Analysis 1: testing PM and ET trends from NTS and ETS
To investigate if and how IKC prevalence could be also related to ongoing changes affecting vegetation activity in consequence of climate change, some metrics were extracted from the abovementioned time series (NTS and ETS). Phenological metrics (PMs) are synthetic descriptors of vegetation activity along its annual growing season. Climate change proved to condition such activity shifting and reshaping past "ordinary" behavior of plants along the year. The following PMs were considered: the start of the growing season (SOS) representing the day of the year (DOY) when phenology is admitted to boost; the end of the season representing the day of the year (DOY) when phenology is admitted to stop; the length of the growing season (LOS) representing the time range (in number of days) separating EOS from SOS; the maximum of NDVI (MAXVI) representing the highest value reached by NDVI during the growing season and proved to be a good predictor of climate change effects on vegetation.
PMs were estimated by TIMESAT 3.3 with STL software [60][61][62] that was specifically developed to enable the monitoring of land surface processes by remotely sensed data. TIMESAT 3.3 with STL [63][64][65] iteratively fits and smooths by mathematical functions the yearly NDVI time-series, finding the best smoothed approximation of the NDVI along the year at pixel level. Once raw data have been approximated by the selected fitting function, PMs can be extracted in correspondence of singular points having a phenological meaning (e.g., EOS, SOS, LOS, MAXVI, etc.) along the local temporal profile.
With reference to ETS and separately for FAV1, FAV2, FAV3, class yearly average ET value was computed and analyzed along the years in the period 2000-2019 at the regional level, to explore eventual trends in ET values, as well.
Since analysis is performed at single pixel and year level, PM estimates are saved as raster layers showing PM spatial distribution at the considered year.

Analysis 2: IKC prevalence vs. PMs/ET
After demonstrating that significant trends could be recognized affecting PMs and ET values from MOD13Q and MOD16A2 datasets, respectively, we tested their potential correlation with IKC Pr at the regional level. As far as IKC/ET comparison was concerned, the yearly cumulative ET was calculated. The a priori hypothesis was that a change in environmental conditions could drive to a change in IKC occurrences (Pr value). It is worth to remind that Pr data were supplied aggregated at the regional level, while PMs and ET measures were mapped at the pixel level over the whole area. Consequently, all deductions refer to general trends that could be possibly improved if more distributed and geolocated data of IKC Pr were available.
IKC Pr was tested against all the computed metrics and modeled by a 2nd order polynomial regression.
priori assumed as favorable land cover classes to test IKC prevalence against to (see Figure 4). They were aggregated into a single macro-class, hereinafter called "FAV" (see Figure 5).

Analysis 1: testing PM and ET trends from NTS and ETS
To investigate if and how IKC prevalence could be also related to ongoing changes affecting vegetation activity in consequence of climate change, some metrics were extracted from the abovementioned time series (NTS and ETS). Phenological metrics (PMs) are synthetic descriptors of vegetation activity along its annual growing season. Climate change proved to condition such activity  (1). Values are reported in Table 4.

Analysis 1: testing PM and ET trends from NTS and ETS
In this work, phenological metrics (SOS, EOS, LOS, MAXVI) were estimated at pixel level by TIMESAT 3.3 processing the whole NTS. A Seasonal Trend decomposition by Loess (STL) was adopted to de-trend NTS pixel profile and removing noise. Seasonal component was refined by Savitzky-Golay filtering to reduce, but not removing, remaining local strong variations. The yearly growing season was recognized with the whole multi annual time series using the sinusoidal harmonics approach.
With reference to de-trended/filtered NTS profiles PMs were extracted using the simple thresholding approach. An arbitrary value of 0.5 was set as reference NDVI value to refer SOS and EOS to. Consequently, SOS was assumed as the DOY when NDVI reached the 0.5 threshold value along the ascending part of the phenological yearly bell; EOS was assumed as the DOY when NDVI reached the 0.5 threshold value along the descending part of the phenological yearly bell. LOS map was computed by grid differencing from EOS and SOS maps. MAXVI was found looking for the highest NDVI value between SOS and EOS.
PMs were mapped over the area as raster layers and spatially averaged with respect to FAV1, FAV2, and FAV3 classes ( Table 3).
Graphs of Figure 6a-c show that average PM values of FAV1, FAV2, FAV3 significantly changed their values in the last 20 years, with a continuous progressive trend observable for all of them. First order polynomials used to model trends where calibrated excluding those PMs and ET estimates/measures whose residuals (computed with respect to a 1st order polynomial) showed a value higher than mean ± 2 times the standard deviation. Interpretation of results is given in the Discussions section of this paper. In terms of strength of changes, with reference to gain values of the estimated 1st order polynomial models (Table 5), it can be observed that: EOS is averagely delaying of about 2.6 d/y and no significant difference can be observed concerning altitude classes.
Differently, SOS appears to averagely anticipate about 2 d/y up to 2000 m a.s.l. (FAV 1 and FAV2) and about 3 d/y at higher altitudes (FAV 3). Consequently, LOS showed to enlarge about 4.7 d/y at lower altitudes (FAV1 and FAV2) and about 6.5 d/y at higher altitudes (FAV3). According to Borgogno et al. [66], potential accuracy of NDVI measurements is about ±0.02; consequently, estimated yearly variations of MAXVI cannot be considered singularly significant. Nevertheless, the cumulated effects along the entire explored period (2000-2019) showed that MAXVI significantly changed, since the accuracy reference value of ±0.02 was largely overcome. MAXVI variations between 2000 and 2019 appear to be positively higher at lower altitudes (about +0.09) while almost stable as altitude increases. These results find strong evidence in different studies in literature [67][68][69][70][71][72][73][74][75][76]. Table 5. Yearly average and cumulated (2000-2019) variations of PMs and ET along the years as estimated by 1st order polynomial regression ( Figure 6). In Table 5, DOY refers to the day along the considered year and ranging between 0 and 365. Positive values indicate a delay or lengthening of season, negative values an advance of season according to the observed trends reported in Figure 6. The first part of the table indicates the yearly gain while the second part (cumulated) represent the overall changes in the observed period analyzed for each PMs from the starting year 2000 till arriving to the end in the 2019. order polynomials used to model trends where calibrated excluding those PMs and ET estimates/measures whose residuals (computed with respect to a 1st order polynomial) showed a value higher than mean ± 2 times the standard deviation. Interpretation of results is given in the Discussions section of this paper. In terms of strength of changes, with reference to gain values of the estimated 1st order polynomial models (Table 5), it can be observed that: EOS is averagely delaying of about 2.6 d/y and no significant difference can be observed concerning altitude classes. Differently, SOS appears to averagely anticipate about 2 d/y up to 2000 m a.s.l. (FAV 1 and FAV2) and about 3 d/y at higher altitudes (FAV 3). Consequently, LOS showed to enlarge about 4.7 d/y at lower altitudes (FAV1 and FAV2) and about 6.5 d/y at higher altitudes (FAV3). According to Borgogno et al. [66], potential accuracy of NDVI measurements is about ±0.02; consequently, estimated yearly variations of MAXVI cannot be considered singularly significant. Nevertheless, the cumulated effects along the entire explored period (2000-2019) showed that MAXVI significantly changed, since the accuracy reference value of ±0.02 was largely overcome. MAXVI variations between 2000 and 2019 appear to be positively higher at lower altitudes (about +0.09) while almost stable as altitude increases. These results find strong evidence in different studies in literature [67][68][69][70][71][72][73][74][75][76].
With reference to ET graphs of Figure 6d, a significant increasing trend can be observed for all FAV classes. According to the values of Table 5, it can be noted that all FAV classes behave similarly; water requirement in a period of 8 days appears to averagely increase of about 0.05 Kg·m −2 (about 0.5%) every year. This determined that, in the period 2000-2019 cumulated increment of ET (every 8 days) is around 1 Kg·m −2 corresponding to a percentage difference in water requirement of about 8%. This could be possibly explained with reference to the previously demonstrated increasing of both biomass production (MAXVI) and enlargement of the growing season, that, consequently make vegetation needing more water yearly. With reference to ET graphs of Figure 6d, a significant increasing trend can be observed for all FAV classes. According to the values of Table 5, it can be noted that all FAV classes behave similarly; water requirement in a period of 8 days appears to averagely increase of about 0.05 kg·m −2 (about 0.5%) every year. This determined that, in the period 2000-2019 cumulated increment of ET (every 8 days) is around 1 kg·m −2 corresponding to a percentage difference in water requirement of about 8%. This could be possibly explained with reference to the previously demonstrated increasing of both biomass production (MAXVI) and enlargement of the growing season, that, consequently make vegetation needing more water yearly.

Analysis 2: IKC prevalence vs. PMs/ET
After demonstrating that some PMs and ET values are currently changing with a significant linear trend, authors tested their correlation with IKC Pr. A 2nd order polynomial was found to well approximate the most of tested relationships. Model calibration was achieved after removing outliers from data. This was obtained excluding all those data that, with respect to the 2nd order polynomial model, showed a percent residual Equation (3) ε > ±100%. After outlier removal, a new calibration was run for the model. Results are reported in Table 6.
As shown in Table 6 and Figure 7, good correlations were found between IKC Pr and LOS, EOS, and ET. A 2nd order polynomial proved to well fit the relationship. Differently, SOS showed a weaker correlation for all the tested situations. attractiveness for chamois. At higher altitudes (FAV3), the relationship still persists, but it becomes weaker, probably due to a lower availability of biomass and a shorter phenological season. It is worth to remind that these results were obtained at the regional level, using highly aggregated data. It is the authors' intent to investigate further to make their deductions more robust and test the capability of generalization of the proposed prediction models in different areas. Moreover, it would be desirable to assess relationships between other diseases and environmental factors to better face future risks possibly related to zoonosis [77][78][79][80] and its dependence on climate change [61] and loss of biodiversity [81].  (d-f) show the most significant relationships that were recognized between IKC prevalence and ET (Cumulative yearly ET) values for FAV1, FAV2, and FAV3, respectively. In all the graphs, the red-orange points represent outliers (see Equation (3)).

Discussions
The functional roles of domestic and wild host populations in infectious keratoconjunctivitis (IKC) epidemiology have been extensively discussed claiming a domestic reservoir for the more susceptible wild hosts; in the most of cases all deductions were based on limited data.  (d-f) show the most significant relationships that were recognized between IKC prevalence and ET (Cumulative yearly ET) values for FAV1, FAV2, and FAV3, respectively. In all the graphs, the red-orange points represent outliers (see Equation (3)).
In particular: ET proved to be a good IKC Pr proxy that appeared to be quite independent from altitude; differently, EOS prediction capability was highly more significant for lower altitudes (FAV1); differently, LOS appeared to represent the main factor conditioning IKC spread at medium-high altitudes (FAV2 and FAV3).
As far as ET is concerned, it is well known that it relates to local micro-climatic conditions. In fact, it is directly impacted by micro-local temperature and humidity that favor vegetation growth and limit soil drought. This occurs mainly at lower altitude (FAV1 and FAV2), determining a higher attractiveness for chamois. At higher altitudes (FAV3), the relationship still persists, but it becomes weaker, probably due to a lower availability of biomass and a shorter phenological season.
It is worth to remind that these results were obtained at the regional level, using highly aggregated data. It is the authors' intent to investigate further to make their deductions more robust and test the capability of generalization of the proposed prediction models in different areas. Moreover, it would be desirable to assess relationships between other diseases and environmental factors to better face future risks possibly related to zoonosis [77][78][79][80] and its dependence on climate change [61] and loss of biodiversity [81].

Discussions
The functional roles of domestic and wild host populations in infectious keratoconjunctivitis (IKC) epidemiology have been extensively discussed claiming a domestic reservoir for the more susceptible wild hosts; in the most of cases all deductions were based on limited data.
With the aim to better assess IKC epidemiology in complex host-pathogen alpine systems, the long-term infectious dynamics and molecular epidemiology of Mycoplasma conjunctivae has been investigated in all host populations from different areas in the Pyrenees and Occidental Alps.
Between the years 2000 and 2019, it was consistently detected in Pyrenean and Alpine chamois (Rupicapra p. pyrenaica) populations, as well as in sheep flocks, and occasionally in mouflon (Ovis aries musimon) from the Pyrenees; statistically associated with ocular clinical signs only in chamois. Chamois populations showed different infection dynamics with low but steady prevalence (4.9%) and significant yearly fluctuations (0.0%-40.0%) between the period 2008-2015 [27,28]. Persistence of specific M. conjunctivae strain clusters in wild host populations is demonstrated for six and nine years. Cross-species transmission between chamois and sheep and chamois and mouflon were also sporadically evidenced. In Switzerland, the chamois affected by IKC was found at altitudes between 550 and 3200 m. The estimated overall mortality was less than 5 per cent, but more than 20 per cent have probably died locally [29]. Host population characteristics and M. conjunctivae strains resulted in different epidemiological scenarios in chamois, ranging from the fading out of the mycoplasma to the epidemic and endemic long-term persistence. These findings highlight the capacity of M. conjunctivae to establish diverse interactions and persist in host populations, also with different transmission conditions. Overall, independent M. conjunctivae sylvatic and domestic cycles occurred at the wildlife-livestock interface in the alpine ecosystems with sheep and chamois as the key host species for each cycle, and mouflon as a spill-over host. Although outbreaks of IKC have been described in Austria, France, Italy, Slovenia, and Switzerland, descriptive studies of the role of environmental patterns and to model the outbreaks on a large scale have often been incomplete, owing to the difficulty of detecting fundamental patterns that can affect IKC spread in chamois that live in remote, inaccessible mountain regions. Under this scenario, the remote sensing techniques and EO data can give certainly a huge hand in the understanding and development of possible forecasting models, as we have tried to do in the present work. With these premises, the present study was intended to explore and propose a method based on free accessible EO data to partially close the above-mentioned knowledge gap.
In Aosta Valley (NW Italy), PMs and ET (as measured from the above mentioned EO data) proved to significantly change their values in the last 20 years, with a continuous progressive trend observable for all of them. In terms of strength of changes, an average delay of EOS was observed by about 2.6 days, independently from the altitude class. SOS proved to averagely anticipate about 2 and 3 days per year at lower (<2000 m) and higher (>2000 m) altitudes, respectively. Consequently, LOS is enlarging by about 4.7 and 6.5 days per year at lower (<2000 m) and higher (>2000 m) altitudes, respectively. While looking at the entire period (2000-2019) MAXVI proved to be significantly changing, showing a positive variation (about +0.09) at lower altitude and no variations at higher one. This can be explained admitting that at lower altitudes, in Aosta Valley, grasslands and pastures are often irrigated. Consequently, farmers can vary water release regimes to face climate change effects (higher temperatures, in particular) with the result of moving forage yields (that NDVI is a predictor of) to higher values.
Differently, where more natural (not managed) systems are located (higher altitudes), the increase of yearly MAXVI can be only related to glacier melting that could compensate the increase of water requirement (as confirmed by the ET analysis) by vegetation. Glaciers are, in fact, dramatically reducing in Aosta Valley. Moreover, another compensating action could come from the surrounding forest areas that have been proved to tolerate summer heatwaves.
With reference to ET, a significant increasing trend was observed, independently from altitude. Eight day water requirement from vegetation appears to averagely increase by about 0.05 kg· m −2 (about 0.5%) every year for a total increase of about 1 kg·m −2 in 20 years (2000-2019), corresponding to a percentage difference in water requirement from vegetation of about 8%. This could be possibly explained by the increasing of biomass production (well represented by MAXVI) and by the enlargement of the growing season, that, consequently make vegetation need more water yearly.
As far as PMs/ET and IKC Pr correlation is concerned, some interesting findings came out. A 2nd order polynomial model was found to well approximate the most of relationships, making possible to support scenarios generation of IKC spreading for forecasting issues.
ET proved to be a good predictor of IKC Pr, with no significant conditioning by altitude. EOS seems to operate good predictions of Pr at lower altitudes, while LOS at medium-high ones. This probably depends on wildlife dynamics, that in autumn and at the beginning of winter, looks for grassland (food) especially at lower altitudes where, in that period, grass is not covered by snow and is wetter. In the case of SOS, a possible interpretation can rely on the fact that, if vegetative season lasting increases, animals descent from mountain to valley can be postponed, thus inducing a higher probability of interaction between potential guests and sick animals at higher altitudes. This certainly can increase also the probability of exposure of animals to disease.
With reference to ET it is mainly related to Pr at lower altitudes where micro-local temperature and humidity favor vegetation growth and limit soil drought, determining a higher attractiveness for chamois. At higher altitudes, ET capability of predicting Pr becomes weaker probably due to a lower availability of biomass and a shorter phenological season.
Authors are conscious that this work just introduces a new way to manage wildlife health problems and cannot be retained conclusive. In the nearer future, more disaggregated investigations should be done, and other areas possibly considered. Nevertheless, the proposed approach is sufficiently innovative in the context of wildlife veterinary and, we hope, could open a new interesting trend to map wildlife diseases and related zoonosis risk associated with the interaction between wild animals and domestic ones. A radical change is expected also by technicians and institutional subjects in their ordinary procedures for recording and managing ground data. In fact, the greatest limit to expand and more focus this research relied on the format of ground data that could be obtained only aggregated at the regional level with no information concerning the specific place where each analyzed animal was found. We invite all involved players to carefully consider the possibility of georeferencing every ground observation that comes to their laboratories. Georeferencing of ground data is at the basis of an effective and reliable spatial based approach like the one here proposed, where EO data (especially if available over a long-time span) play a crucial role. Anyway, this work proved that spatially based forecasting models can be reasonably calibrated for generating maps of risk concerning wildlife diseases and zoonosis spreading in a certain area. Moreover, it showed that relationships between IKC and PMs/ET are probably chancing in terms of strength; in fact, we demonstrated that all the considered predictors are suffering from a significant change possibly related to the ongoing climate change. Consequently, we expect that future approaches should more properly rely on contemporary data spatially distributed in place of aggregated data temporally distributed like the one we processed for this work.

Conclusions
A methodological approach, where geomatics and remote sensing data play a crucial role, can certainly represent a useful "addressing" tool in epidemiology and veterinarian eco-pathology fields of application, especially when studying vector-borne diseases and animal diseases and zoonosis which can affect human beings as well. A major technology transfer, thanks to remote sensing, can still be done in the health sector. A better knowledge of the patho-system dynamics, in particular the relationship between environmental component and disease presence/spreading at a landscape level, can certainly help to aim new studies and think a holistic management of how integrate adaptation, mitigation and prevention. When operative, such approaches could support decision-makers dealing with wildlife and domestic animals management and planning (both at the hunting and pastoral level) and with public health (in the perspective of One Health). The availability of free and global remote data is certainly a valid "systemic" tool for risk analysis that can support ordinary diagnostic techniques, allowing continuous monitoring of the effects that climatic and anthropogenic changes in the Alpine and, in general, mountain and wilderness environments, can cause to animals, biodiversity, and the ecosystems.