Decomposing multiple β‐diversity reveals non‐random assembly of the waterbird communities across anthropogenic subsidence wetlands

Understanding how biotic communities are spatially distributed is essential for biodiversity conservation in human‐modified landscapes. The large‐scale subsidence wetlands generated by underground coal mining in China have been increasingly used by large numbers of waterbirds, which remain understudied. We aimed to explore the role of non‐random assembly in shaping the regional waterbird diversity in non‐breeding seasons, providing insights for the conservation.


| INTRODUC TI ON
How biotic communities and biodiversity are distributed across space or time is a fundamental question in ecology and biodiversity conservation (Alahuhta et al., 2017;Anderson et al., 2010). Biodiversity can be studied at different spatial scales, with α-diversity measuring local diversity and γ-diversity measuring regional diversity (Jost, 2007;Ricotta & Szeidl, 2009). First introduced by Whittaker (1960Whittaker ( , 1972, β-diversity generally measures compositional variations among communities along environmental gradients, linking local (α-diversity) and regional (γ-diversity) biodiversity. Patterns of β-diversity may result from two different processes, that is spatial turnover and nestedness (Baselga, 2010;Soininen et al., 2017). Spatial turnover may be attributed to species replacement among assemblages fuelled by environmental filtering and dispersal capacity allowing species to track suitable environments, while nestedness results from species loss or gain related to ordered extinction-colonization dynamics (Baselga, 2010). Measuring and decomposing β-diversity allows testing hypotheses associated with biodiversity distributions and assembly processes of biotic communities (Baselga, 2010;Si et al., 2016), providing insights for management, conservation and restoration of biodiversity (Devictor et al., 2010;Socolar et al., 2016).
Recently, spatial partitioning of biodiversity has been extended from taxonomic diversity to functional metrics to capture that species have specific functional characteristics (Si et al., 2016;Villéger et al., 2013). Studies of β-diversity only measuring the taxonomic dimension essentially consider species as equally distinct, thus providing little information about functional structure of species assemblages (Villéger et al., 2008). To understand how biotic communities assemble and functionally respond to environmental changes, it is important to take the diverse ecosystem functions performed by species into consideration (Cadotte et al., 2011). For this reason, many metrics have been proposed over the recent decades to measure functional diversity (Cadotte et al., 2011;Mouchet et al., 2010;Petchey & Gaston, 2006). Analogous to taxonomic βdiversity, functional β-diversity measures dissimilarity in functional structure among assemblages and can help disentangle how biotic communities functionally assemble along environmental gradients (Münkemüller et al. 2012;Pavoine & Bonsall, 2010). Likewise, functional β-diversity can be driven by functional turnover due to niche differentiation (functional dissimilarity) among communities, or result from nestedness of functional strategies due to different niche filtering intensity among communities (Bevilacqua & Terlizzi, 2020;Villéger et al., 2013). Therefore, partitioning functional β-diversity into turnover and nestedness components allows direct comparisons among multifaceted β-diversity, and can help uncover the underlying mechanisms for functional dissimilarity between communities (Si et al., 2016;Villéger et al., 2013).
Most β-diversity studies focus on breeding communities or sedentary species with low dispersal capacities (Chen et al., 2020;Villéger et al., 2013;Wu et al., 2017), while few studies partition β-diversity among communities dominated by migratory species (but see Arruda Almeida et al. 2019). The waterbird communities in island-like subsidence wetlands in the North China Plain provide an ideal model to improve our understanding of β-diversity among vagile communities. They also provide an opportunity to understand how terrestrial communities assemble in the worldwide human-modified landscape with intensive anthropogenic disturbances (Ellis, 2011;Muller et al., 2017). Land subsidence resulted from massive underground coal mining in the North China Plain was estimated to be 1 × 10 6 ha in 2011, with an ongoing annual increase of 7 × 10 4 ha (Hu et al. 2014). Due to high groundwater levels and abundant rainfall, more than one third of the subsidence area has been waterlogged, generating hundreds of wetlands. Previous studies have found that these unintended man-made subsidence wetlands are increasingly used by a wide range of waterbird species, especially those migrating along the East Asian-Australasian Flyway . Environmental filtering is suggested to result in functional clustering in these waterbird communities, and the functional α-diversity is highly affected by environmental variables . These assemblages are also significantly nested, which could be attributed to selective extinction processes .
However, no studies have been done to decompose taxonomic and functional β-diversity to compare relative contributions of different components in these assemblages dominated by migrating bird species with high dispersal capacity.
In this study, we applied the framework proposed by Villéger et al. (2013) to the waterbird communities across subsidence wetlands in the North China Plain to quantify the taxonomic and functional β-diversity, and compare their turnover and nestedness components. We expected that (1) taxonomic and functional diversity metrics (α-and β-diversity, as well as their corresponding components) would be positively correlated because more species may hold more functional strategies, (2) functional β-diversity would be higher than taxonomic β-diversity because local environmental filtering may result in wetland-specific functional convergence, and wetlands with various environmental conditions may support communities with different functional strategies , Editor: Zhixin Zhang conditions, changes in the multiple biodiversity metrics should be monitored to inform the ongoing management plans.

K E Y W O R D S
functional β-diversity, nestedness, non-random assembly, subsidence wetland, turnover, waterbirds thus increasing functional dissimilarity among communities, and (3) turnover would dominate taxonomic β-diversity because high movement capacity facilitates different species to inhabit in various environments via species sorting, while nestedness would dominate functional β-diversity due to directional loss of functional strategies filtered by different environmental conditions (Soininen et al., 2017). Furthermore, we simulated assemblages under random processes in null models with the same species richness and pairwise taxonomic β-diversity between communities as the observed, and compared observed functional α-and β-diversity with the simulated counterparts. Given the possible role of environmental filtering in the assembly and the resulting functional convergence , we expected that 4) observed functional α-diversity would be lower than simulated means. Because the different environmental conditions in subsidence wetlands may act as different filters to favour different clusters of functional strategies, we also expected that 5) most of observed pairwise functional β-diversity values would be higher than simulated means. Based on the above analyses, we aimed to understand the patterns of taxonomic and functional β-diversity and the role of non-random process, for example environmental filtering, in shaping the waterbird diversity in the subsidence wetlands, providing insights for the conservation.

| Study area
The Huainan-Huaibei coal mining area (~1.5 × 10 6 ha, 32°44′-33°44′N, 116°02′-117°31′E), located to the north of Huai River in the North China Plain, is one of the 14 largest coal bases in China (Figure 1 in . The average elevation of the area is ~30 m above sea level, and the groundwater level is ~1.5 m. The climate is warm temperate monsoon type with an average annual temperature of 15℃ and an average annual precipitation of 970 mm. The area is dominated by flat agricultural landscape under which there has been a history of more than 100 years of coal mining. The massive underground coal extraction has caused large-scale ground subsidence with an area of 6.2 × 10 4 ha by the end of 2016 .
Nearly half of the subsidence area has been waterlogged to an average water depth of 3-6 m, creating more than one hundred wetlands of different sizes and landscape connectivity. These wetlands were formed in different years, and many of them are still expanding because of the continuing underground coal mining.

| Bird data
We used point-count method to survey waterbirds, that is species of birds that are "ecologically dependent upon wetlands" according to the Ramsar Convention (Gardner & Davidson, 2011), at 55 randomly selected subsidence wetlands approximately every two weeks from September 2016 to April 2017. Bird counts were carried out by two skilled bird observers in clear and sunny days along the boundary of each wetland at 1-6 counting points depending on the wetland area and accessibility. The observation area at each counting point had a radius of 1 km and did not overlap with others to avoid double counting. Waterbirds occurring within each observation area were recorded and identified to species following MacKinnon et al. (2000) within 15 min with the help of binoculars (10 × 42 WB Swarovski) and a telescope (20-60× zoom ATM 80 Swarovski). Birds flying over the observation areas from outside were not recorded. Because of the relatively small wetland area and good vision, bird counting at each wetland was thorough, which is important for quantifying assemblage composition (Si et al., 2018).
Bird observations were categorized into three seasons, that is autumn migration (September-November 2016), wintering season (December 2016 to mid-February 2017) and spring migration (late February-April 2017), with each season including four surveys and each survey covering all the 55 wetlands. To avoid the influence of vagrant species rarely using the wetlands on our analyses, we excluded species with less than three records (2 in autumn, 5 in winter and 4 in spring). Wetlands with less than four species during each season were also excluded because there should be at least four species in each community when defining its functional convex hull in the following. Therefore, a total of 46 wetlands and 51 species of presence-absence data were included in the following analyses. We defined a waterbird community in each wetland as the pool of species recorded in the four surveys during each season. For these 46 wetlands (thus, 46 communities), we calculated pairwise taxonomic and functional β-diversity metrics in the following for all the possible pairs (i.e. 1,035 pairs).

| Functional traits
Consistent with , we selected body mass, main food type, foraging method and substrates as the functional traits of waterbirds (Table 2 in . The selected traits are related to birds' resource use, energy requirements and trophic levels, and are thus commonly used to measure functional diversity (Flynn et al., 2009;Petchey et al., 2007). These trait data were collected from del Hoyo et al. (2004( -2011( ), del Hoyo et al. (1992( -2002 and the BTO database (http://www.bto.org/about birds/ birdf acts). Following Villéger et al. (2013), we used the Sørensen dissimilarity index to measure pairwise taxonomic β-diversity (β sor , Eqn 1 in Appendix S1) as the percentage of dissimilarity in species composition between two communities. The pairwise taxonomic β-diversity was then partitioned into turnover (β sim , Eqn 2 in Appendix S1) and nestedness components (β sne , Eqn 3 in Appendix S1) following Baselga (2012). Baselga (2012) suggests that the mean of pairwise β-diversity would underestimate the regional β-diversity. Therefore, we used multiple-site β-diversity to measure the real overall dissimilarity in species composition among all the studied communities. The multiple-site taxonomic β-diversity (β SOR , Eqn 4 in Appendix S1) was calculated using all the 46 wetlands in each season and was then decomposed into turnover (β SIM , Eqn 5 in Appendix S1) and nestedness components (β SNE , Eqn 6 in Appendix S1) following Baselga (2010, 2012).

| Functional β-diversity and partitioning
To measure functional diversity, we first used the R package FD to calculate functional dissimilarity between all pairs of species using the Gower's distance, which allows a mix of both categorical and continuous traits in our trait data (Gower, 1971). Then, a principal coordinate analysis (PCoA) was performed on the functional distance matrix (Laliberté & Legendre, 2010;Villéger et al., 2008), and the first three axes were used as coordinates to define a multidimensional functional space. Species in the functional space were located by their respective synthetic traits, and the functional richness (functional α-diversity) of a community was defined as the proportion of the convex hull hypervolume filled by all species in each season occupied by one community (Cornwell et al., 2006;Villéger et al., 2008).

| Statistical analysis
We firstly evaluated the quality of the functional space following Maire et al. (2015). Specifically, we used Mantel's test to investigate the correlation between species' Euclidean distances in the threedimensional functional space and the Gower's functional distances to assess the effect of PCoA-induced information loss on the functional distances between species. The following statistical analyses were focused on α-diversity and pairwise β-diversity.
Repeated-measures ANOVA performed by the package nlme and post hoc Tukey's test with a Bonferroni correction was used to compare the differences in taxonomic and functional β-diversity among seasons. The correlation between taxonomic and functional α-diversity was assessed using the Pearson correlation coefficient and significance. Mantel's test was used to investigate correlations between taxonomic and functional β-diversity, as well as between their respective components. To remove any possible effects of spatial autocorrelation, we also used partial Mantel's test to assess the correlations while controlling for the effect of spatial distance between communities. In the Mantel's and partial Mantel's tests, Pearson's correlation coefficients were calculated, with significance tested using randomization procedures (9,999 iterations in this study). The paired t test with a Bonferroni correction was used to compare the taxonomic and functional β-diversity, as well as the respective turnover and nestedness.
To investigate whether the observed functional α-and β-diversity resulted from a random assembly process, we followed Villéger et al. (2013) to compare the observed metrics with randomly simulated means in null models. Firstly, we subsampled species for each wetland, given its species number, from the regional species pool during each season without replacement. For each wetland, 999 null communities were simulated and functional α-diversity was calculated.

F I G U R E 1
The observed and simulated functional α-diversity under random assembly processes in null models maintaining the number of species in each wetland. The 95% confidence intervals of the 999 simulated values for each wetland are shown by the grey bars. The observed functional α-diversity values lower and higher than the simulated are displayed as red and blue circles, respectively. The linear correlations between observed taxonomic and functional α-diversity are shown by blue lines with 95% confidence intervals TA B L E 1 Taxonomic and functional β-diversity and the respective turnover and nestedness components of the waterbird assemblages across the subsidence wetlands in the Huainan- of species unique to community 1 and 2, respectively, while S 12 is the number of species shared by the two communities. For functional metrics, V 1 and V 2 are the volumes of convex hull unique to community 1 and 2, respectively, while V 12 is the volume of the intersection between the two convex hulls

Huaibei coal mining area in China
The t test with a Bonferroni correction was used to compare the observed and the simulated means, which would not be different under random assembly process. If the observed values were significantly lower than the simulated means, we concluded that the assemblages were functionally clustered. In contrast, the assemblages would be overdispersed in the multidimensional functional space.
Secondly, for each of the 1,035 pairs among the 46 surveyed wetlands, we simulated 999 pairs of random communities with the same number of species unique to each community and shared by the paired two communities as the observed pairs. The identities of species were randomly subsampled without replacement from the regional pool.
This procedure maintained the pairwise taxonomic β-diversity and its turnover and nestedness components for each pair of surveyed wetlands. Then, pairwise functional β-diversity and its turnover and nestedness components were calculated for each of the 999 simulations, and its mean value was compared with the observed counterparts using the t test with a Bonferroni correction. The chi-square goodness-of-fit test was used to determine whether the proportion of wetlands with lower observed functional α-diversity or the proportion of wetland pairs with higher pairwise functional β-diversity, turnover and nestedness component was significantly larger than 50%.
Taxonomic and functional β-diversity and their turnover and nestedness components were calculated using the R package betapart (Baselga & Orme, 2012). All analyses were carried out using R 3.6.2 (Appendix S2; R Development Core Team, 2018). Statistical significance level was p < .05, and data are displayed as mean ± standard error.

| Taxonomic and functional richness
We recorded a total of 51 waterbird species (

| Taxonomic and functional β-diversity
Both multiple-site taxonomic and functional β-diversity were similar among the three seasons and were largely contributed by the turnover components ( The pairwise functional β-diversity varied from 0.002 to 1 and was highest in spring with no difference between autumn and winter (F 2, 2068 = 6.16, p = .002). Pairwise functional β-diversity and its nestedness component were significantly higher than the corresponding taxonomic indices, while functional turnover was lower than taxonomic turnover in all the seasons ( Table 1). Functional turnover was significantly lower than functional nestedness (autumn: t = −6.42, p < .001; winter: t = −8.49, p < .001; and spring: t = −5.23, p < .001). The functional β-diversity and the two components were positively correlated with the taxonomic counterparts, no matter controlling for the effect of spatial distance between wetlands ( Table 2 and Figures 2-4).

| Observed and randomly simulated functional diversity
The observed functional α-diversity was significantly lower than the randomly simulated mean in most of the 46 surveyed wetlands in all the three seasons ( Table 3 and Figure 1). The proportion of wetland pairs with higher observed pairwise functional β-diversity than random was significantly larger than 50% in all the three seasons ( Figure 2).
The proportion for functional turnover was larger than 50% in autumn and spring, while the proportion for functional nestedness was not larger than 50% in any seasons ( Table 3, Figures 3 and 4). (NT) . We found that multiple-site taxonomic β-diversity and functional β-diversity were similar and relatively high among the waterbird communities across the subsidence wetlands in the non-breeding seasons. The overall compositional dissimilarity and functional dissimilarity were primarily caused by taxonomic and functional turnover, respectively. The calculation used here measured the real overall dissimilarity across all the studied communities, rather than the mean pairwise dissimilarity used in some studies (Si et al., 2016). Because comparisons between studies should involve the same number of assemblages (Baselga, 2010), we focus on pairwise β-diversity in the following.

F I G U R E 2
The observed and simulated pairwise functional β-diversity under random assembly processes in null models where pairwise taxonomic β-diversity (and the two components) was kept constant. The 95% confidence intervals of the 999 simulated values for each pair of wetlands are shown by the grey bar. The observed functional metrics lower and higher than the simulated are displayed as red and blue circles, respectively. The linear correlations between observed taxonomic and functional β-diversity are shown by blue lines with 95% confidence intervals

| Turnover and nestedness components of β-diversity
Although a previous study found nestedness among the studied waterbird assemblages (Li, Zhao, et al., 2019), turnover contributed more to the pairwise taxonomic dissimilarity, which is consistent with the general pattern found in a variety of taxa and ecosystems (Soininen et al., 2017). In contrast, the contribution of nestedness to functional β-diversity was slightly more than that of turnover ( Table 1). The contrasting pattern of taxonomic and functional components was also found by Villéger et al. (2013) and Si et al. (2016) who suggested that the frequent taxonomic turnover was mainly driven by functionally redundant species. On the one hand, the higher taxonomic turnover indicates that species replacement was prominent among the communities with similar species richness but different taxonomic compositions (Baselga, 2010 (Logez et al., 2010). Because habitat diversity varies significantly among the wetlands , we speculated that the functional nestedness in the metacommunity might arise due to different niche filtering intensity among the wetlands. The higher functional nestedness might also be linked to the selection of functional traits, which may largely determine the functional metrics (Petchey & Gaston, 2006). Compared with categorical variables, continuous variables (e.g. body size) may decrease the probability of higher functional turnover because they will always be present in all the species (e.g. a species cannot have 0 value of body size).

| Observed and randomly simulated diversity
The null models showed that observed functional α-diversity in most wetlands was lower than the randomly simulated mean given the F I G U R E 4 The observed and simulated pairwise functional nestedness under random assembly processes in null models where pairwise taxonomic β-diversity (and the two components) was kept constant. The 95% confidence intervals of the 999 simulated values for each pair of wetlands are shown by the grey bar. The observed functional nestedness lower and higher than the simulated is displayed as red and blue circles, respectively. The linear correlations between observed taxonomic and functional nestedness are shown by blue lines with 95% confidence intervals regional species pool and species richness in each wetland (Figure 1).
The lower functional α-diversity may be resulted from functional convergence due to non-random assembly processes, for example environmental filtering, which simplify functional strategies in a given community (Cadotte & Tucker, 2017;Concepción et al., 2017).
This process is analogous to selective extinction where species unable to tolerate certain environmental conditions will be lost, resulting in co-occurrence of functionally similar species (Si et al., 2016;Wang et al., 2011). Under random assembly processes, pervasive dispersal would facilitate homogeneity of species composition among communities, generating low taxonomic and functional βdiversity (Weinstein et al., 2014). However, in our metacommunity dominated by highly vagile migratory waterbirds, observed functional dissimilarity between most community pairs was higher than the expected in the null models ( Figure 2). The functional dissimilarity among the communities might be strengthened by their different environmental conditions, which can act as various filters to favour species with different traits, resulting in different clusters of functional strategies suitable with certain wetlands (Logez et al., 2010;Si et al., 2016). This pattern provides further evidence for the important role of non-random processes, for example environmental filtering, in the assembly of the waterbird communities. The rationale can also explain why the waterbird communities differed more in functional strategies than in taxonomic compositions in all the three seasons. Given that species differ in their functional strategies, it is reasonable to expect that an increase in taxonomic richness would result in an increase in functional space of a community (Si et al., 2016;Villéger et al., 2013). Despite the positive correlation, environmental filtering impedes the inflation of functional space in communities with increasing number of species, and thus increases pairwise dissimilarity due to environment-specific cluster of species' functional strategies (Concepción et al., 2017). Although most observed functional α-diversities were lower and most β-diversities were higher than the simulated values in null models indicating the dominance of non-random assembly, there were also some higher observed functional α-diversities and lower β-diversities indicating minor roles of random processes in assembly of the waterbird communities in subsidence wetlands (Si et al., 2016).

| Implications for management and conservation
Disentangling the relative contributions of turnover and nestedness to β-diversity patterns has important implications for conservation plans to preserve regional biodiversity (Angeler, 2013). Specifically, if nestedness dominates β-diversity, sites with high α-diversity deserve priority to be protected because complementarity among communities is low. On the contrary, when dissimilarity among communities is mainly caused by turnover, multiple sites should be protected. To protect waterbirds in the subsidence wetlands in the North China Plain, the higher taxonomic turnover implies that most of the wetlands, especially those holding higher functional richness, should be protected. Given the important role of non-random processes, such as environment filtering, in the assembly of the waterbird communities, habitat diversity in these subsidence wetlands should be enhanced to reduce niche filtering intensity among the wetlands, thus providing more niches for various waterbird species.
Furthermore, given the changing environment in such anthropogenic subsidence areas, it is important to monitor changes in multiple biodiversity metrics in the waterbird assemblages to inform the ongoing management plans.

ACK N OWLED G EM ENTS
The study was financially supported by the National Natural Science

CO N FLI C T S O F I NTE R E S T
The authors declare no conflicts of interest.

PEER R E V I E W
The peer review history for this article is available at https://publo ns.com/publo n/10.1111/ddi.13396.

DATA AVA I L A B I L I T Y S TAT E M E N T
All the raw data are provided as Supporting Information.