Local habitat and landscape attributes shape the diversity facets of bird communities in Inner Mongolian grasslands

The loss and fragmentation of natural habitats because of anthropogenic activities are major threats to biodiversity worldwide. In recent decades, vast mosaics of natural and seminatural habitats have been transformed into fragmented agricultural landscapes in Inner Mongolia, China, with potential negative effects on avian diversity. We quantified the effect of local and landscape habitat attributes on the taxonomic, functional, and phylogenetic diversity of bird communities in Inner Mongolian grasslands. We considered eight independent habitat variables obtained by variance inflation factor analysis. We used canonical correspondence analysis to determine how these habitat factors of multiple scales explained variance in species composition. We then fitted Bayesian generalized additive models to analyze the habitat-biodiversity relationships and included a smooth effect of land cover richness to test the intermediate disturbance hypothesis in each model. Our results provided evidence that differences in bird assemblages can be explained, in part, by differences in local and landscape-scale habitat features. The responses of the four diversity indices to these predictors were diverse and scale-dependent. We found species richness and Shannon diversity exhibited similar response, with both being negatively related to bare land percentage while being positively related to plant canopy and impervious surface percentage. Phylogenetic diversity was positively associated with plant richness while negatively associated with forest percentage and impervious surface percentage. We found no statistical evidence for a relationship between functional diversity and any of the variables examined here. Additionally, for the four measures of bird diversity, we did not find any evidence that they peaked at intermediate levels of habitat disturbance. We propose that assessments of regional grassland bird communities should be conducted at multiple scales and that a range of biodiversity metrics are required to better evaluate and inform conservation decision making, especially when the target is preserving not only species but also their evolutionary history and ecological functions. Les caractéristiques d'habitat aux échelles locale et du paysage façonnent la diversité des communautés d'oiseaux dans les prairies de la Mongolie-Intérieure RÉSUMÉ. La perte et la fragmentation des habitats naturels causées par les activités anthropiques représentent des menaces importantes pour la biodiversité dans le monde entier. Au cours des dernières décennies, de vastes mosaïques d'habitats naturels et semi-naturels ont été transformées en paysages agricoles fragmentés en Mongolie-Intérieure (Chine), entrainant possiblement des effets négatifs sur la diversité aviaire. Nous avons quantifié l'effet des caractéristiques d'habitat aux échelles locale et du paysage sur la diversité taxonomique, fonctionnelle et phylogénétique des communautés d'oiseaux dans les prairies de la Mongolie-Intérieure. Nous avons considéré huit variables indépendantes d'habitat obtenues par l'analyse des facteurs d'inflation de la variance. Nous avons réalisé une analyse canonique des correspondances pour déterminer de quelle façon ces variables d'habitat expliquaient la variance de la composition des espèces à plusieurs échelles. Nous avons ensuite bâti des modèles additifs généralisés bayésiens pour analyser les relations habitatbiodiversité et avons utilisé un lissage de la richesse de l'occupation du sol pour tester l'hypothèse de la perturbation intermédiaire pour chaque modèle. Nos résultats ont montré que les différences dans les assemblages d'oiseaux peuvent être expliquées, en partie, par les différences dans les caractéristiques d'habitat aux échelles locale et du paysage. La réponse des quatre indices de diversité à ces variables explicatives était variée et dépendait de l'échelle. Nous avons constaté que la richesse des espèces et la diversité de Shannon présentaient des réponses similaires, les deux étant négativement reliées au pourcentage de sol dénudé tout en étant positivement reliées au pourcentage de couvert végétal et de surface imperméable. La diversité phylogénétique était positivement associée à la richesse des plantes, mais négativement associée au pourcentage de forêts et de surface imperméable. Nous n'avons trouvé aucun lien statistique attestant d'une relation entre la diversité fonctionnelle et l'une des variables examinées. De plus, nous n'avons pas observé que les quatre mesures de la diversité d'oiseaux atteignaient un sommet à des niveaux intermédiaires de perturbation d'habitat. Nous pensons que les évaluations des communautés régionales d'oiseaux de prairie devraient être menées à des échelles multiples et qu'il serait nécessaire d'établir une série de mesures de la biodiversité pour mieux orienter la prise de décision en matière de conservation, en particulier lorsque l'objectif est de non seulement de préserver les espèces, mais aussi leur histoire évolutive et leurs fonctions écologiques.


INTRODUCTION
The conversion of natural habitats into anthropogenic habitats such as croplands or urban areas is one of the primary drivers of species extinction (Pimm andRaven 2000, Jetz et al. 2007). The loss and degradation of suitable habitats in farmed landscapes (Brennan andKuvlesky 2005, Askins et al. 2007) contributed to the fast decline in abundance and distribution of grassland birds across Europe (Donald et al. 2006 and North America (Rosenberg et al. 2019). A similar process of agricultural intensification has also occurred in Inner Mongolia, the main grassland region of northern China, where vast mosaics of natural and seminatural habitats have been transformed into fragmented agricultural landscapes in recent decades (Jiang et al. 2006). For example, formerly dominant grasslands there have become more fragmented, characterized by large blocks of cropland interspersed with smaller, more isolated grassland patches (John et al. 2009). In addition, nearly 90% of the grasslands are now degraded to varying degrees through overgrazing and early mowing of grass fodder fields (Xu et al. 2000), the primary productivity of which is only about 50% of that of the intact grasslands (Jiang et al. 2006). Furthermore, human settlements and planted woodlands (most of them as linear tree belts and windbreaks) fragment the remaining grasslands and create abrupt boundaries that exacerbate edge effects.
Such land cover changes are frequently associated with biodiversity decline. A decline in species richness has been widely observed across biogeographic regions (e.g., Jiguet et al. 2010, Correll et al. 2019). However, species richness alone cannot tell the complete story of an area's biodiversity; it provides only an estimate of the number of unique species present. Functional diversity (FD) was first defined as the value and range of species traits that influence the way ecosystems operate (Buchmann and Roy 2002), and expanded to measure the diversity of traits related to life history and functions (Petchey and Gaston 2002). Now it is widely used to describe the mechanistic link between species, their traits, and how they affect ecosystem resistance, resilience, and functioning (Petchey and Gaston 2006). Land cover change can decrease FD if anthropogenic habitats select species with specific functional traits. Phylogenetic diversity (PD) estimates cumulative evolutionary history by totaling the branch lengths in a community-wide phylogenetic tree (Faith 1992). PD is considered as a better measure of biodiversity, providing a comparable, evolutionary measure not possible with species counts (Miller et al. 2018). Communities with greater PD are predicted to be more resilient to environmental changes and to better preserve unique phylogenetic lineages (Vane-Wright et al. 1991, Srivastava et al. 2012. Land cover changes can greatly alter the degree of phylogenetic relatedness within and among assemblages if anthropogenic habitats select against evolutionarily unique species (Morelli et al. 2016, Sol et al. 2017. Indicators of functional and phylogenetic diversity summarize information on the different functional traits/roles and evolutionary history of the species within a community. Furthermore, areas of high species richness are not always congruent with areas of high FD or PD (Monnet et al. 2014), and therefore quantifying and contrasting all three measures could provide a better understanding of how land cover change affects biodiversity patterns. However, despite recent recommendations in nature conservation optimization encouraging the use of these additional diversity measures (Winter et al. 2013), previous studies investigating the effect of land cover changes on birds in Inner Mongolia focused exclusively on taxonomic diversity (e.g., Wen et al. 2002, Liang et al. 2018).
There is abundant evidence demonstrating that local habitat, e.g., vegetation structure, is critical in determining bird species composition and abundance in farmland (e.g., Andrén 1994, Warren et al. 2005, Sirami et al. 2009). However, it is unknown how well fragmented grassland landscapes function as habitat for birds (Bakker et al. 2002), though evidence suggests that areasensitive grassland birds require large, uninterrupted tracts of treeless grasslands (Andrén 1994, Pickett andSiriwardena 2011). Insights into how such bird species respond to habitat conditions at broader spatial scales would enhance our ability to direct grassland conservation over broad geographic regions and complement what has been learned at local scales (Naugle et al. 1999). Previous studies have revealed that environmental gradients can affect functional diversity but not necessarily in the same direction as it does for taxonomic diversity, depending on the spatial scale considered , Filippi-Codaccioni et al. 2010. To find optimal spatial scales of habitat management, we also need to understand at which scales biodiversity responds to environmental conditions. Understanding how species diversity is maintained has been among the foremost challenges for community ecologists (MacArthur and MacArthur 1961). One famous hypothesis of diversity maintenance under disturbance is the intermediate disturbance hypothesis (IDH). The IDH predicts patterns of maximum diversity under intermediate disturbance regimes, based on a tension between competitively superior species and species that can rapidly colonize following disturbance (Connell 1978, Shea et al. 2004). The IDH has been supported by many studies, such as in terrestrial (Molino and Sabatier 2001), freshwater (Padisák 1993), and marine communities (Sousa 1979). Some empirical studies, however, failed to observe the hump-backed diversity-disturbance relationship predicted by IDH (Shea et al. 2004, Fox 2013, and the ecology of disturbance with relatively mobile birds is likely different but largely unknown (Brawn et al. 2001). For example, Rueda-Hernandez et al. (2015) observed that medium-sized fragments had the highest avian species richness in a fragmented cloud forest landscape in Mexico, contrasting with a former study that found a greater number of species in larger fragments (Martínez-Morales 2005). These results suggest further studies are needed to validate the intermediate disturbance hypothesis. The mosaic of land cover types in Inner Mongolia provides an opportunity to examine the IDH on different diversity measures of bird community. One can logically expect taxonomic, functional, and phylogenetic diversity should be relevant with habitat disturbance gradients because variation in the response of bird composition and different diversity metrics are both likely to be affected by the specific combinations of species traits, acquired through evolutionary processes (Webb et al. 2002).
The aim of this study was to determine the influence of local and landscape attributes on the community composition and species diversity of grassland birds in Inner Mongolia. First, we examined the variation in the composition of species assemblages explained by local and landscape attributes using canonical correspondence analysis. Second, we explored the relationships between these habitat attributes and taxonomic, functional, and phylogenetic diversity of bird communities by fitting Bayesian generalized additive models. We addressed the following questions: (1) How do bird assemblages vary along local and landscape attributes?
(2) What are the independent effects of these habitat attributes on species diversity metrics? (3) Does the IDH apply for bird diversity metrics in our study site in Inner Mongolia? We predict that these avian diversity metrics should peak at intermediate levels of habitat fragmentation, i.e., land cover richness in our case.

Study area
We conducted this study in northeast Inner Mongolia, China (Fig.  1). The study site consists mainly of three vegetation zones: the temperate coniferous and deciduous forests zone, the meadow steppe zone, and the typical steppe zone (Han et al. 2009, Wu et al. 2015. The region has a temperate climate, with mean annual temperature and yearly precipitation ranging between 6 °C and 7 °C and between 300 and 400 mm, respectively (Wu et al. 2015). The study area covers a gradient of land cover ranging from natural grassland to intensive farming areas, the latter dominated by grain products, and includes scattered bare lands, tree plantations, and human settlements. Wheat and maize represent the main grain products in the region. Bare lands are generally areas with little vegetation, composed primarily of exposed soil, sand, gravel, and rock backgrounds, including recently harvested, fallow land, and all other types of land not covered by vegetation such as lake bottoms in dry season. Tree plantations mostly comprise commercial 3-to 10-year-old Populus and Sibirica planted primarily for lumber, fruits, and medicine production. Human settlements are usually small towns or cities holding 10,000-138,000 inhabitants, where buildings, houses, and paved roads increased in density from the rural to the core urban zone.

Bird counts
Bird assemblages were sampled by the point-count method with one visit (Bibby et al. 2000) between 01 May and 15 July 2018, which covers the overall breeding season of grassland birds in Inner Mongolia. Our study comprised 598 sampling points that we roughly assigned to one of four main land cover types: natural grassland (285 sites), farmland (126 sites), village (108 sites), and forest (79 sites). They were randomly distributed, but at least 50 m away from clear-cuts and large obstacles (such as urban or suburban developments and roads) and separated by at least 500 m. Four highly experienced ornithologists participated in the survey. The observers stood at a fixed location for a 10-min period recording all birds seen or heard within a 50-m radius. The 50-m limit was imposed to maximize detectability and decrease potential observer error, which can occur more frequently if attempting to identify cryptic species over long distances. Birds that were flying over were counted only if they were actually using the circular plot, such as for foraging, displaying, etc. Point counts were performed within 4 h after sunrise, avoiding days with fog, steady drizzle, prolonged rain, or winds greater than Beaufort 3 (13-19 km.h -1 ). During the count, the observers mapped the relative position of each individual bird to avoid double counting.

Landscape and local habitat attributes
We characterized local-scale vegetation coverage and structure using four main variables: plant canopy, plant richness, bare ground percentage, and vegetation height. We established three 2 x 2 m quadrats within a 10-m radius of the center of the pointcount station. For each quadrat, we visually estimated the proportion of bare ground (as the percentage of soil that is not covered by vegetation) and the plant canopy (as the percentage of vegetated surface viewed from a distance of four meters and one meter height in front of a Robel pole). We further recorded plant richness as the number of different plant species detected within the quadrat, and vegetation height as the averaged plant height measured at 10 random points along a diagonal line within the quadrat. Values from the three quadrats were averaged for each observation point.
Landscape-scale habitat composition was quantified according to a 30-m resolution raster data derived from the Finer Resolution Observation and Monitoring of Global Land Cover (http://data. ess.tsinghua.edu.cn), in which land uses and land cover were grouped into 10 categories: cropland, forest, grassland, shrubland, wetland, water, tundra, impervious surface, bare land, and snow/ice (Yu et al. 2013). In this scheme, only land covered by crops is included as cropland. Harvested agricultural land and grazed grassland with traces of cultivation are listed under barren land in consideration of their land cover function. Forests are areas that have a distinct canopy texture on Thematic Mapper (TM) images, including natural coniferous or broadleaf trees and single or mixed fruit trees. Impervious surfaces are primarily discriminated based on artificial cover such as asphalts, concrete, sand and stone, bricks, and other cover materials (Yu et al. 2013). For each sampling point, we recorded the proportion of different land cover types (grassland, cropland, bare land, forest, and impervious surface) within a 300-m buffer radius. We used a buffer radius of 300 m to minimize the spatial overlap between two adjacent buffered areas and consequently to avoid high redundancy of landscape attributes for two close point counts.
This distance also offered a compromise among the territory size of most song birds. We assumed that birds recorded within the radius were most closely associated with the specific survey habitats. We obtained edge density, the number of land cover patches and Simpson diversity of land cover types at the same landscape scale. These indices above were calculated using the "landscapemetric" R package (Hesselbarth et al. 2019).

Multifaceted diversity
We measured richness as the number of observed bird species per sampling point. We also calculated alpha diversity with the Shannon index (Hill 1973). We determined the functional diversity by combining the relative bird abundance data with six biological traits describing habitat specialization, nest site, migratory status, diet level, body mass, and clutch size (see Table  A1.1). These traits are generally related to species' habitat selection, resource requirements, or reproduction. For example, body mass is related with home range/territory size. Habitat specialization was counted as the number of habitat types the species is known to occupy (nesting, forage, roost, etc.), which was obtained from https://www.iucnredlist.org/; lower value indicates higher habitat specificity. Averaged body mass (g) was obtained from https://birdsoftheworld.org/, https://avibase.bsceoc.org, andhttps://animaldiversity.org/. Other data (e.g., Wang et al. 2018) were also considered when compiling bird trait information (see Table A1.1).
Based on these functional traits, we computed abundanceweighted functional diversity, using the metric "FD" developed by Petchey and Gaston (2002). First, we calculated Gower's distances between species based on a functional traits matrix, which is a measure of the similarity of species that can cope with mixed trait data (Podani and Schmera 2006). Second, we converted the distance matrix into a dendrogram using the unweighted pair group averaging (UPGMA) clustering algorithm, as it yielded a dendrogram with the highest cophenetic correlation with our original distance matrices and has also been identified to perform best in most cases (Podani and Schmera 2006). The "FD" package ) was used to compute the functional diversity index.
To obtain the phylogenetic tree, we downloaded the global phylogenetic tree of birds from "BirdTree" under the option of "Hackett All Species: a set of 10 000 trees with 9993 OTUs each" (http://birdtree.org; Jetz et al. 2014). A phylogeny for the bird species recorded in this study can be seen in Fig. A1.1. We pruned them with the data subset of bird abundance and traits. The phylogenetic diversity of each site was calculated using the "PD" function in the "picante" R package (Kembel et al. 2010).
Given that FD and PD were both highly correlated with species richness (FD: r = 0.96, p < 0.001; PD: r = 0.84, p < 0.001), null models were thus performed to assess whether the observed PD and FD were significantly different than expected given the observed species richness per se (Mouchet et al. 2010). A null distribution of FD and PD was generated for each site by randomizing community data matrix with the independent swap algorithm (Gotelli 2000) maintaining species occurrence frequency and sample species richness. To test whether the actual FD and PD for each community was significantly higher or lower than the mean of the null distribution, we repeated the randomization process 1000 times and calculated the standardized effective size (SES) of FD and PD values as follows: SES = (X obs -μ null )/ σ null , where X obs is the observed FD or PD value, μ null and σ null are the mean and standard deviation of randomized FD or PD value, respectively. This approach allowed us to determine if sesFD or sesPD are independent from species richness. Negative and positive SES values represent communities that are less or more functionally/phylogenetically diverse than expected given their species richness (Rader et al. 2014).

Statistical analyses
Prior to analyses, we log(x+1)-transformed land cover composition (in percentages) and zero-centered and scaled continuous predictors to improve model fitting and facilitate variable selection and result interpretation. We identified collinearity for all explanatory variables using variance inflation factors (VIF). The VIF is based on the square of the multiple correlation coefficients resulting from regressing a predictor variable against all other predictor variables, the higher the value, the higher the collinearity (Dormann et al. 2012). Much divergence exists in the literature regarding the VIF value to be used as the threshold for collinearity (Cenfetelli and Bassellier 2009). Commonly recommended values are 10 (Hair et al. 1998), 5 (Kline 2015), 3.3 (Petter et al. 2007), and 2 (Kock and Lynn 2012), meaning that a VIF equal to or greater than the threshold value would suggest the existence of collinearity among the variables. Such divergence makes it difficult to derive clear-cut methodological guidelines for researchers, and is in part due to the different contexts in which these values were proposed. We conservatively set the threshold value as two, and used a stepwise approach to exclude all explanatory variables with highest VIF each time. Four variables, i.e., bare ground percentage at the local scale and edge density, cropland percentage, and Simpson diversity of land cover types at the landscape scale were excluded during this procedure. We examined Pearson correlation matrices and reassured that the remaining eight variables were not highly correlated (|Pearson's r| < 0.37).
We used canonical correspondence analysis (CCA) to determine how environmental factors explain variance in species composition. CCA is a constrained ordination combining principles of ordination and regression. CCA allows models to be built using all bird species simultaneously and to test for statistical significance of linear combinations of explanatory variables using a unimodal response model (ter Braak and Smilauer 1998). CCA also allows a visual interpretation (a biplot) of species-environment relationships. Species that occurred on less than four sites were omitted because their distributions are poorly described by ordination techniques (ter Braak 1995). Forward selection and Monte-Carlo permutation test with 1000 permutations were used to find out which canonical axes contributed significantly to explain data variance. To perform CCA, we used "vegan" R package (Oksanen et al. 2013) in R3.6.2 (R Development Core Team 2019).
We checked the Pearson correlation coefficient among diversity measures to see whether the standard effective size of functional diversity (sesFD) and phylogenetic diversity (sesPD) are truly independent of species richness. We then fitted Bayesian generalized additive models implemented in the R package "brms" (Bürkner 2017) to analyze the habitat-diversity relationships. We developed a global model that included main effects for all local and landscape habitat attributes that were retained after VIF analysis, and included a smooth function of land cover richness, a habitat fragmentation measure, to test intermediate disturbance hypothesis. To run the GAM, we used the following formula: Diversity metrics = plant richness + plant canopy + grass height + s(land cover richness) + bare land percentage + grassland percentage + forest percentage + impervious percentage. Specifically, species richness and Shannon's diversity index were modeled using a gamma distribution, while sesFD and sesPD were modeled using a Gaussian distribution. Each model was run on four independent chains with 6000 iterations (warm up 2000 iterations). Chain mixing and convergence were checked by R-hat convergence diagnostic, a metric that compares the between-and within-chain estimates for Bayesian model parameters and other univariate quantities of interest. If chains have not mixed well, i.e, the between-and within-chain estimates do not agree, R-hat is larger than 1.1 (Gelman and Shalizi 2013). We reported the posterior distribution as posterior medians and posterior highest density intervals (HDI). The HDI summarizes the distribution by specifying an interval that spans most of the distribution, say 90% of it, such that every point inside the interval has higher credibility than any point outside the interval. We considered parameters and responses ecologically meaningful if the 90% HDI did not overlap 0. Additionally, we tested spatial autocorrelation of species richness and Shannon diversity using global Moran's I (Table A1.2) in the "lctools" R package (Kalogirou 2016). We further assessed spatial autocorrelation in the posterior means of the residuals of each Bayesian model by typical Q-Q plots ( Fig.  A1.2). We conclude that our model outputs were robust to spatial autocorrelation effects.

Bird surveys
Standardized bird surveys resulted in the detection of 3628 individuals of 92 species across 598-point counts. The species richness of a site ranged from 1 to 9 with an average of 2.5 and a standard deviation of 1.6. Shannon's diversity values ranged from 0.46 to 1.89 with an average of 1.18 (± 0.28). Tree Sparrow Passer montanus displayed the highest occurrence, being detected in 236 of the 598 sampling points (39.5 %), followed by Eurasian Skylark Alauda arvensis (37.5%) and Eurasian Magpie Pica pica (12.3%). All other species had a relative occurrence of less than 10%. The abundance of individual species is reported in Table A1.1. At less than five sites, 44 species were detected, either because they are rare, e.g., not grassland birds, or hardly detectable; hence, we decided to exclude them from subsequent data analyses and estimate the diversity indices from the occurrences and abundances of the remaining 48 species.

Canonical correspondence analysis
The total variance explained by the CCA performed on the eight explanatory variables was 7% for breeding bird assemblages. The first and second axis of the CCA were statistically significant (P < 0.01, Monte Carlo test with 1000 permutations), accounting for 25% and 16% of the constrained variation in the data, respectively. The first ordination axis indicated a gradient from sites with larger bare land percentage to sites with larger grassland percentage, associated with higher plant height, canopy and plant richness. Bird assemblages of grassland-dominated sites were composed of species requiring large areas of open grassy and shrubby habitats, such as Mongolian Lark Melanocorypha mongolica and Jankowski's Bunting Emberiza jankowskii. Species associated with more bare land and land cover richness were Daurian Jackdaw Corvus dauuricus and Red-rumped Swallow Cecropis daurica (Fig. 2). The second axis was mainly a woodland composition gradient from landscapes dominated by forests to landscapes dominated by impervious surface. Carrion Crow Corvus corone and Large-billed Crow Corvus macrorhynchos were associated with large amounts of forest cover, while Common Cuckoo Cuculus canorus, Meadow Bunting Emberiza cioides, and Beijing Hill-warbler Rhopophilus pekinensis typically occurred in impervious surface-dominated mosaics (Fig. 2).

Fig. 2.
Canonical correspondence analysis (CCA) ordination biplot of bird assemblages and habitat attributes. The proximity of species to arrows and their perpendicular distance along an arrow are measures of the relative influence of explanatory variables. Species that distributed closely to centroid were not labeled for clarity.

Diversity metrics
With the null model approach, the functional diversity and phylogenetic diversity are independent from species richness, with a low Pearson correlation coefficient (|r| < 0.15) among them, while species richness and Shannon's diversity are highly correlated (Fig. A1.4). We checked convergence with the potential scale reduction factor, R-hat, which is close to 1, indicating that Table 1. Relationship between avian diversity index and habitat variables, estimated from Bayesian GAM models. "Estimate" is the median of the posterior distribution of each variable. "Std.Error" is the standard error of each parameter estimate, i.e., median absolute deviation. HDI computes the highest density interval for posterior samples, which contains the parameter values of highest probability and that span the 90% most probable values. Any parameter value inside the HDI has higher probability density, i.e., higher credibility, than any parameter value outside the HDI. HDI not containing zero are in bold. Rhat is the potential scale reduction factor or Gelman-Rubin diagnostic and is a measure of how well the chains have converged and ideally should be equal to 1. Fixed and smooth(land cover richness) is the fixed and smooth effect part of land cover richness, respectively. The posterior estimate of smooth(land cover richness) indicates how wiggly the fitted model is, values around 1 indicates the smooth term tend to be close to a linear term. Bayesian GAM models for each diversity metric were well supported by the data (Table 1). Relationships between habitat variables and avian diversity were scale and metric specific (Table  1), and we describe the "significant" responses explicitly. The strongly correlated species richness and Shannon diversity exhibited a similar response pattern; both are negatively related to bare land percentage while positively related with plant canopy and impervious surface percentage, and the parameter estimates of which are nearly identical (Table 1). Phylogenetic diversity was positively associated with plant richness while negatively associated with forest percentage and impervious surface percentage (Table 1). Functional diversity does not relate significantly to any habitat variables used in this study (Table 1).
Additionally, for the four measures of total diversity, we did not observe any significant response to the smooth term of land cover richness, because the 90% highest density interval does include 0. The posterior estimates of the smooth term are less than 2, indicating it tends to be a linear parametric effect (Table 1). We did not find any evidence that diversity metrics peaked at intermediate level of habitat disturbances, i.e., a clear hump-shape curve along the gradient of land cover richness; most of them were flat or monotonic increasing functions (Fig. 3).

DISCUSSION
In the present study, we examined the overall response of bird assemblages to local and landscape-scale habitat features, and identified the key habitat attributes that drive taxonomic, functional, and phylogenetic diversity of bird communities in Inner Mongolia, using a dataset comprising a large number of sampling sites. Several past studies on biodiversity issues either concentrate on single diversity index, e.g., species richness, or solely evaluate the local-scale habitat effects in Inner Mongolia (e.g., Wen et al. 2002, Liang et al. 2018, but see Liang et al. 2019).
Our study advances these works by pointing to clear differences in bird composition along habitat gradients, and revealing that relationships between habitat variables and avian diversity are scale-dependent and metric-specific. Our results suggested a range of diversity indices are likely to be required to better understand how bird assembly is shaped in human-modified Inner Mongolian grasslands. Fig. 3. The relationship between species diversity and land cover richness fitted by Bayesian generalized additive models. Blue line indicates the mean predicated value, and the grey area indicates 90% credible intervals based on the standard errors.
We found that measured habitat attributes are not the only environmental factors structuring bird assemblages given the low amount of variation explained by canonical correspondence analysis. Some biotic factors, e.g., species competition or preypredator relationship, missed in our study or other abiotic variables not captured by these habitat variables might explain this remaining variance. For example, a group of water bird species, e.g., Phalacrocorax carbo and Sterna hirundo, were identified in the upper right of CCA ordination biplot, indicating a potential gradient from upland (left) to lowland (right) or dry to wet along the first ordination axis, which could be affected by precipitation or natural/artificial wetlands that were not recorded in this study. Several grassland-specialist birds, such as Emberiza jankowskii, were associated with features, e.g., higher grass height and grassland percentage, that are typical of native grasslands, while others, such as Cuculus micropterus and Rhopophilus pekinensis were associated with features more characteristic of human settlements. These results are consistent with previous studies (e.g., Han et al. 2020) that have highlighted the importance of natural grasslands for many typical grassland birds. Our results also show that several species, such as Emberiza cioides, could make use of highly modified grasslands. Thus, from a conservation perspective, bird populations in Inner Mongolia would benefit not only from policy strategies focused on the conservation of natural grasslands, but also from management guidelines for agriculture sites.
We found that Shannon diversity and species richness were positively related to increased plant canopy cover, which may reflect higher niche segregation opportunities Rotenberry 1981, Ferger et al. 2014). The negative effect of bare land percentage could be predicted by the classical theory of the species-area relationship (Schoener 1976), that is, lower habitat suitability with increased bare land reduces species abundance, which is supported by abundant empirical evidence (Fahrig 2003). We found that areas with a higher impervious surface percentage were associated with higher levels of bird species richness and Shannon diversity. One possible reason is that some species with generalist niche requirements are less sensitive to human-induced habitat modification and may therefore experience low rates of local extinction across the urban habitat matrix ). For example, species with omnivorous or generalist nesting habits, such as swallows or pigeons, are expected to be positively affected by, or even thrive, in urban environments (McKinney and Lockwood 1999). Another reason could be that areas with higher impervious surface are generally next to human settlements, and could support diverse resources such as garbage or water for bird species in semiarid environments. The presence of impervious surfaces could also correspond to a mosaic of natural habitats and human settlements increasing the overall number of species by allowing habitat specialists and generalists to coexist .
In contrast with taxonomic diversity, phylogenetic diversity significantly decreased with higher impervious surface percentage. Frishkoff et al. (2014) pointed out that two processes were responsible for changes in phylogenetic diversity: species loss and increasing species relatedness. In our case, though the species richness increased with higher impervious surface percentage in urban or suburban areas, they are probably achieved by the incorporation of phylogenetically close species. Our results confirm the expectation that by homogenizing the physical environment to meet the demands of human beings, urban or suburban settlements may be creating filters that determine that only some lineages can persist in these new environments (McKinney 2006). These lineages generally include species-rich, widely distributed families like pigeons, crows, and swallows (Sol et al. 2017). Our results are partly in line with previous studies showing that land cover change, particularly urbanization, causes phylogenetic homogenization (e.g., Liang et al. 2019, Weideman et al. 2020).
Phylogenetic diversity significantly increased with higher plant richness and land cover richness. Niche theory can in part explain the response of PD to plant and land cover richness (Wiens and Graham 2005). Increasing plant richness and land cover richness potentially reflect a greater variety of habitats, i.e., more niches, which in turn should support more species to occur in the same general area. However, the significantly positive link to habitat diversity was only found in phylogenetic diversity instead of species richness or Shannon diversity. Although the mechanisms remain unclear, one possible explanation is that species thriving in mixed habitats tended to be overrepresented in a few distantly related taxa and often lacked close-relatives in the same community; their presence may contribute to explain why diverse land cover had higher phylogenetic diversity than expected by the number of species they sustained (Sol et al. 2017).
In comparison to taxonomic and phylogenetic diversity, functional diversity was not significantly related to any of the local and landscape attributes. Several reasons may explain the weak and insignificant response of the FD indices. First, the data we originally collected for studying taxonomic diversity could be insufficient for capturing FD, as for most empirical datasets, a much higher sampling effort is generally required to reach an adequate precision for FD measurements than for taxonomical diversity measurements (van der Plas et al. 2017). Second, the low variation in Shannon diversity between sites may contribute to a high probability of sensitivity in functional diversity, i.e., making FD sensitive to methodological choices such as distance measure or clustering algorithm (Poos et al. 2009). Although we need to be careful about generalizing from our case study, our results are in line with previous studies emphasizing that human-induced impacts can differentially affect these three biodiversity components (e.g., Devictor et al. 2010, Morelli et al. 2017) and are consistent with recent recommendations in nature conservation optimization in favor of using multiple complementary biodiversity metrics (Winter et al. 2013).
Based on the IDH and diversity-trophic structure hypothesis (DTSH), we were expecting that higher land cover richness will increase resource diversity for associate bird species (Knops et al. 1999). Different diversity metrics, however, either produce a nonsignificant intermediate peak along the gradient of land cover richness, or the trends remained linear, as in Figure 3. Thus, our hypothesis that bird diversity peaks at intermediate levels of habitat disturbance (IDH) was not supported, though there are many other studies showing species richness or abundance follows the prediction of the IDH (reviewed by Mackey and Currie 2001). The use of various indices of diversity as well as the different definitions of ecological disturbance can be large sources of variation in outcomes among studies (Svensson et al. 2012). Our results suggest that consideration of land cover richness alone may give misleading conclusions of its effects on the bird assembly. Several other factors may mediate the effect of land cover richness, such as the interactive effects with grassland percentage or with other unmeasured but important smaller-scale variables, e.g., habitat quality. Furthermore, we observed a marginally positive response of phylogenetic diversity along the land cover richness gradient. The logical or mathematical basis underlying this pattern is still not clear. It appears to be driven by a complex interplay of different environmental predictors. Specific traits of individual species may also play a role in the response of PD to disturbance. However, it is generally difficult to determine what biological attributes the species might be expected to have to provide the essential ingredients for coexistence, and then demonstrating, in the actual community, that such attributes are both present and active in promoting species coexistence (Shea et al. 2004). Additional experiments, however, are required to mechanistically understand the relative importance of these biological attributes.
There are several issues that, at the very least, should be kept in mind. For example, our survey design did not consider well the behavioral differences between species groups, e.g., waterbirds vs passerines, and the point-count methods could be inappropriate for several nocturnal species. Potentially important factors include those that are linked to land use intensity, e.g., road and traffic and human population density, or certain types of microhabitat are not available in our satellite images. It should also be noted that the magnitude of the response to habitat variables was not so strong, thus some other effects may have been missed because of a limited range of variation in our data. For example, the gradients of land cover richness had a small coverage, which may have somewhat impacted our results, so a larger sample size could reduce this variance. Despite these limitations, the data presently available clearly show how local and landscape-scale habitat attributes are influencing multiple biodiversity of birds in Inner Mongolian grasslands.

CONCLUSION
Grassland birds are threatened by habitat loss and fragmentation because of human-dominated lands. Our approach to quantifying local and landscape attributes enabled us to identify the positive effects of plant canopy and plant richness for taxonomic diversity, as well as the negative effects of fragmentation created by the inclusion of impervious surface for phylogenetic diversity of birds in Inner Mongolia. These findings illustrate the diverse response to habitat attributes by grassland birds and the importance of maintaining a local and landscape perspective to conserve avian communities. Our results show that an increase in species richness only can be misinterpreted as a sign of conservation improvement, since a significant loss of phylogenetic diversity may be masked by apparent increases in taxonomic indicators. We thus propose that using species-specific ecological traits and considering potentially different facets of diversity could provide more accurate insights for conservation planning purposes.
Responses to this article can be read online at: https://www.ace-eco.org/issues/responses.php/1745 Figure A1.1 Visualization of the phylogenetic tree for the recorded 92 bird species in Inner Mongolia. Figure A1.2 QQ plots to check the spatial autocorrelation in the posterior means of the residuals of each Bayesian generalized additive model. Figure A1.3 a-e: species richness changed along the different land cover gradients, the blue line indicated a loess-smooth fitted function, and the light blue represents standard errors.
f: boxplots are used to visualize the median, two hinges and all outlying points of the five different land cover variables. Figure A1.4 Pearson correlation coefficient among the four diversity metrics. The bivariate scatter plots with loess smooth lines were drawn below the diagonal, histograms were shown on the diagonal, and the Pearson correlations were reported above the diagonal. This figure was made by the function "pairs.panels" in "psych" R package. Reference: Wang, Y., Si, X., Bennett, P.M., Chen, C., Zeng, D., Zhao, Y., Wu, Y. and Ding, P.
doi:10.1111/ecog.03158 Table A1.2 Results of Moran's I test to examine spatial autocorrelation of species richness and Shannon diversity between sampling sites. The Moran's I statistic ranges from -1 to 1. Values near 0 indicate no spatial autocorrelation (no spatial pattern -random spatial distribution) and values in the interval (0, 1) indicate positive spatial autocorrelation (spatial clusters of similarly low or high values between neighbors and vice versa). The z scores and p values that calculated for resampling and randomization null hypotheses are also listed.