Spatiotemporal Assessment of Habitat Quality in Sicily, Italy

Abstract

We measured the spatiotemporal dynamics of habitat quality (HQ) in Sicily in two different reference years, 2018 and 2050, assuming a business-as-usual scenario. To estimate HQ and related vulnerability, we used the Integrated Valuation of Ecosystem Services and Trade-offs (InVEST) Habitat Quality model and data on land use/land cover provided by the Esri Land Cover 2050 project. We also implemented a Coarse–Filter approach to validate the reliability of HQ measures and detect biodiversity hotspots that require priority conservation. Further, we used spatial statistic tools for identifying clusters or hotspot/coldspot areas and uncovering spatial autocorrelation in HQ values. Finally, we implemented a geographically weighted regression (GWR) model for explaining local variations in the effects on HQ estimates. The findings reveal that HQ in Sicily varies across space and time. The highest HQ values occur in protected areas and forests. In 2018, the average HQ value was higher than it was in 2050. On average, HQ decreased from 0.29 in 2018 to 0.25 in 2050. This slight decline was mainly due to an increase in crop and urbanized areas at the expense of forests, grasslands, and bare lands. We found the existence of a positive spatial autocorrelation in HQ, demonstrating that areas with higher or lower HQ tend to be clustered, and that clusters come into contact randomly more often in 2050 than in 2018, as the overall spatial autocorrelation moved from 0.28 in 2018 to 1.30 in 2050. The estimated GWR model revealed the sign and the significance effect of population density, compass exposure, average temperature, and patch richness on HQ at a local level, and that such effects vary either in space and time or in significance level. Across all variables, the spatial extent of significant effects intensifies, signaling stronger localized influences in 2050. The overall findings of the study provide useful insights for making informed decisions about conservation and land planning and management in Sicily.

1. Introduction

Assessing habitat quality (HQ) is crucial for biodiversity conservation, ecosystem health, sustainable resource management, and climate adaptation. Habitat quality refers to the degree to which ecosystems can support a species or community by providing essential resources, as well as by offering protection from threats [1,2]. High levels of HQ support wildlife species and maintain rich ecosystems, while low levels of HQ lead to wildlife species and habitat degradation [3]. Several key factors affect the level of HQ; for example, food and water availability [4]; shelter and cover [5]; nesting and breeding sites [6]; biodiversity [7]; habitat connectivity [8]; land fragmentation [9]; pollution levels [10]; human impact [11]; and climate conditions [12]. Understanding which factors contribute in which ways to maintaining and enhancing HQ and its variations in space and time is pivotal to monitoring the conditions of ecosystems and productive lands and to forecasting their health into the future. Assessing HQ allows the identification of areas most in need of restoration and determination of which interventions will be most effective for enhancing habitat features; the functionality of ecosystems; and potential sites where new protected areas and ecological connectivity can be established to achieve targets established at the international or national level in line with the 2030 Kunming–Montreal Global Biodiversity Framework [13,14,15,16].

There are different tools and methods for measuring and monitoring HQ based on the spatial distance between habitat and threats caused by human activities. These tools vary depending on habitat type, species of interest, and the goals of the assessment [17]. However, to enhance the accuracy and effectiveness of HQ assessment, a combination of different HQ measuring and monitoring methods is recommended [18,19,20,21,22].

The literature on HQ assessment provides many examples of such tools being used worldwide. Different scales of analysis have been adopted, from local [23] to regional [24], and different ecosystems have been analyzed, for example, lakes [25], rivers [26], urban areas [27], coastal areas [28], forests [29], and mangroves [30].

However, only a limited number of HQ measurements relate to Mediterranean countries [11,31,32,33,34,35], despite the region being worthy of analysis. The Mediterranean Basin is a global biodiversity hotspot, hosting a vast array of unique and diverse species, many of which are found nowhere else [36]. It is an area particularly vulnerable to the effects of climate change, including rising temperatures, altered precipitation patterns, and increased drought [37]. In addition, the Mediterranean Basin is densely populated and experiences significant human activity such as urbanization, agriculture, and tourism, which pose critical threats to the natural ecosystems.

This study contributes to filling this gap in research by measuring the HQ in Sicily (Italy), the largest island in the Mediterranean Sea. We estimate the spatial and temporal patterns of HQ at two time points, 2018 and 2050, assuming a business-as-usual scenario in 2050. To assess HQ, and related vulnerability, we use the InVEST Habitat Quality model, developed under the Natural Capital Project. The InVEST Habitat Quality model allows HQ to be measured by estimating the extent of habitat types in a landscape and their state of degradation by linking land use and land cover (LULC) with habitat threats [38]. Data on LULC in the reference years are provided by the Esri Land Cover 2050 project. Other objectives of the study include the following:

(i)

Validate HQ measurements to identify areas requiring priority conservation and restoration through a coarse–filter approach [39,40,41];

(ii)

Detect hot and cold spots, as well as global and local spatial agglomeration effects through the application of Getis-Ord Gi* statistics, and global and local Moran’s statistics;

(iii)

Examine the impact of socioeconomic, topographic, climatic, and ecological variables on habitat quality (HQ) at both global and local levels, using the Ordinary Least Squares (OLS) and Geographically Weighted Regression (GWR) models [42].

All of these objectives aim to enhance the understanding of HQ dynamics and inform and guide the development of effective conservation and restoration strategies in Sicily.

2. Material and Methods

2.1. Study Context

Sicily (Italy) is the largest island (see Figure 1) in the Mediterranean Sea (37°34′08.1″ N 14°05′09.9″ E), covering approximately 25,500 km² with altitudes ranging from sea level to 3340 m at Mount Etna, the largest active volcano in Europe, and approximately 1000 km of coastline with diverse coastal habitats [43]. The remarkable geological, morphological, and climatic diversity of the island makes Sicily one of the richest biodiversity hotspots in the Mediterranean [44,45]. The dominant land use is agriculture, with agricultural land covering more than half of the island’s surface area. Key crops include wheat, barley, olives, grapes (for wine), almonds, and citrus fruits, particularly lemons and oranges. Although less extensive, pastures are also present, mainly in the island’s interior and hilly areas. Forests are found in mountainous regions, including the Madonie, Nebrodi, and Etna. Common tree species include oak, beech, and pine. Natural grasslands are found mainly in the interior and higher elevations. Currently, approximately 18.95% (ca. 489,500 ha) of Sicily’s landmass is included in areas designated for nature conservation. The regional system of protected area consists of four regional parks (Madonie Park, Nebrodi Park, Etna Park, and Alcantara River Park), 76 regional reserves, and 245 sites of the Natura 2000 network which include 17 priority habitats, and 413 species of vascular flora listed in the Annex I of the Habitats Directive, and three internationally recognized wetlands for migratory birds [46,47,48]. Urban areas are concentrated around major cities such as Palermo, Catania, Messina, and Syracuse. Sicily faces several ecological challenges, including habitat loss due to urbanization and agricultural expansion, water scarcity, fire events, and threats from invasive species. Additionally, climate change poses severe risks to the island’s fragile ecosystems [49,50].

2.2. Material and Methods

2.2.1. InVEST Model

We used the InVEST Habitat Quality model (version 3.11.0) to assess HQ in Sicily. The InVEST Habitat Quality model is the tool most used worldwide to assess HQ. Its popularity is mainly due to the ease of operation, minimal demand for data as input, and the ability to generate thematic end maps. For an updated and complete list of InVEST applications, see [38]. The InVEST Habitat Quality model considers habitat worth and rarity as a proxy for biodiversity and enables the estimation of the extent of a habitat type in a landscape and its state of degradation. The InVEST Habitat Quality model values and maps HQ and related degradation (D), representing the spatial distribution of the level of vulnerability of habitats [38].

The InVEST Habitat Quality model links LULC data with specific threats. Its implementation requires the identification of three elements: (i) key habitats; (ii) threats; and (iii) habitat sensitivity to each threat. HQ is estimated through the following equation [38]:

$$H{Q}_{xj}=H_j\left[1-\left({\displaystyle \frac{D_{xj}^z}{D_{xj}^z+k^z}}\right)\right]$$

(1)

where the following are defined:

-

j represents the key habitat type (with j = 1,…, 5);

-

HQ_xj represents the HQ index of the grid unit x and for the j habitat;

-

$H_j$, which ranges between 0 and 1, represents the habitat suitability score of the j key habitat, with 1 indicating the highest suitability to species;

-

$D_{xj}^z$ is the degree of habitat degradation of grid unit x for the j key habitat;

-

k^z is the semi-saturation constant.

The degree of habitat degradation (D_xj) is measured according to the following equation [38]:

$$D_{xj}={{\displaystyle \sum }}_{r=1}^R{{\displaystyle \sum }}_{y=1}^{Y_r}\left({\displaystyle \frac{w_r}{{{\displaystyle \sum }}_{r=1}^R{w}_r}}\right)r_y{i}_{rxy}{\beta }_x{S}_{jr}$$

(2)

where the following are defined:

-

r is the threat source;

-

R is the number of threat sources;

-

y indexes all grid cells on r’s raster map;

-

Y_r is the set of grid cells on r’s raster map;

-

w_r is the relative weight of each rth threat and indicates the relative destructiveness of a source of degradation relative to all habitats;

-

r_y is a dummy that equals 1 if in the yth grid cell the rth threat is present in the yth grid cell, and 0 otherwise (this information is uploaded in the model through threat raster maps representing the spatial distribution of each threat);

-

i_rxy represents the influence of threat rth on the yth grid in the xth grid habitat, assuming that the degree of threat decreases with the increasing of the distance between the grid and the source of the threat;

-

β_x represents the xth grid vulnerability to threats, which depends on the level of legal, institutional, social, and physical protections. It equals 1 if grid vulnerability is maximum because it is not protected, and 0 otherwise;

-

S_jr is the sensitivity of the key habitat j to the rth threat.

The influence of threat rth on the yth grid in the xth grid habitat (i_rxy) identifies the type of decay which can be calculated in different ways according to the hypothesis of whether the threat affects the habitat quickly or slowly. Two decay relationships can be assumed:

$$\mathrm{Linear}: i_{rxy}=1-\left({\displaystyle \frac{d_{xy}}{d_{r}max}}\right)$$

(3)

$$\mathrm{Exponential}: i_{rxy}=exp\left[-\left({\displaystyle \frac{2.99}{d_{r}max}}\right)d_{xy}\right]$$

(4)

where d_rmax is the max distance (expressed in km) between habitats and threats beyond which the threat does not generate negative effects, and d_xy is the linear distance between grid cells x and y.

In the linear decay function, the effect of the threat on habitat decreases steadily with increasing distance from the disturbance source. In contrast, the exponential decay function describes a process in which the effect of the threat decreases rapidly near the disturbance source and then stabilizes with increasing distance. The choice of a linear or exponential decay model depends on the nature of the threat and the type of habitat under consideration.

To ascertain potential threats and specify decay functions, we reviewed 158 eligible studies published in the period 1980–2023 and identified using the Preferred Reporting Items for Systematic Reviews and Meta-Analyses (PRISMA) framework [51].

Table 1. Parameter values for threats and specification of the decay function.

Threat	Description	Decay Function	drmax	Wr
Intensive agriculture	Arable land in irrigated and non-irrigated areas, vineyards, olive groves, citrus groves	linear	1.0	1.0
Fire events	Fire at anthropogenic risk	exponential	0.5	1.0
Main roads	Primary roads (highways and provincial roads)	linear	0.6	0.5
Urban	Urbanized areas (continuous and discontinuous urban fabric), industrial and commercial areas	exponential	2.0	1.0
Hunting	Areas where hunting is allowed	linear	1.0	0.8
Bare land	Areas which are in natural state	exponential	1.0	0.5
Hydrogeological risk	Areas vulnerable to flooding	exponential	0.5	0.8
Geomorphological risk	Areas vulnerable to landslides	exponential	0.5	0.8

lists the main threats identified in our review and reports the threat values assigned to W_r, d_rmax as well as our hypotheses about the related decay function in the supplementary material, we provided raster maps. Weights to each habitat were assigned on the basis of analysis of the empirical literature. Specifically, we utilized and adapted the expert information gathered by [11] for a study on HQ in the Italian protected areas.

Table 2. Sensitivity table.

ID	Key Habitat	Habitat Weight	Intensive Agricultural	Fire Events	Main Roads	Urban	Hunting	Bare Land	Hydrogeological Risk	Geomorphological Risk
1	Forest	1.0	0.2	0.8	0.3	0.4	0.8	0.8	0.2	0.2
2	Sparse vegetation	0.9	0.5	0.5	0.6	0.8	0.6	0.8	0.4	0.4
3	Wetland	0.9	0.5	0.2	0.4	0.2	0.8	0.3	0.3	0.3
4	Agroecosystem	0.5	0.2	0.2	0.2	0.8	0.8	0.7	0.8	0.8
5	No vegetation area	0.0	0.0	0.0	0.0	0.0	0.0	0.0	0.0	0.0

reports the value of habitat sensitivity assigned to each key habitat category by considering the specific threats (S_jr). For this purpose, we also relied on a comprehensive literature review and considered the judgments provided by local experts.

The InVEST Habitat Quality model’s estimation produces raster maps in which the HQ and D in each grid cell (100 m × 100 m) assume values ranging from 0 to 1. The HQ raster maps were processed using the natural break method in ArcGIS Pro 3.1.0 to identify four new HQ levels: (i) low; (ii) low–medium; (iii) medium–high; and (iv) high (see

Table 3. Reclassification of HQ.

Levels	Reclassification	Value
1	Low	0.00
2	Low–medium	0.01–0.50
3	Medium–high	0.51–0.90
4	High	0.91–1.00

). The same natural break method was used to allocate D original values into five levels: (i) high; (ii) medium–high; (iii) medium; (iv) medium–low; and (v) low (see

Table 4. Reclassification of vulnerability (D).

Levels	Reclassification	Value
1	Low	0.001–0.046
2	Low–medium	0.047–0.104
3	Medium	0.105–0.130
4	Medium–high	0.131–0.161
5	High	0.162–0.231

).

To identify changes in the HQ value at two time in points (2018 and 2050), we overlapped reclassified raster files and used the Raster Calculator tool in ArcGIS Pro 3.1.0 software to obtain a new raster containing HQ values (in 2018, class levels of HQ ranges from 1 to 4, while in 2050, they vary from 10 to 40). Figure 2 presents the procedure used to identify the changes in HQ levels, and

Table 5. Changes in HQ levels from 2018 to 2050.

Label	Level of HQ in 2018	Level of HQ in 2050	Degradation/Improvement in HQ
	0.51–0.90	0.00	Significant degradation
	0.91–1.00	0.00	Significant degradation
	0.01–0.50	0.00	Slight degradation
	0.51–0.90	0.01–0.50	Slight degradation
	0.91–1.00	0.01–0.50	Slight degradation
	0.91–1.00	0.51–0.90	Slight degradation
	0.00	0.00	Stable
	0.01–0.50	0.01–0.50	Stable
	0.51–0.90	0.51–0.90	Stable
	0.91–1.00	0.91–1.00	Stable
	0.00	0.01–0.50	Slight improvement
	0.00	0.51–0.90	Slight improvement
	0.01–0.50	0.51–0.90	Slight improvement
	0.51–0.90	0.91–1.00	Slight improvement
	0.00	0.91–1.00	Significant improvement
	0.01–0.50	0.91–1.00	Significant improvement

reports the size of these changes.

2.2.2. Validation of the HQ Estimates

To validate the reliability of the HQ estimates, mainly with respect to the ability to detect areas and sites requiring priority conservation and restoration actions, we used a coarse–filter approach [39], which is a mixed approach. The filter component examines the integrity of ecosystems and landscapes by focusing on the study of species and communities [40]. The coarse component considers the higher organizational levels as environmental units to preserve the different elements of biodiversity [41]. Each component produces a hotspot analysis. The base map used for the filter component was a raster of species included in the International Union for Conservation of Nature (IUCN) Red List and by BirdLife International. (The underlying species maps come from the IUCN Red List and BirdLife International. Integrated data from the IUCN Red List, World Database of Key Biodiversity Areas and World Database on Protected Areas are available in the Integrated Biodiversity Assessment Tool [52].) This raster maps the value of the biodiversity index (B) which displays the relative importance of each pixel in relation to its aggregate contribution to the distribution of forest-dependent species of mammals, birds, amphibians, and conifers. Data (in .tiff format) were extracted from the online platform Global Forest Watch [52]. The base map used for the coarse component was the LULC Copernicus 2018. Habitat (H) classification and weights assigned to each habitat type were the same as those adopted for HQ assessment. The overlaying of maps produced through the two hotspot analyses allowed the identification of PPAs. Each pixel of the reference rasters corresponds to one hectare of area. The original values range from −3 to +3. Therefore, values were reclassified using ArcGis Pro software 3.3.0 (Raster Calculator tool), assigning 0 to the coldspot areas (where original values range from −3 to −1) and 1 to the hotspot areas (where original values range from +1 to +3). Original values of zero were considered insignificant. In Figure 3 the new overlapped raster, each pixel is assigned a value of 0, 1, or 2. The value 0 indicates areas in which there are no significant hotspots for either raster (H or B). The value 1 indicates areas in which the hotspot area condition is present in only one of the two rasters (H or B). The value 2 indicates areas in which there is a significant high hotspot in both rasters (H and B).

By adding the layer of current protected areas, we calculated the spatial distribution of the following:

-

Protected areas;

-

Unprotected areas with value of 1 in H or B;

-

Unprotected areas with value of 2 in H and B;

-

Areas with degradation level as presented in Table 5.

2.2.3. Transfer Matrix

Dynamics in LULC were analyzed using the so-called transfer matrix (S_ij):

$$S_{ij}=\left[\begin{array}{ccc}{s}_{11}& \dots & s_{1n}\\ \dots & \dots & \dots \\ s_{n1}& \dots & s_{nn}\end{array}\right]$$

(5)

where S_ij represents the transfer areas between the initial (i) and final (j) land use types, accounting for the five (n) different land use types (key habitats) previously listed [17,28,53,54].

2.2.4. Spatial Statistical and Econometric Analysis

To identify clusters or the hotspot/coldspot areas, we relied on the Getis-Ord Gi* statistics:

$${G_i}^{*} = \left({{\displaystyle \sum }}_{j,j\ne i}^n{w}_{ij}{x}_j\right)/({{\displaystyle \sum }}_{j,j\ne i}^n{x}_j)$$

(6)

where n is the number of observations, x_j is the value observed in the region j, and w_ij is the element of spatial weights of a matrix with a diagonal with a non-zero element [55,56].

The Getis-Ord Gi* statistics evaluate whether a particular region has a high concentration of values that are higher (hotspots) or lower (coldspots) in relation to the average values in the surrounding area. If G_i* > 0, then areas are surrounded by relatively high value clusters; if G_i* < 0, then areas are surrounded by relatively low-value clusters. To be a statistically significant hotspot, a feature will have a high value and be surrounded by other features with high values. (For a statistically significant positive Gi* score, the larger the Gi* score, the more intense the clustering of high values (hotspots). For a statistically significant negative Gi* score, the smaller the Gi* score, the more intense the clustering of low values (coldspots).)

To test the existence and intensity of spatial autocorrelation of HQ with respect to spatial location, we estimated either global or local statistics [42,56,57]. Global autocorrelation, which allows the identification of clusters and spatial relationships only for the entire system, was measured by Global Moran’s I statistics:

$$I = {\displaystyle \frac{{{\displaystyle \sum }}_i{{\displaystyle \sum }}_j{w}_{ij}\left(x_i-\overline{x}\right)\left(x_j-\overline{x}\right)}{S^2 {{\displaystyle \sum }}_i{{\displaystyle \sum }}_j{w}_{ij}}}$$

(7)

where x_i is the observations in the region i,$\overline{x}$ is the average of all regions studied, n is the number of regions, and w_ij is part of the spatial matrix W, and

$$S^2={\displaystyle \frac{1}{n}}{{\displaystyle \sum }}_i{\left(x_i-\overline{x}\right)}^2$$

I > 0 means a positive spatial autocorrelation; that is, neighbors tend to have similar pattern; I < 0 indicates negative spatial autocorrelation; that is, neighboring points tend to have different values from each other; I ≈ 0 indicates no spatial autocorrelation; that is, the observed values are randomly distributed.

Local similarity, which refers to features with either high or low values that cluster spatially, was measured through the local Moran’s II statistics:

$$I{I}_i={\displaystyle \frac{\left(x_i-\overline{x}\right){{\displaystyle \sum }}_{i=1}^n{w}_{ij}\left(x_i-\overline{x}\right)}{{{\displaystyle \sum }}_{i=1}^n{\left(x_i-\overline{x}\right)}^2/n}}$$

(8)

where the w_ij elements come from a spatial weight matrix row standardized to one.

Local Moran’s II assumes significantly negative values when point i is surrounded by relatively low values in neighboring points, but it presents significantly positive values when point i is surrounded by relatively high values.

Local Moran’s II statistics and the Getis-Ord Gi* statistics form Local Indices of Spatial Autocorrelation (LISA), by which it is possible to identify spatial patterns of clustering (clusters) and isolated values (outliers) in a spatial data set [58]. The following four types of spatial patterns were identified by LISA:

• High–high-value cluster (HH): area with a high value that is surrounded by other areas with high values (hotspot);
• Low–low-value cluster (LL): area with a low value that is surrounded by other areas with low values (coldspot);
• High–low outlier (HL): area with a high value that is surrounded by areas with low values;
• Low–high outlier (LH): area with a low value that is surrounded by areas with high values.

Getis-Ord, global Moran I, and local Moran II statistics were calculated using ArcGIS Pro 3.1.0.

To investigate how other variables might influence HQ values in each reference year (2018 and 2050), we built and estimated two multivariate linear regression models [59]. First, we estimated, as a preliminary exploratory data analysis, a global linear OLS model, which assumes the stability of parameters ${\beta }_i$ throughout the data:

$$Y={\beta }_0+{\beta }_1 X_1+{\beta }_2 X_2+\dots +{\beta }_k X_k+ϵ$$

(9)

Then, we estimated a local or GWR model, which assumes that coefficients are spatially varying and depend on the presence of homogeneous clusters in space. A GWR model recognizes that the relationships between the dependent variable and the independent variables may differ in different geographic areas [60]. In the GWR, regression models are run locally at each point, weighting nearby observations more than distant ones [61]. The GWR model was specified as follows:

$$H{Q}_i={\beta }_{0 }\left(u_i,v_i\right)+{\beta }_{1 }\left(u_i,v_i\right) X_{1i}+{\beta }_{2 }\left(u_i,v_i\right) X_{2i}+\dots +{\beta }_{k }\left(u_i,v_i\right) X_{ki}+{\epsilon }_i$$

(10)

where $H{Q}_i$ is the HQ in point $i$, $X_{1i}$, $X_{2i}$,…, $X_{ki}$ are the value of the kth factors in point $i, {\beta }_{0 }\left(u_i,v_i\right)$, ${\beta }_{1 }\left(u_i,v_i\right)$,…, ${\beta }_{k }\left(u_i,v_i\right)$ are the value of local coefficients ${\beta }_i$ for given geographical coordinates $u_i,v_i \mathrm{of} \mathrm{point} i$, and ${\epsilon }_i$ is the error terms in point i.

The values of HQ in 2018 and 2050 were selected as dependent variables, and the nine influencing factors were chosen for the OLS regression model. Based on our OLS results, we built and estimated a restricted GWR model using only the following independent variables: population density (PD), compass exposure (CE), annual average temperature (T18), and patch richness (PR).

Regression models were estimated using ArcGIS Pro 3.1.0 and the GWmodel package in R 4.0.5 [62]. In the GWR model, the optimal bandwidth was calculated through the bw.gwr function of the R package GWmodel by assuming a Gaussian distribution for the kernel spatially weighting function (bandwidth is a crucial parameter in GWR; it controls how much neighboring points affect the calculation of coefficients at each location [63]). To facilitate the estimation process, we divided the entire area into a grid of 15,868 hexagons, having a total area of approximately 170 ha, using the Generate Fishnet tool in ArcGIS Pro 3.1.0. A hexagonal grid was preferred over a classic rectangular or square grid because its structure is more effective for detailed spatial analysis and with data distributions that require precision or continuity in space [64].

2.2.5. Data Sources

Data on LULC in 2018 and 2050 come from the Esri Land Cover 2050 project (for more details, see [65]). This project provides data for nine LULC typologies: (i) mostly crop; (ii) grassland-scrub; (iii) mostly deciduous forest; (iv) mostly needleleaf/evergreen; (v) sparse vegetation; (vi) bare area; (vii) swampy; (viii) artificial surface; and (ix) surface water. To simplify operations in the InVEST Habitat Quality model, the above nine typologies were grouped and reduced to the five key habitats, corresponding to the LULC Copernicus 2018 first-level classification [66,67], as listed in

Table 6. Surfaces of key habitats in 2018 and 2050.

Key Habitats		Surface in 2018 (ha)	Surface in 2050 (ha)	Change (ha)
Forest	Mostly deciduous forest	113,408	112,898	510
	Mostly needleleaf	20,866	20,640	226
Shrubs	Sparse vegetation	2660	2660	0
Wetland	Swampy	7676	7676	0
	Surface water	4852	4852	0
Agroecosystem	Mostly crop	1,887,410	1,889,570	−2160
	Grass	377,588	374,703	2885
No vegetation area	Artificial surface	97,775	110,709	−12,934
	bare area	52,013	40,540	11,473

.

Figure 4 presents the spatial distributions of five key habitats in 2018 and 2050.

Table 7. Data sources on threats.

Input Data	2018	2050	Data Type	Source
Land use/Land cover	The layer displays a single global land cover map with a pixel resolution of 300 m	The layer displays a single global land cover map with a pixel resolution of 300 m	Raster	ESRI—Clark labs and European space agency climate change initiative [65] (https://livingatlas.arcgis.com/landcover-2050/, accessed on 21 January 2025).
Intensive agricultural	Areas with intensive agriculture	Areas with intensive agriculture	Raster	ESRI—Clark labs and European space agency climate change initiative [65] (https://livingatlas.arcgis.com/landcover-2050/, accessed on 21 January 2025).
Burned area	Areas having low, medium, or high probability of fire occurrence	Areas having low, medium, or high probability of fire occurrence	Raster	Regione Sicilia [68] (https://data.europa.eu/doi/10.2873/70140, accessed on 21 January 2025)
Hydrogeological risk	Areas having high hydrogeological risk	Areas having low, medium, or high hydrogeological risk	Vectorial	Regione Sicilia [69] (https://www.sitr.regione.sicilia.it/pai/, accessed on 21 January 2025)
Hunting	Area in which hunting is permitted	Area in which hunting is permitted net of potential new protected areas	Raster	Regione Sicilia [68] (https://data.europa.eu/doi/10.2873/70140, accessed on 21 January 2025)
Roads	Main roads	Main roads	Vectorial	Regione Sicilia [68] (https://data.europa.eu/doi/10.2873/70140, accessed on 21 January 2025)
Geomorphological risk	Areas having high geomorphological risk	Areas having low, medium, or high geomorphological risk	Vectorial	Ministero dell’Ambiente [70] (http://www.pcn.minambiente.it/viewer/, accessed on 21 January 2025)
Urban area	Urban areas	Urban areas	Raster	ESRI—Clark labs and European space agency climate change initiative [65] (https://livingatlas.arcgis.com/landcover-2050/, accessed on 21 January 2025).)

reports the sources used to gather data on threats and

Table 8. Independent variables used in the regression models.

Variable	Description	Tool	Source
Nightlight (NTL)	Intensity of the socioeconomic activities and urbanization	Spatial analyst Tools—Zonal—Zonal Statistics as Table in ArcGIS 4.0.1	[71] (https://figshare.com/browse, accessed on 21 January 2025)
Population density (PD)	Population size per square kilometer of land area	Analysis toolbox—Statistics toolset—Summarize Within in ArcGIS 4.0.1	ISTAT [72] (https://demo.istat.it/, accessed on 21 January 2025)
Slope (SL)	Slope (gradient or steepness) from each cell of a raster	Slope (Spatial Analyst) in ArcGIS 4.0.1	Regione Sicilia [73] (https://www.sitr.regione.sicilia.it/geoportale, accessed on 21 January 2025)
Compass exposure (CE)	Values: 0° for north, 90° for east, 180° for south, and 270° for West	Aspect (Spatial Analyst) in ArcGIS 4.0.1	Regione Sicilia [73] (https://www.sitr.regione.sicilia.it/geoportale, accessed on 21 January 2025)
Temperature (T18)	Annual average temperature	Kriging tools in ArcGIS 4.0.1	Regione Sicilia [74] (https://www.regione.sicilia.it/istituzioni/regione/strutture-regionali/presidenza-regione/autorita-bacino-distretto-idrografico-sicilia/siti-tematici/risorse-idriche/report-siccit%C3%A0, accessed on 21 January 2025)
Precipitation (P18)	Average annual rainfall	Kriging tools in ArcGIS 4.0.1	Regione Sicilia [75] (https://www.regione.sicilia.it/istituzioni/regione/strutture-regionali/presidenza-regione/autorita-bacino-distretto-idrografico-sicilia/annali-idrologici, accessed on 21 January 2025)
H_Prime (H_P)	Shannon–Weaver index. It measures diversity within a biological community. The minimum value occurs when there is only one key habitat in the patch, and the maximum value when all key habitats are evenly distributed.	ArcGIS 4.0.1 (ATtILA)—CLC 2018	[76] (https://www.epa.gov/enviroatlas/attila-toolbox—ATtILA, accessed on 21 January 2025)
Level of protection (LP)	It indicates whether a portion of land is protected, taking a value of 0 when the percentage of hectares of protected areas is 0, 1 when the percentage of hectares of protected areas is between 0.1 and 50%, and 2 when the percentage of hectares of protected areas is >50%	Calculation of protected area surfaces within the hexagonal feature	The World Database on Protected Areas, 2020 [77] (https://www.protectedplanet.net/en, accessed on 21 January 2025)
Patch richness (PR)	It describes the variety of habitat types or ecosystems in each geographic area	Command selects by location in ArcGIS 4.0.1	[76] (https://www.epa.gov/enviroatlas/attila-toolbox—ATtILA, accessed on 21 January 2025)

describes the independent variables and sources used for the estimation of the OLS regression model and, partially, for the GWR models.

3. Results and Discussion

Table 9. Transfer matrix.

			LULC 2050
LULC 2018			A		B	C		D		E
			A1	A2	B1	C1	C2	D1	D2	E1	E2
	A	A1						510
		A2						226
	B	B1
	C	C1
		C2
	D	D1
		D2						1424		1461
	E	E1
		E2								11,473

reports values of elements in the transition matrix. Values identify the size of increase or decrease in each key habitat. In 2050, it is expected that crop habitat will have a significant increase of 2160 ha, of which 736 ha will be from forest and 1424 ha will be from grasslands. However, there will be an even more substantial increase in urbanized areas (artificial surface), which will grow by 12,934 ha, of which 1461 ha will be from grasslands and 11,473 will be from bare areas. These projections for LULC changes, which are indicative of increasing anthropogenic pressures, will lead to a reduction in available habitat for wildlife; greater risk of extinction of local plants and animal species; more fragmented landscape; and thus, an overall worsening deterioration of HQ [11,78].

Legend. Letters without numbers identify key habitats reported previously in Table 1.: A (Forest); B (Shrubs); C (Wetland); D (Agroecosystem); E (No vegetation area). Letters with numbers identify other habitats: A1: Mostly deciduous forest; A2: Mostly needle leaf; B1: Sparse vegetation; C1: Swampy; C2: Surface water; D1: Mostly crop; D2: Grass; E1: Artificial surface; B2: Bare area. Figure 5 exhibits the maps of HQ values. The maps do not reveal notable differences in the spatial distribution of HQ values across the two reference years. In general, HQ values correspond to the distribution of LULC habitat types, as already noted by [78]. Higher HQ values prevalently occur in protected areas, and forests that are localized in the northern/northeastern side of Sicily. Medium HQ values occur in crops (cereals) and grasslands which are mostly present in the Sicilian hinterland. As expected, areas of lower HQ value appear in built-up and uncultivated lands. In general, the average HQ value was 0.29 in 2018 (standard deviation = 0.42), while in 2050, it was 0.25 (standard deviation = 0.39). This reduction, which in relative terms is equivalent to 13.79%, is mainly attributable to a strong increase in the size of areas with a low–medium level of HQ (from 514 to 222,127 ha) and to a significant decrease in the size of areas that in 2018 had a medium–high or a high level of HQ (respectively, −26% and −32%) (see Table 9). It is observed that within protected areas, the HQ value was 0.66 in 2018, and this value is projected to decrease to 0.62 by 2050. In the four regional parks, the HQ value in 2018 was slightly higher, at 0.7. These values, which are very similar to those reported by [11] suggest that habitats within the protected areas are generally in good condition. The slight variations observed within the regional system of the protected areas can be attributed to differences in management levels (typically higher in the four regional parks) and the smaller size of the regional reserves and Natura 2000 sites. These smaller protected areas are also affected by reduced spatial and functional connectivity and greater exposure to anthropogenic pressures. The comparison between 2018 and 2050, which anticipates a general decline in HQ even within the regional protected area system, highlights the urgent need for mitigation and adaptation measures to safeguard biodiversity. Protection and restoration actions, along with more careful and sustainable land use planning and management policies—particularly for the preservation of agricultural soils—should also be implemented across the rest of Sicily, where HQ averages are significantly lower, at 0.19 in 2018 and 0.17 in 2050.

Figure 6 presents the average value of HQ for each key habitat type. HQ values in 2018 are systematically higher than the corresponding value in 2050.

Figure 7 presents the map of change in the HQ values at two time points. The map reveals that almost the entire region will retain a stable HQ value during the timeframe.

Table 10. Dimensions of changes in HQ.

Level of HQ	2018		2050		Changes in HQ Level
	Area (in ha)	%	Surface (in ha)	%	Area (in ha)	%
0.00	1759.000	68.39	1759.396	68.41	396	0.00
0.01–0.50	514	0.02	222,641	8.66	222,127	432.15
0.51–0.90	616,103	23.96	456,043	17.73	−16,006	−0.26
0.91–1.00	196,269	7.63	133,454	5.20	−62,815	−0.32
Mean of HQ	0.29		0.25		0.04
Standard Deviation of HQ	0.42		0.39		0.02

details increments/decrements in HQ levels. The results reveal that an area of 1,759,043 ha (95.5%) maintained a stable and low level of HQ, while areas of 21,349 ha and 481,781 ha recorded a medium–high and a high level of HQ, respectively. In addition, the model predicts that 83,218 ha have a high risk of degradation because their HQ level changes from high to medium–high. In contrast, 21,317 ha and 11,620 ha plausibly experience an improvement in the value of HQ, moving from a medium–low to a medium–high level, and from a medium–high to high level, respectively.

Figure 8 illustrates the spatial distribution of vulnerability in 2018 and 2050, while

Table 11. Dimensions of vulnerability class.

Vulnerability Class	2018		2050
	Area (ha)	%	Area (ha)	%
Low	2,265,660	93.72	2,184,699	90.37
Medium–Low	29,806	1.23	1,369,784	5.66
Medium	68,703	2.84	517,404	2.14
Medium–High	35,377	1.46	354,544	1.46
High	17,930	0.74	86,034	0.35

reports portions of regional land falling in each vulnerability class. In 2050, the prediction is that there will be a decrease in land in low-, medium-, and high-vulnerability classes, and an increase in land in the medium–low class vulnerability distributed through-out the region.

The estimated values of the global autocorrelation Moran I statistics indicate the existence of positive spatial autocorrelation in the HQ and denote areas in which higher or lower HQ tended to be clustered. Clusters come into contact randomly more often in 2050 than in 2018. In fact, the Moran’s I values increased from 0.28 in 2018 to 1.30 in 2050, demonstrating that the spatial agglomeration of the HQ improves in 2050.

Figure 9 presents the distribution of clusters and outliers identified through the estimation of the local Moran II statistics. Such distribution reveals several changes across the timeframe. Substantial variations occur particularly in HQ hotspots (areas with high–high levels) and HQ coldspots (areas with low–low levels). In 2018, the distribution of HQ hotspots intercepts a modest percentage (2.94%) compared with the same areas (23%) in 2050. A greater difference is recorded for the HQ coldspots, which in 2018, covers only 4% of the total surface, but in 2050, covers 45% of the total surface. In detail, in 2018, hotspots are primarily concentrated in the northernmost area of Etna Park, partially extending into the territory of Nebrodi Park. Additional hotspots are identified in the Palermo area and the westernmost part of the island near Marsala. Conversely, coldspots are predominantly located in the areas of Custonaci and Alcamo, despite the presence of several Natura 2000 sites (Capo San Vito, Monte Monaco, Zingaro, Faraglioni Scopello, and Monte Sparacio). Smaller coldspot areas are also observed near Termini Imerese in the Palermo province and in parts of the Caltanissetta province, including Special Protection Areas. (SPAs) such as Torre Manfria, Biviere, and Piana di Gela. In 2050, as illustrated in Figure 9, the distribution of clusters appears to be completely different from that in 2018. High–high clusters are primarily concentrated in areas currently occupied by Nebrodi Park, Etna Park, Sicani Mountains, and the Pantalica Reserve. In contrast, low–low clusters are more evenly distributed across Sicily, without any particular concentration in specific areas.

Figure 10 presents the maps of vulnerability clusters and outliers. In 2018, high–high clusters (hotspots) covered 14.3% of the regional area, while in 2050, they will cover 13.4%. Low–low clusters (coldspots) accounted for 45.2% of the area in 2018 and will reduce to 40% in 2050. In 2018, much of the area occupied by low–low clusters fell within major protected areas, except for some more inland and fragmented areas, characterized by high–high clusters. These hotspots are mainly concentrated in northeastern Sicily, in the provinces of Enna and Agrigento, and in some parts of southeastern Sicily. For 2050, the map of vulnerability clusters and outliers shows an opposite configuration to that of 2018. High–high clusters are now concentrated in most protected areas, whereas low–low clusters are mainly found in unprotected areas or, in some cases, in smaller and more fragmented protected areas. This shift underscores the importance of immediate action to address both current protected areas and those vulnerable to degradation, particularly in surrounding regions. Without appropriate management and mitigation strategies, current protected areas may be exposed to negative influences, leading to worsening environmental conditions and greater degradation in the future.

The maps in Figure 11 represent the spatial distribution pattern of hotspots/coldspots under the filter and coarse approaches. Both maps are only for 2018.

Table 12. Distribution of filter and coarse hot/cold spots across confidence levels and values of Getis-Ord Gi* (2018).

Confidence Level	Getis-Ord Gi*	Filter Approach		Coarse Approach
		Surface (ha)	% of Sicilian Surface	Surface (ha)	% of Sicilian Surface
99% Coldspot	−3	9256	35.99	1,534,608	59.67
95% Coldspot	−2	53,991	2.10	18,091	0.70
90% Coldspot	−1	28,249	1.10	9181	0.36
Not significant	0	300,023	11.67	90,325	3.51
90% hotspot	1	2901	1.13	8132	0.32
95% hotspot	2	55,986	2.18	15,302	0.59
99% hotspot	3	973,213	37.84	894,836	34.79

reports further details on the size of hotspot/coldspot areas by distinguishing them according to the value of Getis-Ord Gi* (Gi-Bi) statistics and significance level.

Figure 12 visualizes the map derived by joining the coarse–filter hotspot map with the current PAs network map. The joint map demonstrates that the spatial distribution of hotspots generally overlaps with PAs. The joint coarse–filter approach identifies zones with a high (equals 2) or a medium (equals 1) value of wildlife. In general, 31.57% of the Sicilian terrestrial surface has a high level of biodiversity either in habitats or in plant and animal species, while 34% of the same regional surface has a medium biodiversity value.

The map presented in Figure 13 is obtained by overlaying the map of current protected areas, the vulnerability map, and the map of coarse and filter hotspots. This overlaid map allows the identification of biodiversity hotspots that are currently without legal protection and are also subject to significant degradation. The red polygons correspond to sites with a higher level of vulnerability. The pink polygons correspond to hotspots with the maximum value (2) either for the coarse or the filter approach. The sites (see

Table 13. Distribution of land across degradation level and biodiversity value.

Degradation	Areas with Relevant Biodiversity Value				Areas Without Relevant Biodiversity Value
	Medium		High
	Surface (in ha)	%	Surface (in ha)	%	Surface (in ha)	%
High	8040	5.57	9410	6.52	209	0.14
Medium–High	13,992	9.70	19,285	13.37	806	0.56
Medium–Low	25,565	17.72	36,990	25.64	3118	2.16
Low	12,827	8.89	10,270	7.12	3777	2.62
Total	60,424	41.88	75,955	52.65	7910	5.48

) that intercept 6.52% of the regional terrestrial surface are mainly located near existing protected areas. These results are useful because they suggest new sites to be considered for enlarging current protected areas and places to establish new protected areas by prioritizing highly vulnerable high hotspots. Figure 13 also identifies five areas of special interest for conservation planning, labeled A, B, C, D, and E. Area A is located in the northeastern part of the island, within the province of Messina. It is characterized by the presence of several protected sites, including sites of community interest (SCIs) and special protection areas (SPAs), which encompass zones of high biodiversity value and regions of medium to high vulnerability. One potential action for this area could be the extension of the Nebrodi Park or the creation of a new regional park in the Peloritani Mountains. Area B, situated in the southeastern region of Sicily, includes SCIs, SPAs, nature reserves, and wetlands under the Ramsar Convention. The establishment of the Mt. Iblei National Park is currently under administrative review, following national legislation in 2007, but it has faced delays due to opposition from local political forces and other stakeholders. Our findings could be used to assist policymakers and relevant authorities in defining park boundaries and protection levels. Area C is in the province of Caltanissetta, along the southern coast of Sicily. Like the previous areas, it hosts numerous protected sites, including SCIs, SPAs, and nature reserves, along with districts of high biodiversity value and vulnerable zones. A potential strategy for this locality could involve expanding existing protected zones or creating new ones to encompass areas of greater biodiversity. Area D spans the inland regions of Sicily, particularly in the province of Caltanissetta. It features SCIs, nature reserves, and several high-vulnerability zones. In this case, conservation efforts should promote the creation of larger, interconnected regions that function effectively in maintaining ecosystems. A similar situation exists in Area E, located in the province of Agrigento, which faces comparable challenges in terms of the need for protection and the connection of protected locales.

Table 14. OLS estimates.

	OLS—2018		OLS—2050
	Estimate	t-Value	Estimate	t-Value
Intercept	0.22828	11.65922	−0.11174	−7.39385
NHT	0.000228	1.330921	0.000044	0.328888
PD	0.000398	1.279157	−0.00342	−14.2504
SL	0.00154	3.884854	−0.00024	−0.79765
CE	−0.000046	−0.79362	−0.00024	−5.38672
T18	0.001015	1.4595230	−0.00068	−1.26005
P18	0.000029	2.241114	0.000223	22.643
H_P	−0.00725	−0.82747	0.106683	15.77528
LP	−0.00361	−1.0592	0.20175	76.73853
PR	0.003641	0.896981	0.04277	13.64863
R2	0.001878		0.34516
Adjusted R2	0.001311		0.344788

presents the results of the exploratory OLS regression model. These results reveal striking differences across the two years. The OLS 2050 model outperforms the OLS 2018 model in unexplained variation, estimated by R2, and in the number of statistically significant regressors. In 2018, only the SL and P18 variables have a statically significant influence on HQ. Such effects were positive but very small in size. In 2050, all regressors are individually statistically significantly different from zero, aside from NHT, SL, and T18. The HP, LP, and PR variables work as we had expected; that is, they have a positive effect on HQ, as did T18, whereas PD and CE have a negative effect on HQ.

Table 15. GWR estimates.

	GWR—2018					GWR—2050
	Min	1st_Qu.	Median	3rd_Qu.	Max	Min	1st_Qu.	Median	3rd_Qu.	Max
Intercept	2.64 × 103	2.64 × 103	2.64 × 103	2.64 × 103	0.2639	0.19251396	0.19251396	0.19251396	0.19251396	0.1925
PD	4.24	4.24	4.24	4.24	0.0004	−0.00470832	−0.00470832	−0.00470832	−0.00470832	−0.0047
CE	−3.59 × 10−1	−3.59 × 10−1	−3.59 × 10−1	−3.59 × 10−1	0.0000	−0.00024866	−0.00024866	−0.00024866	−0.00024866	−0.0002
T18	1.20 × 10	1.20 × 10	1.20 × 10	1.20 × 10	0.0012	−0.00023412	−0.00023412	−0.00023412	−0.00023412	−0.0002
PR	3.74 × 10	3.74 × 10	3.74 × 10	3.74 × 10	0.0037	0.08029601	0.08029601	0.08029602	0.08029602	0.0803
Bandwidth: 418 km
R2	0.0004192298					0.04289
AIC	10662					8439

reports the distribution characteristics for each individual parameter of PD, AS, T18, and PR used in the specification of the GWR models. Again, there are striking differences between the models for the two study years. The model goodness-of-fit measures (R²) substantially improve from 2018 to 2050. The superiority of the GWR model in 2050 compared with the GWR in 2018 is further confirmed by the lower value for the Akaike information criterion. Overall, the 2050 GWR model can be considered more satisfactory in the explanation of local variations in the effects of the variables on HQ, number of variables having an effect consistent with expectations, and areas where there is evidence that such variables statistically significantly different from zero at the 10% level.

Table 16. Areas with positive and negative significant effects of GWR variable on HQ.

Variable	2018						2050
	Positive Sign (ha)	Negative Sign (ha)	Positive Sign (%)	Negative Sign (%)	10% Significant Level (ha)	10% Significant Level (%)	Positive Sign (ha)	Negative Sign (ha)	Positive Sign (%)	Negative Sign (%)	10% Significant Level (ha)	10% Significant Level (%)
PD	1,363,879	1,208,138	53	46	360,132	14	0	2,572,017	0	100	2,462,645	95.75
CE	1,129,161	1,442,856	43.9	56.09	254,648	9.9	1,399,564	1,172,453	54.41	45.58	1,284,909	49.96
T18	1,856,290	715,727	72.17	27.82	105,053	4.08	1,336,378	1,235,639	51.95	48.04	367,818	14.3
PR	1,519,845	1,052,172	59.09	40.9	341,181	13.26	2,063,659	508,358	80.23	19.76	2,488,392	96.75

presents, in absolute and relative terms, the sizes of areas in which the HQ is negatively or positively affected by each GWR variable, and the sizes of areas where such effects are statistically significant to at least the 10% level. The identification of these geographical areas relies on estimates of pseudo-t statistics. The analysis of habitat quality (HQ) influencing factors in 2018 and 2050 reveals distinct temporal trends across the selected variables. In 2018, population density (PD) showed a balanced influence, with 53% of the area experiencing positive impacts and 46% experiencing negative impacts. By 2050, the positive influence completely disappears (0%), and the negative influence dominates, affecting 100% of the area. This trend suggests that population growth exerts increasing pressure on habitat quality over time. Furthermore, the area with significant effects expands from 14% (360,132 ha) in 2018 to 95.75% (2.46 million ha) in 2050. Initially, Compass Exposure (CE) exhibited predominantly negative effects in 2018, with 56.09% of the area negatively impacted and 43.9% positively impacted. However, by 2050, a reversal occurs, with positive impacts increasing to 54.41% and negative impacts reducing to 45.58%. This shift indicates a spatially variable response where exposure characteristics increasingly support habitat quality. In 2018, the average annual temperature (T18) had a predominantly positive influence, affecting 72.17% of the area, while 27.82% experienced negative impacts. By 2050, this positive influence diminishes to 51.95%, with negative impacts rising to 48.04%. This trend highlights the growing challenge of temperature fluctuations, likely due to climate change, which increasingly undermines habitat quality. Patch richness (PT) shows a consistently positive trend over time. In 2018, 59.09% of the area experienced positive effects, which increased significantly to 80.23% by 2050. Negative impacts, in contrast, declined from 40.9% to 19.76%. Additionally, the area with significant effects increased from 13.26% (341,181 ha) to 96.75% (2.48 million ha), underscoring the critical role of patch richness in sustaining habitat quality.

Figure 14, Figure 15, Figure 16 and Figure 17 visualize the spatial, temporal, and statistical significance of GWR estimates.

4. Concluding Remarks

This study measured HQ in Sicily (Italy) using the InVEST Habitat Quality model. HQ was evaluated across the area of the island to address the influence of key habitats, and across time to predict to what extent HQ will be affected by changes in LULC. Estimates of HQ were heterogenous. As expected, protected areas and forests had higher HQ values. On average, HQ decreased from 0.29 in 2018 to 0.25 in 2050. This slight decline in HQ reflects projected variation in the LULC. In 2050, growth in crop production and urbanized areas is expected. These projections and the LULC transfer matrix reveal that the increases in these two types of areas will occur at the expense of forest and grassland areas. We validated and merged HQ and vulnerability estimates from the InVEST Habitat Quality model with output from a coarse–filter approach and identified hotspots with high conservation priority. Using spatial statistics, we also find the existence of positive spatial autocorrelation in HQ, demonstrating that areas in which higher or lower HQ tended to be clustered, and that clusters will come into contact more often in 2050 than in 2018. We estimated GWR models by verifying the sign and the significance effect of population density, compass exposure, average temperature, and patch richness. Such variables had different effects on HQ either in space and time or in significance level. Population density is projected to increasingly pressure habitats, while patch richness and compass exposure emerge as more critical positive drivers of habitat quality. Temperature effects become more balanced, reflecting the increasing complexity of climate impacts. Across all variables, the spatial extent of significant effects intensifies, signaling stronger localized influences in 2050.

This study has some limitations to address in future research. First, the selection criteria and data used for implementing the InVEST Habitat Quality model may have suffered because of the subjective process used to assign values to habitat suitability and weights to various threats, as well as of selecting the appropriate decay function for threats. Second, we did not evaluate HQ continuously over time but only at two points (2018 and 2050), and the temporal patterns of HQ in 2050 were based on the business-as-usual hypothesis. Alternative scenarios might be explored, such as the Shared Socioeconomic Pathways (SSPs) identified by the Sixth IPCC Assessment Report [79]. Third, additional key habitats, and alternative high-resolution data sources, such as remote sensing information prepared on cloud computing platforms like Google Earth Engine, could be used for quantifying land use and land cover, as well as other methods, such as Markov chains, should be investigated to project LULC change. Fourth, the way we used the coarse–filter approach may have been too restrictive, and other relevant biodiversity metrics may have been left out of the process. Fifth, the spatial regression models warrant more in-depth investigations to unveil spatial specific distribution effects on HQ. Future GWR analysis should include other variables and use more efficient models such as the beta GWR model, GWR model with heteroscedastic random error, and mixed GWR, in which some parameters are fixed, and others are determined locally.

Despite its limitations, this study offers valuable insights that can inform local policy- and decision makers in promoting transparency and facilitating communication with local stakeholders. These insights are crucial for designing and implementing effective conservation plans, expanding the current protected area system to meet the targets set by the Kunming–Montreal Global Biodiversity Framework [80] and the new EU 2030 Biodiversity Strategy [81], and developing adaptive, sustainable management strategies. Finally, this study provides a solid foundation for further research and exploration in this field.