Characterizing historical transformation trajectories of the forest landscape in Rome's metropolitan area (Italy) for effective planning of sustainability goals

With the aim at developing a landscape dynamics framework for environmental planning and management and testing the effectiveness of protected areas in achieving the 2030 Agenda of the United Nations sustainability goals, we characterized the historical transformation trajectories of forest area changes from 1936 to 2010 in the Metropolitan City of Rome Capital (Italy). Remote sensing‐based products coupled with landscape pattern metrics and fragmentation analysis have been implemented, comparing different historical forest maps. The results show a remarkable forest area gain – from 17.6% to 25.5% – thanks to 68,299 ha of recently established forest. Statistical descriptors showed that the highest relative gain occurred in mountain zones, confirming a wide European forest recovery pattern in marginal areas from past deforestation and overexploitation. Deforestation mainly occurred in the flat and hilly areas where almost 26,000 ha of forests were lost since 1936. In summary, two main forest landscape dynamics were reconstructed: (I) the increase of forest cover fragmentation in the lowland areas; and (II) the rise in the forest area in the interior sectors of the mountain landscape, mainly within protected areas. Restoring the forest ecosystem's bioecological integrity has been highlighted as an urgent action for biodiversity conservation and carbon mitigation. In lowland areas, the study revealed the urgent need to establish new protected areas and rewilding spaces as landscape metrics are relatively below the sustainability targets for healthy forest ecosystems. The proposed framework can be used for testing the effectiveness of environmental planning and management in other forest landscapes to achieve the Agenda 2030 goals.


| INTRODUCTION
Nowadays, few terrestrial ecosystems have remained undisturbed by anthropic activities (Foley, 2005;García-Vega & Newbold, 2020). Therefore, many forest landscapes worldwide reflect ecological and socio-economic history (Caetano-Andrade et al., 2020;Curtis et al., 2018). In Europe, forests have been used for millennia by humans, who have transformed the species composition, structure and spatial patterns of the community (Munteanu et al., 2015), and consequently their natural disturbance regimes and functionality compared to natural ecosystems (Sommerfeld et al., 2018). In the Mediterranean biodiversity hotspot, the current structure and physiognomy of the remaining forests have been altered by recurring degradation and deforestation activities (e.g., Mensing et al., 2018) and tracking the progress in halting habitats and species loss through ecosystem restoration and conservation programmes is an urgent action . In recent times, most forests have been maintained and managed for their wood and non-wood products and for many other functions, such as slope stabilization, to prevent hydrogeological instability (Führer, 2000). Since the beginning of the last century, there has been a natural recovery of forests in many temperate regions worldwide as a consequence of the abandonment of traditional mountain agriculture, driven by socio-economic factors such as immigration in urban areas (Geri et al., 2010;Gibon et al., 2010;Nadal-Romero et al., 2016;Romero-Calcerrada & Perry, 2004;Sitzia et al., 2010;Tasser et al., 2007). In the same period, deforestation activities occurred for farming and urban sprawl in lowland environments, for which characterizing historical forest transformation trajectories is a fundamental step in environmental landscape planning.
It is essential to investigate the landscape changes in forest ecosystems to meet the goals of the 2030 Agenda of the United Nations (www.un.org/sustainabledevelopment/development-agenda) and improve effective planning strategies to increase their resilience (Turner, 2010).
In the global framework of biodiversity conservation and carbon stores, protected areas are recognized as the cornerstones of the global conservation strategy and can reduce forest loss compared to unprotected areas (Wolf et al., 2021). Protected areas can maintain higher biodiversity levels and carbon stocks than neighboring alternative land use areas (Coetzee et al., 2014) and must be seen as an irreplaceable tool to guarantee complex ecosystem functions (Coad et al., 2019). As recognized by the key multilateral environmental agreement aimed at slowing biodiversity decline (UN Agenda 2030UN CBD, 2010), the expansion and effective management of protected areas is needed to mitigate the loss of biodiversity (Watson et al., 2014); thus, describing long-term regional changes in forest cover is critical for effective environmental management planning.
Understanding the dynamics of change and its impact on forest landscape functioning is crucial for biodiversity conservation (García-Vega & Newbold, 2020), and carbon mitigation strategies (Erb et al., 2018) and interesting solutions are being developed thanks to the integration of ecological history and remote sensing data. In this framework, landscape metrics and/or indices are often proposed and evaluated to assess landscape characteristics and to monitor changes in land use (Diaz-Varela et al., 2009;Geri et al., 2010;McGarigal, 2014;Modica et al., 2012;Uuemaa et al., 2009). The survival of threatened species also depends on the landscape dynamics and habitats' spatial configuration (Itani et al., 2020;Palmero-Iniesta et al., 2020;Rocha-Santos et al., 2020;Williams et al., 2020;Zhang, 2020), which change in quality, shape and position. Forest fragmentation can lead to the extinction of some species by isolating populations from each other and creating forest areas that are too small to be functional and maintain viable populations (Kettunen et al., 2007). Fragmentation can cause a change in ecological conditions with consequences on the abundance and distribution of species (May et al., 2019) due to the increase in forest edges and reduction in carbon storage capacity compared to a canopyclosed forest (Brinck et al., 2017).
This study presents a novel approach through the integrated use of historical forest cover data and remote sensing-based prod-  For each forest dynamic (gain, loss and persistence), the values of topographical variables such as elevation and slope were extracted and compared in terms of differences by using univariate statistics and applying the Kruskal-Wallis nonparametric test (Kruskal & Wallis, 1952). Statistical analyses were performed with R statistical software (R Core Team, 2020).

| Forest cover change detection analysis
To detect the changes that occurred in the time interval investigated (1936-2010), we performed a change detection (Singh, 1989) followed by a post-classification comparison approach (Lu et al., 2004;Modica et al., 2017)

| Forest canopy density reference map
To reconstruct the forest canopy dynamics and qualify them from an ecological functional point of view, we used advanced remote sensing-based products that users can suitably customize to distinguish and identify the forest and NF components (Chiarucci & Piovesan, 2020). As the reference layer for forest cover status in 2010, we used the GTCD with a spatial resolution of 30 m Â 30 m, (https://glad.umd.edu/dataset/global-2010-tree-cover-30-m, accessed on 27 January 2021). The GTCD dataset for the year 2010 (Hansen et al., 2013) is per-pixel estimates of percent maximum tree canopy cover, expressed as an integer value percentage (1%-100%).
The dataset was first clipped using the study area boundary. Then, we defined the corresponding forest area using a threshold of percentage tree cover as an area that was larger than 0.49 ha and with a tree cover of more than 10%, consistent with the forest definition of the FAO (FRA, 2018) and resampled at the spatial resolution of 20 m Â 20 m with class intervals every 10%. We performed an accuracy assessment to use the obtained map as a reference layer for forest canopy density. To produce sufficiently precise estimates of the classes' area, the sample size for each mapping class was chosen to ensure that the sample size was large enough (Global Forest Observations Initiative, 2013). Therefore, we calculated an adequate overall sample size for a stratified random sampling distributed among different strata (Cochran, 1977). First, we determined the number of sample units for the study area using the following Equation (1) (Olofsson et al., 2014) n ¼ Where: n is the total sample size, W i is the mapped proportion of the area of class i, and S i is the standard deviation of class i (forest/nonforest). S i was obtained according to the following formula (Equation (2)) Where: U i is the expected user accuracy of class i. For both selected classes (forest/non-forest) we used 0.09 and 0.91 as mapped proportions of the forest class areas for the NF area. We also set an U i of 0.9 and 0.8 for forest and non-forest classes, respectively. The overall sample size resulting from this calculation was 1528 sample points, and we applied stratified random sampling to allocate the samples to each stratum (altitudinal belt). As suggested by FAO (2016) and Congalton & Green (2019), we assigned a minimum size of 100 sample points to each altitudinal belt. We allocated the remaining number of samples proportionally, according to each stratum area (Table 1) Figure S1). Differences in canopy density values were tested applying the Kruskal-Wallis non-parametric test (Kruskal & Wallis, 1952) and the Mann-Whitney U test (Zar, 1996) as a post hoc pairwise comparison method after verifying the normality and homoscedasticity of the data with the Levene test (Levene, 1960) and Shapiro-Wilk test (Shapiro & Bradbury Wilk, 1965).

| Landscape metrics and forest fragmentation analysis
To quantify landscape changes closely linked to ecological processes, we analysed the forest landscape's structure and composition (McGarigal, 2014). To this end, metrics on size, shape and edge were used, selecting a few significant variables from the many available in the literature (Mcgarigal et al., 2002;Uuemaa et al., 2009  and mean fractal dimension (MFRACT). To calculate landscape metrics with the highest accuracy, referring to the input data, we used the vectorial data model and the plugin V-Late (vector-based landscape analysis tools extension) for ArcGIS (Lang & Tiede, 2003). Then, we computed the forest fragmentation using the raster-based forest fragmentation model of Riitters et al. (Riitters et al., 2000;Riitters et al., 2002). In this model, based on image convolution, the input raster must be binary (i.e., all cells are classified as forest = 1 or non forest = 0). (3) perforated (density > 0.6 and densityconnectivity > 0); (4) edge (Pf > 0.6 and densityconnectivity < 0); (5) transitional (0.4 < density < 0.6); and (6) (Fichera et al., 2015).
The edge category concerns those landscape areas between compact forest clusters neighbouring to and compact NF clusters. Considering the scale dependence of spatial metrics and that low data resolution could lead to inaccurate landscape pattern analysis (Wickham & Riitters, 2019), as input data, we used a raster with 10 m Â 10 m of geometrical resolution. Moreover, following Riitters et al. (2000Riitters et al. ( , 2002 and other scholars (Kowe et al., 2020;Li et al., 2011), the square moving window was fixed to 5 Â 5 pixels. into BL ones ( Figure 6).

| Forest canopy density estimation according to the GTCD in 2010
The estimation of the forest canopy density reference map (Table S3) obtained an overall accuracy level of 99.01% for the altitudinal belt I,  (Table S3). The PA for the 'non-forest' classification achieved an overall value of 98% for all three altitudinal belts. In contrast, the UA decreases with altitude but still shows the classification process's reliability.

| Forest canopy density assessment
Forest canopy densities for areas that gained trees show that growth dynamics have led to forests' formation with an average canopy density value of 46% (for the whole study area) ( Figure 7a). As the canopy density percentage's class increases, the forest extension decreases.
The most widespread density percentage class in the cover gain area is 10%. Moving to canopy density up to a threshold value of 60%, the forests' frequency gradually decreases to only 8%. It then increases in frequency again for classes >70% of density values, with the latter class covering 12% of the MCRC. A canopy cover density ratio of 100% is present only on a tiny part of the study area (0.2%; $ 14 ha).
In forest grain areas, BL forest confirms a U shape distribution with a minimum of around 70%-80% of canopy density cover (Figure 7a). In 2010, no coniferous forest reached 100% canopy density (Figure 7a).
The other MF category's frequency shows an increasing trend with increasing canopy density up to 80%, with a 133.33% increase in the distribution (from 21% to 49%) (Figure 7a).
Where persistence dynamics occurred (Figure 7b), the BL forests' frequency increases with canopy density. In contrast, the frequency of MF decreases. Moreover, in this case, coniferous forests do not reach a closed canopy structure (Figure 7b).
Considering the 2010 canopy density of forest cover gain areas within the altitudinal belts, low canopy densities (between 10% and 40%) are mainly distributed in the lowland belt (altitudinal zone I).
Higher canopy density values are mainly distributed between the hilly belt and the mountain belt (zones II and III) ( Figure S4).
The trend is confirmed even considering all MCRC (Figure 8), where the forests in the lowland belts, especially the coastal strip, have a lower canopy density than the forests of the mountain belt. Most deciduous and mixed deciduous forests are intensely coppiced, and these forest types occur in correspondence with sparsely, low-and medium-settled areas. These woodlands, therefore, maintain a lower density value than mountain forests that are managed following high forest systems. Moreover, a higher forest canopy density is a distinctive trait of protected areas (nature parks and Natura 2000 network areas) (Figure 8).

| Forest landscape dynamics
The number of patches shows the quantitative changes in forest cover both in the whole landscape and for different categories (Table 4; Table S5). In 1936, the forest area had minor fragmentation (Table 4) with respect to 2010 (Table S5). The general trend of the increase in the forest area is also linked to growth in the number of patches, as reflected in the MPS metric decrease, observed for all forest categories (Table S5). The MSI showed a general increase in size variability (Table 4) and individual categories ( passed from 3.15% to 6.94% and from 2.40% to 6.42% of the total surface, respectively, causing an increase in the landscape fragmentation metrics (Table S6) (Table 4), which significantly decreased in the four analysed landscape classes. On-the-other-hand, the form of patches in the 2010 landscape is more complex in all landscape classes, particularly for the BL forests (Table S5). with other studies showing forest cover increase in this area (Biasi et al., 2015;Salvati et al., 2017). The depopulation of mountain villages (Falcucci et al., 2007) and the abandonment of agricultural activities in marginal areas are the leading causes of the expansion of forest ecosystems in particular in mountain areas (Geri et al., 2010;Salvati & Sabbi, 2011), also leading to transitions in regimes and values (Gulinck et al., 2018). Conifers and BL forests showed a general increase, with the latter expanding mainly in the mountain belt, confirming the trend observed for most of the Apennines (Malandra et al., 2018).  Figure S8).

T A B L E 4 Summary of landscape metrics analysis for the MCRC
Although the mountain belt represents only 10% of the metropolitan area, forests cover 71% of the zone, with the highest cover gain and persistence rate. This result confirms the rewilding trend that has characterized the Apennine landscape in the past 60 years. One of In the last decades, forest restoration was principally due to a passive landscape rewilding due to the depopulation of mountain areas and ad hoc mountain management policy, mainly through the institution of protected areas. However, small patches of deforestation were detected up to the mountain belt, where new infrastructure and buildings were built primarily for tourist purposes of summer recreational activities and winter skiing ( Figure S8).
A different forest dynamic occurred in the lowland belt. The actual FCI of 14% is well below the 40% of optimal human-modified landscapes for forest biodiversity conservation (Arroyo-Rodríguez et al., 2020). Here, the primary historical process is linked to new forest expansion in marginal areas and deforestation. There are two leading causes for the few forest persistent areas: deforestation for agricultural land reclamation (up to 1950s) and urban expansion (Salvati, 2013). In the flat areas around Rome's city area (campagna Romana), the forest and agricultural lands were the most exposed to increased human pressure driven by urban sprawl (Salvati et al., 2014). In some coastal and flat inland areas, beyond some examples of relict forests characterized by structural and compositional heterogeneity (Pratesi, 2015), woodlands have always been subjected to intensive forestry (e.g., coppices) for providing goods and extensive grazing by domestic animals. Nevertheless, umbrella pine (Pinus pinea Canopy density monitoring and its continuity over time have proved to be a key action to measure forest's multifunctionality. Although the GTCD has a global reach, it has proved to be very reliable on a local scale. Visual interpretation of very high-resolution imagery, employed for data validation, confirms good results in deriving statistics on forest cover over large regions (Schepaschenko et al., 2019). As reported in the literature, photointerpretation can lead to an error that is less than 10% in estimating the extent of forest at a global scale (Bastin, Berrahmouni, et al., 2017;Schepaschenko et al., 2017). Using prior maps with satisfactory accuracy, integrating the accuracy of visual interpretation and performance of classification methods could be implemented in order to extend this approach to broader regions. In addition to the already known advantages of visual interpretation (Schepaschenko et al., 2019), further possibilities to improve forest observation will be given by the rising development of dedicated tools for visual interpretation providing very high-resolution (VHR) satellite imagery that can be viewed by users over many parts of the world (Lesiv et al., 2018).
Canopy density is connected to biomass stocks; high stocking levels ensure water regulation, conservation of biodiversity, carbon sinks and the mitigation of climate change effects. For this reason, maintaining and implementing a continuous forest cover through the application of close to nature forestry will guarantee a better net value of timber, carbon sequestration, production of secondary products, scenic beauty and large numbers of habitat trees, compared to forests managed with higher impacting silvicultural models such as the even-aged coppice system (Peura et al., 2018). Forest canopy density can thus be considered as an indicator of harvesting pressure.
Besides, it is possible to highlight the urgent need for a new environmental policy and management of forest landscapes in the lowland belt to increase their cover as well as the level of naturalness. There is a need for an in-depth revision of planning and management strategies targeted to resolve issues related to anthropogenic pressure (Cosentino et al., 2018)  Several studies confirm that a higher degree of forest cover contributes to mitigating heat waves within the ecosystems (von Arx et al., 2012(von Arx et al., , 2013. In contrast, the loss of cover leads to increased local heat, which aggravates the imbalance between the community's responses and climate change (Zellweger et al., 2020). The climatic processes' analyses are fundamental for understanding the relationship between land use, forest biodiversity and its functioning in this global change era. Moreover, a higher canopy density is connected to a greater mitigation capacity of surface runoff and reduced damage caused by floods during intense rainfall (EEA-European Environment Agency, 2015). Mature forests, ordinarily characterized by deeper soils, guarantee a better regulation of evapotranspiration phenomena, therefore positively affecting the ecosystem's hydrological balance.
Continuous canopy coverage is also linked to a more significant presence and abundance of endangered species such as saproxylic insects (Hardersen et al., 2020;Rossi de Gasperis et al., 2016). For a sustainable future, we need a general and greater availability of mature and old-growth forest habitats for conserving endangered species dependent on ancient trees and deadwood (Peura et al., 2018).
The restoration of old-growth forests remains a relevant goal in these regions because conserving and restoring forest ecosystems with a high naturality level have a priority ranking in environmental planning for a sustainable future (Chiarucci & Piovesan, 2020).
Protected areas can maintain higher biodiversity levels and carbon stocks than neighbouring alternative land use areas (Coetzee et al., 2014) and must be seen as an irreplaceable tool to guarantee complex ecosystem functions (Coad et al., 2019). As recognized by the key multilateral environmental agreement aimed at slowing biodiversity decline (UN Agenda 2030UN CBD, 2010), the expansion and effective management of protected areas is needed to mitigate biodiversity loss (Watson et al., 2014).
The European Union (EU) is preparing a post-2020 global transformation framework centred on biodiversity and forest strategies (Visconti et al., 2019)

DATA AVAILABILITY STATEMENT
Data sharing is not applicable to this article as no new data were created or analyzed in this study.