Influence of urbanization on the avian species-area relationship: insights from the breeding birds of Rome

The species-area relationship (SAR) is one of the most investigated patterns in ecology and conservation biology, yet there is no study testing how different levels of urbanization influence its shape. Here we tested the impact of urbanization on avian SARs along a rural-urban gradient using the breeding birds of Rome (Central Italy). We divided the city into 360 cells of 1 km2. Each cell was classified as rural, suburban or urban using the proportion of impervious surface calculated from remote sensing data. For each of these three landscape categories, we constructed a SAR as a species accumulation curve (Gleason function) using bird species distribution data. SAR intercepts (i.e. the number of species per unit area) decreased from rural to urban areas, which indicates that urbanization depressed the number of species, reflecting the loss of specialized species strictly associated with natural habitats. The slope was highest for the rural curve, indicating that natural landscapes have the highest turnover due to their higher habitat heterogeneity. A higher slope for the urban cells, compared to the suburban ones, can be explained by the presence of green spaces embedded in the built-up matrix which host different avian communities. Previous studies that compared whole cities with natural areas failed to find differences in the respective SARs. Our study, which constructed SARs for different levels of urbanization, indicated significant changes in the SARs along the rural-urban gradient. Further analyses in other cities and taxa will be useful to test how general are our findings.


Introduction
Urban expansion is considered one of the major threats to biodiversity (Czech et al. 2000;McKinney 2002McKinney , 2006McDonald et al. 2013), yet the complexity of the urbanization process led most cities to present a mosaic of habitats that supports the presence of many rich biotic communities (Adler and Tanner 2013;Fattorini 2019).
This mosaic of habitats can be ideally summarized by an environmental gradient, ranging from relatively natural areas, typically at the border of the city, to highly modified urban landscapes in its more interior sectors. Such a gradient is known as the urban-rural gradient and it is frequently used to explore the influence of different levels of urbanization on biotic communities (McKinney 2008;Forman 2014).
Although for most animal groups species richness decreases in response to increasing urbanization intensity, there are several studies reporting a higher species richness at intermediate levels of urbanization (Blair 1996;McKinney 2008;Leveau 2019). A few studies report peaks of plant and invertebrate diversity at the highest levels of urbanization, a situation that, however, has been never observed in vertebrates (McKinney 2008). Moreover, the influence of urbanization on biodiversity is negative when analysed at small scales, but may appear positive at larger scales, because of the admixture of different habitat types that allow the presence of species with different ecological needs (e.g., Pautasso 2007). However, studies comparing species richness along the rural-urban gradient can be biased by the different extent of the areas with different levels of urbanization within a city. For example, the surface of the most urbanized sectors (e.g., the city centre) can be much smaller than those of the less urbanized suburbs. Since larger areas tend to have more species, differences in species richness at different levels of urbanization may reflect their different spatial extent. Comparing species richness of areas with different size is not a trivial task. Simply dividing the observed richness by the extent of the area implies that the species-area relationship (SAR, that is, the increase in species number with increasing area) is linear, whereas the relationship is known to be typically curvilinear (Matthews et al. 2016).
The species-area relationship is one of the best-known ecological patterns (Lomolino et al. 2010) and SAR modelling may be used to compare how species accumulates with area increase in systems with different characteristics. In fact, there are two main types of SARs: those arising from true isolates (ISARs, island species-area relationships, in which the sampled areas are isolates and the function fitted is based on how many species are found in each sampled area) and those arising from the progressive aggregation of smaller sampling areas into larger areas (species accumulation curves, SACs, which present cumulative counts of increased species number with sampling area; see Matthews et al. 2016 for the distinction between these two concepts). Namely, Scheiner (2003) distinguished several types of SARs, according to the underlying sampling schemes: (a) strictly nested quadrats (Type I curves); (b) quadrats arrayed in a contiguous grid (Type II curves); (c) quadrats arrayed in a regular but noncontiguous grid (Type III curves); or (d) areas of varying size, often islands (Type IV curves). While Type IV corresponds to the ISARs, the other types are forms of SACs.
Many mathematical functions have been proposed to model the SAR, but the Arrhenius (1921) power function is the most commonly used since it fits best most empirical data and it is possibly supported by ecological theories (Fattorini 2006;Triantis et al. 2012). The power function is S = CA z , where S is the number of species, A is the area and C and z are fitting parameters. This relationship can be linearized in a log-log space as log(S) = log(C) + zlog(A). In this linear form, the two parameters of the function are easier to be interpreted: z is the slope of the regression equation (i.e. how fast species accumulate with increasing area) and C is the intercept of the regression equation (i.e. the expected number of species per unit area).
Although the Arrhenius power function is frequently found to be the best model for ISARs, it can be not appropriate for SACs (Ugland et al. 2003). For SACs, there is theoretical and empirical support for the use of another model, known as the Gleason (1925) semilogarithmic function (Ugland et al. 2003). The Gleason model is: S = zlogA + C, where C is the intercept of the curve in the arithmetic space (the power function does not have an intercept in the arithmetic space) and z is a direct measure of the initial and overall slope (given the particular log-base used). Because the semi-log model has a real intercept (i.e. in arithmetic space), it may seem to have some advantage over the power model (Lomolino 2001).
In both models, C is the expected number of species per unit area, and therefore can be used to compare species richness in systems with different ranges of area size. By contrast, z is an accumulation parameter, showing how fast species accumulate with area. Several authors considered the z-parameter as a measure of β-diversity (e.g., Lennon et al. 2001;Ricotta et al. 2002;Scheiner 2003Scheiner , 2004Drakare et al. 2006;Polyakova et al. 2016;Sreekar et al. 2018). However, the z-value of a SAC does not quantify compositional heterogeneity, but the rate of change in α-diversity with increasing sampling unit size (Tuomisto 2010). Although z alone does not express β-diversity, larger values of z correspond to faster accumulations of species as area increases, indicating higher rates of species turnover, although limited to species gains.
Recently, Fattorini et al. (2018) used the SAR approach to compare the richness of urban and rural avian communities, finding no differences between the two curves. The two SARs (modelled with the power function) were constructed using Italian towns and natural areas of different sizes, thus following an ISAR approach. However, one could ask how the SAR changes not only between cities and natural areas, but also among different levels of urbanization within the same city. This approach can reveal how urbanization intensity influences SAR parameters and can allow the comparison of species richness for areas with different levels of urbanization which cannot be compared directly because of their different extent. To this end, we used the distribution of breeding birds in the city of Rome (Italy) to construct SARs for different levels of urbanization with a SAC approach. We used a grid system to divide the study area into equally sized cells and assigned each cell to one of the following three categories of urbanization: rural, suburban and urban. Then, we calculated the avian species richness of each cell and, for each of the three categories of urbanization, we modelled the respective SAR by constructing an accumulation curve. This way, we were able to calculate the two curve parameters, C and z, and therefore to test not only how species richness differ among sets of areas with different size under different levels of urbanization, but also how urbanization influences species accumulation with area. Since the Arrhenius model is the most commonly used function to model the SAR, but there is support for the use of the Gleason model in the case of SACs, in this study we used both models and compared the results.
There is a vast literature about how SAR parameters are influenced by environmental characteristics (see, for example, Fattorini et al. 2017). A recent global analysis (Beninde et al. 2015) showed that, although other variables can be important in determining bird richness in urban green spaces (namely vegetation characteristics), area is the most important correlate. A recent global study of urban birds also showed that species richness in urban green spaces is influenced by their extent, but not by the size of the human population of the urban areas . Several studies on avian SARs in urban environments (MacGregor-Fors et al. 2011;Pautasso et al. 2011;Ferenc et al. 2014;Fattorini et al. 2018) compared curves obtained from different cities with those constructed using areas from the surrounding natural matrix, thus using an ISAR approach. To the best of our knowledge, however, our research is the first study that investigates how urbanization intensity influences SAR parameters by comparing SARs constructed for sets of areas with different urbanization within a single city and that uses a SAC approach. A common problem in SAR comparisons is that the interpretation of the function parameters is influenced by the different characteristics of the compared systems. For example, for ISAR studies, Gould (1979) and Fattorini et al. (2017) recommended comparing different archipelagos using the same taxon (to investigate the impact of different ecological characteristics on the SAR parameters of a given group), or different taxa using the same archipelago (to investigate how SAR parameters are influenced by a taxon's biology), since equation parameters are taxon-dependent. In the present research, we compared the SACs for different levels of urbanization for the same group (the breeding birds), thus following the "different archipelagos using the same taxon" approach to investigate the impact of urbanization intensity on the parameters of the avian SAR. The within-city approach is also particularly meaningful to eliminate confounding factors that would be present if levels of urbanization were taken from different cities, which would differ for macroecological aspects, such as climate, regional species pool, age, etc.
Urban-to-rural gradient studies showed that urbanization may influence species richness in different ways (reviewed in McKinney 2002(reviewed in McKinney , 2008. A common observed pattern is a decrease in species richness from rural to suburban and finally to urban areas, which is expected if species are negatively affected by the various forms of habitat alterations and changes in ecosystem functioning produced by urbanization (loss of habitats, fragmentation, alterations in interspecific interactions, pollution, human presence, changes in natural disturbance regimes, etc.). For example, Marzluff's (2001) compilation of bird studies found that 61% showed lower species richness in urban-suburban areas compared with more natural rural areas. However, some studies (e.g. Blair 2001) indicate that species richness tends to be higher in areas with intermediate levels of urbanization (suburban areas) than in rural areas, as a result of the higher environmental heterogeneity (e.g. a higher variety of vegetation types, nesting and resting sites, food resources, etc.) and productivity of these anthropogenic habitats (which is in line with the so-called intermediate disturbance hypothesis; see, for an avian example, Malavasi et al. 2009). Finally, very few studies show an increase in species richness in the transition from moderate to high urbanization for plants and certain invertebrates, possibly including large fractions of alien and very generalist species (McKinney 2002(McKinney , 2008, and no bird study revealed such pattern (Marzluff 2001).
In this study, we tested whether the parameter C (which is a measure of α-diversity, being the number of species per unit area) and z (which expresses the rate of species accumulation) are uniform (null-hypothesis) or vary between sets of areas with different urbanization according to the aforementioned three patterns. Namely, we can formulate the following alternative predictions for deviation of C-values from uniformity: (1) C-values can decrease from the rural to the urban areas, if urbanization affects negatively avian diversity (P 1C ); (2) Cvalue can peak in the suburban areas, if avian diversity is promoted by intermediate disturbance (P 2C ); and (3) urban areas can have a C-value higher than that of the suburban areas, but lower than that of the rural ones, if a higher level of urbanization depresses diversity but, at the same times, allows the presence of more urban exploiters that are less favoured in suburban areas (P 3C ). Similarly, for z-value, we can expect deviating from uniformity, because: (1) z-values can decrease from the rural to the urban areas, if urbanization increases faunal homogenization, thus depressing the velocity of species accumulation rate (P 1z ); (2) z-value can peak in the suburban areas, if an intermediate disturbance generates a mosaic of urban/rural habitats that increases spatial species turnover, and if this turnover exceeds that generated by habitat heterogeneity in more natural areas (P 2z ); (3) urban areas can have a higher z-value than suburban areas, but lower than that of rural areas, if the influence of natural habitat heterogeneity overwhelms the heterogeneity represented by the urban mosaic of habitats (P 3z ).

Avian species distribution
We chose the city of Rome (Fig. 1a) as a study area because of (1) the availability of a bird atlas (Cignini and Zapparoli 1996) reporting species distribution across the city with a regular grid system and (2) its wide extent (about 360 km 2 ), which allowed the possibility of constructing SARs with a large number of units for different levels of urbanization, from densely urbanized areas to natural, well-preserved areas. The atlas used a regular UTM grid of 360 squares of 1 km × 1 km as a censusing scheme (Fig. S1, Supplementary material) and reported, for each breeding bird species, the presence/absence in each cell as recorded between 1989 and 1993 by a team of 67 observers. We digitized all atlas maps and used these data to construct a matrix of species occurrences per grid cell (Table S1, Supplementary material).
Rome is located in Central Italy, near the Tyrrhenian coast. The area shows a typically Mediterranean climate, with mean annual rainfall between 810 mm and 940 mm, and mean annual temperature between 14.8°C and 15.6°C. Elevation in the study area ranges between 15 m and 139 m. Vegetation is mainly represented by fragments of mixed woods of Quercus spp. and bushes, grasslands, and pasturelands. Built up areas include, especially in the city centre, numerous ancient buildings and ruins, in addition to modern constructions.
With a population of about 3 million inhabitants, Rome is the largest city in Italy and one of the largest in Europe. As in other researches (see Fattorini 2011), urban Rome was defined here as the territory of the town encompassed by the great motorway ring that circumscribes an area of about 360 km 2 . Approximately one-half of this area is occupied by ruins, historical villas, archaeological sites, meadows, grasslands, gardens, parks, and suburban uncultivated grounds. The Tevere and Aniene rivers occupy an area of about 25 km 2 .
Rome is one of the oldest cities in the world, continuously inhabited for more than 2500 years. For about two millennia Rome area was largely occupied by a semi-natural landscape, the so-called 'Campagna Romana' (a rural ecomosaic of habitats with various degrees of naturalness that originated by the millenary presence of the city of Rome; Celesti Grapow and Fanelli 1993;Lucchese and Pignatti 2009), with built-up areas covering a very small proportion of current city area. Until 1860-1870, when Rome was annexed to the Italian Kingdom, the city was about 14 km 2 wide (with just one third of the town occupied by built-up areas) and had 200,000 inhabitants. The surrounding 'Campagna Romana' (3500 km 2 wide, 4 inhabitants/ km 2 ) was characterized by a very pristine environmental condition, because economic and social conditions, coupled with the presence of malaria, prevented an effective colonization. Most of this territory was occupied by wild grazing and wheat areas. Thus, urban expansion in the twentieth century took place on an area that had been practically uninhabited. Modern urbanization proceeded in an uncontrolled fashion, with relicts of Campagna Romana embedded in the urban tissue, sometimes as green enclaves completely surrounded by built-up areas. A few larger green spaces extend outside the urban territory, being in direct connection with the rural landscape.
Overall, Rome hosts 75 species of breeding birds, which is one of the highest values of species richness recorded in Italian cities (Fattorini et al. 2018). This value is however lower than those recorded from other European cities of comparable size (Ferenc et al. 2014), which can be explained by the southern location of Rome. In Europe, the number of breeding birds in urban areas increases with latitude (i.e. northern cities have more species than southern ones), possibly because higher climatic variability promotes generalist species (which are expected to be favoured in urban areas) (Ferenc et al. 2014(Ferenc et al. , 2019.

Landscape analysis
Urban-to-rural gradient studies typically identify three main levels of urbanization: urban, suburban, rural. However, there is no consensus about how landscapes along the gradient should be classified into these (or other) categories Clergeau et al. 2006;Magura et al. 2008;Iannella et al. 2016;Tzortzakaki et al. 2019). As the percentage of impervious surface (i.e., the cover of soils with impervious materials such as concrete, metal, glass, tarmac, and plastic) is the most obvious feature that influences ecosystem functioning and biodiversity patterns in cities (McKinney 2002;Strohbach et al. 2019;Tzortzakaki et al. 2019), we adopted this measure to quantify the level of urbanization in the study area.
To calculate the percentage of impervious surface in each recording cell, we used a Landsat 5 image taken between 1 and 16 June 1996. We selected this image from a set of Landsat 5 photos available for the period 1990-2000 as it was the most synchronous with the period of bird censusing and showed a minimum amount of cloud cover. Moreover, the selected photo was taken at the start of the summer, when vegetation is more clearly detectable through spectral signatures. Then, we applied the Semi-Automatic Classification Plugin (SCP) algorithm (Congedo 2016) to classify each pixel of this image into one of the following four landscape classes: vegetation (including woods, natural grasslands, and vegetative crops), water (including rivers and lakes), built-up (including buildings and surfaces covered with concrete, asphalt, or similar material) and bare soil (including exposed soil and non-vegetative crops). To create a training file for the automatic classification on the basis of spectral signatures we generated 100 random points on the map, for each point we created a "region of interest" (RoI) and manually assigned each region to one of the abovementioned classes. Then, we merged the spectral signature of the RoI of the same class. To help the manual identification of RoIs, we used a high resolution (5 m) photo (available at: http://wms.pcn. minambiente.it/ogc?map=/ms_ogc/WMS_v1.3/raster/IGM_ 25000.map) of the same area taken in the same period (the LandSat 5 image has a 30 m resolution).
As additional help, we also visualized the Landsat 5 photo using the RGB = NIR Red Green scheme as false-colour composites, which correspond to Landsat 5 bands code 4-3-2, because in this way vegetation appears in red (due its high reflectance in the NIR) and can be hence readily detected (Wessman et al. 1988;Elvidge 1990;Meng et al. 2009). We also used the RGB = Blue NIR Shortwave infrared scheme, which corresponds to Landsat 5 bands code 3-4-7, because it allows a clear distinction between vegetation (that appears in different shades of green) and constructions (shades of purple) (Congedo 2016).
To asses if the four classes could be correctly identified by the automatic classification on the basis of their spectral signature, we applied various distance measures: Jeffries-Matusita, Spectral angle, Euclidean distance, and Bray-Curtis (Congedo 2016). Spectral signatures were statistically different for all these distance measures. However, we obtained the best separations with the Jeffrey-Matusita distance in combination with the maximum likelihood algorithm (see Kavzoglu and Mather 2002; Ahmad and Quegan 2012), which was then selected for the automatic identification. Finally, we assessed the accuracy of the obtained classification by choosing 100 random points on the Landsat photo and checking if they were correctly classified using the same higher resolution map used for the construction of the training set (Congalton and Green 2019). We found that 96% of pixels were classified properly, which indicates that the classification done by the algorithm was accurate. We aggregated the three categories "vegetation", "water", and "bare soils" as "natural", whereas "built-up" was considered as impervious surface. Then, we calculated the proportion of impervious surface to total area (Fig. 1b).
The emerging pattern (Fig. 1c) indicates a somewhat radial urban-rural gradient. To divide the urban-rural gradient on the basis of the relative extent of the impervious surfaces, McKinney (2002McKinney ( , 2008 proposed the following threshold values: >50% of impervious surface: urban core; 20-50%: suburbia; <20%: urban fringe/rural. We adopted the following slightly different thresholds as more suitable for our study area: >60% of impervious surface: urban; 20-60%: suburban; <20%: rural (Fig. 1d). In fact, we evaluated various alternative threshold values and compared visually the resulting spatial patterns. The selected threshold scheme appeared best suited to reflect the landscape structure specific of the study area, as also emerging by vegetation analysis (Fratarcangeli et al. 2019).

SAR construction and comparison
We constructed species-area relationships (SARs) separately for each urbanization level. SARs were constructed as species accumulation curves (SACs), since our data correspond to Scheiner's (2003) Type II of SARs. A SAC is a graph recording the cumulative number of species recorded in a particular environment as a function of the cumulative sampling effort. With this approach we modelled the increase in recovery of new species as a function of the number of cells that are progressively added within each of three levels of urbanization. In constructing an accumulation curve, the order in which samples are included influences the way in which species are added, and hence the resulting species number. To remove the influence of the order of cell accumulation on species richness, we calculate the average number of species for each level of accumulation using the analytical Mao tau sample-based procedure as implemented in EstimateS 9.1.0 software (Colwell 2019), which analytically calculates the species richness as a function of number of samples. This provided, for each area size, the number of breeding species independently from the order of accumulation. Then, we modelled the resulting SARs using the Arrhenius power function model its linearized form (log-log, base 10) and the Gleason semilogarithmic model (log 10 ). Since the R 2 computed from log-transformed variables holds only for the linearized space obtained by the transformation, we fitted the Arrhenius models also using a non-linear algorithm with the R package sars .
We used analysis of covariance (ANCOVA) to test for differences in C and z values between urban, suburban, and rural SARs (Fattorini et al. 2017). ANCOVAs were done in R 3.0 (R Core Team 2015). Significance was set at α = 0.05.

Results
Out of the 360 cells included in the study area, 104 were classified as urban, 157 as suburban, and 99 as rural. Overall, 75 bird species were recorded from the study area. The rural cells hosted, on the whole, 73 species (97.3% of the total avifauna); the rural cells hosted 69 species (92.0%) and the urban cells hosted 62 species (82.7%). In all cases, the Gleason function modelled the SARs better than the Arrhenius function (Table 1, Fig. 2; Tables S2 and 3, Figs. S2 and S3, Supplementary material). Thus, we restricted our comparisons to the curves obtained with the Gleason model.
The SARs modelled using the Gleason function for the three levels of urbanization (Fig. 2) showed an overall difference in their z-values (ANCOVA: F 2 , 3 5 4 = 23.10, p < 0.0000001). Post-hoc Tuckey tests indicated significant differences between slopes in all pairwise comparisons (p < 0.001). In summary, the rural curve had a slope (z ≈ 18) significantly higher than those of the suburban and urban curves, and the urban curve had a slope (z ≈ 17) higher than that of the suburban one (z ≈ 16).
As regards the intercepts, they should be properly compared only for curves that have equal slopes. However, the intercepts resulted significantly different globally (F 2,356 = 3734, p < 0.0000001) and in all pairwise comparisons (Posthoc Tuckey tests, p < 0.001 in all cases), decreasing in the order: rural (C ≈ 37) > suburban (C ≈ 34) > urban (C ≈ 29). However, C-values overestimated species richness when compared with the average number of species for one sample unit in rural, suburban and urban sectors (33.12, 27.59, 24.88, respectively).

Discussion
Our study indicated that the avian species-area relationships (SARs) obtained as species-accumulation curves (SACs) were Table 1 Ordinary least squares regressions for the Gleason model: S = z × log(A) + C, applied to the avian species-area relationship at three levels of urbanization in the city of Rome. S = number of species, A = area, R 2 = goodness of fit, t = Student's t, p = probability. C and z are fitted parameters best described by the Gleason model, with SAR intercepts (i.e. the number of species per unit area) decreasing from rural to urban areas. This is an indication that urbanization depressed the number of species. At the same time, the slope was highest for the rural curve, indicating that natural landscapes have the highest turnover, but the suburban curve had a slope lower than that of the urban one. Although the Arrhenius power function is frequently indicated to be the best fitting model for the SAR, there is some consensus that the SAR can in fact follow different models, and the best fitting function can be found only empirically (Triantis et al. 2012). A possible cause of confusion is that SARs can be constructed in different ways, and while the Arrhenius model might be the most appropriate for island SARs (ISARs), it might not be the best for SACs. We found that avian SARs in the city of Rome are better fitted by the Gleason model than the Arrhenius one. This may be a consequence of the fact that our SARs are SACs, not ISARs, which is in accordance with the idea that species accumulation curves are well approximated by the Gleason model (Ugland et al. 2003).
Differently from previous studies, which found no substantial differences between urban and rural SARs (Pautasso et al. 2011;Ferenc et al. 2014;Fattorini et al. 2018), our results suggest that urbanization intensity affects avian SARs. A possible explanation of this lack of differences in previous studies might be that urban SARs were constructed using the areas of whole cities, which included also peripherical rural areas (so all the three categories urban-suburban-rural). By contrast, we constructed separate SARs for different levels of within-city urbanization, and we found that a higher level of urbanization tends to impact negatively on the avian diversity by reducing the number of species per unit area (i.e., lowering the C parameter) in the order: rural > suburban > urban. This is consistent with our hypothesis H 1c . However, C-values overestimated species richness per unit area by 11% for the rural cells, 20% for the suburban cells, and 14% for the urban cells, respectively. Notably, the most overestimated species richness per unit area is that of the suburban cells, whose SAR was that with the lowest goodness-of-fit. This overestimation might be due to the presence of a threshold effect. It is possible that, below a certain area size, the number of species drops with a higher velocity because of the strong influence of units with very few species due to undersampling or unusually low environmental heterogeneity. However, the use of SARs indicated clearly that avian species richness diminished along the rural-urban gradient, which is consistent with the fact that, for most taxa, species richness decreases along the rural-urban gradient (for reviews see McKinney 2008; Lepczyk et al. 2017). Fig. 2 Avian species-area relationships at three levels of urbanization in the city of Rome (blue circles: rural; green triangles: suburban; yellow squares: urban) modelled by the Gleason function, S = z × log(A) + C, where S is the number of species, A is area, and C and z are fitted parameters. Dashed lines are ordinary least squares regression lines. Regression parameters are given in Table 1. (For interpretation of the references to colour in this figure legend, the reader is referred to the web version of this article) For insular ecosystems, a decreasing trend in the C parameter has been associated with a decrease in the island carrying capacity of compared islands and an increase in their isolation (Fattorini et al. 2017). This explanation may be also applied to our case, because more urbanized sectors are expected to have a lower carrying capacity than rural ones (due to the lower vegetation cover, and hence productivity and biomass) and are more isolated from the avian species pools of the exurban natural areas.
As regards the slope (z), the highest value was found for the rural SAR. Disentangling the effects of the intercept and the slope in a regression line can be difficult, because they are strongly and negatively correlated: raising the slope by pivoting on the joint mean, the intercept drops (Gould 1979). However, in our case, the highest z-value was found for the rural SAR, which also showed the largest C-value. In this case, the increased slope cannot be viewed as influenced by a lowered intercept, which suggests a true higher rate of species accumulation. Thus, we can conclude that the avian SAR of the rural sectors is characterized by a larger number of species per unit area and a higher species accumulation, compared to the suburban and urban sectors. This pattern can be explained by (1) the loss of the most specialized species strictly associated with natural habitats (which leads to a decrease in C-values), and (2) the predominance, in the suburban and urban areas, of synanthropic and generalist species (Faeth et al. 2011), which results in higher faunal homogenization (McKinney 2006), and hence lower species accumulation rates (z-values). Sol et al. (2020), in a global study on bird functional diversity along urbanization gradients, found that the loss of functional diversity in moderately and highly urbanized environments relative to surrounding natural and rural habitats reflects a decrease in both species richness and abundance evenness. It is possible that a reduction in functional diversity, due to the loss of functionally redundant species, is also implied in the reduced species richness observed in the urban sectors of our study area. However, whereas the C-value of the urban cells was significantly lower than that of the suburban ones, thus indicting that increasing urbanization progressively reduces richness, the slope of the urban curve was slightly higher than that of the suburban one, which indicates a higher species accumulation rate in the urban sectors than in the suburban ones. This result is consistent with our hypothesis H 3z and can be explained by the presence of green spaces (e.g. urban parks, gardens, historical villas, etc.) within cells that, because of the prevalence of impervious surfaces, are classified as urban. These small green spaces, strongly isolated by the surrounding built-up matrix, may act as islands hosting different avian communities, thus increasing the species accumulation rate. Previous studies in the study area found that bird richness was correlated with the size of urban parks (Arca et al. 2012) and that old villas (large-sized green areas embedded in the built-up matrix of the inner sectors of the city) may have more species than peripherical fragments (Vignoli et al. 2013). In particular, the presence of old villas might increase substantially the number of species because of vertical layering of vegetation and the structural ageing of their vegetation that facilitates the occurrence of woodland (specialized) species when compared to the spatial (horizontal) heterogeneity occurring in the Campagna Romana (where old woods are rare).
We would like to stress that, in our study system, our minimum sampled area corresponded to the C-value. Thus, one could compare values of species richness per unit area even simply using the estimated values for one quadrat as obtained by the sample rarefaction curves (Mao tau). However, in other cases, the number of species richness per unit area might be below the minimum sampled area; for example, the unit area might be 1 km 2 , but even the smallest sampling unit might be larger than 1 km 2 (a common situation in ISARs, when even the smallest island is larger than 1 km 2 ). In these cases, it is necessary the use of C-values as an extrapolation to obtain the number of species per unit area. The closeness between our Cvalues and the rarefied values per one km 2 indicates that Cvalues are good estimators of species richness per unit area.

Conclusions
We found that the avian species-area relationship has different patterns according to the level of urbanization. C-values (the number of species per unit area) decrease with increasing urbanization. This indicates that urbanization depresses avian richness. Urbanization decreases species accumulation rates from rural to suburban areas, but the urban areas had a higher z than the suburban ones, possibly because of the presence of sparse green spaces in inner city-areas. Our approach to express levels of urbanization as the proportion of impervious surface from remote sensing imagery is simple and can be easily replicated in virtually any urban area. Thus, we hope that our work will promote similar studies in other cities and taxa to investigate how general are the patterns outlined in this paper.
Data Availability Data are provided as Supplementary material.

Compliance with ethical standards
Conflicts of interest/competing interests The authors have no relevant financial or non-financial interests to disclose.
Ethics approval Not applicable.
Consent to participate Not applicable.

Consent for publication Not applicable.
Code availability Not applicable.
Open Access This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if changes were made. The images or other third party material in this article are included in the article's Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article's Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this licence, visit http://creativecommons.org/licenses/by/4.0/.