Investigating sacred natural sites and protected areas for forest area changes in Italy

Abstract

Forests will be critical to mitigate the effects of climate and global changes. Therefore, knowledge on the drivers of forest area changes are important. Although the drivers of deforestation are well known, drivers of afforestation are almost unexplored. Moreover, protected areas (PAs) effectively decrease deforestation, but other types of area‐based conservation measures exist. Among these, sacred natural sites (SNS) deliver positive conservation outcomes while making up an extensive “shadow network” of conservation. However, little is known on the capacity of SNS to regulate land‐use changes. Here, we explored the role of SNS and PAs as drivers of forest loss and forest gain in Italy between 1936 and 2018. We performed a descriptive analysis and modeled forest gain and forest loss by means of spatial binomial generalized linear models with residual autocovariates. The main drivers of forest area changes were geographical position and elevation, nonetheless SNS and PAs significantly decreased forest loss and increased forest gain. Although the negative relationship between SNS and forest loss is a desirable outcome, the positive relationship with forest gain is concerning because it could point to abandonment of cultural landscapes with consequent loss of open habitats. We suggest a legal recognition of SNS and an active ecological monitoring and planning to help maintain their positive role in biodiversity conservation. As a novel conservation planning approach, SNS can be used as stepping stones between PAs increasing connectivity and also to conserve small habitat patches threatened by human activities.

INTRODUCTION

Forests play a pivotal role in lessening the harmful impacts of climate change and, more generally, of global environmental change (Bonan, 2008; FAO, 2018). Besides representing the primary carbon stock on land and driving a large portion of carbon fluxes throughout the biosphere, they host a major portion of terrestrial biodiversity (Brockerhoff et al., 2017; Gibson et al., 2011; Lindenmayer, 2009). Moreover, they provide additional ecosystem services such as soil formation and protection, climate and flood regulation, slope stabilization, nutrient cycling, water and air purification and supply a wide variety of products (e.g., timber, fibers, food, and medicines). Forests have also been crucial to the livelihoods of local communities and maintained important cultural and spiritual values for long periods of time (Cooper et al., 2016; Stara et al., 2015). Hence, their status and dynamics are critical for human well-being and a sustainable future (Chiarucci & Piovesan, 2020; Ellison et al., 2017; FAO, 2018).

When observed over short timeframes, forest area and cover appear to be quite static, whereas major fluctuations, typically due to natural dynamics, can be observed at longer time scales, decades or more. The reduction of forest area and human disturbance in natural and semi-natural landscapes can lead to losses of biodiversity and erosion of ecosystem functioning and services. The global reduction of natural forests is a major challenge for the near future (Chiarucci & Piovesan, 2020) and several large-scale reforestation projects aim to counteract it worldwide (see e.g., Cao, 2008; Vadell et al., 2016; Wang et al., 2019; Yao et al., 2019). However, forest gain may also cause concerns in cultural landscapes, because afforestation can result in a reduction, rather than enhancement, of species adapted to open spaces. These phenomena mainly occur in mountainous and rural areas leading to the loss of semi-natural open ecosystems following the abandonment of agroforestal practices (Agnoletti et al., 2019; Amici et al., 2015). Afforestation by tree planting may be unsuccessful, especially when using species unadapted to local climate (Rackham, 2006), or may have negative impacts on regional water availability (Xiao & Xiao, 2019). On the other hand, secondary successions can lead to a reduction of forest fragmentation and to a general restoration of ecological functionality, but also to the formation of wooded patches with low or moderate natural values. At larger scales, rewilding landscapes can create sinks for CO₂, mitigating climate change and conserve biodiversity (Moomaw et al., 2019; Navarro & Pereira, 2015).

Monitoring changes in forest area, as well as identifying their spatial correlates, is therefore relevant for the sound management of biodiversity and ecosystem services. Although the drivers of deforestation and forest degradation have been extensively investigated, and are typically identified with commodity production, silvicultural practices, shifting agriculture, wildfires, and urbanization (Curtis et al., 2018; DeFries et al., 2010), few studies have explored the spatial correlates of afforestation. Clement et al. (2009) found that in Northern Vietnam, afforestation was positively correlated with proximity or presence of wood industry and distance from major roads, whereas it was negatively correlated with housing allocation. Upton et al. (2014) identified a number of physical, economic, and policy drivers for afforestation in Ireland, with proportions both of different soil types and private forests being the most relevant predictors. In the surroundings of Siena (Tuscany, Italy), forest gain was associated with higher elevations and steeper slopes, as a direct consequence of rural exodus and landscapes going wild (Geri et al., 2010). A similar pattern was recently reported for the province of Rome (Solano et al., 2021).

Protected areas (PAs) are the cornerstone of biological conservation and are usually expected to positively affect forest distribution and naturalness, thanks to conservation-oriented management. Andam et al. (2008) found that without legal protection 10% of forests within PAs in Costa Rica would have been cut between 1960 and 1997. Without the establishment of PA deforestation in China would have increased by up to 50% between 2000 and 2015 (Yang et al., 2019). PAs in Sumatra not only had lower deforestation rates than unprotected areas but also a similar pattern was observed even for adjacent areas, up to 10 km from PA borders (Gaveau et al., 2009). However, it must be noted that PAs are often located in remote or inaccessible areas where human impact tends to be lower, and when this pattern is not accounted for in the models, it can produce biased results (Joppa & Pfaff, 2009; Pfaff et al., 2015). PAs in the Amazon that are closer to roads and cities have a higher effectiveness on forest conservation (Pfaff et al., 2015), while Clemente et al. (2020) reported similar results for the Brazilian Cerrado. Besides their spatial location, other elements can influence the effect of PAs on deforestation and maintaining forest habitat integrity (Cropper et al., 2001; Herrera et al., 2019; Jones et al., 2018; Leberger et al., 2020; Miteva et al., 2019).

Besides PAs, additional non-formal area-based conservation measures have also been in place, sometimes much longer than PAs. Among these, sacred natural sites (SNS) are probably the best known. SNS are defined as “areas of land or water holding special spiritual significance for people and communities” (Wild & McLeod, 2008) and are often regarded as the oldest form of habitat protection in human history. They can be viewed as forming a shadow network of unofficial PAs which conserves biodiversity, provides ecosystem services, and potentially strengthens the “official” protected network (Avtzis et al., 2018; Dudley et al., 2009; Frascaroli et al., 2019; Zannini et al., 2021). SNS are often found in association with ancient forest remnants (see e.g., Cardelús et al., 2013; Shakeri et al., 2021) and have been shown to locally halt or reduce deforestation (Campbell, 2004, 2005). However, large-scale studies accounting for both the effects of PAs and SNS in controlling forest area changes are currently lacking.

In the present study, we assess whether PAs and SNS are important factors regulating forest loss and forest gain processes by means of descriptive analyses and modeling. We take Italy as case study, using nationwide data on both forest changes in a time frame of 82 years (1936–2018) and the spatial occurrence of SNS. While focusing on the role of SNS and PAs in respect to forest area changes, the present work also provides insights on other correlates of forest area changes at national level.

METHODS Study area

Italy is an elongated peninsula located in the middle of the Mediterranean Basin. It stretches from approximately 8.2°E to 18.5°E in longitude and from 35.5°N to 47.1°N in latitude. Its area is about 300,000 km², one fifth of which is made up of islands. It has a rugged topography, characterized by two main mountain chains, the Alps in the north, and the Apennines extending along the peninsula. The largest plain areas are in the North, the Po flood plain and the Venetian-Friulian flood plains, whereas other minor plains are scattered across the peninsula and the islands. The population is unevenly distributed, being concentrated in the urban areas along the coasts and in the northern plains. This pattern is the result of a rural exodus that started at the end of the XIX century and dramatically accelerated after World War II, resulting in the abandonment of large portions of interior areas of the Alps and the Apennines.

Forests and woodlands in Italy cover about 33% of land area (FAO 2020) and range from oro-boreal forests dominated by conifers to Mediterranean sclerophyll forests, passing through temperate broad-leaved forests (Mucina et al., 2016).

Data—forest area changes

To analyze forest gain and forest loss for the period 1936–2018, we produced a map of forest area changes (Figure 1a) following the approach proposed by Camarretta et al. (2018). We used the Italian Kingdom Forest Map 1936 by (IKFM) and the Corine Land Cover 2018 (CLC). We retained all the forest polygons without considering the forest type (see Camarretta et al., 2018 for a list of used land covers). Layers were then converted to raster format at 500 m resolution and overlaid, to produce a map of forest area changes. The value of each pixel was determined as the combination of two binary states, that is, forest or no-forest, one for each period, resulting in 4 possible outcomes: No forest (no-forest in 1936 and no-forest in 2018), Forest loss (forest in 1936 and no-forest in 2018), Forest gain (no-forest in 1936 and forest in 2018), and Forest persistence (forest in 1936 and forest in 2018). In addition, we added inland waters and wetlands to identify areas unsuitable for the subsequent analysis, where no forest area changes could be observed as they inherently could not host forests.

FIGURE 1. Forest cover changes in Italy for the period 1936–2018. “No forest” are areas that were not covered by forests in the 1936 nor in the 2018; “Forest persistence” are areas that were covered by forests in both periods; “Forest loss” are areas that were covered by forests in the 1936 but were not in 2018; “Forest gain” are areas that were not covered by forests in 1936 but were in 2018. a) Map of forest cover changes; b) percentage distribution of forest cover changes in the whole country (Italy) and in Protected Areas (PA) and Sacred Natural Sites (SNS) alone

Data—protection status

PAs in Italy mostly consist of two networks. One is the network of national parks and local reserves that are part of IUCN's Official List of Protected Areas (OLPA), with the first national park in the country having been established in 1922. The other network is the European conservation scheme Natura 2000, which was set up in the EU countries in the 1990s, following the implementation of both the Birds Directive and the Habitats Directive. These two networks overlap greatly, with Natura 2000 containing most of OLPA.

In addition, Italy is dotted by a number of SNS, characterized by a wide spectrum of environmental conditions and cultural uses (Frascaroli et al., 2016). Recently, Frascaroli et al. (2019) assembled a geo-referenced list of SNS across Italy and compared their spatial and landscape features with those of Italian PAs, highlighting the complementarity of the SNS and PA networks.

Protected areas

We considered as PAs all the areas belonging to the Natura 2000 network or OLPA or both, independently of the year of establishment. We aggregated the two data sets in a single layer, covering 62,739.5 km² of land area. We treated all the PAs as a unique category of protection status, alternative to SNS. Even though the PA network includes different categories, we did not take them into account as we wanted to focus on the differences between unofficial and official area-based conservation measures, rather than discriminate between different protection regimes within PAs.

Sacred natural sites

We used the data set assembled by Frascaroli et al. (2019), consisting of 2332 SNS throughout the whole Italian territory. These SNS are mainly shrines, hermitages, abbeys, monasteries and cloisters, and they range in size from small chapels of few m² to large monastic estates spanning over many hectares. In addition, management in these SNS varies largely, some SNS are abandoned or seldom managed, mainly during annual ceremonies and related pilgrimages; at the other end, there are permanently inhabited and cultivated SNS. Moreover, all the 2332 SNS are Catholic, even though some dates back to pre-Christian ages (whereas the most recent ones were funded in the 20th century). Although Frascaroli et al. (2019) noted that the distribution of SNS in the data set is biased toward certain geographical areas (i.e., Central and NW Italy), we assumed the data set to be sufficiently robust for our analysis by virtue of its large size. Following Frascaroli et al. (2019), we considered all the pixels in a circular buffer of 1250 m radius around each SNS as its area of influence, resulting in 10,413 km² of land area. The 1250 m radius is an arbitrary size that in many cases does not properly reflect the real area of influence, as it can be smaller or bigger. However, it enabled us to model forest area changes, without retrieving site-specific information.

Data—other explanatory variables

In addition to Protection Status, we extracted a number of potentially relevant variables to model forest gain and forest loss. A first set of variables consisted in geographical and topographical variables, that is, Longitude, Latitude, Elevation, Slope, and Terrain Ruggedness Index (TRI; Riley et al., 1999). Moreover, we considered anthropogenic factors other than Protection Status, that is, Distance to Nearest Road and Population Density (number of people per square kilometer). To compare model coefficients, we standardized explanatory variables as Z-scores by subtracting from each observation the variable's mean and dividing by the variable's standard deviation.

Data analysis—descriptive statistics

To first explore the frequencies of different forest area changes (and lack of change) across Italy and the different protection statuses, we calculated a contingency table with protection status versus the four categories of our forest area change map.

Data analysis—modeling

To investigate forest gain and loss, we modeled the two processes separately, by means of binomial generalized linear models (GLMs). First, we removed current inland waters as we presumed they were mostly already present in 1936. We also removed pixels at elevations above 1850 m a.s.l. as we assumed it as a reasonable treeline in the region. We converted the explanatory variables to raster format, aligning the resulting layers with the map of forest area changes. To avoid multicollinearity, we checked correlation among all pairs of explanatory variables. We found that Elevation, Slope, and TRI were all highly correlated (Pearson's ρ > 0.7). We, thus, retained only Elevation, assuming that it adds specific information on the geographical setting, compared to the other two. Then, we split our data set into two subsets: areas that were covered by forest in 1936 (919,856 observations) and areas that were not (236,345 observations). In this way, observations for each data set could have only two outcomes: no change and change. Notably, observations from the first data set (forested areas in 1936) could only be Forest persistence (i.e., no change) and Forest loss (i.e., deforestation). Conversely, observations from the second data set (unforested areas in 1936) could only be No forest (i.e., no change) and Forest gain (i.e., afforestation). Then, we ran some preliminary binomial GLMs to explore the validity of our approach, but binned plots showed that binary response variables were unsuited for our models. Hence, we transformed the response variables by applying a moving window of 2.5 km, thus calculating the percentage of Forest gain and Forest loss within the moving window. To reduce spatial auto-correlation (SAC) of the residuals, we added Longitude and Latitude quadratic and cubic terms to the models (i.e., trend surface analysis with polynomial regression; Fletcher & Fortin, 2018). In addition, as response variables showed an unimodal trend with respect to Elevation, we also used its quadratic term. Therefore, we checked for SAC in the residuals of the models by calculating Moran's I. As we still found the residuals to be significantly spatially autocorrelated, we tried to remove SAC by adding a residuals auto-covariate term (RAC) to each model. Finally, we calculated the explained deviance (D²) of both the final models and the models without RAC, we checked for SAC in the residuals of the final models by calculating Moran's I and calculated variable importance. All analyses were performed with R 3.6.3 (R Core Team, 2020). Data import and preparation were performed with sf (Pebesma, 2018), raster (Hijmans, 2020), magrittr (Bache & Wickham, 2014), dplyr (Wickham et al., 2020) and rgrass7 (Bivand, 2019) packages. D² and Moran's I were calculated by means of modEvA (Barbosa et al., 2013) and spdep (Bivand & Wong, 2018) packages, respectively. Variable importance was calculated with caret (Kuhn, 2021) package. We used ggplot2 (Wickham, 2016), ggspatial (Dunnington, Dunnington & Thorne, 2020) and RStoolbox (Leutner et al., 2019) packages to produce graphical outputs.

RESULTS Descriptive statistics

Between 1936 and 2018, Italy experienced forest gain in about 14.6% of its territory, with even larger values in PAs and SNS (21.7% and 18.8% of their areas, respectively; Figure 1b). Forest loss occurred in about 4.6% of the whole territory, with moderately larger values in PAs and SNS (6% and 5.2%, respectively; Figure 1b).

Modeling—forest gain

The model explained most of the observed variability (D² = 0.89). However, this was largely due to the addition of the RAC (D² of the model without RAC = 0.38). After adding the RAC, SAC decreased significantly (without RAC: Moran's I = 0.931, CI = [0.93, 0.932]; with RAC: Moran's I = 0.663, CI = [0.661, 0.664]). However, the residuals of the model still showed highly significant SAC (p-value <2e-16). All coefficients, except for the quadratic term of Longitude, were significant (Figure 2). Moreover, the size of the coefficients, as well as variable importance, varied greatly, with Longitude, Latitude, and Elevation terms being much greater than the others. Moreover, RAC, Elevation, and Northing were also the most important variables (Figure 3a). In particular, forest gain was strongly associated with Northern and Western locations and with mid-altitude elevations. Protection Status was associated with increased forest gain, with PAs having a slightly greater effect than SNS. Forest gain increased further from roads and in sparsely populated areas (Figure 4).

FIGURE 2. Standardized coefficients estimates for residual auto-covariate binomial generalized linear models (GLMs) modeling Forest gain and Forest loss in Italy for the period 1936–2018. SNS means sacred natural sites, and PAs means protected areas. Solid dots represent statistically significant coefficients (p [less than] .05), and empty dots represent statistically nonsignificant coefficients, whiskers r 95% confidence intervals. X-axis is log10 transformed. Residual auto-covariate estimate is not shown

FIGURE 3. Variable importance plots of (a) auto-covariate binomial generalized linear model (GLM) of Forest gain; (b) auto-covariate binomial GLM of Forest loss

FIGURE 4. Marginal effect plots of Forest gain residual auto-covariate binomial generalized linear model. SNS means sacred natural sites, and PAs means protected areas. Residual auto-covariate estimate is not shown. Continuous covariates were back-transformed from normalization before plotting. Gray areas and error bars represent 95% confidence intervals

Modeling—forest loss

The model explained a large amount of observed variability (D² = 0.78) but, similarly to the other model, this was mostly due to the addition of the RAC (D² of the model without RAC = 0.22). Still, SAC decreased significantly after adding the RAC (without RAC: Moran's I = 0.891, CI = [0.888, 0.894]; with RAC: Moran's I = 0.602, CI = [0.599, 0.604]), but we were not able to fully remove SAC from the residuals (p-value <2e-16). Population Density was the only non-significant predictor (Figure 2). Again, Longitude, Latitude, and Elevation terms showed the largest absolute values, but with different trends than for forest gain (Figure 5). In addition, Longitude, Latitude, and Elevation were the most important variables along with RAC (Figure 3b). Forest loss was more strongly associated with Southern and Western locations and low- and high-altitude elevations. Forest loss increased in densely populated areas and close to roads. Protection Status was associated with a decrease in forest loss, with PAs having a slightly larger effect than SNS.

FIGURE 5. Marginal effect plots of Forest loss residual auto-covariate binomial generalized linear model. SNS means sacred natural sites, and PAs means protected areas. Residual auto-covariate estimate is not shown. Continuous covariates were back-transformed from normalization before plotting. Gray areas and error bars represent 95% confidence intervals

DISCUSSION

Our results, based on calculating forest gain and loss at the national scale in the time window of about a century, indicated that the main trend was toward a net increase of forest area throughout Italy (Antrop, 2004; Camarretta et al., 2018; Falcucci et al., 2007). Rural abandonment since the 1960s and the subsequent afforestation of lands formerly dedicated to grazing and small-scale agriculture are commonly indicated as the main drivers of this tendency (Camarretta et al., 2018). This is quite clearly reflected in the land use dynamics of SNS, which in Italy are largely associated with traditional rural landscapes (Frascaroli et al., 2019) and had a relatively higher increase of forest area than Italy as a whole (18.8% vs. 15.5%) in the period considered. At the same time, both SNS and especially PAs had relatively higher rates of gross forest loss.

However, our model of forest loss contradicts this last observation. Indeed, after accounting for other geographical and anthropogenic factors through multiple regression, SNS and PAs in Italy underwent significantly less forest loss than the whole territory between 1936 and 2018. Although this pattern was already known for PAs (see e.g., Andam et al., 2008; Gaveau et al., 2009; Pfaff et al., 2015), this has not been tested for SNS, but see Cardelús et al., 2013. Despite being significant, the association of SNS and PAs with reduced deforestation is small compared to the geographical variables and the RAC. One of the reasons could be that a large portion of the deforestation pattern of SNS and PAs is masked by their geographical location, and especially Elevation. Notably, SNS in Italy are mostly found at mid elevation and PAs at mid-high elevations (Frascaroli et al., 2019), which in turn are the elevation ranges where forest loss was at minimum. In certain cases, the loss of active governance due to the abandonment of hilly and mountain areas may have led to a marginalization of SNS and their role in cultural landscapes. As for the weak association between PAs and forest loss, it could be due to the temporal scale of this study, as the vast majority of PAs in Italy were established after 1936 and their effects on forest conservation may be delayed.

Our results also support the hypothesis that some form of land protection promotes forest gain. However, similarly to what we observed for forest loss, the effect of Protection Status was significant, but relatively small. This again could partly be a consequence of Elevation explaining portions of variability, which could have been explained by Protection Status, as SNS and PAs are mostly found at elevations where forest gain is at its maximum. The positive association between Protection Status and forest gain is not necessarily a desirable phenomenon, especially where cultural landscapes have great relevance for biodiversity and ecosystem services (Blondel, 2006). Indeed, the loss of open habitats such as traditional pastures or semi-natural grasslands can lead to landscape transformation toward closed forests. Moreover, particular attention to forest gain should be posed in areas where open habitats are rare or particularly relevant for biodiversity such as mountain areas or islands, as well as agricultural landscapes (Deák et al., 2020; Tsiakiris et al., 2009).

Besides Protection Status, the models explained most of the observed variability after the addition of the RAC. The great difference in explained variability between the RAC and the other covariates is most likely due to the large number of zeros in the response variables, which, in turn, reflects the intrinsically clustered arrangement of forests and anthropogenic land uses. This clustered arrangement probably increased the importance of geographical variables as well (i.e., Longitude, Latitude, and Elevation), which were the most important predictors after the RAC.

Intriguingly, the effects of the anthropogenic variables, Protection Status, Distance from Nearest Road, and Population Density, were rather weak and explained only a small portion of variability. Nevertheless, the relationships were almost always significant and their directions unsurprising. Indeed, unprotected, less remote, and densely populated areas suffered more forest loss, whereas protected, remote and sparsely populated areas were associated with forest gain. However, the patterns observed for Protection Status could raise some concern. On one hand, active protection of these areas may have prevented deforestation. On the other hand, the increased forest gain seems to hint at the fact that the reason may have been land abandonment. In mountain sites, accidental rewilding (i.e., rewilding as a consequence of land abandonment) is leading to landscapes dominated by forests reducing the typical patchiness of cultural landscapes. In addition, forest gain is likely related to the loss of secondary grasslands and pastures, which may be relevant for the conservation of priority EU open habitats (Natura 2000) in this biogeographical context. However, afforestation likely reduces forest fragmentation, potentially restoring functional forest ecosystems. Moreover, rewilded forest landscapes are more resistant to climate change, thus allowing the conservation of threatened species and ecosystem services such as C0₂ removal, soil protection and slope stabilization (Moomaw et al., 2019; Navarro & Pereira, 2015).

Three main areas of possible policy interventions emerge from our findings. The first is the official acknowledgment of the conservation role of SNS. This role has been already documented (see Dudley et al., 2010 for a review), especially with regard to the biological composition of local areas. Our study demonstrates the importance of SNS in maintaining forest area at a national scale and over a broad time frame. As evidence grows on the contribution of SNS to nature conservation at different spatial and temporal scales, it also becomes increasingly important that this contribution is adequately recognized and ratified. Inclusion in PAs can be an easy way to grant such ratification, although it might not be the most effective approach (Mallarach & Verschuuren, 2019). Fortunately, new policy instruments are gaining traction, such as the definition of “Other Effective Area-based Conservation Measures” (OECMs, see e.g., Dudley et al., 2018), which could be more suited for SNS (IUCN WCPA, 2019). Moving toward a more formal recognition of the conservation role of SNS at regional, national, or EU scales, whether by adopting existing frameworks or elaborating new ones, should be prioritized (Frascaroli et al., 2019).

The second indication is the need to investigate the role of management at SNS for biodiversity conservation. Indeed, we identified a clear trend toward the afforestation of SNS. As in case of other cultural landscapes, this can lead to the loss of rare open habitat patches typical of semi-natural ecosystems. Notably, SNS can be associated with a variety of livelihood practices (in Italy, for example, pastoralism and agroforestry), that have contributed to shaping and maintaining heterogeneous multifunctional landscapes where domestic and natural biodiversity coexist (land sharing landscapes). Many of these agrarian practices are rapidly disappearing due to societal changes and the economic marginality of such activities. The present study has delineated landscape transformation areas providing a mapped knowledge base for monitoring and evaluating trade-offs in ecosystem services between cultural and natural landscapes in a changing society. Further efforts should ideally rest on interdisciplinary research linking historical ecology, traditional ecological knowledge (Molnár & Babai, 2021) and conservation biology aiming to document and understand past and current management practices at SNS and their role in nature conservation. However, it should be emphasized that human populations have been the main actors in historical landscape changes (Mensing et al., 2020). Therefore, the current mountain landscape going wild is an indicator of a transforming society in search of a more balanced relationship that combines natural resource use and biodiversity conservation, balancing land sharing and sparing at different spatial scales (Ekroos et al., 2016).

Finally, our study offers broader insights that go beyond PAs and SNS alone. Indeed, we found that the effect of SNS and even PAs on reducing deforestation has ultimately been marginal, the registered effect is more due to geographical and topographical factors. In other words, isolation and difficult access have been more important forms of protection than official or cultural human norms and regulations (Cervellini et al., 2017). For these reasons, the construction of new forest roads that could favor the exploitation of strategic areas for forest biodiversity conservation should be resisted (Ibisch et al., 2016), by establishing roadless areas (Kati et al., 2022). Future efforts to broaden and integrate conservation and restoration of ecological networks should primarily focus on more easily accessible areas, as these historically have been more at risk. In addition, this would help close a well-known gap in the spatial cover of PAs (Joppa & Pfaff, 2009), which is highly skewed toward isolated high-elevation areas at the expense of lowland and densely populated ones.

ARTICLE IMPACT STATEMENT

Sacred natural sites and protected areas play an important role in forest conservation and afforestation processes throughout Italy.

ACKNOWLEDGMENTS

The research was supported by the Ministry of Education, University and Research, Italy (MIUR) through PRIN project 2015P8524C "Biodiversity and ecosystem services in sacred natural sites (BIOESSaNS)." Piero Zannini has been supported by LifeWatch Italy through the project LifeWatchPLUS (CIR-01_00028).

CONFLICT OF INTEREST

None.

AUTHORS' CONTRIBUTIONS

PZ, FF, JMH, and AC conceived the idea and developed the methodology. PZ performed the analysis. PZ lead the writing. All the authors commented and edited previous versions of the manuscript.

DATA AVAILABILITY STATEMENT

Data and scripts are available at 10.5281/zenodo.6412324.