Assembly processes of waterbird communities across subsidence wetlands in China: A functional and phylogenetic approach

Although assembly processes have been studied in a wide range of taxa, determining assembly rules remains controversial, particularly in assemblages consisted of species with strong dispersal capacities. Moreover, few studies focused on communities in recently human‐created habitats. We tested two prevailing but opposing hypotheses, environmental filtering and limiting similarity, in waterbird communities across subsidence wetlands created by underground coal mining in China, with an aim to better understand assembly processes in communities composed of highly mobile species in human‐dominated landscape.


| INTRODUC TI ON
How biotic communities are assembled from regional species pool is one of the basic questions in ecology (Cornwell & Ackerly, 2009).
In spite of other possible assembly processes (Weiher et al., 2011), two prevailing but opposing mechanisms are environmental filtering (Cornwell, Schwilk, & Ackerly, 2006) and limiting similarity (Chesson, 2000;MacArthur & Levins, 1967). Environmental conditions may act as filters selecting species with similar traits to survive the same local environment, resulting in a larger ecological similarity and closer phylogenetic relatedness among coexisting species than would be expected by chance (i.e., functional or phylogenetic clustering; Mouchet, Villeger, Mason, & Mouillot, 2010). Alternatively, niche theory highlights the importance of interspecific competition in assembly processes, and thus, species with similar ecological characteristics or phylogenetic relatedness are prevented to stably coexist (i.e., functional or phylogenetic overdispersion; MacArthur & Levins, 1967;Mouchet et al., 2010).
In the progress of exploring assembly mechanisms in biotic communities, multifaceted metrics have been increasingly advocated to quantify assemblage structures, that is, functional (FD) and phylogenetic diversity (PD), in addition to traditional measures of biodiversity such as species richness (SR) and taxonomic diversity (TD; Gómez, Bravo, Brumfield, Tello, & Cadena, 2010;Mouchet et al., 2010;Tucker et al., 2017). When using SR and TD to measure biodiversity, species are considered equally distinct, providing little information about ecosystem functions and phylogenetic relatedness of the co-occurring species (Sodhi & Ehrlich, 2010;Tscharntke et al., 2012). However, FD and PD consider species-specific functional traits and evolutionary histories, respectively, and may therefore lead to a better understanding of mechanisms driving assembly processes. FD measures variability in species traits linked to ecosystem functions and thus should provide insights into how species coexist along gradients of environmental variables (Barbaro, Giffard, Charbonnier, van Halder, & Brockerhoff, 2014;Mouchet et al., 2010).
PD measures evolutionary differences among species in a community (Srivastava, Cadotte, MacDonald, Marushia, & Mirotchnick, 2012), bridging gaps between ecological processes with long-term evolutionary potential of species (Sandel, 2018). These metrics are complementary, quantifying multiple dimensions of biodiversity, and therefore facilitate a comprehensive understanding of assembly processes (Pavoine & Bonsall, 2011).
Although assembly processes have been studied in a wide range of taxa from bacteria to plants and animals, the topic is still controversial (Pavoine & Bonsall, 2011;Stevens, Cox, Strauss, & Willig, 2003;Zak, Willig, Moorhead, & Wildman, 1994). Most of the studied communities comprise relatively sedentary species with low dispersal capacities. However, few studies have focused on assemblages dominated by taxa with high dispersal capacities, such as migratory waterbirds, that may assemble differently from sedentary taxa (Mendez et al., 2012;Wiens, 1992). Waterbirds with various biological and ecological characteristics perform a wide array of ecosystem functions, such as nutrient cycling, ecosystem engineering, and dispersal of seed, pathogens and invertebrates (Green & Elmberg, 2014;Sekercioglu, 2006). They are also one of the most conspicuous taxa sensitive to wetland degradations and are often used as bioindicators of environmental changes (Green & Elmberg, 2014). Therefore, waterbirds are an ideal species group to explore assembly processes and assess effects of habitat changes on biotic communities. Moreover, in the context of global losses and degradations of natural wetlands (Butchart et al., 2010;Kar, 2013), waterbirds have been increasingly found to use artificial wetlands as complementary habitats, such as paddy fields, aquaculture ponds and water reservoirs (Elphick, 2015;Navedo et al., 2012;Petchey, Evans, Fishburn, & Gaston, 2007). An understanding of waterbird responses to environmental variables and assembly processes in human-dominated landscape can increase the success of management and conservation actions of waterbirds (Thompson & Starzomski, 2007;Tscharntke et al., 2012).
Subsidence wetlands created by large-scale underground mining are unintended man-made wetlands, overlooked by the conservation community, but they harbour a large variety of flora and fauna species (Lewin, Spyra, Krodkiewska, & Strzelec, 2015;Nawrot, Kirk, & Elliott-Smith, 2003;Townsend et al., 2009). China is rich in coal resources, and the annual production is approximately 50% of the worldwide total (Dong, Samsonov, Yin, Yao, & Xu, 2015). The massive underground coal mining has resulted in extensive land subsidence, estimated to be 1 × 10 6 ha as of 2011, with an annual increase of 7 × 10 4 ha (Hu et al., 2014). Because of high groundwater levels and abundant rainfall, subsidence land has been flooded in the North China Plain, one of the most important coal production areas in China. The newly created artificial wetlands are still expanding and have been found to be used by a wide range of waterbird species, especially by long-distance migrants from the East Asian-Australasian importance of environmental filtering and habitat variables in structuring assemblages dominated by species with high dispersal capacities and suggests that increasing habitat diversity and reducing disturbances will contribute to waterbird conservation in this human-dominated landscape.

K E Y W O R D S
assembly process, environmental filtering, functional diversity, phylogenetic diversity, subsidence wetland, waterbird assemblage Flyway (Li, Zhao, & Wang, 2019). These artificial wetlands have clear boundaries and vary in spatial complexity, providing good chances to explore assembly processes of communities dominated by migratory species in a human-modified landscape.
In this study, we quantified taxonomic, functional and phylogenetic diversity of waterbird communities in subsidence wetlands in the North China Plain. We tested for a seasonal effect and for independent contributions of multi-scale habitat variables to variations in the multifaceted biodiversity metric data. Because the waterbird communities are dominated by migratory species that can rapidly respond to environmental changes (Li et al., 2019), we hypothesized that the assemblages should predominantly be structured by environmental filtering. We expected stronger effects of habitat variables on functional structures on which environmental filtering directly operates. We also expected a shift from functional or phylogenetic clustering to overdispersion in the assemblages with increasing habitat diversity, as heterogeneous habitats provide more ecological niches to be partitioned. With these analyses, we contribute to a better understanding of assembly processes in biotic communities, exemplified by waterbird communities dominated by migratory species in these expanding man-made wetlands.

| Study area
The study was carried out in the Huainan-Huaibei coal mining area, which occupies an area of ~1.5 × 10 6 ha in the North China Plain producing ~4.1% of the national output (Hu et al., 2014). The massive underground mining has created large-scale ground deformation and land subsidence. The total subsidence area was more than F I G U R E 1 Subsidence wetlands surveyed for waterbirds in the Huainan-Huaibei coal mining area in China 30,000 ha in 2010, with an annual expansion of over 2,000 ha (Xie, Zhang, Yi, & Yan, 2013). Nearly half of the subsidence area has been flooded because of high groundwater levels and abundant rainfall, creating hundreds of independent wetlands of different sizes. These subsidence wetlands were created in different years, and many of them are still expanding due to the ongoing underground mining.
They differ in a wide range of environmental attributes, which may influence biotic communities formed after waterlogging. These man-made wetlands are located within the East Asian-Australasian Flyway and have been found to harbour a wide array of migratory waterbird species (Li et al., 2019).

| Bird data
From September 2016 to July 2017, we carried out a total of 13 point-count waterbird surveys, each covering the same 55 subsidence wetlands. The total sampled area was 6,226 ha, accounting for ~40% of the subsidence wetlands in this region. These wetlands were selected randomly, representing a wide range of environmental conditions. Depending on the wetland area and accessibility, we placed 1-6 counting points along its boundary with unobstructed views of each wetland. We defined the areas within a radius of 1 km at each counting point as observation areas, and these areas were not overlapping to avoid double counting. The relatively small areas and clear boundaries of the wetlands guaranteed good detection of waterbirds, which is important when quantifying assemblage structures (Si et al., 2018).
Each field survey was completed within three clear and calm days, and the "look-see" total counting method (Delany, 2005) was employed by the same two experienced bird observers (CL & SY) to record waterbird species and abundance at each counting point.
Waterbirds occurring within each observation area were identified to species level following Jetz, Thomas, Joy, Hartmann, and Mooers (2012), within 15 min, with the help of binoculars (10 × 42 WB Swarovski) and a telescope (20-60 × zoom ATM 80 Swarovski). We defined a waterbird species as a species that is "ecologically dependent upon wetlands" according to the Ramsar Convention (Gardner & Davidson, 2011). We only recorded waterbirds that used the sampled wetlands, ignoring those flying over. Species with less than three records during all surveys were excluded from the following analyses.
According to the migration chronology of the waterbirds, we divided the field surveys into four periods: autumn migration (September-November 2016), wintering season (December 2016 to mid-February 2017), spring migration (late February-April 2017) and summer breeding season (May-August 2017). The season division allowed us to better understand the variations in the multiple diversity metrics under influence of assembly processes that vary among seasons. During each nonbreeding period, there were four surveys with breaks of minimally two weeks in between. These surveys were pooled in the analyses by summing the abundances of each species recorded in each wetland.
We carried out only one survey in the summer breeding season (in July 2017). Therefore, we had four species × sites matrices, one in each period. Here, we considered a community as a pool of waterbird species co-occurring in a given wetland in each period.

| Habitat variables
We measured eight habitat variables at both local scale and landscape scale in each sampled wetland (Table 1) where p i is the proportion of the wetland area occupied by the ith of n habitat types, that is, open water, aquatic vegetation and ground (Simpson, 1949). We also calculated Wetzel's (1975) shape index (SW) to quantify the irregularity of each wetland, as the amount of deviation from a perfect circle. There were no significant changes in habitat variables during field surveys, and thus, the same variables × sites matrix was used in all analyses.

| Functional traits
Functional diversity may depend on arbitrary decisions on the selected traits, and there is no consensus on the selection of species characteristics to capture all functions of species. To facilitate meaningful comparisons among studies, the same suite of traits is suggested to be selected (Petchey & Gaston, 2006). We followed Petchey et al. (2007) and Jia, Wang, Zhang, Cao, and Fox (2018) to select four traits to measure the functional diversity of waterbird assemblages, comprising one continuous and three categorical attributes (

| Biodiversity metrics
For each wetland during each period, we defined species richness (SR) as the total number of species and calculated biodiversity metrics, in the R packages ape and picante (Kembel et al., 2010;Paradis & Schliep, 2018), representing the divergence dimension for taxonomic (TD), functional (FD) and phylogenetic diversity (PD). TD was quantified using the Simpson index (Simpson, 1949) under the mathematical framework of Rao's quadratic entropy (Rao, 1982), which is defined as.
where d ij is the distance between species i and j in a community containing a total of S species. The distances were weighted by the relative abundances of each species (p i and p j backbone (Hackett et al., 2008;Jetz et al., 2012), and summarized these trees to generate a 50% majority rule consensus tree using SumTrees (Sukumaran & Holder, 2010). Based on the functional dendrogram and phylogenetic tree, respectively, we calculated two met- pairwise functional or phylogenetic distances (branch lengths on the functional or phylogenetic dendrogram) among co-occurring species, representing an overall divergence of the community. MNTD quantifies the mean functional or phylogenetic distance between nearest neighbours, and it can describe the degree of terminal clustering among co-occurring species (Webb, 2000).
To examine whether co-occurring species in communities were clustered or overdispersed on the functional dendrogram or phylogenetic tree, we compared the observed MPD and MNTD in each wetland during each period with the corresponding mean value generated by null models by shuffling the tip labels of the trees (Swenson, 2014). In this process, 999 communities with species richness and occurrence frequencies equal to the observed communities were randomly generated using an independent swap algorithm, and metrics were calculated and averaged (Gotelli & Entsminger, 2001).
The species pool used in the simulation process was defined as all species recorded during the given period.

| Phylogenetic signal
The use of PD to estimate ecosystem functioning of species depends highly on the strength of phylogenetic signals in functional traits (Srivastava et al., 2012). Phylogenetic signals measure the statistical dependence among species' traits associated with their phylogenetic relatedness, and significant phylogenetic signals indicate higher similarity than expected by chance among closely related species (Revell, Harmon, & Collar, 2008). We used the statistic D (Fritz & Purvis, 2010) to measure phylogenetic signals in categorical traits, with lower D values indicating a more conserved trait evolution (i.e., stronger phylogenetic signal). If D approaches 0, the trait is distributed as expected under the Brownian motion model of evolution (i.e., conserved trait evolution), and D < 0 suggests a highly clustered trait. A value of D = 1 or > 1 indicates that the trait is randomly distributed (i.e., no signal) or overdispersed on the phylogenetic tree.
We used the Pagel's λ (Freckleton, Harvey, & Pagel, 2002;Pagel, 1999) to measure phylogenetic signal in body mass, a continuous trait. A value of λ = 0 indicates no correlation between evolution in body mass and phylogeny. If body mass evolves according to Brownian motion along a phylogeny, λ will be equal to 1. A value of λ between 0 and 1 indicates that body mass has evolved according to a process in which the effect of phylogeny is weaker than in the Brownian model (Pagel, 1999

| Statistical analyses
The package nlme was used to perform repeated-measures ANOVA to test for seasonal differences in SR, TD, MPD and MNTD (Pinheiro, Bates, DebRoy, & Sarkar, 2018). Post hoc Tukey's test with a Bonferroni correction was used to compare metrics among the four periods. We used one-sample t test to determine whether NRI and NTI were significantly different from 0, which would be expected by chance. We employed a hierarchical partitioning method (Chevan & Sutherland, 1991;Mac Nally, 2002) to quantify independent contributions of habitat variables towards variations in NRI and NTI. In hierarchical partitioning, models containing all possible combinations of explanatory factors were considered in a hierarchy and goodness of fit (R 2 ) was calculated for each model. A randomization procedure with 1,000 iterations was performed to test the significance of the effect of each variable.
The package hier.part was used to perform the hierarchical partitioning (Mac Nally, 2002;Nally & Walsh, 2004); all analyses were carried out in R 3.3.4.

| RE SULTS
During the 13 field surveys, we recorded a total of 62 waterbird species across all sampled wetlands. After discarding species with less than three records, we obtained 51 species (i.e., 48 species in autumn, 44 in winter, 47 in spring and 18 in summer) which were used in the following analyses. Phylogenetic signals were significantly related to body mass and most categorical traits, indicating strong phylogenetic niche conservatism (Table 2).
Species richness (SR) was on average highest in autumn and lowest in summer with no difference between winter and spring (F 3,150 = 55.67, p < 0.001). Taxonomic diversity (TD) was on average higher in autumn but no differences were found among other periods (F 3,150 = 6.00, p < 0.001). Phylogenetic MPD was on average highest in autumn and lowest in summer, but there were no significant differences between autumn and spring, between winter and spring, and between winter and summer (F 3,150 = 8.94, p < 0.001).
Phylogenetic MNTD was highest in summer with no differences among other periods (F 3,150 = 17.30, p < 0.001). Functional MPD was highest in autumn and spring, and lowest in summer (F 3,150 = 14.77, p < 0.001), while functional MNTD was higher in winter and spring than those in autumn and summer (F 3,150 = 8.69, p < 0.001; Figure 2).
Phylogenetic NRI was greater than 0 during all four periods.
Phylogenetic NTI was greater than 0 in winter and spring, but not different from 0 during the other two periods. Functional NRI was not different from 0 in spring, but greater than 0 during the other periods. Functional NTI was greater than 0 in summer, but not different from 0 during the other three periods (Figure 3).
(2) Standardized effect size = M null − M obs ∕SD null F I G U R E 2 Seasonal differences in multiple biodiversity metrics of waterbird communities in subsidence wetlands in the Huainan-Huaibei coal mining area in China. Boxplots with the same letter indicate no significant difference as determined by post hoc Tukey's tests after the repeated-measures ANOVA. MPD: the mean pairwise distance; MNTD: the mean nearest taxon distance; SR: species richness; TD: taxonomic diversity The proportions of variation explained by the variables with significant effects ranged from 10.4% to 53.4%. The selected habitat variables had stronger effects on the standardized effect sizes of functional diversity than on those of phylogenetic diversity (Table 3).
We only found a negative effect of habitat diversity (HD) on phylogenetic NRI in autumn, and a negative effect of wetland age (AG) on phylogenetic NTI in winter. Wetland age (AG) was a positive predictor for functional NRI and NTI in all periods. Functional NTI in winter was also positively affected by wetland area (AW), area of open water (AO) and distance to human settlements (DH), but negatively with area of aquatic vegetation (AA) in summer. As habitat diversity (HD) increased, functional NRI and NTI decreased in all periods except in winter. Total area of wetlands within a 5-km buffer area (TA) was a negative predictor for functional NRI in winter, but a positive predictor for functional NTI in summer. There was no effect of wetland shape (SW) on any biodiversity metrics.
F I G U R E 3 Standardized effect sizes of functional and phylogenetic diversity and their 95% confidence intervals (along with p-values of one-sample t tests) of waterbird communities in subsidence wetlands in the Huainan-Huaibei coal mining area in China. Open triangle: phylogenetically nearest relative index (phylogenetic NRI); filled triangle: phylogenetically nearest taxon index (phylogenetic NTI); open circle: functionally nearest relative index (functional NRI); filled circle: functionally nearest taxon index (functional NTI) TA B L E 3 The independent contribution (%) of each habitat variable to variations of the multifaceted diversity metrics of waterbird communities in subsidence wetlands in the Huainan-Huaibei coal mining area in China Notes. Negative correlations are indicated by a minus sign, and significant relationships (p < 0.05) are in boldface. AG: wetland age; DH: shortest Euclidian distance from the boundary of each wetland to the nearest human settlement occupying an area >50 ha; AW: area of each wetland; AO: area of open water in each wetland; AA: area of aquatic vegetation in each wetland; HD: habitat diversity within each wetland; TA: total area of wetlands (>1 ha) within a 5-km buffer area surrounding each wetland; SW: Wetzel's (1975) shape index. NRI: nearest relative index; NTI: nearest taxon index.

| D ISCUSS I ON
We found that species diversity of waterbird communities in subsidence wetlands was on average highest during the autumn migration period (Figure 2). In the context of natural wetland loss along the East Asian-Australasian Flyway (Kirby et al., 2008), the newly created subsidence wetlands provide important complementary staging habitats for the migratory waterbirds, particularly during their southward migrations when there are more birds due to the post-breeding population increase (Li et al., 2019). Seasonal patterns of functional and phylogenetic diversity measured by MPD were similar to SR and TD, indicating that overall functional and phylogenetic divergence of the communities increased with species enrichment in autumn. However, functional and phylogenetic MNTD were rather low in autumn, and were quite different from the other metrics, providing evidence of clustering along the functional or phylogenetic dendrogram for the cooccurring species in our sampled wetlands. MPD and MNTD are two divergence metrics of biodiversity, measuring the overall divergence and terminal clustering, respectively, and may differently respond to species enrichment (Webb, 2000). The asynchrony of different biodiversity facets found in our study and in previous studies (Che et al., 2018;Monnet et al., 2014) justifies measurement of both MPD and MNTD when quantifying divergence in species assemblages.
The positive values of NRI and NTI indicated that these waterbird communities were composed of species with similar functional traits and close phylogenetic relatedness, suggesting the predominant importance of environmental filtering in the assembly processes ( Figure 3). Under this hypothesis, habitat variables act like filters to allow species with similar ecological traits to coexist, resulting in phylogenetic clustering with high niche conservatism (Webb et al., 2002). Our results contrasted with the expectation of the limiting similarity hypothesis that interspecific competition prevents coexistence of species with similar ecological niches or close relatedness (Chesson, 2000;MacArthur & Levins, 1967). Because of high dispersal capacities and rapid responses to environmental changes, waterbird communities should also be less likely determined by stochastic process or "priority effect," whereby pioneer species largely influence the community development (Fukami, 2015;Lok, Overdijk, & Piersma, 2013). The relative importance of environmental filtering and limiting similarity in assembly processes is hypothesized to be scale-dependent, exemplified by shifts from functional or phylogenetic clustering to evenness in some mammal and plant assemblages with decreasing geographical scales (Bryant et al., 2008;Cardillo, Gittleman, & Purvis, 2008;Cavender-Bares, Kozak, Fine, & Kembel, 2009). Contrary to this hypothesis, we found functional and phylogenetic clustering at local scales, like those found in Neotropical forest birds (Gómez et al., 2010). This may be attributed to the fact that the many waterbird species largely rely on wetlands, resulting in strong waterbird-habitat associations (Cao & Fox, 2009), and the effect of environmental filtering may therefore be stronger than interspecific interactions, even at local scales. Besides, migration phenology differs among species and their temporary interactions are not stable in a given wetland compared to those among resident taxa, particularly during migration periods, making limiting similarity less important in assembly processes (Kirby et al., 2008).
The assembly processes of these waterbird communities were mediated by habitat variables with different effects among seasons (Table 3). In general, the effects were stronger on the standardized effect sizes of functional metrics than on those of phylogenetic dimensions, suggesting that environment filtering on phenotypic traits were more prominent. Unlike the regional species pool, which is determined by broadscale and long-term colonization-extinction processes (Emerson & Gillespie, 2008;Ricklefs, 2006), species present in local communities vary in space and time, and select habitats in response to local ecological conditions. It is more prominent in assemblages dominated by migratory waterbirds, due to their large dispersal capacities and rapid responses to habitat changes. The environmental effects were also stronger on functional NTI than on functional NRI. Ecological processes may have very different effects on functional divergence depending on how they are measured. In contrast to functional NRI capturing the overall dissimilarity of the taxa, functional NTI is more sensitive to the distribution of lineages close to the tips of trees (Swenson, 2014). Therefore, NRI might be more influenced by factors that operated in the distant past, while NTI patterns are likely to reflect more recent changes (Che et al., 2018). This pattern provides further evidence for the importance of environmental filtering over a variety of temporal scales in assembly processes of assemblages comprised of mobile species.  (Li et al., 2019), providing limited opportunities for resource partitioning. Similar to heterogeneous habitats, wetlands with better landscape connectivity can attract more waterbird species (Che et al., 2018), resulting in functional overdispersion, particularly in winter. However, waterbirds may rely more on local environment during the breeding season when communities were functionally clustered (Wiens, 1992). Furthermore, the waterbird communities were more likely to be functionally clustered in wetlands far away from human settlements, indicating increased functional redundancy under lower disturbance levels. On the contrary, lower functional redundancy in more disturbed environments, such as reported in other empirical studies (Laliberté et al., 2010;Luck, Carter, & Smallbone, 2013), may increase vulnerability of communities to future disturbances (Naeem, 1998).
In conclusion, we found functional and phylogenetic clustering in waterbird communities in subsidence wetlands in the Huainan-Huaibei coal mining area, indicating the importance of environmental filtering in structuring assemblages dominated by species with high dispersal capacities. Habitat variables played important roles in the assembly processes, with stronger effects on functional diversity, which is more likely determined by ecological processes.
The patterns in multiple biodiversity metrics were not synchronous, reflected by different responses to habitat variables and seasons.
The assemblages somewhat shifted from functional clustering to overdispersion with increasing habitat diversity, which underlines the importance of environmental heterogeneity in maintaining functional and phylogenetic diversity. The subsidence wetlands provide complementary habitats for waterbirds migrating along the East Asian-Australasian Flyway; however, these waterbirds are facing a variety of threats that should be further studied (Li et al., 2019). In order to better manage and protect waterbirds in these man-made wetlands, we suggested (a) increasing habitat diversity and improving landscape connectivity within the wetland network, (b) reducing human disturbances and (c) systematically monitoring these waterbird communities along with the expansion of subsidence wetlands, with more emphasis on their early stages.

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