A global consistent positive effect of urban green area size on bird richness

Abstract

Background

Although the species-urban green area relationship (SARu) has been analyzed worldwide, the global consistency of its parameters, such as the fit and the slope of models, remains unexplored. Moreover, the SARu can be explained by 20 different models. Therefore, our objective was to evaluate which models provide a better explanation of SARus and, focusing on the power model, to evaluate the global heterogeneity in its fit and slope.

Methods

We tested the performance of multiple statistical models in accounting for the way in which species richness increases with area, and examined whether variability in model form was associated with various methodological and environmental factors. Focusing on the power model, we analyzed the global heterogeneity in the fit and slope of the models through a meta-analysis.

Results

Among 20 analyzed models, the linear model provided the best fit to the most datasets, was the top ranked model according to our efficiency criterion, and was the top overall ranked model. The Kobayashi and power models were the second and third overall ranked models, respectively. The number of green areas and the minimum number of species within a green area were the only significant variables explaining the variation in model form and performance, accounting for less than 10% of the variation. Based on the power model, there was a consistent overall fit (r² = 0.50) and positive slope of 0.20 for the species richness increase with area worldwide.

Conclusions

The good fit of the linear model to our SARu datasets contrasts with the non-linear SAR frequently found in true and non-urban habitat island systems; however, this finding may be a result of the small sample size of many SARu datasets. The overall power model slope of 0.20 suggests low levels of isolation among urban green patches, or alternatively that habitat specialist and area sensitive species have already been extirpated from urban green areas.

Background

Urban expansion is a major driver of global environmental change (Grimm et al. 2008). Since 1950, rural populations in the majority of countries have grown slowly whereas urban populations have grown rapidly (United Nations 2014). In fact, it is expected that 66% of the world’s population will live in urban areas by 2050 (United Nations 2014). Thus, a major sustainable development challenge is to create environmentally healthy cities for the large number of people who will reside in urban areas in the near future. The provision of green areas within cities is a key issue in addressing this challenge. Green areas in cities are essential for biodiversity conservation because they constitute refuges that enable certain native taxa to persist in urban environments (Grimm et al. 2008; Sukhdev 2013; Nielsen et al. 2014). Thus, conserving and providing additional green areas is an efficient method to improve both human wellbeing and biodiversity conservation in cities (Miller 2002; Chivian and Bernstein 2004; Dunn et al. 2006; Fuller et al. 2007; Stott et al. 2015).

Urban green areas are generally surrounded by built-up areas, and hence can be regarded as habitat-islands—relatively isolated habitats in which biological interactions between green areas are reduced by matrices dominated by impervious surfaces (Tjørve and Turner 2009; Szlavecz et al. 2011; Matthews 2015; Fattorini et al. 2018a). As urban green areas can thus be viewed through the lens of island theory, it is expected that various island macroecological and biogeographical patterns should be observed in these systems (Matthews 2015; Fattorini et al. 2018a). For example, birds in urban green areas have been shown to follow the general pattern of increasing species richness with increasing island area, i.e. they are characterized by positive species-area relationships (‘SAR’; Fernández-Juricic and Jokimäki 2001; Beninde et al. 2015). However, the factors influencing the way in which species richness accumulates with increasing island area in urban systems (the species-urban green area relationship, hereafter referred to as the SARu) remains little explored. For example, although the power model, which describes a convex increase of the number of species with area (see example in Fig. 1b), has often been assumed and/or shown to provide the best fit to SAR data (Connor and McCoy 1979; Drakare et al. 2006; Watling and Donnelly 2006; Dengler 2009; Triantis et al. 2012; Hanski et al. 2013; Matthews et al. 2016a; but see also Tjørve 2009; Scheiner et al. 2011), a critical examination of the performance of different SAR models to SARu data has rarely been undertaken (but see Natuhara and Imai 1999), and a global scale analysis is lacking. More generally, whilst urban green areas can be regarded as habitat islands, and despite their relatively small size rendering them amenable study systems for testing island theory, few urban ecologists contextualize their work in terms of island biogeography (see Fattorini et al. 2018a, b).

A common aim of many SAR meta-analyses has been to analyze how and why the slope of the power SAR model (z) varies across datasets (Drakare et al. 2006, Triantis et al. 2012, Matthews et al. 2016a), in a way to understand the rate of species loss with habitat loss and fragmentation. Moreover, many factors have been shown to play significant roles in determining the way species richness declines with area at the global scale (He and Legendre 1995; Tjørve 2003; Matthews et al. 2016a). For example, sigmoid SAR models may inform threshold size of green areas for species richness conservation, whereas positive linear models imply a continuous benefit of increasing green areas for species richness. Although the importance of area in driving green area species richness in cities has recently been discussed at the global scale (Beninde et al. 2015), it is still unknown to what degree SARu characteristics (e.g. general form, z values from power model fits) vary across cities (see examples in Fig. 1), and what variables, including both biological (such as the traits of the taxon under study) and environmental (such as the degree of connectivity amongst green areas) variables, drive this variation.

Fig. 1

Examples of different species-urban green area models fitted with some of our datasets. a Linear model fitted with data obtained from Marseille, France, r² = 0.52 (Lizée et al. 2011); b Power model fitted to data from Seoul, South Korea, r² = 0.75 (Park and Lee 2000); c Kobayashi model fitted to data from Rennes, France, r² = 0.53 (Croci et al. 2008); and d Weibull 3 model fitted to data from Santiago de Chile, Chile, r² = 0.43 (Urquiza and Mella 2002)

Here, we compare the fit of 20 competing SAR models to SARu data from 47 cities across five continents. Model performance was evaluated using three criteria: the number of times a model provided the best relative fit to a dataset, model generality and model efficiency (outlined below). We also examine whether the z value of the commonly used power SAR model varies systematically across cities. We tested three primary hypotheses:

• 1. Both ground nesting bird species and migrant bird species have been shown to be more sensitive to urban green area size than other types of birds (Park and Lee 2000; Zhou and Chu 2012). As such, we predict that the functional composition of local bird species assemblages will influence the form of the SARu, resulting in steeper positive slopes (larger z values) in assemblages with a greater proportion of ground-nesting species and/or migrants.

• 2. Datasets containing a large proportion of species with high dispersal ability, indicated by their body size, may result in shallower SARu slopes (smaller z values) due to the exchange of species among green patches, leading to a reduction in the importance of patch area (Hubbell 2001; Drakare et al. 2006).

• 3. Geographical variables have also been shown to affect the form of the SARu, and the SAR more generally. For example, the degree of isolation between green areas is expected to influence the form of the SARu. Thus, we hypothesize that a set of less isolated urban green areas may lower the slope and fit of the SARu by facilitating the movement of birds between green areas and thus reducing the importance of patch area (Fattorini et al. 2018a).

Moreover, we examine the influence of a variety of variables related to the dataset, including species traits, environmental conditions and methodological constraints, on model performance, and the form of the SARu across datasets (Triantis et al. 2012; Matthews et al. 2015). The search for consistent SARu patterns, along with the identification of factors associated with variability in SARu form across cities, will provide useful information (such as improvement in our ability to predict the number of extinctions resulting from urbanization and habitat loss and fragmentation) for the development and conservation of green spaces in cities; a major conservation issue due to global urban expansion and intensification.

Methods

Data collection

We searched scientific articles, theses and unpublished reports in Google Scholar, Scopus (1960‒2013), Biological Abstracts (2002‒2011) and JSTOR (1913‒1999) during 2013. More recently published articles (i.e. during 2014 and 2015) found without a systematic search were also included. The inclusion of theses and unpublished reports was carried out to correct for any potential bias in effect sizes, because published studies tend to have significantly larger effect sizes than studies in the grey literature (Hopewell et al. 2007). We used the keywords [avian OR bird*] AND [green OR park OR cemeter* OR remnant* OR golf] AND urban, in English as well as in Spanish. Finally, we searched for additional articles in reviews of urban birds (Fernández-Juricic and Jokimäki 2001; Chace and Walsh 2006; Garden et al. 2006; Ortega-Álvarez and MacGregor-Fors 2011).

We defined urban green areas as patches dominated by vegetation and surrounded by a matrix of built-up land uses (Matthews et al. 2016a). This definition included urban parks, remnants of native vegetation, golf courses and cemeteries. When full species lists were not reported in the paper, we contacted the authors directly for the relevant data. In addition, for studies focused on a certain bird assemblage (e.g., forest specialists or migrant birds) or that excluded some species, such as the rock dove (Columba livia), the study authors were asked to provide the total number of species in each green area; otherwise, those studies were not considered in the analysis. However, studies that excluded waterfowl were considered for analysis because their low number in urban green areas probably did not affect significantly SARu parameters (Murgui 2007). Only those studies considering more than three green areas were included to enable analysis of variance calculations in a subsequent meta-analysis (see below). We were careful to include each dataset only once to avoid pseudo replication; thus, datasets used in two or more articles were included only once and, whenever data within a study were available for several times of the year (e.g. papers that analyzed seasonal dynamics of communities), we selected data only from the breeding season because they were the most commonly reported. For articles that did not provide species richness data, we contacted the source paper authors for the relevant data, or we took the number of species from graphics within the source papers using DataThief (Tummers 2006).

Whenever possible, we explored the lack of independence between green areas by using maps reported in each study or locating them using Google Earth when the name of the green areas was provided. Green areas separated by less than 100 m were considered non-independent, because it was likely that the same individual bird would use both green areas (Blair 1996; Fernández-Juricic 2001). In such cases, the green area with the highest bird richness was selected for the analysis and the other areas were discarded.

Biological and environmental variables

For each of the collected datasets, we obtained the following biological and environmental variables that could influence SAR parameters: (1) the proportion of migratory species; (2) the proportion of bird species nesting at low or medium height, on shrubs or on the ground (hereafter ground nesting); (3) the mean length of all species; (4) the biome; (5) latitude of the city; (6) year of city foundation; (7) human population size of the city at the study time; (8) minimum distance to other urban green areas (hereafter called Isolation 1); and (9) minimum distance to non-urban areas (hereafter called Isolation 2). Non-urban areas are different land use types out of the city area that do not contain buildings or other impervious surfaces. Variables 1‒3 were calculated using information provided in each study; otherwise, we found the information in handbooks and field guides (Mitchell 1957; Yamashina 1961; Peterson et al. 1973; Wild Bird Society of Japan 1982; Hilty et al. 1986; National Geographic Society 1999; Del Hoyo 1994–2011; Kazmierczak and van Perlo 2000; Hume 2002; Hilty 2002; de la Peña 2010). Information about species body size was obtained from the Handbook of the Birds of the World (Del Hoyo 1994–2011), and we took the mean body size of all species in each study. Mean body size was a proxy of the dispersal ability of each species (Paradis et al. 1998; Sutherland et al. 2000). To obtain the biome classification of each study, we used the georeferenced map of Olson et al. (2001). Information on human population size, year of foundation of the city and latitude were obtained from the source articles or, in cases where these data were not provided, from Wikipedia. Year of foundation was an indicator of the time since birds species colonized cities and has shown to be positively related to their densities in urban areas (Møller et al. 2012).

Isolation measures (variables 8 and 9) were obtained using information in the source papers. First, where possible, we located each green area on Google Earth by using the information supplied in each article (maps, name of green areas); otherwise, we contacted the source paper authors for this information. Then, for each urban green area within a study, the minimum distance to another urban green area of at least 1 ha was obtained using functions in Google Earth Pro and the mean distance for all green areas was calculated (Isolation 1 metric). Isolation 1 was only calculated for studies in which we obtained information on distance for at least 50% of the urban green areas. The minimum distance from all green areas within a study to the nearest non-urban area (Isolation 2) was obtained from maps provided in the source paper, or when the relevant data were not available, from Google Earth. Datasets without isolation data were still included to enable analysis regarding the importance of the other methodological and biological variables.

We also obtained biological data in regards to the number of species in each dataset: the minimum number of species on a green urban area in a dataset (S_min), the maximum number of species (S_max) and the range of species richness values (S_scale).

Methodological variables

For each of the collected datasets, we also derived the following methodological variables (Matthews et al. 2016a): (1) the number of green areas (N_i); (2) bird survey sampling method; (3) the season when bird surveys were carried out; (4) the smallest green area (A_min); (5) the largest green area (A_max); (6) the range of green area sizes (A_scale; calculated as A_max/A_min); and (7) additional characteristics of the study’s sampling method. The bird survey sampling method was classified as one of: transects, point counts, both, or intensive search. The ‘both’ classification related to studies in which point counts were used in small green areas, and transects in large green areas. The intensive search category related to studies that included multiple survey methods within individual green areas (e.g. transects and mist nets). The season the bird survey was undertaken was classified as one of: breeding season, non-breeding season, and both seasons. Finally, we made a classification of studies according to the following three aspects of the sampling methodology: (1) whether studies specified the sampling method technique; (2) the number of times each green area was sampled; and (3) whether the amount of sampling units increased with the size of the green area. Based on these classifications, studies were classified into three types: (1) stringent—studies that specified the sampling method, each green area was sampled four or more times, and the number of visits increased with green area size; (2) intermediate—studies that specified the sampling method but for which one of the other two conditions of stringent studies were not met; and (3) lax—studies that did not specify the sampling methods and for which both the other two conditions were not met.

Model comparison

Our analyses were focused on island species-area relationships (cf. Matthews et al. 2016a), that is, SAR/SARus in which each data point was an individual green urban area. For datasets in which we obtained the area and species richness of each individual green area, we compared the fit of twenty SAR models using an information theoretic approach (Burnham and Anderson 2002) and the methodology outlined in Triantis et al. (2012); see also Matthews et al. (2016a). With the exception of the linear model, the 19 SAR models were fitted using the ‘sars’ R package (Matthews et al. 2019) which contains functions to non-linear model fitting (see Triantis et al. 2012). The linear model was fitted using standard linear regression in R (R Core Team 2019). A given model fit was considered to be satisfactory if model residuals were normally distributed (using a Shapiro normality test) and homoscedastic (using Pearson correlations), and the optimization algorithm converged.

For each dataset, we compared model fits using the Akaike information criterion corrected for small sample size (AIC_c; Burnham and Anderson 2002). We considered the smallest AIC_c value to represent the best model and all models within < 2 ΔAIC_c of the best model were considered as having similar support (Burnham and Anderson 2002). From these AIC_c values we also calculated the AIC_c weights (wAIC_c) and combined them to form a model selection profile.

We determined the observed shape (linear, convex or sigmoid) of the best fitting model using the algorithm in Triantis et al. (2012; note that the observed shape can occasionally be different to the general shape of the model). We determined if an observed fit was asymptotic using the method of Triantis et al. (2012), whereby, for asymptotic models, we analyzed the fitted parameters to check whether the estimated asymptote was within the range of the data.

Considering the model fits across all datasets, we calculated model generality, efficiency and overall model rank. Generality was calculated as the proportion of datasets for which a model provided a satisfactory fit, whilst efficiency represented the average wAIC_c for all datasets in which a model provided a satisfactory fit. According to Triantis et al. (2012), the overall model rank was derived by standardizing generality and efficiency values [(criterion value − mean criterion value)/standard deviation] and summing the resultant values. Based on the findings of Matthews et al. (2016a), as a final check to determine whether the overall ranks of the linear model and the power model were sensitive to the number of islands in a dataset, we undertook the model comparison using all datasets with seven or more islands and computed the model ranks. We then repeated this procedure using all datasets with eight or more islands, and so on, calculating the models ranks each time (see Matthews et al. 2016a).

Factors explaining variation in SARu patterns

Factors explaining variation in the model selection profile

Following Matthews et al. (2016a), we used constrained analysis of principal coordinates (CAP; Bray–Curtis dissimilarity, 9999 permutations; Anderson and Willis 2003) to assess whether any of our predictor variables explained variation in the model selection profile. CAP is a redundancy analysis of the results of Metric (Classical) Multidimensional Scaling, and it allows for the use of non-Euclidean distances, in this case Bray–Curtis dissimilarity. If Euclidean distances are used, CAP analysis produces the same results as a redundancy analysis (see Anderson and Willis 2003, for a full outline of the method). This was done twice, once considering the best fitting model, and once considering the best model shape (see Matthews et al. 2016a for further details).

For certain datasets, we were not able to obtain estimates for all of the explanatory variables. Thus, we ran the aforementioned CAP analyses twice: once considering all datasets, but excluding the variables for which we did not have complete coverage, and once using all predictor variables and the subset of datasets for which we had data for all variables (n = 21; note that this involved re-running the model selection with these 21 datasets as the AIC_c values needed to be re-calculated). To determine which variables to include in the model, we tested for variable normality using histograms, and multicollinearity using variance inflation factors (VIFs). Several variables were log transformed (a constant, 0.1, was added to variables including zero values) to induce normality, and certain variables were removed due to multicollinearity. The year of foundation was transformed by first taking the absolute value of the oldest year of foundation in our dataset (a constant of 1 was added), and then adding this value to each year of foundation value. In summary, the following predictor variables were used for the all dataset analyses: city population (log transformed), year of foundation (squared), latitude (squared), number of green urban areas (N_i, log transformed), the range of green area sizes (A_scale, log transformed), the minimum number of bird species in green areas (S_min, log transformed) and the range of bird species richness (S_scale, log transformed). All VIFs were under 2. For the subset of datasets analyses, we used the same variables minus S_min, with the addition of both isolation measures (log transformed), the proportion of migrant species (log transformed), proportion of ground nesting species and mean body size (log transformed). Non-significant terms were removed using a backwards selection procedure in order to derive a minimum adequate model.

Factors influencing the fit and slope of the power SARu model

As the power model is the most widely used model in SAR studies, we also focused on the parameters and fit of the power model on its own. The power model is given by the equation S = cA^z, where S is the number of species, A is the area, and c and z are free parameters. In this paper, we used the logarithmically transformed version of the model, which describes a linear increase in logS with logA, with a given slope (z) and intercept (c) (Rosenzweig 1985; Hubbell 2001).

We used the statistical technique of meta-analysis (Hedges and Olkin 1985; Rosenberg et al. 2000) to test for heterogeneity in the fit (r, correlation coefficient) and the slope parameter of the power model (z), and to determine whether any of our biological and methodological variables were related to these properties. We focused on z as c has been shown to be dependent on sample size, making comparisons across datasets problematic (Rosenzweig 1985; Drakare et al. 2006). Using a statistical meta-analytic approach, we estimated the overall magnitude of the power model parameters, giving us a measure of the ‘effect size’. We obtained a weighted overall effect size by taking into account the sampling variance of each study; this is inversely proportional to the effect size of each SARu fit and slope (Rosenberg et al. 2000).

The effect size of the slope was the z value, and its variance was obtained by squaring the standard error (Drakare et al. 2006). We obtained the z and r values from the published papers. Otherwise, we obtained data of bird species richness and green area size by using information in tables of the published papers, data from figures using Datathief or contacting authors of the papers. Then, we used the log transformed values of bird species richness and area in a simple linear regression model. The fit of the power model (r) was obtained by the square root of the regression coefficient. We used the Fischer-Z-transformed correlation coefficients (r_Z) as an effect size for the goodness-of-fit of the association between species richness and area (Rosenberg et al. 2000), where r_Z = 0.5 ln[(1 + r)/(1 − r)]. For r_Z, the variance was calculated from the number of observations (N) as 1/(N − 3) (Borenstein et al. 2009).

Effect sizes were estimated using random-effect models, which are an appropriate approach in ecological studies (Gurevitch and Hedges 1999). For both the z values and r_Z we calculated the global effect size. The significance of the effect sizes was obtained by using 999 randomizations of the effect sizes to calculate confidence intervals and evaluating the non-overlapping of bias corrected confidence intervals (95% CI). To evaluate heterogeneity in effect sizes among case studies, the Q test was employed. Q tests were used to analyze the global heterogeneity among studies and the difference among studies related to categorical variables (Borenstein et al. 2009). The relationships between effect size and the continuous variables were analyzed by meta-regression, which is a weighted linear regression (Borenstein et al. 2009); the significance of the meta-regression was obtained using 999 randomizations.

As it has been shown that studies with significant results are more likely to be published and, therefore, to be included in any meta-analysis, there is a chance of publication bias in our analysis (Borenstein et al. 2009). We estimated the fail-safe number (Rosenthal 1979) that provides the number of non-significant cases needed to turn the results of a given meta-analysis non-significant. A high fail-safe number relative to the number of case studies suggests the absence of publication bias (Rosenberg 2005). A robust fail-safe number should be higher than 5n + 10, where n is the number of case studies in the meta-analysis. All meta-analyses were undertaken using Metawin 2.0 (Rosenberg et al. 2000).

Results

Of the 101 articles we found that analyzed bird communities in urban green areas, relevant information was found in 49 studies. These studies were published between 1976 and 2015 (mean = 2004). Of these 49 studies, 44 were conducted solely in urban parks, four in urban parks and golf courses or cemeteries, and one study was conducted only in cemeteries. A total of 17 studies (35%) were from Europe, 14 (29%) from South America, nine (18%) from North America, seven (14%) from Asia and two (4%) from Oceania (Fig. 2). Of the 49 studies, the sampling methodology employed was classified as stringent in 22 (45%), intermediate in 17 (35%), and lax in 10 (20%) (Additional file 1: Table S1). For 46 studies, we obtained data on the area and species richness of each green area in the study (see below).

Fig. 2

Location of the case studies included in the meta-analysis

Multi-model comparison and the best SARu model

Considering the 46 datasets with area and species richness data, we were able to perform model comparisons using 37 datasets. Eight of the remaining nine datasets were excluded as they had fewer than seven islands (seven is the minimum required to calculated AIC_c with these 20 models; see Triantis et al. 2012) and one dataset was excluded as none of the twenty models provided a satisfactory fit to the data according to our criteria. Overall, the linear model provided the best fit to the most datasets (n = 10; Fig. 3a). In regards to the proportion of datasets in which a given model provided a satisfactory fit according to our criteria (i.e. our generality metric), the weibull4 model was the best ranked model (Fig. 3b), followed by a number of the other more complex models. However, the top 13 models all had very similar generality scores (Fig. 3b), indicating that most models provided satisfactory fits to the datasets. With respect to the efficiency metric, the linear model performed the best, followed by the power and Kobayashi models (Fig. 3c). Considering the overall model rank, the linear model was ranked first, followed by the Kobayashi and power models (Table 1). However, our analysis of model rank sensitivity to sample size indicated that the overall rank of the linear model substantially decreased as datasets with fewer islands were iteratively removed (Fig. 4); in contrast, the rank of the power model (Fig. 4) and the Kobayashi model was relatively insensitive to variations in the number of islands in a dataset.

Fig. 3

The performance of 20 species-area relationship models fitted to 37 urban green area datasets. Performance was measured in three ways: a the proportion of studies for which a model provided the best fit; b the proportion of studies for which a given model provided a satisfactory fit (generality); and c the average AICc weight for datasets that fitted satisfactorily a given model (efficiency). See Table 1 for model names

Table 1 The twenty SAR models compared in the model selection, including the number of parameters, general model shape, whether or not the model contains an asymptote, and the overall rank obtained based on the best fitting model, model generality and efficiency (see “Methods” section)

Fig. 4

Change in model performance with the number of urban green areas (islands) in a dataset. Results for two models are displayed: a linear and b power models. Model ranks were first determined for all datasets with seven or more islands, and then for all datasets with eight or more islands, and so on, iteratively up to datasets with 20 or more islands

With respect to the observed model shape, the average wAIC_c of convex model shapes across datasets was 0.77, compared to 0.12 for sigmoidal models and 0.20 for linear models. The presence of asymptotic model fits was relatively rare (average wAIC_c for non-asymptotic model fits = 0.68; average wAIC_c for asymptotic model fits = 0.14).

Factors explaining variation in the model selection and shape profiles

All datasets

When the model selection profiles of all datasets using a subset of predictor variables (n = 37) were used as the response variable, and 0.05 was used as the significance threshold for predictor variables in the CAP analysis, the only significant predictor variable was the number of islands (F = 5.6; P < 0.001), and this explained only 8.5% of the variation in the data. The performance of the linear model was influenced by the number of green areas and the minimum number of species analyzed in each study, whereas the performance of more complex models such as Beta-P and Extended Power 1 were positively related to the number of green areas in each dataset (Fig. 5a). When model shape (i.e. vectors of wAICc summed across models for each shape) was used as the response variable, the only significant predictor was the minimum species richness observed in a dataset (F = 2.0; P = 0.03), but again, this only explained a small amount of variation in the data (5.4%). The performance of the linear model was positively related to the minimum number of species in each dataset, whereas convex models were negatively related to that variable (Fig. 5b).

Fig. 5

Constrained analysis of principal coordinates (CAP) analysis plots showing the effect of a the number of islands (logA_n) in a dataset on the AIC_c model selection profiles, and b the minimum number of species (logS_min) on an island in a dataset on the model shape profiles. CAP analysis is a redundancy analysis of the results of metric multidimensional scaling. The CAP analysis was undertaken to assess how much of the model selection profiles and model shape profiles were explained by our predictor variables (see “Methods”). Non-significant predictors were removed using a backwards selection procedure. Bray–Curtis dissimilarities were used in both cases

Subset of datasets including all predictor variables

Re-running the model selection using the subset of datasets that included data for all predictor variables (n = 21) resulted in the same qualitative results; that is, the linear model still provided the best fit to the most datasets and had the highest efficiency value (results not shown). Interestingly, when the model selection profiles of the subset of datasets that included data for all predictor variables (n = 21) were used as the response variable, with a 0.05 significance threshold, the CAP model selection analysis did not find any significant variables for either the model selection profile or the model shape profile. That is, the final models only contained one variable that was not significant.

Factors explaining variation in the fit and slope of the power SARu model

We obtained 49 r_Z values and 45 slope values. The mean global fit of the power model (Pearson correlation coefficient) was 0.71 (CI 0.65, 0.75), explaining 50% of the variance in the data. There was no significant heterogeneity in r_Z values among studies (Q = 49.34, df = 48, P = 0.42; see Additional file 2: Fig. S1a), suggesting consistency of the fit at the global scale. Therefore, regardless of differences among cities related to size or geographic location, there was a similar rate of species increase with area. The fail-safe number was higher than the number of studies (5063 vs 255, respectively), indicating that there was no publication bias among studies. We found no significant effects of the moderator variables on the fit of the SARu (Table 2).

Table 2 Results of a weighted mixed-model meta-analysis assessing how the fit r_Z of the power model species-area relationship varies between urban green areas in different cities, in relation to continuous and categorical moderator variables

The z value (slope) showed significant consistency among studies (Q = 47.10, df = 45, P = 0.39; see Additional file 1: Fig S1b) and the mean global value was 0.20 (CI 0.17, 0.23), which is quite close to the canonical value of 0.25. The fail-safe number was again higher than the number of studies (3186 vs 255, respectively), suggesting the absence of publication bias among studies. The slope values were unaffected by the moderator variables (Table 3).

Table 3 Results of a weighted mixed-model meta-analysis assessing how the slope (z) of the power species-area relationship model varies between urban green areas in different cities, in relation to continuous and categorical moderator variables

Discussion

Our results show that increasing species richness with increasing area is a consistent pattern in urban green areas across cities worldwide. Beyond differences among cities related to size or geographic location, there was a similar rate of species increase with area. Although most of the SAR models provided reasonable fits to our SARu datasets, the linear model was the top ranked model overall, indicating a constant increase of species with green area size. However, the performance of the linear model was influenced by the number of green areas and the minimum number of species analyzed in each study. Arguably, the power and Kobayashi models were stronger models predicting the increase of species richness with area, indicating a decelerated increase of species at the highest green area sizes. Put another way, our results suggest that species loss due to habitat fragmentation is faster at the smallest green area sizes.

Unlike previous reviews (Magle et al. 2012; Beninde et al. 2015) and global scale studies (Pautasso et al. 2011) on biodiversity in urban areas, our analysis included a relatively high proportion of studies conducted in South America. The inclusion of Spanish key words in our search is probably the main factor determining the high proportion of South American studies. The lack of studies we found from Africa (see also Nielsen et al. 2014), a continent that has experienced substantial urban expansion in recent decades (Seto et al. 2011), is remarkable and highlights either a need for future data collection in Africa, or for meta-analyses to include key words in searches in a wider range of languages, or both.

Global variation of SARu model form

Our study showed that the way in which species richness increases with green area size was related to the number of green areas in each study. Species richness increased constantly with area in those case studies that contained the least number of green areas, as previously found by Matthews et al. (2016a). This may somewhat reflected a scale-dependency of SARs (i.e. the greater the number of green areas, the larger the spatial extent covered by the system of green areas within a city), although the amount of variation in the model selection profile explained by this variable was low. In the datasets with fewer islands, the addition of the intercept in the linear model may improve its performance relative to the power model which fixes the intercept to zero (zero area, zero species). In contrast, the performance of the Kobayashi and power models, the next two highest ranked models, was not influenced by the number of islands. Both the Kobayashi and power models have a convex shape, lack of an asymptote and are simpler than most of the other models examined, such as the rational model and the Lomolino model (Guilhaumon et al. 2010; Matthews et al. 2016a). Thus, in general our analyses indicate that simpler models outperform more complex models, i.e. those with higher numbers of parameters, when fitted to SARu data; this finding is consistent with previous quantitative comparisons of SAR form (Connor and McCoy 1979; Triantis et al. 2012; Matthews et al. 2016a). Moreover, the convex shape of these models and the lack of asymptote suggest that species loss due to green area loss is the fastest at the smallest green area sizes, whereas species increase with area is the slowest at the greatest green area sizes.

Several methodological constraints are known to influence the form of the SAR. For example, different functions or shapes of SAR may exhibit scale-dependency (Drakare et al. 2006; Triantis et al. 2012), and the range of island areas can influence the suitability of different models to fit the SAR (He and Legendre 1995; Tjørve 2003; Matthews et al. 2016a). Studies that cover relatively small and intermediate island sizes are expected to follow exponential and power model shapes, respectively; whereas studies that span small to very large areas are expected to be better fit by sigmoidal models (He and Legendre 1995; see also Connor and McCoy 1979, for discussion). Nonetheless, our results showed no significant effect of the range of green area sizes on model type or shape, although we found little support for sigmoidal models in general in our SARu datasets. A possible explanation is that the majority of SARu datasets do not contain islands of the size predicted to generate sigmoidal SAR patterns (He and Legendre 1995; Matthews et al. 2015). In addition, the range of green area sizes within individual datasets may not be large enough to exert an effect on SAR forms.

Consistency in the fit and slope of the power model fitted to SARu

The average slope of the SARu power model, obtained after standardization for differences in sample size using a meta-analysis approach (global z = 0.20), was consistent across cities worldwide. This value is towards the lower end of the range of values theoretically predicted for island archipelagos (z = 0.20‒0.35: Preston 1962; MacArthur and Wilson 1967; see Connor and McCoy 1979, for review), although our z value approaches the upper value for terrestrial or intraprovincial SARs (Rosenzweig 1985). However, it is relatively low compared with that observed in terrestrial habitats for several taxa, including birds (z = 0.25 for independent islands; see Fig. 1 in Drakare et al. 2006), and in urban datasets for multiple taxa (z = 0.27, Matthews et al. 2016a). This suggests that, whilst urban green areas can be considered as isolates (Szlavecz et al. 2011; Matthews 2015; Fattorini et al. 2018a), as far as birds are concerned, urban green areas are distinguishable from other true and habitat island systems due to the relatively gradual rate of increase in species richness with area. This could potentially be the result of low turnover rates of bird species among patches within systems of urban green areas, which is consistent with the species-fragmented area relationship perspective (Hanski et al. 2013), but as we looked at island SARs as opposed to accumulation curves (Matthews et al. 2016b), further tests are needed to evaluate this hypothesis.

Our results indicate that, on a global scale, the slope of the power model fitted to SARu data was significantly positive and consistent across different urban environments. This suggests both that the size of green areas is a fundamental determinant of bird species richness in cities worldwide, and that the habitat structure of green areas and the urban matrix across cities are likely similar, which results in a similar response of bird species to increasing patch area. Alternatively, bird communities inhabiting urban green areas in cities may have similar functional roles, such as feeding on seeds and fruits on midstorey or understory vegetation (La Sorte et al. 2018), which result in filtered assemblages that have little functional roles in common with native bird communities (Evans et al. 2018; Vaccaro et al. 2019). Nonetheless, we found that the size of urban green areas explained about 50% of the mean local variation in bird richness, and the slope and fit values were non-significant for approximately a third of the datasets (Additional file 2: Fig. S1a, b, see Additional file 3: Appendix), which is likely due to other predictor variables not considered in this study. Future studies are needed to assess whether other predictor variables representing variation in habitat structure (related to nesting site availability) and primary productivity (related to food resources), such as vegetation volume (Chavez-Almonacid 2014), percent cover of native or non-native vegetation (Urquiza and Mella 2002; Garaffa et al. 2009), NDVI (Bino et al. 2008), habitat diversity (Faggi and Perepelizin 2006), the distance to other green areas (Batllori and Uribe 1998; Urquiza and Mella 2002; MacGregor-Fors and Ortega-Álvarez 2011) or the building cover surrounding green areas (Leveau and Leveau 2016), are also significant determinants of SARu parameters.

Several quantitative syntheses of SARs have found systematic variation in the parameters of the power model as a function of, amongst other things, the sampling scheme employed, spatial scale, types of organisms, habitats or ecosystems studied, and matrix type (Rahbek 1997; Drakare et al. 2006; Matthews et al. 2016a). However, contrary to our predictions, the different variables we examined had a negligible influence on the slope of the power model. For example, there was no significant effect of bird traits, such as migratory and nesting behaviors. A possible explanation for this finding is that migratory and ground nesting species are responding to other environmental factors, such as the distance to rural areas (MacGregor-Fors et al. 2010), the proportion of green areas in surrounding parks (Husté and Boulinier 2011; Leveau and Leveau 2016), or the presence of nest predators (Jokimäki and Huhta 2000). Our results also showed a lack of a significant effect of connectivity on the slope and fit of the SARus, in agreement with Beninde et al. (2015) who also found no significant effect of stepping stones on bird diversity in urban areas.

Conclusions

The consistency of the SARu, in combination with the low average slope of the power model we have observed, may ultimately be associated with the occurrence of many synanthropic generalist species that respond to urbanization by increasing their abundance to the detriment of habitat specialist species (Faeth et al. 2011 and references therein), resulting in biotic homogenization (McKinney 2006). The elevated abundance of such urban adapted species along with the loss of habitat specialists in urban environments (Shochat et al. 2006; Faeth et al. 2011) may promote a departure from lognormality and equilibrium conditions predicted by the power SAR model for geographical islands (Preston 1962; Connor and McCoy 1979), thus favoring a linear SARu. Given that green areas can be composed of non-native vegetation which may negatively affect bird habitat specialists (Munyenyembe et al. 1989; Burghardt et al. 2009), we recommend conserving or restoring native vegetation communities in urban green areas. Moreover, further studies are needed for a comprehensive assessment of the influence of species-abundance distributions on the shape of the SARu. Alternatively, the good performance of the linear model in our model comparison may simply be due to various characteristics of the datasets analyzed, such as the number of green areas sampled, or the minimum number of species in a set of green areas. Overall, we suggest that management actions to increase the size of green areas would have positive effects on the conservation of urban biodiversity worldwide.