1 Introduction

Despite the widespread presence of human pressures, cities can host a great diversity of habitats, plant communities as well as rare and threatened species (Schwartz et al. 2006). This richness can be so high that sometimes exceeds that of the surrounding rural areas (Wania et al. 2006; McKinney 2008; Faeth et al. 2011). The reason for this richness is relatively poorly understood but it seems partly related to the heterogeneity of the urban environment (Sukopp 2004; Kowarik 2011) which with a wide range of urban green areas types, such as gardens, backyards, parks and natural and seminatural vegetation patches, allows the coexistence of a great variety of species and plant communities. Urban parks usually are the most heterogeneous and most diverse areas within the cities (Hadidian et al. 1997; Nielsen et al. 2014); their conservation is increasingly important in urban planning (Secretariat of the Convention on Biological Diversity 2012). To enhance biodiversity, ecosystem's integrity and functionality and at the same time to grant the ecosystem services they provide, appropriate management strategies are required (Gaston et al. 2005; Shwartz et al. 2013). Managers of protected areas are constantly faced with the need to ensure both the protection of nature and the public fruition (Rossi et al. 2009). In urban context, where natural areas are mostly used for recreational activities and human pressure can be considerable, the protection of biodiversity could be an even more difficult challenge (Minuto et al. 2006; Hamberg et al. 2008). The heavy frequentation for recreational purposes is a known source of threat for wildlife conservation (Cole 1995) and management choices not backed by scientific basis, such as monitoring or assessments of the potential impact of visitors, may have dramatic consequences on ecosystems (Minuto et al. 2006; Shwartz et al. 2014). Scientific and ecological knowledge is therefore an essential prerequisite for urban planners’ decisions.

The necessity of a knowledge-based management practice is particularly evident in the case of the network of protected areas established in 1998 (Cazzola 2005) in the city of Rome. This system includes 14 reserves with different types of ecosystems, degree of human impact and management requirements. Some of them are embedded within the urban matrix and represent remnants of seminatural areas that maintain the original features of the so called “Campagna Romana”, the typical hilly landscape of extensively managed grasslands intermixed with forest patches, that surrounds the city (Celesti-Grapow and Fanelli 1993). The last remaining fragments of seminatural vegetation within the city greatly contributes to the high level of Rome diversity (Celesti-Grapow et al. 2013) and assessing the effectiveness of the protected areas through specific monitoring programmes is paramount to guarantee their protection.

Thanks to the availability of floristic data collected before the establishment of the park (Bianco 1994) and by repeating the survey after about 30 years, we analyzed the effect of the change in management after protection on the flora of the “Pineto Urban Regional Park”, one of the most important areas of the network. In turn, results were used to suggest guidelines aimed at the conservation of the high floristic diversity (677 species) and historical values of the park.

2 Methods

2.1 Study area

Located in the north-western district of the city, in a highly urbanized area at about 2 km from the Vatican City, the “Pineto Urban regional Park” (Fig. 1) is one of the richest floristic areas in Rome (Celesti-Grapow 1995; Fratarcangeli et al. 2019) and for this reason it was widely investigated by botanists in the past (Montelucci 1953; De Lillis et al. 1986; Bianco 1994; Bianco et al. 2003).

Fig. 1
figure 1

a Location of the “Pineto Urban Regional Park” within the city of Rome. b Study area with the main geomorphological features and the location of the 28 one-hectare plots surveyed during the floristic investigation

Full size image

The area represents one of the last fragments of the typical “Campagna Romana” landscape escaped from urbanization and until about 25 years ago it was grazed by sheep. The presence of sheep dogs discouraged the frequentation of the park and grazing hindered the growth of woody vegetation. After the establishment of the park in 1987, the area has experienced an abrupt change in management with the abandonment of grazing activities and the increase of frequentation for recreational purposes.

Characterized by a typical Mediterranean climate (Blasi 1994), due to the geomorphological heterogeneity (Funiciello and Giordano 2008) (Fig. 1) the park presents a diversified vegetation (Bianco 1994; Fanelli 2002; Bianco et al. 2003). Pliocene marine clays at the bottom of the valleys, impeding water infiltration, allow the formation of ponds and wetlands. In top of the hills, a mixture of Cistus salviifolius garrigues, Tuberaria guttata grasslands and cork oak (Quercus suber) woodlands can be found. On the west of the park dry Mediterranean meadows characterize the tuffaceous plateau originated as a result of volcanic activity.

2.2 Data

We compared the floristic censuses conducted in 1991 (Bianco 1994) and in 2018 in the same 28 one-hectare permanent plots (Fig. 1) randomly distributed within the park and selected during the first campaign. Taxonomy follows Bartolucci et al. (2018) and Galasso et al. (2018) for native and alien species, respectively.

As a predictor of floristic change, we considered the minimum distance between plots and the park boundary. Distances were measured both along the paths during the field campaigns and as the crow flies on the QGIS software (QGIS Development Team 2020). The former measure broadly reflects the accessibility of each site by the citizens; the latter measure accounts for other kinds of disturbances associated with the urban landscape such as pollution or an increase in the urban heat island.

2.3 Statistical analysis

All statistical analyses were run with the software R (R Core Team 2019). Basic statistics on species richness and presence of aliens were computed for each sampling year and differences between the two surveys were tested for significance using the Wilcoxon test (Wilcoxon 1945). To analyze the temporal changes of the community structure, Rank-Abundance diagrams (MacArthur 1957) were built with “BiodiversityR” package (Kindt and Coe 2005). In Rank-Abundance diagrams species are ranked from the commonest to rarest on the x-axis and abundance is plotted on the y-axis. We derived species abundance for the whole area as the sum of their presence in single plots.

To quantify the temporal change of floristic heterogeneity, β-diversity was calculated for each sampling year as mean of plot-to-plot dissimilarity values (Whittaker 1972; Legendre et al. 2005). Dissimilarity was computed with the Jaccard index (Jaccard 1900, 1901) using the “vegdist” function from the R package “vegan” (Oksanen et al. 2019). The difference between the two years was tested for significance using the “BetaDispersion 2.0” function developed by Bacaro et al. (2012, 2013). The procedure performs an ANOVA test on differences in mean plot-to-plot dissimilarities between groups (i.e. years) and, by randomization without replacement of the within-group dissimilarities among groups, ensures the selection of an adequate null model.

We calculated the turnover between 1991 and 2018 for the total flora (Jaccard) with the Jaccard index. To quantify changes in floristic composition a “Temporal Beta‐diversity Index” (“TBI”) (Legendre 2019) was performed for each sample plot using the R package “adespatial” (Dray et al. 2019). The index represents the dissimilarity value computed between two presence-absence vectors representing the same plot in 1991 and in 2018. Dissimilarities, once again, were computed with the Jaccard index. To discover significant changes, the TBI index was tested for significance with a permutational approach on the “TBI” R function and p-values were adjusted with the Holm procedure. Moreover, the TBI value was partitioned into its components of “Gain” and “Loss”, which inform about the directionality of the change in each plot (Legendre 2019).

The floristic similarity between plots both on a temporal and spatial scale was also explored with the multivariate-ordination analysis non-metric Multidimensional Scaling (NMDS), performed with the “isoMDS” function on “MASS” R package (Venables and Ripley 2002).

Finally, we regressed the species richness of each plot (richness 1991, richness 2018), the number of aliens (aliens 1991, aliens 2018) and the indices of temporal diversity (TBI, Gain and Loss) against the distance from the border in order to understand how this variable affects the floristic richness and composition as well as the temporal change in β-diversity. For each response variable we created two distinct models placing as a single predictor the distance calculated along the paths (path) and the distance as the crow flies (crow flies). We therefore obtained 14 different models: richness 1991 vs. path, richness 2018 vs. path, richness 1991 vs. crow flies, richness 2018 vs. crow flies, TBI vs. path, Gain vs. path, Loss vs. path, TBI vs. crow flies, Gain vs. crow flies, Loss vs. crow flies. Poisson and quasi-Poisson Generalized Linear Models (GLM), suitable for presence-absence data, were fitted on number of species and aliens, while Binomial Generalized Linear Models (GLM) were performed on TBI and on its terms of decomposition (Gain, Loss). Models assumptions were verified testing for residuals normality.

3 Results

The total number of species collected in 28 1-hectare plots decreased from 560 in 1991 to 536 in 2018, while the alien species increased from 20 to 35. The average number of species per plot increased from 113 to 123 while the mean number of aliens remained essentially the same (8 species in 1991 and 7 in 2018). The Wilcoxon test did not detect significant temporal changes in terms of both floristic richness (p-value = 0.175) and presence of aliens (p-value = 0.848).

The pattern of Rank-Abundance diagrams (Fig. 2) shows a shift from a nearly broken-stick distribution to an essentially geometric model, suggesting the transition from a higher evenness community (1991) to an assemblage with greater species dominance (2018).

Fig. 2
figure 2

Rank-abundance diagrams: dashed line represents 1991 sample; solid line represents 2018 sample

Full size image

The Betadispersion 2.0 function (Table 1) detects a significant loss in \(\upbeta \)-diversity (p-value < 0.05) which drops from 0.840 in 1991 to 0.717 in 2018.

Table 1 Betadispersion 2.0 function, ANOVA output
Full size table

The total flora shows a moderate turnover (Jaccard = 0.38). Instead, TBI analysis shows that a high turnover per plot, on average of 80% (TBI = 0.79), occurred in the area over 27 years. The TBI value (Fig. 3) ranges from 0.65 and 0.91 but no plot changed its floristic composition in an exceptional way (p-values > 0.5; Appendix Table 3). Overall, the temporal change in floristic heterogeneity is more connected to the process of species gain (mean Gain = 0.415) than to a process of species loss (mean Loss = 0.376). The directionality of the change in terms of species gains or losses (Appendix Table 3) is summarized in the map in Fig. 3.

Fig. 3
figure 3

Temporal Beta diversity (TBI) values in each plot. The colour of plots summarizes the directionality of their change

Full size image

The NMDS ordination technique (stress = 21.61%; Fig. 4) shows a higher floristic homogeneity of the 2018-sample plots compared to the same squares investigated 27 years earlier. All plots from the second census are very close to each other, while the ones investigated in 1991 show a greater dispersion within the bidimensional space. Moreover, the NMDS procedure seems to detect a clear-cut temporal distinction along its first dimension (Dim1) as most of plots investigated in 1991 are on the left of the graph while all those investigated during the second survey are located on the right side.

Fig. 4
figure 4

NMDS ordination plot: triangles represent sample plots investigated in 1991, dots represent the same squares investigated in 2018

Full size image

A significant relationship emerged in 2018 between the distance from the boundary calculated along the paths and the number of species (richness 2018 vs. path; Table 2; Fig. 5). The model richness 2018 vs. path returns a moderate value of percentage of explained deviance (21.8%) and this relationship emerges only if the distance is calculated along the paths and not when it is calculated as the crow flies (richness 2018 vs. crow flies). All the other models (Appendix Table 4) are not significant. If an outlier is removed in the model richness 1991 vs. path the predictor becomes marginally significant but with a very shallow slope (Appendix Table 4).

Table 2 GLM models output: species 1991 vs. path and species 2018 vs. path
Full size table
Fig. 5
figure 5

Relation between plot species richness and the distance from the boundary of the park in 1991 (left) and in 2018 (right)

Full size image

The models that correlate the distance from the boundary with the temporal \(\upbeta \)-diversity and its components Gain and Loss, are not significant (Appendix Table 4) suggesting the lack of a spatial structure of these variables. Moreover, these models return negative values of adjusted R2 meaning that the distance from the boundary is not able to predict any variation on the response variable.

4 Discussion

Biotic homogenization, one of the main threats to biodiversity (Kühn and Klotz 2006; McKinney 2006), is defined as the gradual increase of communities’ similarity over time (Rahel 2002; Olden and Poff 2003) and/or between different sites (Devictor et al. 2008), usually as a consequence of the spread of common species and the simultaneous decrease of the rare ones (Mckinney and Lockwood 1999). Researches on homogenization are generally focused on the dynamics related to the spread of the alien species (Schwartz et al. 2006; Dar and Reshi 2014) but the complex phenomenon of biotic homogenization can also be the result of habitat alteration (Rahel 2002) and of the spread of generalist species and the extirpation of rarer species (Olden and Poff 2003; McKinney 2006). Temporal analysis conducted on species pools from the same sites at different time are powerful tools to analyze the process of homogenization (Dar and Reshi 2014).

Our results seem to indicate a process of homogenization in the time laps studied. The significant decrease of β-diversity (p-value < 0.05) of the whole area, which drops form a value of 0.839 in 1991 to 0.717 in 2018, is perhaps the strongest evidence of the decline of floristic heterogeneity. This trend is further highlighted by the ordination analysis (NMDS) in which plots of 2018 are considerably more concentrated than the 1991 samples. Notwithstanding the decline of \(\upbeta \)-diversity we found an increase of α-diversity (113 species per plot in 1991; 123 species per plot in 2018) and the simultaneous decrease in γ-diversity (560 species in 1991, 536 in 2018). Biotic homogenization is not always accompanied by a loss in plant richness (Naaf and Wulf 2010) and similar contrasting trends have been found in several studies (Smart et al. 2006; Jurasinski and Kreyling 2007).

This homogenization is also not related to an increase in the abundance of alien species. They increase in number from 20 species in 1991 to 35 in 2018, but their frequency remained stable: 8 species per plot in 1991, 7 species per plot in 2018. Conversely, we observed a significant turnover, as expressed by the temporal pattern of β-diversity. Turnover is moderate at the level of the total flora (Jaccard = 0.38) and is very high on average in single plots (TBI = 0.79). Species have changed their pattern of distribution in the period analyzed and Rank-abundance diagram (Fig. 2) suggests that changes in dominance relationships occurred. The pattern of 2018 curve is characterized by high dominance of a few species, resulting in the so-called geometric model, while the community of 1991 was more equitable (higher evenness), returning shapes corresponding to the log-normal and broken-stick distributions.

Changes of the community structure can sometimes be related to changes in the regime of disturbance (Battisti et al. 2016). Such a change appears if we consider the distribution of species richness in relationship with distance from the border of the park (Fig. 5). In 1991 it was not possible to detect a distinction between the margin of the park and the core areas, whereas this is evident in 2018 and is expressed by the significant increase of the slope of the curve richness 2018 vs. path (Fig. 5). The relationship becomes nonsignificant if the distance is expressed not along the paths but as the crow flies (richness 2018 vs. crow flies). Our results suggest that the main driver of the ecosystem change is an increase in frequentation of the park through the paths. The park, scarcely frequented in 1991 with the exception perhaps of a very narrow fringe near the main roads, is now experiencing a greater influence of the anthropogenic disturbance in the more accessible areas near the borders.

Frequentation increased not only because of the increase of the population of the city, but because of management choices in the mid-90s, in particular the removal of a sheep herd that used to graze in the park. The traditional management practice discouraged frequentation and maintained the vegetation in a stable, diverse condition, whereas the management promoting a recreational use and accessibility by the people tend to increase the species dominance and the floristic homogenization (Figs. 2 and 4). It is well known that human trampling can have deleterious effects on plant communities by changing their composition, diminishing the vegetation cover, homogenizing them and threatening their conservation (Cole 1995; Hamberg et al. 2008; Rossi et al. 2009). The higher degree of homogenization in habitats subject to greater anthropic disturbance has also been found by Lososová et al. (2012) in several European cities.

In addition to our findings, the removal of sheep herd has led to a floristic ruderalization and an increase of soil nitrogen highlighted in a previous study (Bianco et al. 2003) and it seems to have significantly affected also the insect communities (Carpaneto et al. 2005).

5 Conclusions

The conservation of biodiversity and ecosystem resilience in urban protected areas requires appropriate management (Gaston et al. 2005; Shwartz et al. 2013). A possible approach can entail the maintenance of a traditional landscape resulting from the interaction between extensive human activities and the natural processes of the ecosystems (Berkes et al. 2000; Middleton 2013). The few relics of the “Campagna Romana” landscape are an example of how a moderate human pressure can increase the number of available niches and decrease competition, allowing a great number of species to thrive. Among many other practices, grazing promotes a high environmental heterogeneity (Schwabe et al. 2013) and species diversity (Aavik et al. 2008; Dostálek and Frantík 2008) and sometimes can prevent the increase of ruderal species (Schwabe et al. 2013).

To halt the ongoing homogenization process that is characterizing the Pineto park, the reintroduction of grazing and a limitation of the fruition of the park in outer areas close to the urban matrix seem to be optimal and sustainable management options. Managed grazing, already adopted in other urban parks, can enhance the ecosystem services provided by the green urban areas, for instance dairy products which are already sold in the “Appia Antica Regional Park” in Rome. Besides the maintenance of a high level of environmental diversity, the reintroduction of a moderate grazing could therefore promote the active rediscovery of a rural environment that is becoming increasingly rare in urban context, providing a different, more sustainable way to enjoy nature.