Study Area

Fieldwork was carried out in Merida (approximately 20.9° N, 89.6° W, 15 m a.s.l.), the main city of the Yucatan Peninsula (YP), southeastern Mexico with more than 1.1 million people (Consejo Nacional de Población, 2015). YP is one of the Mexican biogeographic regions with highest levels of bird species richness (Navarro-Sigüenza et al., 2014), and it holds considerable importance for wintering and transient Nearctic-Neotropical species (Calmé et al., 2015). The area was originally covered by seasonally dry tropical forest characterized by a dry season that may last between 7 and 8 months (Torrescano-Valle & Folan, 2015). Currently, more than 25% of the native species that conform the urban flora belong to the Fabaceae, Euphorbiaceae, and Poaceae families, while common exotic species include Flamboyant (Delonix regia), Golden rain tree (Cassia fistula), and Indian almond (Terminalia catappa) (Peraza-Contreras, 2011).

A total of 22 green spaces were surveyed within Merida Municipality limits (Figure 1). Herein, we use green spaces to refer to urban open spaces dominated by trees and shrubs that are used by humans for several purposes such as recreation, exercise, education, or others, and to which access can be unrestricted or restricted. For this work, we surveyed botanical gardens, public parks, an archaeological site, a zoo, a sport complex, and a reforested area within an industrial plant (Figure 2). Green space size ranged from 0.5 to 39.1 ha and distance to the nearest native vegetation patch ranged from 101 to 5685 m. Some of the surveyed sites harbor artificial waterbodies and were included in the sample to acknowledge the importance of these habitats for local bird diversity given that superficial waterbodies are scarce in YP (Torrescano-Valle & Folan, 2015).

Figure 1.

Left: Study Area (Represented as the Red Dot on the Mexico Inset Map) and Overview of the 22 Surveyed Green Spaces Within Merida Municipality, Yucatan, Mexico. Right: Zoom to the Shadowed Area on Left.

Figure 2.

Surveyed Green Spaces in the City of Merida Encompass a Wide Variety of Habitats.

Bird Surveys

Bird surveys were conducted during the fall of 2016 from September 27 to November 15. Birds were surveyed using 5-min fixed-width point counts (25 m) separated by at least 150 m. We chose a 25-m radius to increase the probability that all individuals would be detected in all the surveyed habitats (Hutto et al., 1986). The number of sampling points ranged from one to eight and corresponded to green space size. Counts were made during the first 4 hr after sunrise under suitable weather conditions (Ralph et al., 1996). Each point was visited thrice and bird data were collected by a single observer (R. N. D.). All detected birds were included in the avifauna description except Chaetura vauxi and Stelgidopteryx serripennis, since these species were flying over the plots and hence were unlikely to be using the habitat within the plot (Gates, 1997). Statistical analyses were performed with landbirds only since we did not measure the main habitat features that influence the distribution of aquatic birds (Rosa et al., 2003). Raptors were not included in statistical analysis because the count method is not suitable for estimating their numbers (Fuller & Mosher, 1981). For each sampling point, we pooled data from all three visits to get cumulative lists of detected species and generate a species presence/absence table. To determine whether our survey effort was enough to provide a representative sample of the bird community in the time surveyed, we computed the nonparametric incidence-based estimator Jacknife 1 (González-Oreja et al., 2010) using EstimateS ver. 9.1.0 (Colwell, 2013).

Habitat Characterization

We measured nine variables in the field to evaluate the habitat using 25-m-radius circular plots centered on each bird sampling point (Table 1). We considered five classes of environmental variables that could influence the distribution of birds in green spaces. Vegetation composition was evaluated by native tree species richness, exotic tree species richness, native shrub species richness, and exotic shrub species richness. Vegetation structure was assessed by maximum tree height, maximum basal area, maximum bush height, and maximum bush basal area. To quantify the extent of the urban infrastructure in green spaces, we counted the number of poles.

Table 1.

Descriptive Statistics of Environmental Variables Recorded in Green Spaces of Merida, Yucatan.

Variables’ names, abbreviation, description, and the corresponding predictor set are shown. SE = standard error.

Vegetation coverage within each green space was estimated using the Soil Adjusted Vegetation Index (SAVI), which is a modified version of the Normalized Vegetation Index (NDVI) (Huete, 1988). NDVI and SAVI are strongly correlated to several vegetation parameters including vegetation density and percent green vegetation cover (Huete, 1988; Purevdorj et al., 1998). In this study, we employed SAVI because it minimizes errors due to soil substrate optical properties (Huete, 1988). Using a Copernicus Sentinel-2 satellite image and an open source geographic information system (QGIS Development Team, 2020), we computed SAVI using the formula:

NIR = near-infrared band;

red = red band;

L = 0.5.

The quality of images acquired during the bird sampling season was low due to cloud cover so we used a satellite image acquired on January 25, 2017 (cloud cover percentage = 0.0%). The resulting SAVI layer had a 10-m spatial resolution. SAVI values were classified in 12 classes ranging from one (built environment) to 11 (very dense vegetation) and zero values represent water. High-resolution Google Maps images were used to assign SAVI classes that represent vegetation, considering from very scattered vegetation to highly dense vegetation. Once pixels were classified in 1 of these 12 classes, we obtained the total number of pixels for each class. Then, total counts were used to estimate vegetation cover.

Green space area can affect the probability of occupation of species in different fashion (Roberts & King, 2017), so we estimated green space size to introduce it in the models. Similarly, isolation of green spaces can influence the composition structure of urban bird communities (Charre et al., 2013; Fernandez-Juricic, 2002). In this study, we used two alternative approaches to quantify green space isolation. For the first one, we calculated the Euclidean distance from each green space to the closest continuous native vegetation patch. The second approach considered the fact that vegetation cover is not homogeneously distributed through the city and that vegetation cover adjacent to green spaces can influence bird richness and abundance (Shanahan et al., 2011). Therefore, we calculated the extent of vegetation cover in a 100-m width buffer around each green space. Estimations of vegetation cover adjacent to green spaces did not include estimations of vegetation cover within each green space. For this purpose, we employed the same SAVI classes used to identify vegetation cover within green spaces, but in the case of vegetation cover adjacent to green spaces, vegetation-related SAVI classes were grouped in two broad categories: scattered and dense vegetation.

Bird Trait Data

All species in our database were characterized based on three functional traits: diet, foraging strata, and body mass. Diet was expressed as the percentage use of each of the 10 food items categories considered (Table 2). Percentages of diet composition for each species sum to 100. Foraging strata trait was expressed as the estimated percentage use of each one of the seven strata considered. The values of the seven strata sum to 100. Functional traits were sourced from Wilman et al. (2014) and del Hoyo et al. (2018). In addition, residence status in the area was determined based on digital species distribution maps (BirdLife International, 2018) and considering four categories: residents, wintering, transient, and wintering/transient species.

Table 2.

Bird Traits Used in This Study.

Data Analyses

Here, we investigate the effect of potentially influential environmental factors on the functional composition of avian assemblages of urban green spaces through RLQ and fourth-corner tests, which allow to analyze trait-environment relationships (Dray et al., 2014). RLQ is a three-table ordination aimed to identify the main co-structures between an environmental table (R) and a trait table (Q) with the constriction of a species table (L). On the other hand, the fourth-corner approach quantifies and tests the significance of bivariate associations between traits and environmental variables (Dray et al., 2014). RLQ combines three separate ordinations which summarize the main structures of each table. In this way, RLQ relate species traits and environmental variables considering a sites-by-species table (ter Braak et al., 2012). In this study, we employed a binary species table (presence or absence). For the independent ordinations, we followed Borcard et al. (2018) to choose the ordination method based on the type of the variables and to assign rows and columns weights. To test the significance of the association between the environmental and trait tables, several permutation models have been proposed (Thioulouse et al., 2018; see details of permutation models in Borcard et al., 2018). We used a single global test that consists of two independent models whose null hypotheses are species compositions in the sites are not related to environmental conditions of the sites (Model 2) and species distribute according to their environmental preferences but irrespective of their traits (Model 4) (Borcard et al., 2018). The maximum p value of both permutation tests becomes the overall p value to attain a correct Type 1 error (ter Braak et al., 2012).

Measured environmental variables describe broad habitat features, so we grouped them in distinctive sets: vegetation composition (4), vegetation structure (4), vegetation cover (6), and urban infrastructure (1). The fifth set, patch extra descriptors (4) includes green space size and isolation measures. We performed an RLQ analysis including all the predictor sets (that total nineteen variables). In addition, we performed alternative RLQ analyses excluding one of the sets in each ordination since some variables can introduce noise and reduce significance. These RLQ ordinations with reduced environmental tables were obtained to explore the contribution of broad habitat characteristics in determining the functional composition of these bird assemblages. The sum of the correlation L metric of the first two RLQ axes was used as an indicator of the goodness-of-fit of the RLQ ordination (Bernhardt-Romermann et al., 2008). This metric serves to compare the correlation between the trait-based species scores and the environmental-based site scores generated by the RLQ ordination and the correlation of the sites and species scores of the separate ordination of the species table (Thioulouse et al., 2018). The ordinations with better fit are reported and used in the species grouping (further details afterwards), and displayed in the plots.

The fourth-corner analysis tests the relationships between species traits and environmental variables, one at each time (Thioulouse et al., 2018). In this study, we employed jointly fourth-corner analysis and output RLQ axes, which can be interpreted as either environmental gradients or trait syndromes (Dray et al., 2014). Since multiple tests are performed, p values need to be corrected. We set to 9,999 as the number of permutations and used the false discovery rate method to adjust p values in order to avoid Type 1 error (see details in Thioulouse et al., 2018). We decided to report results for which p value < .10 to increase the power of the test given the small sample size (Zar, 2014) and in order to detect likely associations.

To distinguish groups of species that share traits and respond in similar ways to environmental characteristics, we grouped species based on their resulting RLQ scores using Ward’s hierarchical clustering. To determine the optimal number of groups, we considered Calinski–Harabasz index (Borcard et al., 2018). Only those species that occurred in three or more green spaces were included in the analyses to reduce the disproportionately large effects of rare species (Legendre & Gallagher, 2001). We investigated patterns of response of two avian categories: resident and wintering species, so we performed RLQ analysis separately for these subsets of species. The experimental unit in all the multivariate analyses was the green space; hence, data from the sampling points were pooled for each green space. All statistical analyses were performed with R (R Core Team, 2018), applying functions from vegan package (Oksanen et al., 2017), and ade4 package (Dray & Dufour, 2007) and found in Kleyer et al. (2012).

Avifauna of Green Spaces in Merida

A total of 87 species from 16 orders and 32 families were detected in green spaces of Merida (Online Appendix 1). Species richness per green space ranged from 9 to 43. Results of the richness estimation indicate that our sample captured 78.4% of the species present in the surveyed space and time. There was a clear predominance of resident birds over migratory ones (Table 3). The most abundant species were Great-tailed Grackle (Quiscalus mexicanus) and White-winged Dove (Zenaida asiatica). Six species were widespread occurring in more than 80% of the green spaces, while 34 species such as Masked Tityra (Tityra semifasciata) or Vermilion Flycatcher (Pyrocephalus rubinus) were detected in just one green space. Aquatic birds, such as Blue-winged Teal (Spatula discors) or Little Blue Heron (Egretta caerulea), were present in only three green spaces and comprised 17 species including migratory ones. We recorded three species endemic to the YP (Calmé et al., 2015): Yucatan Woodpecker (Melanerpes pygmaeus), Yucatan Jay (Cyanocorax yucatanicus), and Orange Oriole (Icterus auratus) and two exotic species in the area: Rock Pigeon (Columba livia) and Eurasian Collared-Dove (Streptopelia decaocto).

Table 3.

Species Richness, Abundance, and Estimated Biomass of Species of Different Categories Recorded in Green Spaces of Merida, Yucatan.

RLQ Ordination for Birds of Urban Green Spaces

The correlation L score indicated that the optimal environmental set for resident and wintering species were similar (Table 4). The environmental table used in the RLQ analysis of resident birds included the vegetation composition, vegetation structure, urban infrastructure, and patch extra descriptor sets. For wintering species, the environmental table consisted of the same descriptor sets except for urban infrastructure. The first two axes accounted for most of the variability explained by the separate ordinations (92% and 81% for resident and wintering species, respectively), so the covariance between environmental variables and species traits was well described by RLQ analysis. Considering the first RLQ axis, variability was better captured for the environmental table (95% for resident birds, 92% for wintering birds) than for the trait table (63% for resident birds, 52% for wintering birds).

Table 4.

Summary of RLQ Ordinations for Resident and Wintering Species With Different Sets of Explanatory Variables.

Ranking of RLQ ordinations was based on correlation L metric. Only the five best ranked ordinations for each group of species are listed. Variance captured by separate ordinations represents the maximum value, to which variance captured by corresponding RLQ axes is compared (coinertia R and coinertia Q). Percentage of variability captured by RLQ Axes 1 and 2 is shown together with p values obtained for each ordination. Abbreviation for sets of variables are (numbers in brackets correspond to the number of variables in the set): comp, vegetation composition (4); stru, vegetation structure (4); cove, vegetation cover (6); infr, urban infrastructure (1); patc, patch extra descriptors (4).

Drivers of Functional Composition: Resident Birds

RLQ Axis 1 extracted 87.5% of the co-inertia, and it defined a gradient of green space isolation, driven mainly by decreases in the amount of dense vegetation and increases in the distance to native vegetation patches (Figure 3A). The number of plant species, especially native species of both trees and shrubs, proved to be another relevant driver along RLQ axis, and it showed an opposite association with isolation (Figure 3B). The ordinations indicated that large-bodied species were common in sites far from native vegetation patches and scarcerlysurrounded by dense vegetation, while small-bodied species were present in green spaces that occupy an intermediate position along this axis. Furthermore, canopy-forager species were present in green spaces surrounded by relatively large amounts of dense vegetation and where native flora predominated (Figure 3C). This includes species such as Scrub Euphonia (Euphonia affinis) or Green Jay (Cyanocorax yncas) (Figure 3D). On the contrary, ground-foraging species whose diet is dominated by seed and plant material became more common in sites of reduced native flora richness. RLQ Axis 2 separated green spaces based on the number of exotic tree species, but it only extracted 5.1% of the co-inertia, so we do not consider it for further discussion.

Figure 3.

RLQ Analysis Results for Resident and Wintering Bird Species (Left and Right Columns, Respectively). Plots show the ordination of (A) surveyed sites, (B) environmental descriptors, (C) species’ traits, and (D) species. Environmental descriptors and traits significantly associated to Axis 1 are shown with purple, while marginally significant associations are shown with orange labels. Bird species are grouped based on the output RLQ scores, using Ward’s hierarchical clustering. Codes for environmental variables are shown in Table 1.

Global testing showed that the link between traits and environment was significant for resident birds (Model 2 p value: 0.0001; Model 4 p value: 0.03), and this finding was supported by fourth-corner tests on RLQ axes. The joint approach of RLQ analysis and fourth-corner tests indicated that numerous environmental factors of green spaces were related to the functional composition of bird communities (Figure 3B). In particular, the amount of dense vegetation adjacent to green space showed the highest influence on trait composition. Regarding species’ traits, the association of body mass and aerial-foraging strategy with Axis 1 was marginally significant (Figure 3C). There was no evidence for significant associations for RLQ Axis 2. Results from the joint approach of RLQ and fourth-corner analysis are summarized in the ordinations plots (Figure 3B and C).

Drivers of Functional Composition: Wintering Birds

RLQ Axis 1 (69.3% of co-inertia) separated less isolated green spaces harboring a large number of shrub species, both native and exotic ones, from those green spaces with impoverished shrub richness, but with an elevated number of exotic tree species and far from native vegetation patches (Figure 3A). Body size was related to this axis together with the aerial-foraging strategy. Other traits, such as understorey-foraging and midhigh-foraging were not clearly associated with any particular characteristic of green spaces.

Global RLQ testing did not support a significant association between traits and environment (Model 2 p value: 0.12; Model 4 p value: 0.76). However, fourth-corner tests indicated few marginally significant associations (Figure 3B and C): dense vegetation extent and native shrubs richness were the characteristics of green spaces associated with trait composition, whereas body mass and the aerial-foraging strategy revealed an association with the environmental gradient defined by RLQ Axis 1.

Classification of Species

Differences between subsequent values of Calinksi–Harabasz criterion suggested the clustering of resident species in three groups. The ordination plot showed that Groups 2 and 3 clearly occupied opposite extremes of the environmental gradient extracted by Axis 1. Group 1 was the most numerous and included small-bodied species. This group comprised species that feed predominantly on the understorey and midhigh vegetation strata with a diet consisting mainly of invertebrates. On the contrary, Group 2 included large-bodied birds that feed on invertebrates and seeds. Species belonging to this group tend to forage on the ground and use to a much lesser extent other strata. The third group includes medium-sized species with a more diversified use of food items and foraging strata.

Regarding wintering birds, species were clustered in three groups. Species classified in Groups 1 and 3 were small-bodied birds. Although species in Group 1 were mainly insectivores and used the understorey and midhigh strata to forage, species belonging to Group 3 were characterized by the use of the ground, understorey, and canopy strata to forage. Group 2 comprised Summer Tanager (Piranga rubra) only, a large-bodied species with a wider use of foraging strata and a predominantly insectivore diet.