Consistency in bird community assembly over medium-term along rural-urban gradients in Argentina

The analysis of bird community assembly rules is fundamental to understand which mechanisms determine the composition of bird species in urban areas. However, the long-term variation of community assembly rules has not been analyzed yet. The objectives of this study are (1) to analyze the variation of community assembly rules along rural-urban gradients of three cities in central Argentina and (2) to compare the patterns of community assembly between two periods separated by 6 years. Bird surveys were performed along transects in urban, suburban, and rural habitats during 2011 and 2017. Departures from null models that took into account differences in species richness (standardized effect size, SES) were calculated for functional and phylogenetic diversities. A total of 57 species were recorded. Bird species richness was higher in suburban than in urban and rural habitats. SES of functional diversity increased over the years and was significantly lower in urban habitats than in rural habitats, showing a pattern of functional clustering in the most urbanized areas and functional randomness in rural ones. Phylogenetic diversity was higher in both suburban and urban habitats than rural ones, and the phylogenetic clustering in rural bird assemblages changed to randomness in suburban and urban habitats. Bird communities in urban habitats were phylogenetically random and functionally clustered, evidencing environmental filtering by urbanization. In contrast, bird communities in rural areas tended to be phylogenetically clustered, evidencing that certain clades are adapted to rural areas. The processes structuring bird communities along rural-urban gradients were consistent between the 2 years compared.


Introduction
Urban areas are expanding continuously, causing declines of bird diversity mainly through habitat loss and fragmentation, pollution, and altered resource flows (Grimm et al. 2008;Aronson et al. 2014;Shanahan et al. 2014;Beninde et al. 2015). This diversity decline is reflected in the taxonomic, functional, and phylogenetic diversity (Devictor et al. 2010;Morelli et al. 2017;Liu et al. 2019;Leveau and Leveau 2020). Taxonomic diversity is characterized by the number of species in a site (Cadotte and Tucker 2018). Functional diversity reflects the range and variety of resources that species obtain (Díaz and Cabido 2001), whereas phylogenetic diversity is related to the phylogenetic relationships among species . Birds in urban areas are mainly habitat and food generalists, so urban bird communities show reduced functional diversity (Schütz and Schulze 2015;Oliveira Hagen et al. 2017;Palacio 2020). On the other hand, a reduced phylogenetic diversity in urban areas has been found due to generalist species being closely related to each other Sol et al. 2017;La Sorte et al. 2018).
The analysis of community assembly rules is fundamental to understand urban bird communities because it allows predicting how species from the regional pool are assembled in urban areas (Magura et al. 2018). Two aspects during the last decade have allowed to analyze community assembly rules along urbanization gradients (Moreno-Contreras et al. 2019, Weideman et al., 2020: (1) the availability of information regarding phylogenetic relationships among coexisting bird species and their functional characteristics , Wilman et al. 2014 and (2) the possibility to compare the observed phylogenetic and functional diversity against null models that take into account differences in species richness (Swenson 2014).
The contrast of observed functional diversity values with those obtained with null models allows suggesting departures from random, or stochastic, processes to deterministic processes molding bird communities, such as environmental filtering or competition (Petchey et al. 2007). Observed functional diversity lower than expected by chance, hereafter functional clustering, suggests a mechanism of environmental filtering from the regional pool, allowing the coexistence of species with similar traits (Petchey et al. 2007). A higher observed diversity than expected by chance, hereafter functional overdispersion, suggests species competition because coexisting species are functionally more different than expected, thus avoiding competition for resources (Petchey et al. 2007). On the other hand, the phylogenetic diversity departures from null models show patterns of phylogenetic clustering or overdispersion (Webb et al. 2002), where coexisting species are closely related or more distant than expected, respectively.
The integration of functional and phylogenetic departures from null models results in four scenarios (Webb et al. 2002;Emerson and Gillespie 2008): (1) functional and phylogenetic clustering of communities, which suggests that similar traits among coexisting species are phylogenetically conserved; (2) functional clustering with phylogenetic overdispersion or randomness of communities, which suggests that similar traits of coexisting species resulted from convergent evolution; (3) functional overdispersion with phylogenetic clustering or randomness of communities, which suggests that different traits resulted from convergent evolution among coexisting species; and (4) functional and phylogenetic overdispersion, which suggest that different traits among coexisting species resulted from conserved traits.
Analyses of bird community assembly rules along urbanization gradients are scarce and have shown that urban bird communities are phylogenetically and functionally random (Weideman et al., 2020), or functionally clustered with phylogenetic randomness, suggesting a mechanism of environmental filtering (Moreno-Contreras et al. 2019). On the other hand, the temporal variation of community assembly rules along urbanization gradients has seldom been assessed.
Interannual variations in climatic conditions that affect resource availability and habitat structure for birds and species colonization or extinctions can affect bird species richness and community assembly patterns (Li et al. 2015, Rota et al. 2017, Weyland et al. 2019. However, bird communities in urban areas are expected to be more stable than in rural areas (Suhonen et al. 2009, Leveau et al. 2015, presumably because of the less variable habitat structure and food availability in urban areas (Shochat et al. 2006, Leveau 2018. This is underlined by the fact that the effects of hurricanes, tornados, and fire on vegetation structure are reduced or suppressed in urban areas (Argañaraz et al. 2015;Kingfield and De Beurs 2017;Pickens et al., 2017). Activities such as irrigation, maintaining vegetation, or using fertilizers and pesticides may buffer the interannual variation of primary productivity in urban areas (Shochat et al. 2004;Leong and Roderick 2015) and the interannual variation of habitat structure in suburban ones (Guthrie 1974, Lepczyk et al. 2004. The aims of this study were (1) to analyze the variation of community assembly rules along rural-urban gradients of three cities in central Argentina and (2) to compare the patterns of community assembly between two periods separated by 6 years. A functional clustering with phylogenetic randomness was expected in bird communities of urban areas of both years. However, long-term changes in community assembly are expected in rural areas as a response to climatic and habitat structure changes, species extinction, and colonization.

Study area
The study was carried out in three cities of the Pampas Region, central Argentina ( . The three cities were in the same region, separated by a maximum distance of 59 km; therefore, the influence of latitude or climate was negligible. The landscape surrounding the cities is composed mainly by grazing land, cropland, and exotic tree plantations. The climate is temperate, with mean annual precipitation of 923.6 mm and mean annual temperature of 14°C (data from the National Meteorological Service, www.smn.org.ar).

Bird surveys
Three habitat types were selected for bird surveys: (1) rural, composed of crops and pastures; (2) suburban, composed of detached houses with yards; and (3) urban dominated by impervious surfaces (see Leveau, 2019 for details). Five 100 × 50-m transects separated by at least 200 m were placed in each habitat (see also Rajashekara and Venkatesha 2015). Transects were visited once during each December 2011 and 2017. Surveys were carried out in the first 4 h after dawn on days without rain or strong winds, and all birds seen or heard were counted, except those flying over the top of buildings or trees, or below that height but showing no feeding activity.

Sample coverage
Because transects visited once each year may not be enough to detect most of the species present, sampling completeness was estimated for each transect each year. The sample coverage was estimated with the function iNEXT of the iNEXT package (Hsieh et al. 2016) and is the proportion of the total number of individuals that belong to the species detected in the transect. The sample coverage was used as an explanatory variable in models.

Environmental and climatic changes between periods
In Mar del Plata, 2017 had a higher annual mean temperature than 2011 (15.17°C versus 13.92°C, respectively), and the annual precipitation was also higher (1405.8 mm in 2017 vs. 1103 mm in 2011). Moreover, whereas the habitat structure of urban and suburban habitats remained stable, most of the transects in the rural habitat had a replacement of crops by abandoned fields.

Taxonomic, functional, and phylogenetic diversity
The taxonomic alpha diversity of each transect during each year was determined as the bird species richness, which was the total number of species in each visit. The functional alpha diversity for each transect during each year was measured using a matrix of bird traits as columns, and bird species as rows was constructed. Functional traits were related to the way species obtain resources and were obtained from the EltonTraits 1.0 database (Wilman et al. 2014), which provides information about diet, foraging strata, and body mass (g) (see Supplementary Table S1). The database contains information of percent usage of each diet and foraging strata category. Additionally, information about clutch size and residency status was obtained from the Handbook of the Birds of the World (HBW) online (https://www.hbw. com/). With the matrix of species and traits, we calculated the species dissimilarities using the Gower distance because our data traits had continuous as well as binary data, using the function gowdis of the FD package in R (Laliberté et al. 2014, R Core Team 2017. Because functional diversity is generally positively related to species richness Gaston 2002, Weideman et al. 2020), the function ses.mpd of the package picante was used (Kembel et al. 2010). This function calculates the standardized effect size (SES) of mean pairwise distances in communities (hereafter, SES of functional diversity), taking into account differences in species richness with the null model "richness" using 999 randomizations. SES of functional diversity was based on species presence/absence data.
The phylogenetic alpha diversity was measured for each transect during each year, using a dated phylogeny of all species in this study from the BirdTree database  (Ericsson et al. 2006). A 50% majority-rule consensus tree was constructed using the software TreeAnnotator (Rambaut andDrummond 2002-2015). The tree was used in R with the functions read.nexus and as.phylo of the ape package (Paradis et al. 2004). This function ses.mpd was used to calculate the standardized effect size (SES) of mean pairwise distances in communities (hereafter, SES of phylogenetic diversity), taking into account differences in species richness with the null model "richness" using 999 randomizations. SES of phylogenetic diversity was based on species presence/absence data.

Statistical analysis
The differences of alpha diversity in relation to sample coverage, habitats, years, and the interaction between habitat and years were analyzed through generalized linear mixed models (GLMM). Sample coverage was a continuous variable, whereas habitat and year were categorical variables. Therefore, the relationships between alpha diversity, habitat, and years were analyzed depending on a fixed value of sample coverage (James et al. 2013). Bird species richness variation was analyzed using a negative binomial distribution of errors due to overdispersion with the glmer.nb function of the lme4 package (Bates et al. 2018). SES of functional and phylogenetic diversity variations were analyzed using a Gaussian distribution of errors with the function lme of the nlme package (Pinheiro et al. 2013). A random structure with transects nested within city was used to control for spatial and temporal pseudoreplication. Models were constructed by backward elimination of non-significant variables (P > 0.05) from the full model using the anova function. Final models were contrasted against null models using a likelihood ratio test (LRT test) (P < 0.05). Differences of means between habitats were explored with Tukey tests using the function glht of the multcomp package (Hothorn et al. 2016). Multicollinearity was explored with the function vif of the car package (Fox et al. 2012), and no multicollinearity was detected between predictor variables. Plots were made with the plot_model function of the sjPlot package (Lüdecke et al. 2018).

Results
A total of 57 species were recorded (Table 1). The House Sparrow (Passer domesticus) and the Rock Dove (Columba livia) were common in urban habitats. The Picazuro Pigeon (Patagioenas picazuro), the Eared Dove (Zenaida auriculata), and the Red-bellied Thrush (Turdus rufiventris) were common in suburban habitats, whereas the Rufous-collared Sparrow (Zonotrichia capensis), the Grassland Yellow-finch (Sicalis luteola), and the Monk Parakeet (Myiopsitta monachus) were common in rural habitats. Most of the species were more common during 2017 (Table 1).
The mean sample coverage of transects was 0.80, ranging from 0.18 to 1.00. Therefore, there was a considerable heterogeneity in the sample coverage between transects. The sample coverage was included in the models explaining the relationship between the response variables and the habitat type and years.
Bird species richness varied between habitats and in relation to sample coverage (LRT = 43.32, P < 0.001; Table 2). Bird species richness was higher in suburban habitats than in urban and rural ones (Tukey tests, P < 0.05) and decreased with increasing sample coverage, suggesting that transects with high species richness were undersampled (Fig. 2a, b).
SES of functional diversity varied in relation to habitats and years (LRT = 21.94, P < 0.001; Table 2). SES of functional diversity was significantly lower in urban habitats than in rural ones (Tukey test, P < 0.05; Fig. 2c), showing a pattern of functional clustering in the most urbanized areas and functional randomness in rural habitats. Moreover, SES of functional diversity increased between the first and the second census, showing functional clustering in 2011 and functional randomness in 2017 (Fig. 2d).
SES of phylogenetic diversity varied between habitats and in relation to sample coverage (LRT = 10.39, P = 0.006; Table 2). SES of phylogenetic diversity was higher in urban and suburban than in rural habitats (Tukey tests, P < 0.05; Table 2; Fig. 2e), showing phylogenetic randomness in urban and suburban habitats and a phylogenetic clustering in rural ones. Moreover, SES of phylogenetic diversity decreased with sample coverage, suggesting that transects with bird communities characterized by high phylogenetic diversity were undersampled (Fig. 2f).

Discussion
Although the mean sample coverage of transects was high, there was a variation among transects that was significantly related to bird species richness and phylogenetic diversity. Therefore, the results obtained must be interpreted with caution.
Community assembly rules varied along the gradients, suggesting environmental filtering and phylogenetic randomness in urban areas. Conversely, rural areas were associated with functional randomness and phylogenetic clustering. Moreover, although functional diversity changed between years, these changes were similar along rural-urban gradients, suggesting consistent community assembly rules. Bird species richness decreased with increasing sample coverage of each transect, showing that those transects with the highest species richness values had the lowest estimated sampling completeness. Therefore, this relationship suggests that transects with the highest species richness were undersampled, and some species were not detected. In spite of this, consistent changes in species richness between habitats were found.
Bird species richness was the highest in suburban areas, and this pattern may be related to a higher habitat complexity and more resources in suburban areas compared to urban and rural ones (Guthrie 1974, Vale and Vale 1976, Blair 1996, Leveau 2019). On the other hand, urban areas were dominated by a few exotic species that nest in buildings and have an omnivorous diet, such as the House Sparrow and the Rock Dove, and this pattern has been documented by several authors worldwide (Blair 1996;Clergeau et al. 1998;Conole and Kirkpatrick 2011;van Rensburg et al. 2009;Silva et al. 2016;Escobar-Ibáñez et al. 2020).  The SES of functional diversity was the lowest in urban areas showing functional clustering. This suggested that urban bird assemblages were molded by environmental filtering, allowing the permanence of bird species with certain functional traits, such as omnivorous diet and resident status. This result agrees with those found by other authors at local scales in the USA (Evans et al., 2018), Germany (Schütz and Schulze 2015), South Africa (Weideman et al. 2020), but see Moreno-Contreras et al. (2019) in Mexico. However, studies conducted at larger spatial scales found the opposite pattern (Oliveira Hagen et al. 2017), suggesting that urban areas can provide a variety of resources for birds, increasing their functional diversity. The contrasting patterns found between studies suggest that the relationship between functional diversity and urbanization varies with the spatial scale that is considered in each study.
The SES of functional diversity was higher in 2017 than in 2011, showing a pattern of functional clustering in 2011 and randomness in 2017. This pattern could be related to more favorable climatic conditions in 2017. This year was warmer and wetter than 2011, possibly Fig. 2 Bird species richness variation between a habitat types and b in relation to sample coverage; SES of functional diversity variation between c habitat types and d years; and SES of phylogenetic diversity variation between e habitat types and f in relation to sample coverage along ruralurban gradients in central Argentina. SES, standardized effect size; values higher than zero indicate a greater diversity than expected by the null model, whereas values lower than zero indicate lower diversity than expected by the null model. Points are means, and vertical lines and grey areas are confident intervals resulting in higher food availability for birds (Weyland et al. 2019). In addition, rural areas had a replacement of crops by abandoned fields in 2017, which probably resulted in increased habitat heterogeneity and resources availability for birds (Orłowski 2005, Herzon et al. 2014). These interannual increases in resources and habitat heterogeneity probably allowed a change from functional clustering to functional randomness where the variety of traits is proportional to the observed species richness.
On the other hand, the increase in functional diversity between years may be related to the colonization by species that are functionally different to the other coexisting species. For example, the Barn Swallow (Hirundo rustica) is a migrant species that typically breeds in North America during the boreal summer and migrates to winter in South America. However, this species started to breed in the coast of central Argentina in the early 1980s, and currently, it is expanding its breeding range (Martínez, 1983, Idoeta et al. 2011, Segura 2017. In this study, the species was recorded in rural areas during 2011, whereas it was recorded in the three habitats during 2017. According to Wilman et al. (2014) and Brown and Brown (2020), the Barn Swallow feeds on aerial invertebrates which contrasts with the dominant functional types in urban and suburban areas, comprised mainly of groundfeeding and omnivorous (Leveau 2013).
The SES of phylogenetic diversity was lower than random in rural habitats, which suggests that bird species adaptation to rural areas is phylogenetically preserved. Rural areas were particularly inhabited by seed-eater species of the Emberizidae family, such as Grassland Yellow-finch (Sicalis luteola), the Doublecollared Seedeater (Sporophila caerulescens), and the Rufous-collared Sparrow (Zonotrichia capensis), and also by migratory and insectivorous species of the Tyrannidae family, such as the Tropical Kingbird (Tyrannus melancholicus) and the Fork-tailed Flycatcher (Tyrannus savana) (Table 1). This pattern agrees with those obtained by Weideman et al. (2020), who found more clustered phylogenetic structures in rural than urban areas. However, Moreno-Contreras et al. (2019) found a random phylogenetic structure in rural and urban green habitats. Differences between studies may be related to differences in the species pools and the habitats compared. For example, Moreno-Contreras et al.
(2019) compared bird communities in urban green areas with agricultural and natural areas of the Chihuahuan Desert.

Conclusions
Bird communities in rural areas tended to be phylogenetically clustered and functionally random, suggesting that phylogenetic history influenced species presence in rural areas (Fig. 3). In contrast, bird communities in urban areas tended to be phylogenetically random and functionally clustered, evidencing environmental filtering. Moreover, due to species with similar traits but belonging to different clades, a possible process of functional convergence was related to urbanization ( Fig. 3) (Webb et al. 2002). Along the rural-urban gradients analyzed, there was a significant increase of bird functional diversity after a 6-year gap in data collection, probably related to increases in precipitation and temperature between years. However, the processes structuring bird communities probably were stable in the long term. Multiscale studies are needed to get insights into bird community assembly rules of urban areas.
Abbreviations GLMM: Generalized linear mixed models; LRT: Likelihood ratio test; SES: Standardized effect size
Additional file 1: Table S1. Functional traits for the species recorded during December 2011 and 2017 along three urban-rural gradients of central Argentina. Clutch: mean clutch size; Inv: invertebrates; Vend: mammals and birds; Vect: reptiles, snakes, amphibians, salamanders; Vfish: fish; Fig. 3 Summary of the results obtained in the study. Species inhabiting rural areas (black symbols) are phylogenetically clustered, and they are functionally diverse (different symbols). On the other hand, species inhabiting urban areas (black symbols) belong to different clades but have similar traits (same symbols) Vunk: vertebrates-general or unknown; Scav: scavenge, garbage, offal, carcasses, trawlers, carrion; Nect: nectar; PlantO: plant material; Watbelowsurf: foraging below the water surfaces; Wataroundsurf: foraging on or just (<5 inches) below water surface; Understory: Foraging below 2 m in understory in forest, forest edges, bushes or shrubs; Midhigh: Foraging in mid to high levels in trees or high bushes (2 m upward), but below canopy