Relative importance of deterministic and stochastic processes for beta diversity of bird assemblages in Yunnan, China

. Evaluating the relative importance of deterministic and stochastic processes underlying taxonomic and functional beta diversity is crucial in community ecology, because it can reveal the dominant processes of community assembly. However, studies of bird communities remain rare and of limited spatial extents. In this study, we described the taxonomic and functional beta diversity patterns of 32 passerine bird assemblages of Yunnan Province, China. We constructed null models based on observed species beta diversity and used multiple regressions on distance matrices to evaluate the relative contributions of deterministic and stochastic processes to passerine bird assemblage dissimilarity. Our results showed signi ﬁ cant geographic distance decay in taxonomic and functional similarity, with passerine bird assemblages located in the northwest and southwest of the province having higher functional beta diversity values than expected. Environmental distance and geographic distance explained a similar amount taxonomic beta diversity, but environmental distance explained much more functional beta diversity. Our results suggest that both deterministic and stochastic processes drive taxonomic beta diversity, but that deterministic processes, particularly environmental ﬁ ltering, play a dominant role in driving functional beta diversity of passerine bird assemblages at sub-national scale.


INTRODUCTION
Deterministic and stochastic processes are thought to be important in governing the structure of natural communities. Deterministic processes include niche-based environmental filtering and competition, while stochastic (neutral) processes include dispersal limitation and ecological drift (Hubbell 2001, Levine andHilleRisLambers 2009). It is generally accepted that both types of process are at work simultaneously in most communities, but it is difficult to reach consensus on their relative contribution to community structure (Legendre et al. 2009, Qian andRicklefs 2012). In fact, many specific processes of assemblage or coexistence are usually capable of producing observed biodiversity patterns, and ecologists must apply indirect approaches to identifying dominant processes as a result (Kraft and Ackerly 2010, Swenson et al. 2012, Chun and Lee 2017.
Beta diversity was originally proposed to represent the spatial compositional turnover between natural assemblages in order to bridge the gap between local (alpha) and regional (gamma) measures of diversity (Whittaker 1960). Given the limitations of taxonomic beta diversity measures for capturing meaningful ecological differences between species, recent research has shifted in emphasis toward trait-based (functional) patterns of diversity, which might better capture community assembly across environmental gradients (Pavoine and Bonsall 2011, Münkemüller et al. 2012, García-Giróna et al. 2019. Taxonomic and functional diversity should in general be positively correlated, but if two communities inhabit similar environments, then they may develop similar trait values even if geographically isolated from one another, resulting in high taxonomic but low functional beta diversity. Conversely, two communities may have both high taxonomic and functional beta diversity if their habitats are connected but in different environmental conditions (Smith and Wilson 2002, Fukami et al. 2005, Dalmolin et al. 2019. Integrating these two aspects of beta diversity could therefore help to achieve better understanding of the dominant processes in community assembly (Siefert et al. 2013, Dalmolin et al. 2019. Most community assembly studies that utilize these different aspects of beta diversity have focused on plant communities (e.g., Swenson et al. 2011, Siefert et al. 2013, Wang et al. 2019). Studies of animal communities have used taxonomic beta diversity almost exclusively, with functional beta diversity evaluated only over limited spatial and temporal extents (e.g., for ants, Bishop et al. 2015;amphibians, Dalmolin et al. 2019;fish, Villéger et al. 2012;birds, Monnet et al. 2014, Si et al. 2016, Jarzyna and Jetz 2018, Fluck et al. 2020and mammals, Penone et al. 2016). As a globally distributed vertebrate group, bird assemblages have taxonomic and functional beta diversity that is dependent on different processes at different spatial and temporal scales (Jarzyna and Jetz 2018), as with other taxa (Swenson et al. 2011, Chun andLee 2017). However, few studies have investigated regional-scale taxonomic and functional beta diversity in bird assemblages to date, leaving a substantial gap in understanding how community structure is established (but see Weinstein et al. 2014 for hummingbirds).
Yunnan (21°09 0 and 29°15 0 N, 97°32 0 and 106°12 0 E) is located in a transitional zone between tropical southeast Asia and temperate east Asia, positioned at a sutural zone between Gondwana and Laurasia (Jin 2002, Metcalfe 2006. The region forms a major part of the Indo-Burma biodiversity hotspot (Myers et al. 2000) and was ranked as the richest province for birds in China, with more than 840 bird species recorded there-about 65.5% of the total number in China in just 4.1% of the country's land area (Ji 2006). Its unique geological history, as well as diverse topography and climate, makes the region an ideal experimental system to investigate the regional-scale drivers of the biodiversity dissimilarity between bird assemblages.
In this study, we explore the relative contributions of deterministic and stochastic processes, represented by environmental filtering and dispersal limitation, in driving taxonomic and functional beta diversity of passerine bird (Aves: Passeriformes) assemblages of Yunnan at the province scale and county resolution. We assessed taxonomic and functional beta diversity patterns between assemblages, constructed null distributions to test whether the functional dissimilarity was higher or lower than random expectation given the observed species beta diversity, and assessed spatial and environmental explanations using multiple regressions on distance matrices (MRM). We focused on the following questions: (1) whether the taxonomic and functional dissimilarity increases with increasing geographic distance (i.e., shows distance decay); (2) whether the functional dissimilarity between assemblages from extreme habitats (e.g., high latitude, high elevation) and other habitats are significantly higher than expected (indicating deterministic structuring); and (3) whether spatial or climatic explanations of taxonomic and functional beta diversity are significantly different. Specifically, we hypothesized that taxonomic and functional dissimilarity would increase with increasing geographic distance based on the general pattern of distance decay of biodiversity (Soininen et al. 2007). Further, because bioclimatic variables that co-vary with latitude and elevation are important factors for diversity patterns of terrestrial vertebrates (Buckley and Jetz 2008), we expected that functional dissimilarity between the mountainous northwest region and the others would be v www.esajournals.org significantly higher than under a null model. Finally, we hypothesized that both deterministic and stochastic processes work simultaneously, as most former research has revealed, but that deterministic environmental processes would make a larger contribution to functional beta diversity than stochastic spatial processes. This was expected to follow from the constraining effect of the environment on bird functional traits, relating to resource requirements and foraging, while space has less importance due to the study taxa having high mobility.

Bird species composition
We based our analysis mainly on passerine bird composition at province scale drawn from the authoritative book, The avifauna of Yunnan, Part II (Passeriformes) (Yang and Yang 2004), which records 508 passerine species in the province. We excluded bird species with a migratory status of winter resident/passing bird/straggler or a distribution status of rarely seen to retain stable bird assemblages. We also excluded counties with less than 90 recorded passerine bird species in consideration of unbalanced recording. In addition, we excluded counties with areas less than 2500 km 2 and more than 10,000 km 2 to minimize bias in assessment of species and functional dissimilarity between assemblages. This left a data set of 389 passerine species distributions across 32 counties, including 103 species with some subspecies that were only partially resident (Appendix S1: Table S1). There was no significant relationship between species richness and county area in this final data set (Pearson's correlation, r = 0.20, P = 0.265). The bird taxonomy and nomenclature follow the BirdLife taxonomic checklist v8.0 (Bird-Life International 2015).

Functional traits
Morphological traits were initially used to illustrate species interactions within natural assemblages and the ecological space the assemblages occupied (Ricklefs and O'Rourke 1975, Travis and Ricklefs 1983, Dehling et al. 2014, Silva et al. 2016). However, morphological traits are unable to fully capture ecological niches represented by, for example, foraging characteristics in birds (Miles andRicklefs 1984, Pigot et al. 2016). In this study, we therefore combined morphological and dietary traits to assess the functional diversity of passerine bird assemblages. The four morphological traits used were body mass, beak length, wing length, and tarsus length. Body mass is considered to be the most informative trait of animals related to resource and/or energy requirements and metabolism (Speakman 2005). Bill length is related to prey size and handling while wing length and tarsus length have strong relationships to movement capacity and foraging behavior (Miles andRicklefs 1984, Luck et al. 2012). Data on the four morphological traits were compiled from Yang and Yang (2004), with supplementary data taken from Zhao (2001) and Dunning (2007; Appendix S1: Table S2). Average values for adult males and females were calculated separately if both were available and also for species with more than one recorded subspecies.
Dietary functional traits relate to ecosystem function and niche utilization (Greenberg et al. 2000(Greenberg et al. , Şekercioglu et al. 2004. In a global data set produced by Wilman et al. (2014), dietary traits were classified into 10 different food categories, which were not mutually exclusive, based on the percentage usage of each resource type. In this study, the categories vend (mammals and birds), vect (reptiles, snakes, amphibians, salamanders), vfish (fish), vunk (vertebrates-general or unknown), and scav (scavenge, garbage, offal, carcasses, trawlers, carrion) were combined to one category as vers, in consideration of the few passerine species that use these food types. Our final database was comprised of four morphological traits: body mass, beak length, wing length and tarsus length, and six dietary traits: inv (invertebrates-general), vers, fruit (Fruit, drupes), nect (nectar, pollen, plant exudates, gums), seed (seed, maize, nuts, spores, wheat, grains), and plant (other plant material; Appendix S1: Table S2). All functional trait data were normalized to ranges between 0 and 1 with min-max normalization prior to further analysis.
Previous studies demonstrated dependence between certain bird functional traits. For example, body mass usually showed a positive relationship to other aspects of body size such as wing length (Balen 1967, Pigot et al. 2016, while beak length can have important implications for diet (Nebel et al. 2005). In consideration of v www.esajournals.org potential correlations between the 10 functional traits, and therefore redundant data after Pearson correlation analysis (Appendix S2: Table S1), we ran principal component analyses (PCA) using the rda function in package vegan (Oksanen et al. 2019) in R (R Core Team. 2020). We kept the components needed to explain at least 90% of the variance in the data. As a result, the first five components which together accounted for 94% of the variation in functional traits were kept and used in the following analysis (Appendix S2: Table S2).

Environmental data
Climate is believed to be important to bird distributions both directly and indirectly, especially at regional scales (Parmesan andYohe 2003, Spooner et al. 2018). We chose temperature and precipitation which usually are associated with the diversity and large-scale distributions of terrestrial vertebrates (Buckley and Jetz 2008;Qian and Ricklefs 2012). Data for annual mean temperature (bio1), temperature seasonality (standard deviation × 100; bio4), mean temperature of warmest quarter (bio10), mean temperature of coldest quarter (bio11), annual precipitation (bio12), and precipitation seasonality (coefficient of variation: mean/SD × 100; bio15) were obtained from WorldClim version 2.0 with 2.5 arc min spatial resolution (approximately 5 km; the target temporal range of the data set was 1970-2000, Fick and Hijmans 2017). Annual mean temperature, temperature seasonality, annual precipitation, and precipitation seasonality are indicators of energy and water availability and seasonality; mean temperature of the warmest quarter and mean temperature of the coldest quarter represent energy limitations (Qian and Ricklefs 2012). The raster function of the R (R Core Team 2020) package raster (Hijmans 2020) was used to read the data of the six bioclimatic variables, and polygon data describing the boundaries of counties in Yunnan Province were read by the shapefile function of the same package. Then, mean values of all variables for the 32 counties were extracted with the extract function of the package raster. Each variable was again min-max normalized. As described before, we ran principal component analyses (PCA) after Pearson correlation analysis on the six bioclimatic variables (Appendix S2: Table S3 and Table S4). The first two principal components explained 93% of the variation and were used in the analysis.

Data analysis
Taxonomic and functional beta diversity between assemblages.-Taxonomic beta diversity was calculated as the Bray-Curtis distance using species composition data (Magurran 1988). Geographic distance between each pair of counties was calculated as the shortest line between their centers taking the curvature of the earth into account (Appendix S1: Table S4). The geographic distance was used to check whether the taxonomic and functional beta diversity increase with geographic distance (distance decay; Nekola andWhite 1999, Soininen et al. 2007).
Many indices to assess the functional beta diversity among communities have been proposed (e.g., Anderson et al. 2006, Ricotta and Burrascano 2008, Villéger et al. 2011. We chose two widely used presence-absence weighted functional beta diversity metrics for the bird occurrence data in this study: the mean pairwise functional distance (D pw ) and the mean nearest functional neighbor distance (D nn ). They represent two main classes of functional beta diversity and are believed to be unrelated Burrascano 2009, Swenson 2011). Both metrics were calculated based on a trait distance matrix (distance-based) describing the functional dissimilarity between species. The mean pairwise functional distance (D pw ) is a basal metric of beta diversity, calculated as the average distance between all species in one community and all species in a second community. The mean nearest functional neighbor distance (D nn ) is a terminal metric measuring the mean distance between each species in one community and their functionally closest species in a second community (Swenson 2011). The selected first five PCA axes of the functional traits were used to construct the Euclidean trait distance matrix.
The mean pairwise functional distance (D pw ) and the mean nearest functional neighbor distance (D nn ) can be quantified as follows: where n k1 represents the number of species in community k 1 ; n k2 represents the number of species in community k 2 ; δ ik2 is the mean pairwise functional trait distance between species i in community k 1 to all species in community k 2 and δ jk1 is the mean pairwise functional trait distance between species j in community k 2 to all species in community k 1 ; minδ ik2 is the nearest neighbor functional traits distance between species i in community k 1 to all species in community k 2 and minδ jk1 is the mean pairwise functional trait distance between species j in community k 2 to all species in community k 1 . We excluded conspecifics between two assemblages in calculating the nearest neighbor functional traits distance.

Null model of functional beta diversity
To test whether the functional dissimilarity was higher or lower than random expectation given the observed species beta diversity, a null distribution was constructed for each pair of assemblages by shuffling species names on the functional trait distance matrix 999 times, while maintaining the observed spatial distribution of species (Swenson 2014). A standardized effect size (SES) was quantified for functional beta diversity as follows (Gotelli and Graves 1996): where X obs is the observed functional beta diversity value (D pw or D nn ) between two communities, X null is the mean value of the null distribution, and SD(X null ) is the standard deviation of the null distribution. Positive values indicated a higher observed functional dissimilarity than expected while negative values indicated a lower dissimilarity than expected. The average value of the SES for each assemblage against all others (n = 31) was used to check the patterns of regional dissimilarity in functional trait composition.

Spatial and climate explanation of taxonomic and functional beta diversity
To distinguish the effects of spatial (geographic) distance and environmental (climate) dissimilarity at county resolution, MRM were used (Lichstein 2007). Multiple regressions on distance matrices is an extension of partial Mantel analysis which involves a multiple regression of a response matrix on any number of explanatory matrices. We used MRM to partition variance in Bray-Curtis distance, D pw and D nn , between unique contributions of spatial and environmental distances, an effect shared between them, and unexplained variance. The response matrix provides Bray-Curtis distances, D pw and D nn , and two explanatory matrices gave geographic distances between each pair of counties and climate Euclidean distances constructed from the first two climatic PC axes.
The analyses were performed with the vegdist function of the package vegan (Oksanen et al. 2019), comdist and comdistnt functions of the package picante (Kembel et al. 2010), abind function of the package abind (Plate and Heiberger 2016), and lm function in R (R Core Team 2020).

Beta diversity patterns and the standardized effect sizes of functional beta diversity metrics
Passerine bird species richness (SR) of the 32 assemblages ranged from 92 to 197 (Appendix S1: Table S1). Taxonomic beta diversity was positively correlated with geographic distance (distance decay; P < 0.01, Fig. 1a). Both the mean pairwise functional distance (D pw ) and the mean nearest functional neighbor distance (D nn ) between assemblages increased moderately with geographic distance (P < 0.01, Fig. 1b and Fig. 1  c). Meanwhile, the D pw and D nn showed significantly positive correlations (Mantel's r = 0.61). Taxonomic beta diversity was significantly related to D pw (Mantel's r = 0.52) and D nn (Mantel's r = 0.57). Environmental dissimilarity also increased with geographic distance (P < 0.01, Fig. 1d).
With construction of null models maintaining the observed taxonomic beta diversity while randomizing functional traits, values of the SES of the mean D pw and the mean D nn showed similar patterns (Fig. 2). Most assemblages showed lower observed mean D pw and mean D nn than the null distribution (84% and 72%, respectively), indicating more functional similarity than expected. All the counties with higher-thanv www.esajournals.org expected mean D pw and eight out of nine counties with higher-than-expected mean D nn were located in the northwest and southwest region of the province (Fig. 2), demonstrating that these assemblages were more functionally dissimilar to others than expected.

Spatial and environmental explanation of taxonomic and functional beta diversity
The MRM results showed that 65% of the variation in taxonomic dissimilarity could be explained by space and environment, with the variation partitioned into pure space (13%), shared space/environment components (38%), and pure environment components (14%; Appendix S2: Table S5 and Fig. 3). In total, space and environment could only explain 17% and 25%, respectively, for D pw and D nn , however, substantially less than the amount explained for taxonomic dissimilarity. The functional dissimilarity explained by space (pure space plus shared components) was 3% and 12%, respectively, for D pw and D nn , whereas the variation explained by the environment (pure environment plus shared components) was 17% and 24%, respectively, for D pw and D nn . In the space + environment model for functional dissimilarity, environment remained significant, but space was no longer   2. Results of null model analyses of functional beta diversity. The numbers are the mean standardized effect sizes of (a) the mean pairwise functional distance (SES D pw ) and (b) the mean nearest functional neighbor distance (SES D nn ) for each focal assemblage to others. Positive values indicate higher functional beta diversity than expected given the observed species beta diversity; negative values indicate lower functional beta diversity than expected. SES values bigger than 1.96 indicate significantly higher functional beta diversity than expected; SES values smaller than −1.96 indicate significantly lower functional beta diversity than expected given the observed species beta diversity. The elevation was extracted from the CGIAR-CSI SRTM data set version 4.1 (Jarvis et al. 2008). The map was generated using ArcGIS 10.2 (ESRI 2014). significant after controlling for the environment (Appendix S2: Table S5).

DISCUSSION
We described the taxonomic and functional beta diversity of passerine birds across Yunnan Province at county resolution and sought to evaluate the relative contribution of deterministic and stochastic processes to these patterns. Our results showed a strong increase of taxonomic dissimilarity and a moderate increase of functional dissimilarity with geographic distance (Fig. 1). The observed functional beta diversity of most assemblages was lower than would be expected by chance, but assemblages located in the northwest and southwest region showed higher values (Fig. 2), providing evidence of underlying deterministic environmental filtering. Environmental distance alone was also found to explain much more functional beta diversity than geographic distance alone (Fig. 3), suggesting that different underlying processes govern taxonomic and functional beta diversity of passerine birds. Overall, both deterministic and stochastic processes drive taxonomic beta diversity simultaneously, while deterministic processes, particularly environmental filtering, play a dominant role in driving functional beta diversity at the scales examined in our study. The combination of different aspects of biodiversity should help us to gain a more detailed insight into biodiversity patterns and processes. Our results support the view that biological dissimilarity usually increases with geographical distance, according to our first question and hypothesis. Former research revealed that increasing dissimilarity in biodiversity with increasing geographical distance could result from either deterministic environmental or neutral assembly processes (Legendre et al. 2005, Soininen et al. 2007. We find a different balance of these affecting taxonomic and functional beta diversity. This may indicate that spatial dispersal limitation affects taxonomic beta diversity more than functional beta diversity (Chase and Myers 2011). The relative importance of deterministic and stochastic processes underlying biodiversity was believed to depend upon the spatial scale of analysis (Mod et al. 2020). To our knowledge, this study is the first to analyze bird community assemblages using functional beta diversity at administrative scale (province and counties). The results therefore provide a general description of taxonomic and functional dissimilarity between areas in the region, and the relative importance of environment and space in generating these patterns.
We found that the geographic patterns of the mean SESs of both indices of functional beta diversity were similar (Fig. 2), assemblages in the central region of Yunnan were more similar in terms of functional traits than expected, whereas assemblages in the northwest and southwest of the province were more dissimilar to the others than expected, implying deterministic structuring (Fukami et al. 2005). As a mountainous region, Yunnan has a terraced terrain from high mountain summits in the northwest to low valleys in the south (Wu 1987). The latitude of Yunnan ranges from 28°N to 21°N, leading to a temperature and precipitation gradient from north to south. The high latitude and elevations of the northwest (Fig. 2) produce lower temperatures and less precipitation. In contrast, the southwest region is at low latitude and elevation, with higher temperature and more precipitation. In addition, the flora of the northwest has been affected mostly by the uplift of the Himalayas, while the flora of southern Yunnan has evolved with the influence of mainly tropical Asian elements (Hua 2012). Typical vegetation types of northwest Yunnan are coniferous forest, alpine shrub, and alpine meadow, while tropical rain forest and evergreen broadleaved forest dominate in the southwest. As a result, the northwest and southwest counties of Yunnan have extremely different environmental conditions than other counties. The climate and vegetation conditions are known to affect passerine bird assemblages; for instance, species such as Fringillidae prefer seeds that are mostly found in the north region of Yunnan, whereas species such as Dicaeidae and Nectariniidae forage nectar found in the south (Yang and Yang 2004, Appendix S1: Table S1). The fact that these differences link to clear differences in functional passerine bird beta diversity indicates that environmental filtering plays a key role in structuring these communities. For a given level of species turnover, functional beta diversity will be higher than expected if the environment changes rapidly and lower than expected if the environment remains constant in space or time (as known from plant assemblages; Swenson et al. 2011, Siefert et al. 2013. We found that functional beta diversity of passerine bird assemblages was higher than expected in the northwest and southwest region of the province, but lower than expected for the counties located in the central region. This conforms to the common pattern found in plant assemblages and confirms the dominance of deterministic processes. The importance of nichebased ecological process in functional turnover of animal assemblages has also been reported (Weinstein et al. 2014, Hill et al. 2019. The final question addressed by this research is whether there are differences between the spatial and climatic effects on taxonomic and functional beta diversity. Our finding that geographic distance alone explained much more variation in taxonomic than in functional beta diversity (13% vs. 0% and 1%; Fig. 3) supports the view that taxonomic composition turnover is more strongly influenced by neutral processes than functional composition turnover (Swenson et al. 2012). It is unsurprising that neutral processes, especially dispersal limitation, should have a more important role in taxonomic beta diversity at the province scale than in functional beta diversity, because passerine bird species turnover could be increased by geographic distance even while functional traits converge (Smith and Wilson 2002). Qian and v www.esajournals.org Ricklefs (2012) found that both deterministic and stochastic processes drive taxonomic species turnover at global and regional scales in four classes of terrestrial vertebrates including birds, and our study also supports this view for passerine bird taxonomic turnover at the province scale. It is also notable that the variation unexplained by environmental and spatial factors was more than 75% for functional beta diversity and 35% for taxonomic beta diversity, suggesting impacts of other factors not considered here. For example, it is difficult to evaluate species interactions at county resolution, with competition between species within assemblages usually detected at the local scales at which interactions occur (Cavender-Bares et al. 2009). Geographic barriers, especially high mountain ranges in this region, could also be associated with spatial species turnover and functional trait turnover (Buckley and Jetz 2008), but are not analyzed here.
With the development of beta diversity, dozens of measures of taxonomic and functional beta diversity have been proposed (e.g., Whittaker 1972, Chao et al. 2005, Villéger et al. 2011. In this study, we did not partition beta diversity into its spatial turnover and nestedness-resultant components. Instead, we used a traditional species beta diversity Bray-Curtis distance and two popular distance-based measurements of functional beta diversity (Baselga 2010). It is likely that the functional turnover component dominates overall functional beta diversity in this study system, but additional studies will be needed to confirm and quantify this (Bishop et al. 2015, Si et al. 2016.
Global change and disturbance pose important threats that are likely to alter the distributions of living species (Brawn et al. 2001, La Sorte et al. 2017. For the past few years, conservation ecologists have paid close attention to the functional diversity of natural communities in addition to their taxonomic diversity, as functional diversity is likely to be more sensitive to environmental change and anthropogenic influence (Sitters et al. 2016, Anjos et al. 2019, Vaccaro et al. 2019. Here, the environment explained more variation in functional dissimilarity than did space, giving some insight into the environmental dependence of functional trait turnover of passerine bird assemblages, and hence their vulnerability to environmental change.

CONCLUSION
Our results suggest that both deterministic and stochastic processes drive taxonomic beta diversity, whereas deterministic processes, particularly environmental filtering, play a dominant role in driving functional beta diversity of passerine bird assemblages at cross-county scale. Investigations that make use of taxonomic and functional beta diversity are likely to achieve better identification of dominant processes underlying community assembly. The lack of current studies that make use of these measures in bird communities at regional scales is a clear research gap, and one that may obscure the key role of environmental change in determining functional diversity.