A global consistent positive effect of urban green area size on bird richness
Avian Research volume 10, Article number: 30 (2019)
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.
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.
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 (r2 = 0.50) and positive slope of 0.20 for the species richness increase with area worldwide.
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.
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.
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:
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.
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).
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.
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 (Smin), the maximum number of species (Smax) and the range of species richness values (Sscale).
For each of the collected datasets, we also derived the following methodological variables (Matthews et al. 2016a): (1) the number of green areas (Ni); (2) bird survey sampling method; (3) the season when bird surveys were carried out; (4) the smallest green area (Amin); (5) the largest green area (Amax); (6) the range of green area sizes (Ascale; calculated as Amax/Amin); 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.
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 (AICc; Burnham and Anderson 2002). We considered the smallest AICc value to represent the best model and all models within < 2 ΔAICc of the best model were considered as having similar support (Burnham and Anderson 2002). From these AICc values we also calculated the AICc weights (wAICc) 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 wAICc 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 AICc 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 (Ni, log transformed), the range of green area sizes (Ascale, log transformed), the minimum number of bird species in green areas (Smin, log transformed) and the range of bird species richness (Sscale, log transformed). All VIFs were under 2. For the subset of datasets analyses, we used the same variables minus Smin, 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 = cAz, 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 (rZ) as an effect size for the goodness-of-fit of the association between species richness and area (Rosenberg et al. 2000), where rZ = 0.5 ln[(1 + r)/(1 − r)]. For rZ, 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 rZ 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).
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).
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 AICc 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.
With respect to the observed model shape, the average wAICc 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 wAICc for non-asymptotic model fits = 0.68; average wAICc for asymptotic model fits = 0.14).
Factors explaining variation in the model selection and shape profiles
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).
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 rZ 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 rZ 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).
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).
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.
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.
Anderson MJ, Willis TJ. Canonical analysis of principal coordinates: a useful method of constrained ordination for ecology. Ecology. 2003;84:511–25.
Batllori X, Uribe F. Aves nidificantes de los jardines de Barcelona. Misc Zool. 1998;12:283–93.
Beninde J, Veith M, Hochkirch A. Biodiversity in cities needs space: a meta-analysis of factors determining intra-urban biodiversity variation. Ecol Lett. 2015;18:581–92.
Bino G, Levin N, Darawshi S, Van Der Hal N, Reich-Solomon A, Kark S. Accurate prediction of bird species richness patterns in an urban environment using Landsat-derived NDVI and spectral unmixing. Int J Remote Sens. 2008;29:3675–700.
Blair RB. Land use and avian species diversity along an urban gradient. Ecol Appl. 1996;6:506–19.
Borenstein MH, Higgins LV, Rothstein JPT. Introduction to meta-analysis. Chichester: Wiley; 2009.
Burghardt KT, Tallamy DW, Gregory Shriver W. Impact of native plants on bird and butterfly biodiversity in suburban landscapes. Conserv Biol. 2009;23:219–24.
Burnham KP, Anderson DR. Model selection and multimodel inference: a practical information-theoretic approach. New York: Springer Science & Business Media; 2002.
Chace JF, Walsh JJ. Urban effects on native avifauna: a review. Landsc Urban Plan. 2006;74:46–69.
Chavez-Almonacid CA. Relación entre la avifauna, la vegetación y las construcciones en plazas y parques de la ciudad de Valdivia. Tesis de licenciatura: Universidad Austral de Chile, Valdivia; 2014.
Chivian E, Bernstein AS. Embedded in nature: human health and biodiversity. Environ Health Perspect. 2004;112:A12.
Connor EF, McCoy ED. The statistics and biology of the species-area relationship. Am Nat. 1979;113:791–833.
Croci S, Butet A, Georges A, Aguejdad R, Clergeau P. Small urban woodlands as biodiversity conservation hot-spot: a multi-taxon approach. Landsc Ecol. 2008;23:1171–86.
De la Peña M. Nidos de aves argentinas. Santa Fe: Universidad Nacional del Litoral; 2010.
Del Hoyo J, Elliott A, Christie D (1994–2011) Handbook of the birds of the world. Barcelona: Lynx editions
Dengler J. Which function describes the species–area relationship best? A review and empirical evaluation. J Biogeogr. 2009;36:728–44.
Drakare S, Lennon JJ, Hillebrand H. The imprint of the geographical, evolutionary and ecological context on species–area relationships. Ecol Lett. 2006;9:215–27.
Dunn RR, Gavin MC, Sanchez MC, Solomon JN. The pigeon paradox: dependence of global conservation on urban nature. Conserv Biol. 2006;20:1814–6.
Evans BS, Reitsma R, Hurlbert AH, Marra PP. Environmental filtering of avian communities along a rural-to-urban gradient in Greater Washington, DC, USA. Ecosphere. 2018;9:2402.
Faeth SH, Bang C, Saari S. Urban biodiversity: patterns and mechanisms. Ann NY Acad Sci. 2011;1223:69–81.
Faggi A, Perepelizin P. Riqueza de aves a lo largo de un gradiente de urbanización en la ciudad de Buenos Aires. Revista del Museo Argentino de Ciencias Naturales nueva serie. 2006;8:289–97.
Fattorini S, Mantoni C, De Simoni L, Galassi D. Island biogeography of insect conservation in urban green spaces. Environ Conserv. 2018a;45:1–10.
Fattorini S, Lin G, Mantoni C. Avian species–area relationships indicate that towns are not different from natural areas. Environ Conserv. 2018b;45:419–24.
Fernández-Juricic E. Avian spatial segregation at edges and interiors of urban parks in Madrid, Spain. Biodivers Conserv. 2001;10:1303–16.
Fernández-Juricic E, Jokimäki J. A habitat island approach to conserving birds in urban landscapes: case studies from southern and northern Europe. Biodivers Conserv. 2001;10:2023–43.
Fuller RA, Irvine KN, Devine-Wright P, Warren PH, Gaston KJ. Psychological benefits of greenspace increase with biodiversity. Biol Lett. 2007;3:390–4.
Garaffa PI, Filloy J, Bellocq MI. Bird community responses along urban-rural gradients: does town size matter? Landsc Urban Plan. 2009;90:33–41.
Garden J, Mcalpine C, Peterson ANN, Jones D, Possingham H. Review of the ecology of Australian urban fauna: a focus on spatially explicit processes. Austral Ecol. 2006;31:126–48.
Grimm NB, Faeth SH, Golubiewski NE, Redman CL, Wu J, Bai X, Briggs JM. Global change and the ecology of cities. Science. 2008;319:756–60.
Gurevitch J, Hedges LV. Statistical issues in ecological meta-analyses. Ecology. 1999;80:1142–9.
Guilhaumon F, Mouillot D, Gimenez O. mmSAR: an R-package for multimodel species–area relationship inference. Ecography. 2010;33:420–4.
Hanski I, Zurita GA, Bellocq MI, Rybicki J. The species-fragmented area relationship. P Natl Acad Sci USA. 2013;110:12715–20.
He F, Legendre P. On species-area relations. Am Nat. 1995;148:719–37.
Hedges L, Olkin I. Statistical models for meta-analysis. New York: Academic Press; 1985.
Hilty SL. Birds of Venezuela. New Jersey: Princeton University Press; 2002.
Hilty SL, Brown WL, Brown B. A guide to the birds of Colombia. New Jersey: Princeton University Press; 1986.
Hopewell S, McDonald S, Clarke M, Egger M. Grey literature in meta-analyses of randomized trials of health care interventions. Cochrane Database Syst Rev. 2007. https://doi.org/10.1002/14651858.MR000010.pub3.
Hubbell SP. The unified neutral theory of biodiversity and biogeography. California: Princeton University Press; 2001.
Hume R. Complete birds of Britain and Europe. London: Dorling Kindersley; 2002.
Husté A, Boulinier T. Determinants of bird community composition on patches in the suburbs of Paris, France. Biol Conserv. 2011;144:243–52.
Jokimäki J, Huhta E. Artificial nest predation and abundance of birds along an urban gradient. Condor. 2000;102:838–47.
Kazmierczak K, van Perlo B. A field guide to the birds of India, Sri Lanka, Pakistan, Nepal, Bhutan, Bangladesh, and the Maldives. New Delhi: Om Book Service; 2000.
La Sorte FA, Lepczyk CA, Aronson MF, Goddard MA, Hedblom M, Katti M, MacGregor-Fors I, Mörtberg U, Nilon CH, Warren PS, Williams NS. The phylogenetic and functional diversity of regional breeding bird assemblages is reduced and constricted through urbanization. Divers Distrib. 2018;24:928–38.
Leveau LM, Leveau CM. Does urbanization affect the seasonal dynamics of bird communities in urban parks? Urban Ecosyst. 2016;19:631–47.
Lizée MH, Mauffrey JF, Tatoni T, Deschamps-Cottin M. Monitoring urban environments on the basis of biological traits. Ecol Indicat. 2011;11:353–361.
MacArthur RH, Wilson EO. The theory of island biogeography. Monographs in Population Biology, vol. 1. New Jersey: Princeton University Press; 1967.
MacGregor-Fors I, Ortega-Álvarez R. Fading from the forest: bird community shifts related to urban park site-specific and landscape traits. Urban For Urban Green. 2011;10:239–46.
MacGregor-Fors I, Morales-Pérez L, Schondube JE. Migrating to the city: responses of neotropical migrant bird communities to urbanization. Condor. 2010;112:711–7.
Magle SB, Hunt VM, Vernon M, Crooks KR. Urban wildlife research: past, present, and future. Biol Conserv. 2012;155:23–32.
Matthews TJ. Analysing and modelling the impact of habitat fragmentation on species diversity: a macroecological perspective. Front Biogeogr. 2015;7:60–8.
Matthews TJ, Guilhaumon F, Triantis KA, Borregaard MK, Whittaker RJ. On the form of species–area relationships in habitat islands and true islands. Global Ecol Biogeogr. 2016a;25:847–58.
Matthews TJ, Triantis KA, Rigal F, Borregaard MK, Guilhaumon F, Whittaker RJ. Island species–area relationships and species accumulation curves are not equivalent: an analysis of habitat island datasets. Global Ecol Biogeogr. 2016b;25:607–18.
Matthews TJ, Triantis K, Whittaker RJ, Guilhaumon F. sars: an R package for fitting, evaluating and comparing species–area relationship models. Ecography. 2019;42:1446–55.
McKinney ML. Urbanization as a major cause of biotic homogenization. Biol Conserv. 2006;127:247–60.
Miller JR. Hobbs RJ Conservation where people live and work. Conserv Biol. 2002;16:330–7.
Mitchell MH. Observations on birds of southeastern Brazil. Toronto: University of Toronto Press; 1957.
Møller AP, Diaz M, Flensted-Jensen E, Grim T, Ibáñez-Álamo JD, Jokimäki J, Mänd R, Markó G, Tryjanowski P. High urban population density of birds reflects their timing of urbanization. Oecologia. 2012;170:867–75.
Munyenyembe F, Harris J, Hone J, Nix H. Determinants of bird populations in an urban area. Aust J Ecol. 1989;14:549–57.
Murgui E. Effects of seasonality on the species–area relationship: a case study with birds in urban parks. Global Ecol Biogeogr. 2007;16:319–29.
National Geographic Society (US). Field guide to the birds of North America. New York: National Geographic Society; 1999.
Natuhara Y, Imai C. Prediction of species richness of breeding birds by landscape-level factors of urban woods in Osaka Prefecture, Japan. Biodivers Conserv. 1999;8:239–53.
Nielsen AB, van den Bosch M, Maruthaveeran S, van den Bosch CK. Species richness in urban parks and its drivers: a review of empirical evidence. Urban Ecosyst. 2014;17:305–27.
Olson DM, Dinerstein E, Wikramanayake ED, Burgess ND, Powell GVN, Underwood EC, D’amico JA, Itoua I, Strand HE, Morrison JC, Loucks CJ, Allnutt TF, Ricketts TH, Kura Y, Lamoreux JF, Wettengel WW, Hedao P, Kassem KR. Terrestrial ecoregions of the world: a new map of life on earth. Bioscience. 2001;51:933–8.
Ortega-Álvarez R, MacGregor-Fors I. Dusting-off the file: a review of knowledge on urban ornithology in Latin America. Landsc Urban Plan. 2011;101:1–10.
Paradis E, Baillie SR, Sutherland WJ, Gregory RD. Patterns of natal and breeding dispersal in birds. J Anim Ecol. 1998;67:518–36.
Park CR, Lee WS. Relationship between species composition and area in breeding birds of urban woods in Seoul, Korea. Landsc Urban Plan. 2000;51:29–36.
Pautasso M, Böhning-Gaese K, Clergeau P, et al. Global macroecology of bird assemblages in urbanized and semi-natural ecosystems. Global Ecol Biogeogr. 2011;20:426–36.
Peterson R, Mountfort G, Hollom PAD, Díaz G. Guía de campo de las aves de España y demás países de Europa. Barcelona: Omega; 1973.
Preston FW. The canonical distribution of commonness and rarity: part I. Ecology. 1962;43:185–215.
Rahbek C. The relationship among area, elevation, and regional species richness in neotropical birds. Am Nat. 1997;149:875–902.
R Core Team. R: A language and environment for statistical computing. R Foundation for Statistical Computing, Vienna, Austria, 2019. http://www.R-project.org/.
Rosenberg MS. The file-drawer problem revisited: a general weighted method for calculating fail-safe numbers in meta-analysis. Evolution. 2005;59:464–8.
Rosenberg MS, Adams DC, Gurevitch J. MetaWin: statistical software for meta-analysis. Sunderland: Sinauer Associates; 2000.
Rosenthal R. The file drawer problem and tolerance for null results. Psychol Bull. 1979;86:638.
Rosenzweig ML. Species diversity in space and time. Cambridge: Cambridge University Press; 1985.
Scheiner SM, Chiarucci A, Fox GA, Helmus MR, McGlinn DJ, Willig MR. The underpinnings of the relationship of species richness with space and time. Ecol Monogr. 2011;81:195–213.
Seto KC, Fragkias M, Güneralp B, Reilly MK. A meta-analysis of global urban land expansion. PLoS ONE. 2011;6:e23777.
Shochat E, Warren PS, Faeth SH, McIntyre NE, Hope D. From patterns to emerging processes in mechanistic urban ecology. Trends Ecol Evol. 2006;21:186–91.
Stott I, Soga M, Inger R, Gaston KJ. Land sparing is crucial for urban ecosystem services. Front Ecol Environ. 2015;13:387–93.
Szlavecz K, Warren P, Pickett S. Biodiversity on the urban landscape. In: Concotta RP, Gorenflo LJ, editors. Human populations, its influences on biological diversity. Ecological studies, vol. 214. Berlin: Springer-Verlag; 2011.
Sukhdev P. Foreword. In: Elmqvist T, Fragkias M, Goodness J, Güneralp B, Marcotullio PJ, McDonald RI, Parnell S, Schewenius M, Sendstad M, Seto KC, Wilkinson C, editors. Urbanization, biodiversity and ecosystems services: Challenges and opportunities. Dordrecht: Springer; 2013.
Sutherland GD, Harestad AS, Price K, Lertzman KP. Scaling of natal dispersal distances in terrestrial birds and mammals. Conserv ecol. 2000. https://doi.org/10.5751/ES-00184-040116.
Tjørve E. Shapes and functions of species–area curves: a review of possible models. J Biogeogr. 2003;30:827–35.
Tjørve E. Shapes and functions of species–area curves (II): a review of new models and parameterizations. J Biogeogr. 2009;36:1435–45.
Tjørve E, Turner WR. The importance of samples and isolates for species–area relationships. Ecography. 2009;32:391–400.
Triantis KA, Guilhaumon F, Whittaker RJ. The island species–area relationship: biology and statistics. J Biogeogr. 2012;39:215–31.
Tummers B. Data Thief III (v. 1.1). 2006. http://www.datathief.org/. Accessed 25 Mar 2017.
United Nations. World urbanization prospects: the 2014 revision. Highlights (ST/ESA/SER.A/352). 2014.
Urquiza A, Mella JE. Riqueza y diversidad de aves en parques de Santiago durante el período estival. Boletín Chileno de Ornitología. 2002;9:12–21.
Vaccaro AS, Filloy J, Bellocq MI. What land use better preserves taxonomic and functional diversity of birds in a grassland biome? Avian Conserv Ecol. 2019;14:1.
Watling JI, Donnelly MA. Fragments as islands: a synthesis of faunal responses to habitat patchiness. Conserv Biol. 2006;20:1016–25.
Wild Bird Society of Japan. A field guide to the birds of Japan. Tokyo: Kodansha International Limited; 1982.
Yamashina Y. Birds in Japan: a field guide. Tokyo: Tokyo news Limited; 1961. p. 1961.
Zhou D, Chu LM. How would size, age, human disturbance, and vegetation structure affect bird communities of urban parks in different seasons? J Ornithol. 2012;153:1101–12.
We greatly appreciate the help with QGis software provided by Carolina Ramos. LML thanks the support of Lucia Gonzalez Salinas. We are very thankful to Carolina Berget, Andrés Seijas, Karen Ikin, Michael O´Neal Campbell and Pilar Carbó-Ramírez, who generously gave data for the analysis. The comments of Simone Fattorini and three anonymous reviewers greatly improved earlier versions of the manuscript. The research was funded by the Consejo Nacional de Investigaciones Científicas y Técnicas and the Universidad de Buenos Aires (Argentina). The authors dedicate this piece of work to our co-author M.I. Bellocq, dearest colleague and friend who will always be in our hearts and memories.
There was no any funding to carry out this research work.
Ethics approval and consent to participate
Consent for publication
The authors declare that they have no competing interests.
M. Isabel Bellocq—Deceased on 9 July 2019
List of studies analysed with their environmental and bird trait information and fit and slope values of each study.
Forest plot for the fit (a) and the slope value (b) of the SAR in urban green areas.
List of papers used in the meta-analysis.
About this article
Cite this article
Leveau, L.M., Ruggiero, A., Matthews, T.J. et al. A global consistent positive effect of urban green area size on bird richness. Avian Res 10, 30 (2019). https://doi.org/10.1186/s40657-019-0168-3
- Conservation macroecology
- Habitat islands
- Species-area relationships
- Species-urban green area relationship
- Urban parks