Land-use change undermines the stability of avian functional diversity

Survey data

To assess impacts of land-use change on bird diversity, we began by collating surveys from the PREDICTS database, a repository of species occurrence and abundance data sampled across multiple land-use types⁷⁰. We removed 30 datasets because they lacked abundance data (n = 12) or were incomplete (n = 18), usually because sampling was limited to particular guilds or methods, such as camera traps (Supplementary Information). This process produced a baseline of 72 datasets that we augmented by conducting a systematic literature review using Web of Science to identify further published bird surveys that targeted land-use gradients (Supplementary Information). After contacting authors for data, we received 29 suitable datasets that we added to the PREDICTS database. Our sample in this study contains data from these 29 surveys, along with 5 additional datasets released in the latest version of PREDICTS⁷¹. Finally, to improve sampling in megadiverse regions, we integrated further independent datasets generated by intensive surveys in Bornean⁷² and Amazonian rainforests⁷³. Geographical location and sources for all published surveys used in our analyses are presented in Fig. 2, Supplementary Table 1 and Supplementary Data 1.

Most survey data in the PREDICTS database are organized as a hierarchy. Survey sites are nested in study blocks, and blocks are nested in study landscapes. Study blocks are spatially segregated but not always temporally defined⁷⁰. We found 16 datasets in PREDICTS that contain surveys sampled in different years or seasons; therefore, we subdivided the data into separate study blocks partitioned by location and time of survey. To ensure consistency, we also collapsed 17 surveys in PREDICTS into 7 studies by combining all data extracted from the same original source publications. We then partitioned these seven studies into distinct study blocks representing geographical and temporal subsets. After restructuring, our final dataset consisted of 98 study landscapes (Fig. 2), each with numerous survey sites clustered into study blocks. Subdividing our data in this way enabled us to account for spatial, annual and seasonal effects across studies by including study block as a random effect in our models (Supplementary Information).

We converted survey data into species assemblages to enable comparisons across land-use types. In some study landscapes, species assemblages reflect the total number of species identified in a study site, usually pooled across a series of transects or point-counts conducted at intervals between dawn and midday. Other published studies focused at finer resolution, sometimes defining each point-count as a separate survey site. Single point-counts generally undersample species richness, and neighbouring sites may be very close together, which caused problems for our analyses. In these cases, to facilitate calculation of functional metrics and to minimize the risk of pseudo-replication, we grouped species into larger assemblages by aggregating survey sites with similar land-use types in the same study block (Supplementary Information). Species assemblages were therefore defined as all species encountered in a restricted, spatially and temporally segregated area, largely confined to the same land-use type. We do not use the term community because we do not have direct data confirming species interactions⁷⁴.

Land-use classification

To classify land-use types for each survey site, we used PREDICTS data to estimate the predominant type and stage of vegetation and the intensity of human use. First, we assigned sites to one of six vegetation classes: primary vegetation, secondary vegetation, plantation forests, pasture, cropland and urban. We then classified lightly or intensively used primary vegetation sites as disturbed primary vegetation, whereas minimal-use primary vegetation sites were classified as pristine primary vegetation (a proxy for undisturbed natural vegetation). Based on previous analyses showing reduced avian FD in intensely urbanized areas⁵⁵, we also split minimal-use urban sites from sites with more intensive urbanization. Finally, to account for the effects of vegetation structure at different successional stages⁷⁵, we partitioned secondary vegetation according to age class (mature, intermediate, young, indeterminate; Supplementary Information). Indeterminate age secondary vegetation was removed from our dataset.

Our final dataset consisted of 98 study landscapes distributed across 6 continents (Fig. 2a and Supplementary Table 1), representing a total of 1,281 avian assemblages in 10 distinct land-use types: pristine primary vegetation (n = 177); disturbed primary vegetation (n = 281); mature secondary vegetation (n = 44); intermediate age secondary vegetation (n = 77); young secondary vegetation (n = 86); plantation forest (n = 218), pasture (n = 184); cropland (n = 107); and urban, including both minimal-use (n = 46) and intense-use (n = 61) urban landscapes.

Functional trait data

Species traits can provide information about sensitivity to perturbations (response traits) and the impacts of species presence or absence on ecological function (effect traits). In both cases, data availability is often patchy for major taxonomic groups at a global scale⁷⁶. We obtained morphometric measurements for all 3,696 species reported in our study landscapes from the AVONET trait database⁴⁴. Species means were compiled for seven traits: beak length (culmen), beak length (tip-to-nares distance), beak depth, beak width, tail length, tarsus length and wing length (Supplementary Data 1). These traits have been shown to predict a range of key ecological niche axes, including diet and foraging strategy¹⁶ (Supplementary Table 3). We also included data on the hand–wing index (HWI), a metric of wing elongation that predicts aerial lifestyle and dispersal distance in birds⁷⁷. HWI is widely used as a proxy for dispersal ability⁷⁸. Species mean values for all traits used in our study were calculated from an average of 11 individuals per species (41,515 individual birds measured in total).

Avian morphological traits are often intercorrelated because of an underlying association with body size¹⁶. Accordingly, all traits in our dataset were strongly correlated with the body size axis (Extended Data Fig. 2a), apart from HWI (R = 0.22). Following previous studies^40,79, we removed the association with body size through a two-step PCA, which reduced our seven linear morphometric traits into three niche axes related to ecological functions (Extended Data Fig. 2b). We performed two separate PCAs on trophic traits (related to beak morphology) and locomotory traits using all species in our dataset. In both cases, the first principal component (PC) was strongly correlated with body size; therefore, we used the second PC to represent the dominant axis of variation, which is effectively independent from body size (Extended Data Fig. 2c). We then performed a third PCA on the first PC scores from both the trophic and locomotory PCAs, taking the resultant first PC to represent the body size axis.

We use this metric of body size because it correlates strongly with body mass while also reflecting trophic niche differences. As the body size axis is extracted from linear measurements of beak, wing, tail and tarsus, it more closely reflects trophic and locomotory niches than conventional body-mass estimates. For example, our estimates of body size will distinguish between hummingbirds with equal body mass but different beak lengths, thereby reflecting associations with different foraging niches linked to pollination. Finally, we supplemented these three derived trait axes (trophic, locomotory and size) with a fourth morphological trait axis consisting of the log-transformed HWI (related to dispersal ability).

Taxonomic matching

Bird species names and classifications vary over time and between different taxonomic treatments. This is problematical for global datasets based on published field surveys because different authors use a variety of taxonomic approaches, including English or local names. We converted all species names into a single taxonomy using published cross-walks⁴⁴ and verified taxonomic assignments with geographical range maps⁸⁰ (Supplementary Information). This enabled accurate alignment with species trait data. Some taxa reported in survey data were impossible to assign directly to species because they were only identified to the genus level. Deleting these taxa would result in missing data, which can reduce the accuracy of FD estimates⁸¹. Instead, we created pseudo-species representative of the genus. Given that avian life history and morphological traits tend to be highly conserved within genera⁴⁴, we assigned trait values to pseudo-species by averaging the trait values of all congeners potentially occurring at the locality. To generate trait data for averaging, we used geographical range maps to provide a list of all members of the focal genus with geographical distributions overlapping the site location (Supplementary Information). We synthesized data for 73 pseudo-species in 133 of our 1,281 study assemblages.

Dietary data

Birds mediate a wide range of ecological processes and services depending on their trophic interactions, including seed dispersal by frugivores and pest control by invertivores^82,83. The morphological trait dataset used in this study was strongly correlated with avian diets and associated foraging behaviours^16,43. However, the connection between morphology and dietary niche was noisy and weak in some taxonomic groups (Supplementary Information). Therefore, we also included standard diet classifications in functional metric calculations. We used published estimates of the proportion of species diets across nine major resource types: herbivore (aquatic), herbivore (terrestrial), nectarivore, granivore, frugivore, invertivore, vertivore (aquatic), vertivore (terrestrial) and scavenger. The data were extracted from a previous publication¹⁶ and were primarily based on the EltonTraits dataset⁸⁴ with extensive updates and reorganization based on subsequent literature.

To define dietary groups for analyses, we used published data that classified species into trophic guilds according to their primary food source, with any species obtaining >60% of their diet from a single food type defined as a trophic specialist⁴⁴. Species that obtained resources more equally across different food types were classed as omnivores^16,85. Our main analyses focused on trophic guilds rather than omnivores because bird species with more specialized diets have a higher certainty of contributing to particular ecological roles and services⁴⁰. However, generalists are often abundant, which suggests that they may contribute substantially at the population level to ecological processes such as seed dispersal and insect predation. Thus, non-specialist omnivores may help to stabilize ecosystems by providing additional redundancy. To assess whether our results were sensitive to trophic guild classification and inclusion of generalists, we generated broader dietary groupings containing all species that obtained >25% of their diet from a single food source. We then repeated our main analyses on these expanded groups (Extended Data Fig. 6 and Supplementary Information).

Calculating FD and redundancy

To calculate functional metrics, we created trait probability densities (TPDs) using the TPD package in R⁸⁶. The TPD approach uses species mean trait-values and intraspecific trait variation to calculate probabilistic hypervolumes in which the potential position and extent of occupancy for each species can be predicted along multiple trait axes. By using axes of trait variation to define ecological niche axes, TPD hypervolumes represent a Hutchinsonian niche⁴⁵. We constructed species hypervolumes on the basis of their diet proportion data across nine major resource types and their distribution along four derived morphological trait axes (locomotory, trophic, dispersal and size). We then estimated FD for each assemblage as the total cumulative volume occupied in trait space by all species in the assemblage. Functional redundancy was calculated as the proportion of this total volume shared by multiple species, weighted by the relative abundance of species occupying the same regions of trait space⁴⁷ (Supplementary Information).

To create TPDs, we first calculated distance matrices using the R package gawdis, which is specifically designed to combine compositional data (such as our proportional diet data) into a single axis of variation⁸⁷ (Supplementary Information). We calculated distance matrices for each of our 98 study landscapes using diet and morphological data for all species present. Following previously described methods⁸⁸, we back-transformed our distance matrices into three-dimensional coordinates representing the relative position of each species in functional trait space (Supplementary Information).

Calculation of TPDs requires the estimated position of species means in trait space and the square-root of intraspecific variability. The latter is required to generate a probability kernel around each species-mean position. This intraspecific variation kernel (IV kernel) is taken as the niche of each species, represented in functional trait space. Ideally, the IV kernel is estimated by directly comparing measurements of many conspecific individuals and calculating the standard deviation across each dimension in functional trait space. However, as bird diet data are only available as species-mean estimates, we were unable to use this method for our main analysis. Thus, following previous studies^47,89,90, we approximated the IV kernel using a bandwidth estimator, which calculates an equally sized density kernel around each species mean based on the distances between co-occurring species (Supplementary Information). As the volume covered by the IV kernel can vary between assemblages, we calculated the dimensions of the kernel for each assemblage and took the square-root of the mean dimensions as our common IV kernel for each species (Supplementary Information). For each of our 98 study landscapes, we calculated a landscape-level TPD using the TPDsMean function, which generates a TPD using the species-mean positions in each assemblage, alongside the IV kernel dimensions. The functional trait space of each TPD was divided into 125,000 equally sized grid cells, with a value reflecting the likelihood of occupancy.

For each assemblage in each study landscape, we calculated FD and functional redundancy by running the in-built REND and redundancy functions from the TPD package across respective landscape-level TPDs⁴⁷. Both functions filter the landscape-level TPDs for the species that occur in each assemblage. FD was then calculated as the number of grid cells with a likelihood of occupancy >0. Redundancy was calculated as the average number of species that could be removed from each grid cell without reducing the total occupied area of functional trait space⁴⁷ (Supplementary Information). By setting a minimal threshold for occupancy (>0), the FD value reflects the volume of occupied trait space independent of species abundance. Conversely, functional redundancy calculations are shaped by abundance and highly sensitive to the number of species that share a similar area of trait space (Supplementary Information).

To analyse the effect of land-use change on specific ecological roles, we calculated functional metrics across five species subsets to assess how functional trait structure changed within different trophic guilds in each assemblage. We focused on the entire assemblage (n = 1,281) and all generalist species (n = 1,281), as well as three key dietary guilds sampled across all land-use types: granivores (n = 1,271), frugivores (n = 944) and invertivores (n = 1,274). Sample sizes varied between guilds as TPDs cannot be created when three or fewer members of the guild are recorded in a given study landscape.

We did not standardize functional metrics by the number of species present and instead allowed both FD and redundancy to correlate with species richness. We assumed that each additional species added to an assemblage will either increase the number of trophic processes performed or increase the probability that multiple species deliver a particular function, or both⁵. We therefore allowed increased species richness to drive greater FD and redundancy values in larger assemblages. To examine how this decision influenced our results, we also conducted a supplementary analysis in which we modelled the relationship between functional resistance and a species richness standardized redundancy metric (Extended Data Fig. 8 and Supplementary Information).

To address whether our method of estimating the dimensions of our IV kernel affected our conclusions, we recreated our TPDs using intraspecific variation calculated from multiple measurements of conspecific individuals, using morphological measurements extracted from AVONET⁴⁴ (Supplementary Information). The results of the sensitivity analyses based on these revised TPDs were similar, which showed the same general patterns of decline in FD or redundancy with land-use change (Extended Data Fig. 4).

As overall assemblage redundancy reflects the total amount of shared trait space in the assemblage, declines in redundancy are driven by either species losses or reduced amount of niche overlap per species (that is, niche differentiation). To decipher whether the effect of redundancy losses on functional resistance was driven by declines in species richness or reduced niche overlap per species, we also calculated the relative redundancy for each assemblage using in-built functions in the TPD package (Extended Data Fig. 8).

Calculating sensitivity scores

To calculate functional vulnerability and functional resistance metrics, we began by scoring sensitivity to disturbance, which reflects the likelihood that a particular species would be removed from the assemblage by future perturbations. For each species in each assemblage, we generated two forms of sensitivity score: (1) trait-based and (2) rarity-based. Our main trait-based sensitivity score was based on four general response traits associated with extinction risk: geographical range size, body size, diet specialism and dispersal limitation^91,92 (Supplementary Information). This is a broad bandwidth score that uses traits associated with any form of disturbance (for example, fire, storms, drought, habitat loss, pollution or human exploitation) to reflect uncertainty regarding the source and severity of future perturbations. To provide a more explicit test of sensitivity associated with a known stressor, we calculated a secondary trait-based score using response traits associated with sensitivity to climate change: higher elevational distributions, narrower temperature niche, longer generations and dispersal limitation⁹³. To calculate our rarity-based sensitivity scores, we extracted the inverse (that is, negative) abundance of each species in the assemblage based on the assumption that rarer species are more likely to be removed from an environment by population fluctuations⁹⁴. All sensitivity scores were scaled by their standard deviation and centred to have a mean of zero.

Functional vulnerability

Previous studies have associated functional stability with the distribution of unique functional traits⁹⁵ and have linked assemblage vulnerability to the variation in disturbance sensitivity of species traits⁶⁰. We devised a new metric of functional vulnerability to encapsulate both these concepts. To estimate the relationship between trait redundancy and species sensitivity, we quantified the amount of redundancy provided by each species. Species redundancy was estimated as the change in assemblage redundancy after removal of the focal species. To calculate this value, we separately removed each species from the full assemblage and recalculated the assemblage redundancy. We then calculated two functional vulnerability scores as the covariance between species redundancy and either trait-based or rarity-based sensitivity scores using Spearman’s rank correlation coefficient. We reversed the direction of the covariance so that a high positive correlation indicates that the most sensitive species in the assemblage tend to be the least redundant (that is, most unique).

Functional resistance

To estimate functional resistance, we ran species extinction simulations using sequential removals and quantified the associated decline in FD. The sensitivity of species to environmental change is strongly influenced by their functional response traits, with certain trait combinations predicting the likelihood of local extinction^36,41. Moreover, species abundances are often a strong indicator of local extinction risk, as rarer species are typically more sensitive to environmental perturbations^96,97. Therefore, for each assemblage, we ran two extinction simulation scenarios by sequentially removing each species according to both sensitivity scores (trait-based and rarity-based).

Following previous methods^47,98,99,100, we plotted extinction curves to quantify how FD declines as the proportion of species occurring in the original assemblage is reduced to zero. Using a standard approach, we then measured the AUC as an estimate of functional resistance^99,100. If an assemblage maintains constant FD when species are removed, the AUC remains large, which indicates high levels of functional resistance. Conversely, when extinctions drive declines in FD, the AUC decreases (Extended Data Fig. 1 and Supplementary Fig. 1). For each extinction curve, the AUC was measured using the MESS package in R¹⁰¹.

As we were specifically interested in quantifying the pace of FD decline, we standardized the extinction curves by scaling the FD values between 1 and 0, where 1 equals the FD calculated for the full assemblage before any species was removed. This standardization prevented the magnitude of the FD value from influencing the AUC, which ensured that AUC values reflect variation in the shape of the extinction curve only⁹⁸. For each assemblage, we calculated our functional resistance metric across the same five subsets of species as in the preceding FD and redundancy analyses: all species, all dietary generalists, granivores, frugivores and invertivores.

One drawback of the AUC approach is its sensitivity to the order in which species are lost from an assemblage. For example, when a single morphologically unique species is lost before other more redundant species, this causes a steep initial decline in FD. Therefore, following the same methods used to generate functional vulnerability values, we also estimated functional resistance as the half-life (t_1/2) of each extinction curve (Supplementary Information).

In our main simulations, we implemented extinction scenarios based on a set of traits associated with general disturbance. To assess the robustness of our results to methods and choice of traits, we conducted a series of sensitivity analyses. First, we calculated extinction curves generated using trait-based sensitivity scores based on a different combination of response traits specifically related to climate change tolerance. Second, we re-ran analyses under passive (probability-weighted) species loss scenarios, in which 0–2 species were removed at each time step, and the probability of a species being removed at each step was equal to its sensitivity score (Supplementary Information). Under this procedure, highly tolerant species were allowed to remain in the assemblage longer, and assemblages with a high proportion of tolerant species did not necessarily lose species at each time step, thereby increasing the AUC. Results of these sensitivity analyses were similar to our main analyses (Extended Data Fig. 7).

Statistical analyses

To assess how the distribution of key functional traits were affected by land-use change, we subdivided assemblages (n = 1,281) into four categories: (1) primary vegetation, (2) secondary vegetation, (3) agriculture (including plantation forests) and (4) urban. We then constructed a set of univariate linear mixed-effects models with land-use as a single predictor variable and nine separate response variables covering all response and effect traits analysed in this study (see Supplementary Information for the rationale). Models were assessed by comparing coefficient estimates for each land-use type against primary vegetation. When coefficient estimates were positive, land-use change was inferred to filter species with low functional trait values.

To assess the impact of land-use on functional structure, we first constructed three univariate linear mixed-effects models. In all three models, land-use was the sole predictor variable and each model analysed the effect of land-use change on FD, functional redundancy or functional vulnerability. Each model was conducted across all species and separately across four different data subsets related to dietary guild (generalists, granivores, frugivores and invertivores), and two subsets related to climatic region: tropical study landscapes and non-tropical study landscapes. For these models, land-use was split into ten distinct categories: pristine-primary vegetation, disturbed primary vegetation, mature secondary vegetation, intermediate age secondary vegetation, young secondary vegetation, plantation forests, pasture, cropland, minimal urban and intense urban.

To address how land-use change affects the functional resistance of bird assemblages, we ran six additional sets of univariate mixed-effects models. These analyses modelled the effects of land-use change on the functional resistance of each assemblage under both of our main species-loss scenarios (general trait-based AUC and rarity-based AUC) and our four alternative functional stability (general trait-based t_1/2, rarity-based t_1/2, climate trait-based AUC and passive AUC). For these final analyses, we also split land-use into four categories: (1) primary vegetation, (2) secondary vegetation, (3) agriculture (including plantation forests) and (4) urban. In line with FD and redundancy analyses, we conducted our functional resistance models separately across all species and our four dietary guild subsets (generalists, granivores, frugivores and invertivores).

All models were conducted using the lme4 package in R¹⁰². We added study landscape to account for among-study differences in sampling methods and study block to account for confounding variables related to temporal and geographical distinctions in each study. Models were interpreted by comparing the change in estimated regression coefficients for each land-use type using the functional metric value calculated for the least disturbed land-use category as our reference. Our hierarchical modelling approach did not require all land uses to be present in each study landscape to produce standardized regression coefficients. For most models, our least disturbed category was termed pristine primary vegetation, which accounted for 177 (13.8%) of our survey sites and was present in 50 out of our 98 study landscapes (51%). Although no primary vegetation landscape is entirely unaffected by human disturbance, this term is used to differentiate between more heavily disturbed vegetation types. For our functional stability models, we combined all primary vegetation sites into a single land-use type and used this broader classification as our least-disturbed land-use category.

To account for seasonal effects, we ensured that assemblages in the same study block were surveyed in the same season. Therefore, by including study block as a random effect, we removed detection biases arising from variation in sampling season. Moreover, incomplete sampling owing to undetected rare or cryptic species was reduced by aggregating survey sites into larger species assemblages to increase overall sampling depth (Supplementary Information).

We assessed all models for normality of residuals and heteroskedasticity and did not find that any of our models violated the assumptions of a linear model. Owing to the hierarchical structure of our data, it was difficult to incorporate covariance structures that account for spatial autocorrelation between local survey sites in the same study landscape into our global models. However, following a previously described method¹¹, we assessed spatial autocorrelation in study landscapes and study blocks separately using Moran’s I and found that it did not affect our results (Extended Data Fig. 9).

Reporting summary

Further information on research design is available in the Nature Portfolio Reporting Summary linked to this article.