Taxonomic and functional components of avian metacommunity structure along an urban gradient

Identifying biological processes that structure natural communities has long interested ecologists. Community structure may be determined by various processes, including differential responses of species to environmental characteristics, regional-level spatial influences such as dispersal, or stochasticity generated from ecological drift. Few studies have used the metacommunity paradigm (interacting communities linked by dispersal) to investigate avian community composition along an urban gradient, yet such a theoretical construct may provide insights into species turnover even in unnatural settings such as rural to urban gradients. We measured the influence of spatial and environmental characteristics on two aspects of avian community structure across a gradient of urbanization: 1) taxonomic composition and 2) functional richness based on diet, foraging strategies, nesting locations and morphology. We also measured the relationship between species traits and environmental variables with an RLQ-fourth corner analysis. Together, environmental and spatial processes were significantly related to taxonomic structure and functional richness, but spatial variables accounted for more variation than environmental variables. Fine spatial scales were positively correlated with insectivorous birds and negatively correlated with body and wing size. Urbanization was positively correlated with birds that forage at the canopy level, while emergent wetlands were negatively correlated with birds that nested in cliffs and frugivorous birds. Functional richness and urbanization were significantly related to fine spatial variables. Spatial and environmental factors played an important role in taxonomic and functional structure in avian metacommunity structure. This study highlights the importance of studying multiple aspects of biodiversity, such as taxonomic and functional dimensions, especially when examining effects of complementary spatial and environmental processes.


Introduction
Community ecology aims to understand the primary mechanisms influencing species abundance and composition at the local level [1].Metacommunity theory differentiates local and regional processes that influence community structure [2] and can be simplified into four frameworks: neutral theory, patch-dynamics, species-sorting and mass-effects [2].Initiated by Hubbell [3], neutral theory assumes differences among species regarding their niches are nonexistent or at least unimportant to community structure and focuses on how randomly fluctuating demographic processes and dispersal limitation influences diversity of local communities [2].The patch-dynamics framework focuses on how tradeoffs, for example in competitive and dispersal abilities, influence temporal dynamics of communities in environmentally homogenous landscapes [2].The species-sorting framework focuses on species response to a heterogeneous landscape, whereby patterns of presence and absence of species reflect selection for suitable habitats [2,4].Building upon the metapopulation framework of sources and sinks, the mass-effects framework focuses on how high rates of dispersal of individuals of multiple species into less suitable habitats facilitates coexistence at the local level [2,4].
Environmental filtering can strongly influence species coexistence [4].However, environmental characteristics often explain no more than 50% of the variation in taxonomic diversity at the local level [4].While spatial processes or stochasticity often strongly characterize taxonomic composition of communities [4], environmental filtering can have strong influences on related functional characteristics [4][5][6].Functional traits are important characteristics that are relevant to performance and are expected to have a strong association with environmental gradients [4].By examining a variety of functional traits within a large species pool, it may be possible to gain a deeper understanding of regional species distribution and abundance [4].
Combining functional characteristics with metacommunity models also has the potential to provide understanding of effects of anthropogenic modifications on the contemporary biota [4].Avian community composition has been frequently examined as a response to urban gradients at local scales [7,8].Nonetheless, only a few studies have used a metacommunity framework to investigate relative contributions of local environmental and regional spatial processes on avian community composition.Local (e.g., resources, patch quality or species interactions) and regional (e.g., connectivity) processes are hypothesized to influence urban community composition [7,9].This suggests that environmental filtering, habitat selection or dispersal may affect urban community structure [7,9] to an equal or greater extent than natural communities.Furthermore, metacommunity theory can be applied to gain insights into species turnover from rural to urban settings, as well as the traits that are related to environmental filtering [7,10].
To estimate metacommunity structure, we examined taxonomic (species identity) and functional characteristics (morphology, diet, foraging strategy, nesting location, and native/invasive status) of birds distributed throughout Texas.We examined trait-environmental relationship with RLQ and fourth-corner analyses [11][12][13].Both analyses are dependent on three matrices, L (species x site matrix), R (environment x site), and Q (species x trait matrix), but provide different perspectives on the structure of communities [13].While RLQ analyses examine associations among traits and environmental characteristics from a multivariate perspective [13].Fourth-corner analysis examines species pairwise comparisons of traits and environmental variables [13].Because both approaches are complimentary, Dray et al. proposed combining both approaches into one [13].
Increasing urbanization can adversely affect species and can cause strong habitat filtering [9,14].In addition, predation on nest sites is often higher in urban areas compared to natural environments [15].Species with nests close to or on the ground can be more impacted by predation within urban environments [15][16][17].Furthermore, there is unequal access to food resources for different dietary groups.Insect populations are lower in cities due to increased use of insecticides [18].However, increases in seed dispensed by humans cause urban areas to be rich in resources for granivorous bird species [18].Urbanization increases local extinctions of species and populations [19].For all these reasons, urban communities tend to be primarily composed of introduced invasive species or highly adaptable native species [20].Due to the multiple factors mentioned above, we predict that: 1) birds that are granivorous, nest off the ground, or are invasive will be more associated with urban environments, 2) birds that are native, insectivorous, and ground nesters will be more associated with natural environments, and 3) species richness will be lower in areas of high urbanization and will increase in more natural or rural areas.

Study Area
This study focuses on avian communities in Texas, the state with the greatest number of recorded bird species (647) within the United States [21].Texas is highly urbanized.Houston (2,303,482), San Antonio (1,492,510), and Dallas (1,317,929), represent metroplexes with three of the ten highest human populations in the United States.Besides being highly urbanized in parts of the state, Texas is composed of multiple ecoregions.Eastern Texas has a gradient of pine forests to prairies with wetlands in the south, while central Texas has a gradient of cross timbers to open grasslands.Continuing westward there is an increase in mesquite, hill country, canyons and deserts.Precipitation increases from west to east, and temperature increases from north to south [21].

Bird Data -L matrix
We collected data from the Breeding Bird Survey (https://www.pwrc.usgs.gov/bbs/results/)and eBird (http://ebird.org/ebird/data/download) from 1 May to 31 August 2013 to 2017.The Breeding Bird Survey was developed by the US Geological Survey's (USGS) Patuxent Wildlife Research Center and Environment Canada's Canadian Wildlife Service to monitor North American bird populations.The eBird database is a citizen science project whose data has demonstrated to be effective for description of patterns of diversity at multiple spatial scales and to be comparable to other more standardized data sets [23].To limit the inherent biases within eBird data, we followed methods developed by Callaghan et al. and Sullivan et al. [24][25][26].We removed the following kinds of checklists: (1) those that did not report all observed species, (2) were duplicates from multiple observers who participated in the same sampling event, or (3) situations where the observer traveled greater than 5km or covered more than 500 ha.We kept checklists (4) that had a duration between 5 to 240 minutes and (5) followed traveling, random or stationary protocols.We focused on observations of Passeriformes (215 species), Columbiformes (9 species) and Psittaciformes (2 species).We included orders Columbiformes and Psittaciformes to discern how environmental and spatial processes affect invasive species.Observations from avian surveys were plotted in ArcGIS (ESRI, Redlands, California, USA).
We created a grid across Texas with cells of 40 km 2 .In the center of each square, we created a buffer with 10 km radius and an area of 314 km 2 .We expect the buffer to be greater than the home range of all the studied species.Home range size is related to body size in animals [27,28].Studies on species home ranges are limited, but many passerines within this study have a home range of less than 9 km 2 [29][30][31].Some of the larger species within the orders Passeriformes, Columbiformes, and Psittaciformes have home ranges ranging from 0.005 km 2 (feral rock pigeon, Columba livia) [32] to 325 km 2 (common grackle, Quiscalus quiscula) [33].
The buffer used here encompassed the site of multiple species home ranges, suggesting that multiple populations are present within each site.Sites were larger when compared to other studies and this increase can increase the probability of species detection [35].The 20 km distance between each site is greater than the home range of each studied species, which limits the possibility of a single individual being included in more than one site, thereby enhancing independence of data points.We extracted Birding Bird Survey and eBird GPS points that intersected the buffers to make up the communities.To ensure that we limited analyses to those communities that were well sampled, we removed communities with less than 50 individuals or those that did not exhibit an asymptote in species richness based on a rarefaction curve.

Trait Data -Q matrix
For ecological traits, we collected information on diet and percent foraging strategies from Wilman et al. [35] and morphological measurements and nesting strategies from Oberholser [36,37] and Ricklefs [38].We collected data on native status from Oberholser [36,37], to which we defined invasive species as non-indigenous birds, whose distribution expanded due to human facilitation.Because dietary variation depends on location and season, we identified species dietary guild using the item that was the highest proportion of their diet, as done by Wilman et al. [35].In cases where morphological measurements for females and males were collected, we averaged means between the sexes.Furthermore, we used presence (1) or absence (0) for nesting strategies (Table 1).Average body size, dietary characteristics (diet, foraging strategy and bill length) and nest type are related to environmental characteristics of niches [39][40][41].Wing lengths of birds are related to energetic costs of flight and facilitate movements across fragmented locations [42].We characterized trait richness using Rao's quadratic entropy [43][44][45].Trait information was not available for all species; therefore, we preformed these analyses on 189 species, 17 fewer than the taxonomic analyses.Many traditional methods examine trait diversity by organizing species into groups, rather than assessing quantitative species characteristics [45].We measured functional trait richness for diet, foraging type, nesting locations, bill length, wing length and body mass by using the R package "SYNCSA" [46].Furthermore, we calculated species richness of invasive species and native species for each community.

Spatial and Environmental Data -R matrix
We collected land cover data from the 2016 USGS National Land Cover Dataset that had a resolution of 30 m 2 [47].Using ArcGIS, we aggregated the land cover types within each buffer and then expressed each land cover type as a percentage.Moreover, we extracted the averages of precipitation and temperature from May to August of 2013 to 2017 for each site using PRISM that had a resolution of 4 km 2 .To obtain spatial coordinates, we extracted projected Universal Transverse Mercator coordinates (WGS 84, Universal Transverse Mercator, zone 14) for the center of each site via ArcGIS.

Statistical Analyses
To examine metacommunity structure, we created derived environmental variables with a principal components analysis (PCA) and derived spatial variables with principal coordinates of neighborhood matrices (PCNM).For environmental variables, we performed a PCA with a correlation matrix to produce a reduced set of important environmental variables that were orthogonal.We used the broken stick method to determine which PCs were significant (Peres-

Neto et al. 2003).
To construct PCNM's, we followed the protocol of Borcard and Legendre [48] using packages "vegan" [49] and adespatial [50] in R 4.0.3by: 1) creating a Euclidean distance matrix, 2) computing principal coordinates from a truncated distance matrix with a defined threshold of 4 times the nearest neighbor sampling distance as suggested by Borcard and Legendre [48], 3) tested significance of the global model with a canonical correspondence analysis (CCA) with species occurrences as the dependent matrix and all PCNMs as the independent matrix and, 4) after a significant global model, we assessed significance with forward selection based on a CCA and retained only eigenvectors with positive eigenvalues that were significant.
Variation partitioning of a constrained ordination, such as a CCA, partitions variation among multiple groups of explanatory variables [51] and allows for examination of unique variation related to a particular explanatory group (i.e., environment) after controlling for shared variation with other explanatory groups (i.e., space) [51].We examined relationships among spatial factors, environmental factors and species composition with a partial CCA to determine: 1) species composition uniquely related to environmental characteristics based on significant PCs, 2) species composition uniquely related to spatial factors based on significant PCNMs, and 3) species composition related to spatially structured environmental characteristics [52,53] in Canoco 5 [54].The dependent matrix was comprised of species presence/absence for 79 communities, and independent matrices were five significant spatial PCNMs and three significant environmental PCs.Since multiple independent variables are being used, we used adjusted coefficient of determination (R 2 ) to estimate effect size.We used Monte Carlo permutation (999 permutations) to determine significance of unique variation accounted for by environmental and spatial variables.Following a similar procedure as for taxonomic composition using a partial constrained analysis, we examined unique relationships of trait richness with environmental and spatial factors with a partial redundancy analysis (RDA) [52,53] in Canoco 5 [54].
We examined the relationship between species traits and environmental variables with fourth-corner and RLQ analyses: R represented the site by environment matrix, L represented the site by species matrix, and Q represented the species by trait matrix [51].Matrices Q and R were coupled using an ordination to create linear combinations that are then linked to species presence/absence in matrix L [55].The fourth-corner tests the relationship between traits and environmental variables for each species separately, while an RLQ analysis is a multivariate analysis of the three matrices [13] that examines all species simultaneously.Combining both RLQ and fourth-corner analyses, we tested for global significance by examining two permutation models: Model 2 and Model 4. Model 2 permutes sites to determine if the relationship between species and environment was significant [56].Model 4 permuted species occurrences to examine if the relationship between species and traits was significant [56].P-values are estimated by comparing the observed value with statistics of the permutation distribution [13].This method has been suggested by Dray and Legendre to limit type 1 error [56].RLQ and fourth-corner analyses were performed using the R package "ade4" [57].

Bird Data -L matrix
We obtained data on species composition of 79 well-sampled communities.There were 205 species present within communities, with mourning doves (Zenaida macroura) being present at the most locations (97% of sites,

Spatial and Environmental Data -R matrix
Most common landcover type among the communities was shrubland with an approximate average area of 94.34 km 2 ± 24.03.This is followed by pastures (47.38 km 2 ± 14.79), grassland (35.37 km 2 ± 12.27), croplands (34.76 km 2 ± 13.84), and urbanized areas (31.04 km 2 ± 12.26).Principal components analysis on the 15 land cover types and average precipitation and temperature yielded 3 significant PC's (Table 2) that together account for 55.19% of the variation in environmental characteristics.The first PC accounted for 24.38% of the variation among sites regarding environmental variation and was interpreted as a gradient from shrub and grassland to more developed, urban areas (Table 2).The second PC accounted for 17.12% of the variation and encompassed a gradient of high urbanization to pastoral lands.

RLQ and Fourth Corner Analyses
There was a significant global relationship between species occurrences and environmental variables (Model 2: Observed value = 0.25, P < 0.001, Fig 3a).Furthermore, the relationship between species occurrences and traits, while preserving the link between species and environmental variables, was also significant (Model 4: Observed value = 0.25, P < 0.001,

Fig 3b
).This indicated that there was a significant multivariate pattern between traits and environmental variables.Spatial and environmental variables accounted for 86.29% of the variation in trait variables, with the first RLQ axis accounting for 74.93% of the variation.Fourth-corner analysis indicated that 8 out of 176 bivariate associations were significant (Fig 3c).Environmental principal component 1 was positively associated with foraging at canopy level (r = 0.09, P = 0.038).Environmental principal component 3 was negatively associated with a frugivorous diet (r = -0.06,P = 0.038) and nesting on cliffs or buildings (r = -0.07,P = 0.036).For spatial variables, PCNM 4 was positively related to nesting in shrubs (r = 0.07, P = 0.036), PCNM 14 and PCNM 30 were positively related to insectivorous diet (r = 0.10, P = 0.038; r = 0.10, P = 0.038), and PCNM 30 was negatively related to body mass (r = -0.070,P = 0.036) and wing length (r = -0.09,P = 0.038).
There were significant associations between RLQ axes, traits and explanatory variables when combining fourth-corner and RLQ approaches (Table 4).The first RLQ axis was negatively correlated with environmental PC 1, environmental PC 3, PCNM 9, PCNM 14 and PCNM 30.For trait variables, RLQ axis 1 was negatively associated with an insectivorous diet, foraging at mid-high and canopy level, and being native.In contrast, the first axis of the RLQ was positively related to granivorous diets, body mass, wing length, nesting in shrubs, and being invasive.Axis 2 was significantly related to PCNM 4 but no other variables.Communities in areas of increased urbanization and emergent wetlands were mostly associated with insectivorous, native species, and species that forage at mid-high to canopy level.Individuals

Discussion
Urbanization, land-use and climate were not primary driving forces of taxonomic and functional metacommunity structure.While significant, spatial and environmental processes accounted for less than 25% of variation in taxonomic characteristics or trait richness.Spatial characteristics accounted for more unique variation than environmental characteristics with respect to taxonomic components and trait diversity.Spatial and environmental variables also accounted for a large portion of the variation in RLQ matrices.However, environmental variables, specifically urbanization, were not as influential in community structure as we originally predicted.Canopy foraging was the only trait significantly related to increased urbanization.
In a meta-analysis of metacommunities conducted by Cottenie [53], the majority of studies demonstrated that environmental factors were the main driving force of community structure, indicating that species-sorting was the most common framework explaining metacommunity structure.The second most prominent structure was a combination of spatial and environmental variables indicating a combination of mass-effects and species-sorting [53].
Of the 158 studies, the four avian metacommunity studies examined exhibited species-sorting or mass-effects [53].Other avian metacommunity studies have also detailed the importance of environmental variables (species-sorting framework) on community structure [58][59][60][61][62]. Similar to our study, stochasticity was a major characteristic of metacommunity structure in a humanmodified landscape in France, but environmental filtering and dispersal were significant contributors [62].
Our study demonstrated that environmental filtering played a significant role in taxonomic and functional components of community structure.The gradient of emergent wetlands to forests played a major role in functional richness and species presence.Mesic conditions within wetlands can contribute to an increase in food and water availability, being related to increased avian richness and abundance compared to uplands [63,64].Richness in bill lengths, wing lengths and foraging strategies had the highest association with wetlands and Passerine dietary niches were related to their bill size, locomotion and foraging behavior in passerines [65].Bill length has been demonstrated to be longer for passerines in wetlands than in species from drier lands [66].Increased bill length can facilitate greater dietary breadth by increasing ability to reach insects in crevices [65,67].
There was a negative relationship between occurring in wetlands and having a frugivorous diet or nesting on buildings/cliffs.While longer bills may be more beneficial for birds in wetland areas, wider, well-developed, and sometimes shorter bills provide advantages to frugivores [68].Shorter bills allow birds to better manipulate fruit [68].Moreover, birds that nest on buildings or cliff sides tend to occur within cities especially in regions with little topographic relief such as the high plains of Texas.Canyons and mountains are in the western arid regions of Texas and serve as valuable nesting locations for cliff dwelling birds.Buildings within cities can provide pseudo-cliff structures for nesting, which invasive species such as feral pigeons take advantage of when nesting in cities [69].
We predicted that urbanization would influence species composition and functional traits.
However, the gradient of shrubland to urbanization was not significantly related to species composition but was correlated with invasive species and nesting richness.Urbanization greatly influences biodiversity via multiple synergistic mechanisms [9].Although these mechanisms can impact multiple native species, invasive birds thrive in cities and can outcompete similar native species [70][71][72].Adaptations, such as behavioral flexibility [10,73,74], ecological generalism [75,76], certain life history traits such as brood value (how many broods over a lifetime) [76] and human tolerance [74,77] have all been attributed to success of invasive species in urban environments.
Native species and urbanization had a negative association with RLQ axis 1 (Fig. 5).
Urbanized areas are characterized by decreased biodiversity compared to surrounding natural habitat [9].However, there are many attributes that can increase biodiversity in cities such as higher productivity and resource availability [9].Increased primary productivity via urban parks can alleviate the negative effects of urbanization [9] or the harsh environment of arid cities [78].
Parks can increase richness by contributing a wide variety of food and nesting sites [79].Ground nesting birds often are more prevalent in parks, while cavity and tree nesters are often more prevalent in allotment garden (e.g., community garden) [80].Canopy foragers are often prevalent in parks [81], explaining the positive correlation between urbanization and canopy foragers.
While cities are notoriously known for decreasing biodiversity, parks within cities may mitigate this effect [80].Due to the increase in food, water, and temperature (heat-island effect), birds may be able to thrive in many cities [9,81].Therefore, examining characteristics within cities may be key to understanding how urbanization influences species.
Principal coordinates of neighborhood matrices represent eigenvector decompositions of spatial scales and sites [83].Following Borcard and Legendre classification, PCNMs for this study can be classified into three different spatial groups: broad (PCNM 1, 4 and 9), intermediate (PCNM 14) and fine scale (PCNM 30) [83].For the multivariate analysis of traits, moderate to fine scale spatial variables were correlated with nesting in shrubs, body mass, wing length and insectivory.Similar to variation partitioning, correlations are potentially rooted in environmental variables that are related to finer spatial scales [84] or biotic processes (e.g., competition) [83].
My study focused on broad scale landscape variables.However, habitat quality (e.g., food availability) can influence avian distributions [85].Because of high mobility and larger species having larger home ranges [27,28], the negative correlation of the fine scale PCNM 30 with body mass and wing length may be related to dispersal ability.
Future studies would benefit from more censuses in the western region of Texas and from considering more environmental variables.Most of the communities come from the eastern region, which is probably due to human population density being greater in the east.Moreover, vegetation complexity and quality are influential to avian distribution and the addition of these variables to the model might improve understanding of metacommunity structure.
Overall, this study provides insight into how spatial and environmental variables play an influential role in avian community structure.Surprisingly, urbanization, or other measured environmental variables, were not the central driving force in metacommunity structure.Spatial variables play a pivotal role in avian communities, whether they be due to dispersal, environmental variables that are spatially structured, or biotic interactions.

Fig 1 .
Fig 1. Species Commonality in Texas.The commonality of avian species throughout the 79 communities during the summer of 2013 to 2017.

Fig 2 .
Fig 2. Relationship between spatial and environmental variables and taxonomy and

Fig
Fig 2d).Spatial variables PCNM 14 and PCNM 30 were negatively correlated with the first axis

Table 1 .
Traits gathered for species in Texas avian metacommunities in the summer of 2013 to 2017.

Table 2 .
Eigenvalues of significant principal components that are derived from land cover types, precipitation, and temperature characteristics for Texas avian metacommunities in the summer of 2013 to 2017.Bold eigenvalues indicate land cover types that assist in the gradient description

Table 3 .
Correlations of environmental and spatial to axes of a Canonical CorrespondenceAnalysis of species occurrences and Redundancy Analysis of trait richness for birds in Texas during the summer of 2013 to 2017.Degrees of freedom for both analyses were 77.Bold = Pvalue less than 0.05.Degrees of freedom for all correlations equal 77.

Table 4 .
RLQ Axes correlations with trait and explanatory variables.Bold = P-value less than