A global analysis of avian island diversity– area relationships in the Anthropocene

Research on island species-area relationships (ISAR) has expanded to incorporate functional (IFDAR) and phylogenetic (IPDAR) diversity. However, relative to the ISAR, we know little about IFDARs and IPDARs, and lack synthetic global analyses of variation in form of these three categories of island diversity-area relationship (IDAR). Here, we undertake the first comparative evaluation of IDARs at the global scale using 51 avian archipelagic data sets representing true and habitat islands. Using null models, we explore how richness-corrected functional and phylogenetic diversity scale with island area. We also provide the largest global assessment of the impacts of species introductions and extinctions on the IDAR. Results show that increasing richness with area is the primary driver of the (non-richness corrected) IPDAR and IFDAR for many data sets. However, for several archipelagos, richness-corrected functional and phylogenetic diversity changes linearly with island area, suggesting that the dominant community assembly processes shift along the island area gradient. We also find that archipelagos with the steepest ISARs exhibit the biggest differences in slope between IDARs, indicating increased functional and phylogenetic redundancy on larger islands in these archipelagos. In several cases introduced species seem to have 're-calibrated' the IDARs such that they resemble the historic period prior to recent extinctions.

Functional diversity (FD) measures the combination of functional traits expressed by a set of species in a community and can provide a link between species composition and ecosystem function (Petchey & Gaston, 2006). Phylogenetic diversity (PD) incorporates the evolutionary relationships among species in an assemblage and measures the amount of evolutionary history those species represent (Faith, 1992). Collectively, the ISAR, IFDAR and IPDAR have been termed island diversity-area relationships (herein IDARs) and together their analysis aids in generating a more comprehensive understanding of the mechanisms driving the scaling of diversity (Ding et al., 2013;Leclerc et al., 2022;Mazel & Thuiller, 2021;Wang et al., 2023). However, in comparison with the ISAR, we know relatively little about IFDARs and IPDARs, and we lack synthetic comparative global analyses of variation in the form of these three categories of IDAR.
A wide range of metrics have been proposed for measuring FD and PD. To construct IFDARs and IPDARs that compare easily with standard ISARs, FD and PD are often expressed as metrics that sum the branch lengths (e.g. of a functional dendrogram or phylogenetic tree) connecting all species co-occurring on an island (Dias et al., 2020;Mammola et al., 2021;Morlon et al., 2011;Si et al., 2022). While the use of tree/dendrogram-based FD and PD metrics ensures the ISAR, IFDAR and IPDAR are comparable, such metrics are generally correlated with species richness. For this reason, the calculation of FD and PD using tree metrics is often combined with a null model to generate (standardized) effect sizes (ES) that are independent of richness (Mazel & Thuiller, 2021;Tucker et al., 2017). In addition, the analysis of ES values has been argued to provide insights into the community assembly processes involved (e.g. neutral dynamics vs. competition vs. habitat filtering) and how these may change with island area Mazel & Thuiller, 2021;Münkemüller et al., 2020;Schrader et al., 2021). Herein, we refer to the emergent FD.ES and PD.ES patterns (random, overdispersed, clustered) as assembly patterns, and the potential mechanisms underlying these patterns (neutral dynamics, competition, habitat filtering) as assembly processes. However, we know very little about how FD.ES and PD.ES values scale with island area (rather than across continuous scales; see Kraft & Ackerly, 2010), and previous authors have called for a greater focus on scaling patterns in order to better understand community assembly processes on islands (Dias et al., 2020;Leibold & Chase, 2018;Münkemüller et al., 2020;Schrader et al., 2021;Si et al., 2022;Zhao et al., 2020).
Analysing variation in IDARs among archipelagos can emphasize the (i) form/shape of the relationship (e.g. Mazel et al., 2014) and (ii) slope of the curve. The former is important as different relationship forms (e.g. asymptotic vs. non-asymptotic or convex vs. sigmoidal) have different theoretical and conservation implications (Lomolino, 2000;Triantis et al., 2012). The latter tends to be undertaken using the power model, of the form S = c*A z , where S and A are richness and area, respectively, and c and z are fitted parameters Rosenzweig, 1995). Several studies have tested for systematic variation in ISAR slopes (e.g. Matthews et al., 2021;Matthews, Rigal, et al., 2019;Rosenzweig, 1995;Triantis et al., 2012). However, there have been no comparable analyses of variation in the zvalues (slopes) of IFDARs and IPDARs.
Many island systems have been particularly affected by extinctions and the introduction of non-native species (herein 'introduced species') (Blackburn et al., 2016;Boyer & Jetz, 2014;Hume, 2017;Matthews et al., 2022;Matthews & Triantis, 2021;Whittaker & Fernández-Palacios, 2007). Recent work on the impacts of humans on island biogeographic patterns has illustrated how the exclusion of extinct species and the inclusion of introduced species can affect the form of ISARs Cardoso et al., 2010;Helmus et al., 2014), but how these decisions affect other types of IDAR is less understood .
Here, we undertake the first comparative synthetic evaluation of IDARs at the global scale using a collection of 51 avian archipelago data sets representing different island types (true and habitat), encompassing 1051 individual islands and 2111 species. True islands are those surrounded by water (i.e. oceanic, continental-shelf, continental fragments and lake islands), while habitat islands are those surrounded by contrasting terrestrial matrices (e.g. forest fragments surrounded by pasture; Matthews, 2021). True island data sets were further split into volcanic oceanic archipelagos, a subset of true island data sets comprising archipelagos of mainly volcanic origin never connected to continental land masses (all currently isolated from the mainland by >100 km), and other true island archipelagos (e.g. continental-shelf islands, inland islands). For all bird species (extant [native and introduced] and extinct), we collected nine continuous trait measurements. In combination with phylogenetic data, we constructed the ISAR, IFDAR and IPDAR for all data sets. We used null models to generate FD.ES and PD.ES values and explore how these scale with island area. We also provide the largest global assessment of the impacts of species introductions and extinctions on the IDAR, thus furthering our understanding of the 'island biogeography of the Anthropocene' (Helmus et al., 2014). Figure 1 provides an overview of the methodological F I G U R E 1 The analytical workflow followed, linked to the four primary research questions. Top row: we used the full presence-absence matrix for a given data set (used to calculate island species richness, S), alongside a functional dendrogram (used to calculate island functional diversity, FD) and a phylogeny (used to calculate island phylogenetic diversity, PD). Second row: we fit a set of 20 diversity-area relationship (IDAR) models to the island area, SR, FD and PD data (here the coloured lines represent the fit of the power model; blue = ISAR, red = IPDAR, and yellow = IFDAR), and assessed variation in the slope of the power model. Third row: for each data set, we used a null model to calculate island FD and PD effect sizes (ES) independent of richness. For each data set, we fitted a linear model to the log 10 (Area)-ES relationships, comparing it with an intercept model. We used the ES values to test the association between FD.ES significance and PD.ES significance at the island level (+ = significantly positive ES value; − = significantly negative ES value; NS, non-significant ES value). Bottom row: for subsets of data sets, we ran the analyses three times-the historic fauna including extinct species (A), and the current fauna excluding (B) or including (C) introduced species. framework employed. We used this framework to answer four primary questions: Q1: Do richness, FD and PD scale with area in different ways (i.e. do different models provide the best fit to the ISAR, IFDAR and IPDAR)?
Q2: Does the power model slope differ between the ISAR, IFDAR and IPDAR for a given archipelago, and what are the archipelago characteristics that determine such variation?
Q3: To what extent does island functional and phylogenetic community assembly depart from random expectation and do assembly processes vary with island area in a systematic fashion?
Q4: To what extent does the inclusion or exclusion of extinct and introduced species affect different IDAR properties?

Theoretical expectations
For each of the four primary questions above, we developed a theoretical expectation based on previous research on IDARs: Q1: We expect asymptotic models to provide relatively better fits to IFDAR and IPDAR data, compared with ISAR data, due to the previously reported finding of increasing functional and phylogenetic similarity between species (often interpreted as redundancy) with increasing area (e.g. Dias et al., 2020;Mazel et al., 2014).
Q2: For the same reason as in Q1, we expect IFDAR and IPDAR power model slopes to be systematically less steep than ISARs.
Q3: The Equilibrium Theory of Island Biogeography (MacArthur & Wilson, 1967) assumes in its simplest form that species are functionally equivalent and thus represents a null model of island assembly. By extension, there should be no relationship between richnesscorrected FD and PD (ES) values and island area (Ross et al., 2019;Si et al., 2017), the increase in FD and PD with area being simply a function of richness. However, MacArthur and Wilson (1967) recognized entirely random assembly to be simplistic and subsequent work suggests that the relative importance of different traits and assembly processes could potentially vary along the island area gradient, thus influencing the scaling of FD.ES and PD.ES values with island area.
In theory, community assembly may depart from random towards either clustering or overdispersion of traits. Considering true islands, small islands tend to contain a limited and simpler array of habitat types and more extreme abiotic conditions (Chen et al., 2020;Ross et al., 2019;Sfenthourakis & Triantis, 2009). As a result, only a subset of closely related taxa with specific traits are adapted to these conditions and can persist (Liu et al., 2020;Schrader et al., 2021;Si et al., 2017). This should lead to a degree of functional and phylogenetic clustering on small true islands, consistent with some recent empirical analyses (e.g. Matthews et al., 2020;Ross et al., 2019;Schrader et al., 2021;Si et al., 2017;Zhao et al., 2020).
Conversely, larger true islands will often support a broader range of habitats and potential niches (Whittaker & Fernández-Palacios, 2007), allowing a wider set of species to be able to colonize and persist, leading to neutral or overdispersed patterns . Should it be general that island assembly patterns shift from clustering to random/overdispersion along the area gradient, we should then expect a positive relationship between FD.ES and PD.ES and area for true islands. A similar logic applies to habitat islands, where high habitat heterogeneity in large fragments (e.g. due to topographical variation or the presence of environmental gradients; dos Anjos et al., 2022) can support a broader range of bird guilds (e.g. Willrich et al., 2019). However, we predict less consistent patterns for habitat islands in general, which tend to be much noisier systems (Matthews, 2021).
Q4: Regarding the inclusion of extinct species in oceanic true island data sets, we predict that IDAR slope will increase from the historic period to the current period. This prediction is based on the conceptual model of Franklin & Steadman (2008; see also Steadman, 2006) that was developed in the context of land birds on tropical oceanic islands, whereby, within an archipelago, most species are predicted to have occurred on each high elevation island above a minimum size prior to human colonization, and contemporary positive ISARs are mostly the result of species being harder to drive to extinction on larger islands (e.g. due to larger population sizes and more refugia). We predict the slope of contemporary IDARs should increase with the addition of introduced species, as larger islands are known to experience more introductions . We also predict that extinctions and introductions will have dampened the theoretically expected slope of the ESarea relationships (Q3). This is because extinction and introduction are typically non-random processes, involving species with particular traits (e.g. large body size in regard to extinct species) and from certain taxonomic groups (Boyer, 2008;Fromm & Meiri, 2021;Matthews et al., 2022), which together would act to reduce FD.ES and PD.ES values (i.e. reduce overdispersion and increase clustering), particularly on the larger islands.

Data collection
We sourced true and habitat island bird data sets from the literature. For most data sets, we used previous synthetic ISAR analyses (e.g. Matthews et al., 2021;Triantis et al., 2012) to locate potential data sets and returned to the source papers (and subsequent papers by the source | 969 paper authors) to obtain the species lists for each island. True island data sets were also supplemented using Baiser et al. (2017) and Sin et al. (2022), and habitat islands using Chase et al. (2019). For the former, we updated some of the data sets using a range of literature sources (see Appendix S1). For a number of true island cases (the Ryukyus Islands, the Azores, Canaries, New Zealand) we created new data sets through comprehensive literature and database searches (Appendix S1). For inclusion, data sets needed to contain at least seven islands (to enable the calculation of AIC c , discussed below) and possess an accessible bird species list for each island. An exception was made for a Hawaiian data set (Baiser et al., 2017) which only had six islands, as its extreme isolation means it has particular value in representing isolated oceanic archipelagos. Note that two data sets sourced from Baiser et al.
(Society Islands and Cook Islands) classify atolls (collections of small islets connected by sand banks) as individual islands (Appendix S1). All island areas were converted to km 2 . For the analyses, it was necessary to impose a criterion of a minimum of one species on an island, leading to the removal of a small number of islands with zero species (these were only present in a handful of data sets).
As a first step, for each data set and using either data provided by the source paper authors or using species range maps provided by the IUCN Red List (IUCN, 2021), we classified all species as native or introduced to that archipelago (or region for habitat islands). We then excluded from analysis all introduced species (but for some data sets created alternative versions with introduced species included; see below). Otherwise, we used the data sets as originally published in the source papers, meaning that the exact types of species included varied slightly between data sets due to the decisions of the original source paper authors (e.g. including / excluding marine and nocturnal species). However, to roughly standardize the data sets, we also created an alternative version of each by removing the marine, coastal, wetland and riverine species to produce a land birds only version, for which we re-ran the analyses (see Appendix S2 for details). This standardization process involved removing two data sets when analysing just land birds as it resulted in several islands in these archipelagos having no or very few species (Appendix S2). We removed extinct species (when present) from the data sets, but also created alternative versions of certain data sets with extinct species included (discussed below). For each data set, we formatted all species names, including extinct species (see Appendix S2), to match the nomenclature in the phylogenies provided by Jetz et al. (2012) (see Appendix S1).

Data set characteristics
For each data set (archipelago), and using only the islands/species present in the data set, we recorded a number of variables predicted to affect IDAR form (see Matthews et al., 2021;Matthews, Rigal, et al., 2019;Triantis et al., 2012): (1) number of islands (Ni), (2) the ratio between the area of the largest and smallest island (AreaScale), (3) archipelago species richness (Gamma), (4) total archipelago land area (ArchArea), (5) annual mean temperature and (6) maximum island elevation. FD and PD Gamma were calculated as the total FD or PD of an archipelago. For each true island data set, we also calculated (7) isolation (distance) from the mainland and (8) intra-archipelago isolation (MeanDist). Appendix S2 details how these variables were sourced and calculated.

Functional traits
For functional traits, we sourced data for all of the world's 9993 species (BirdTree taxonomy) from the AVONET trait data set (Tobias et al., 2022), allowing us to build a functional space using all of the world's birds and ensure distances between species in functional space represented the best estimates of the true distances. We used eight continuous morphological measurements: (1) total beak length (from the tip to the skull), (2) beak length to the nares, (3) beak width and (4) depth (at the nares), (5) wing length, (6) secondary length, (7) tail length and (8) tarsus length. These measurements have been shown to provide accurate information on the functional role and trophic status of birds at the global scale (Pigot et al., 2020). We also sourced body mass estimates (g) for each species from AVONET (Tobias et al., 2022). The four kiwi species (Apteryx) represent extreme outliers in terms of the wing length, secondary and tail length traits (e.g. for wing length, the kiwis had values 267 times smaller than the species with the next smallest wing length). To avoid these species affecting the functional space to an extreme degree (which occurred even when logtransforming the traits), for these three traits, we replaced the trait values for the four kiwi species with the mean values across all extant species excluding the kiwis. This approach was preferred to the option of simply removing the kiwis, as one of our analysed data sets comprised islands in New Zealand.
Four of the extinct species in our data sets were also in BirdTree and AVONET. For the remaining 154 extinct species in our data sets, we sourced data for the same set of traits (described below). Our final trait data set comprised 10,147 species. All nine traits were log-transformed and then scaled to have a mean of zero and unit variance.
Because the eight morphological traits are correlated with body mass, we also re-ran the analyses using bodysize corrected traits, generated by running eight simple linear regressions with body mass as the predictor and a given morphological trait as the response (both log-transformed). Here, the scaled residuals from each model were then used as the new trait along with logtransformed and scaled body mass.

Calculating FD, PD and effect sizes
We used FD and PD metrics based on summing branch lengths to ensure our diversity metrics shared the same mathematical framework and are thus directly comparable (i.e. they incorporate the sum of the differences in diversity accumulated between species; Dias et al., 2020;Mammola et al., 2021;Tucker et al., 2017). In addition, the use of trees allowed us to include islands with few species (e.g. one or two), which is not possible with FD metrics such as convex hulls when multiple traits are used (Jarzyna et al., 2021;Petchey & Gaston, 2006). For FD, we built a global dendrogram comprising all 10,147 species. A Euclidean distance matrix was generated using all species and the nine trait axes. We then transformed this distance matrix into a dendrogram using the agglomerative hierarchical clustering method UPGMA (Petchey & Gaston, 2006). We checked the dendrogram quality using the tree.quality function in the 'BAT' R package (Cardoso et al., 2015). The values for our dendrograms were relatively high (0.70 and 0.95 for the dendrogram using the uncorrected and bodysize corrected traits, respectively; one corresponding to maximum quality of the functional representation). For each island, we used the global dendrogram to calculate Petchey & Gaston's (2006) FD metric (including the tree root) using the 'picante' R package (Kembel et al., 2010).
For PD, we based our analyses on the BirdTree phylogenetic trees from Jetz et al. (2012) using the Ericson backbone tree with 9993 species. We obtained a posterior distribution of 3000 trees from BirdTree and created a maximum clade credibility tree (node heights = median heights) including all bird species, using the TreeAnnotator program (v1.10.4, Drummond & Rambaut, 2007). The resultant consensus tree had a small number of negative branch lengths which we resolved by converting negative branch lengths to zero, while shortening only the two branches immediately below by the same absolute amount to ensure the tree remained ultrametric and there were no polytomies (we have added this functionality to the 'BAT' R package; tree.zero function). The PD values generated using the original consensus tree and the consensus tree with the negative branches removed were highly correlated (Pearson's r = 0.999). The 154 extinct species not in BirdTree were grafted on to this consensus tree (detailed below). We used this global maximum clade credibility tree to calculate Faith's PD metric (including the tree root; Faith, 1992) for all islands in a data set as outlined for FD. As a sensitivity check, we re-ran the analyses using a randomly selected tree from the 3000 (grafting the extinct species onto this selected tree).
As both FD and PD can be correlated with species richness, to calculate standardized FD and PD values we created a variant of the 'taxa. labels' null model (999 iterations) and the ses.pd function in the 'picante' R package. This null model worked by only shuffling the names of species found in a given data set on the global tree/ dendrogram (i.e. the null model, for a given data set, uses the archipelagic species pool, not the global species pool, but does not prune the tree). We did this to ensure a consistent tree (i.e. the global tree) was used for calculating FD/PD across data sets, given that pruning the tree was found to affect DAR slopes in a small number of cases (full details provided in Appendix S2).
Generally, standardized values of FD and PD are calculated using standardized effect sizes (SES). However, SES assume a normal distribution of null values, an assumption that is often violated, particularly where some samples contain most, or very few, of the species in the pool. Thus, we instead used the effect size (ES) approach used in Matthews et al. (2020 ). This works by calculating the empirical probability (P) that the observed value is less than expected using the formula: where null is the vector of null distribution values, obs is the observed value and n is the number of null model iterations (here n = 999). This probability was then probit transformed to obtain the ES value (see Appendix S2 for further details). This process was done using both FD and PD, resulting in FD.ES and PD.ES values for each island in each data set. Positive ES values >1.96 were considered to represent cases of significant functional/phylogenetic overdispersion, and negative ES values < −1.96 were taken to represent significant clustering. Non-significant ES values (−1.96 < ES <1.96) were considered to represent random community structure.

IDAR multimodel comparison
For each data set, we fitted twenty SAR models (see Table S2 in Appendix S2) to our three diversity variables (species richness, FD and PD) using least squares nonlinear regression and the 'sars' R package . These models represent a range of curve shapes (linear, convex-upward, sigmoidal), number of model parameters (2-4) and properties (asymptotic and non-asymptotic) (see Triantis et al., 2012 for a review).
We designed a grid search method for selecting model starting parameter values to be used in the non-linear regressions; this method has now been added to the 'sars' package (version 1.3.5; available from CRAN) (see Appendix S2 for details). For each twenty-model set, models were compared and ranked using Akaike's information criterion corrected for small sample size (AIC c ) (Burnham & Anderson, 2002). As the denominator in AIC c must not be negative, for the Hawaii data set (with only six islands) it was necessary to exclude the two four-parameter models from the model set. For each data set and diversity metric, we stored the model ranks, and a multi-model curve was constructed using the AIC c weights from all converged model fits (see Matthews,  . In each case, we also extracted the z-value and c-value from the non-linear power model fit; this model converged in all cases. To ensure the same models were fitted in all cases, we did not remove model fits based on residual assumption checks (e.g. normality). However, given that the least square parameter estimates equal the maximum likelihood estimates only under the assumption of normal errors with constant variance (see discussion in the 'sars' package vignette), and given we are using AIC c , we re-ran the model selection including checks for both these properties. As an additional sensitivity test, we also re-ran the various power model z-value analyses using the z-value from the log 10 -log 10 (linear) power model.
The majority of these twenty models were originally chosen to fit SAR data Triantis et al., 2012) based on the expected shape of the SAR (e.g. convex-upward or sigmoidal). However, the work on island FD.ES-area and PD.ES-area relationships to date (e.g. Diaz et al., 2020;Matthews et al., 2020;Schrader et al., 2021) has shown that they are not well characterized by such model shapes and thus it is not necessarily appropriate to fit the same set of models to these data. Instead, based on our theoretical expectations, we compared the fit of a linear regression model with an intercept-only null model in semi-log space (i.e. log area but not richness) using AIC c . We used a semi-log transformation (log 10 ) as there is no a priori reason to logtransform ES values and it has previously been shown to be an effective method for assessing ES-area relationships Schrader et al., 2021).

Exploratory modelling of IDAR slope variation
First, we tested whether there were significant differences between the z-values of the ISAR, IFDAR and IPDAR. As the IDAR z-values within a data set were not independent, we compared the z-values between data sets using a generalized linear mixed effects model (beta family and logit link; fit using restricted maximum likelihood) and the 'glmmTMB' R package (Brooks et al., 2017), with diversity type (i.e. richness, FD, PD) as a categorical fixed effect, and the data set as a random effect. We used the same approach to compare the slopes from the FD.ES and PD.ES-area relationships, except here we used the Gaussian family as many of these slopes were negative.
Second, we assessed what archipelago characteristics drove variation between data sets in the (i) z-value of the ISAR, IFDAR and IPDAR, and (ii) slope of the ES-area relationships. Following Marx et al. (2017), we undertook an exploratory modelling approach using Pearson's correlations, which were preferred over regression analyses due to the relatively small size of our data set. The z-values of one of the ISAR, IFDAR or IPDAR were correlated against each of our archipelago-level predictor variables in turn, using log-transformations when necessary to meet assumptions. As the IFDAR and IPDAR z-values were correlated with the ISAR z-value, we also ran a series of partial correlations with these two variables, allowing us to control for the ISAR z-value. We ran the modelling using all data sets (i.e. true and habitat islands) and using as predictor variables: Ni, Gamma (or FD / PD Gamma), AreaScale, ArchArea, (maximum) elevation, temperature and the power model c-value. We then re-ran the modelling using just the true island data sets and adding in as predictors both MeanDist and isolation from the mainland. We then re-ran these correlation tests but instead used the slope of the linear model fitted to the FD.ES-area and PD.ES-area relationships as variables.

The effect of including extinct and introduced species on diversity scaling relationships
For ten true island data sets (Canaries, Cook Islands, Hawaii, Lesser Antilles, Marianas, Society Islands, Cape Verde, New Zealand, Azores, and Ryukyu Islands), there were a relatively large number of species introduced to each archipelago (ranging from 11 to 60% [median = 19%] of the total contemporary archipelago bird fauna). For these data sets, we also created alternative versions representing the current faunas with introduced species included (we only considered currently established introduced species; see Appendix S1). For the first eight of those, we were also able to build data sets representing the historic fauna, that is, the island composition around 1500 CE, including extinct species and extirpated extant species. For five data sets where (coarse) data were available (Hawaii, Marianas, Cook Islands, New Zealand, Canaries), we also built data sets representing the prehistoric fauna (i.e. prior to human colonization of the islands; including all species known to have gone extinct in the last ~125,000 years) excluding marine species (Appendix S2). For the Marianas and Cook Island prehistoric data sets, we removed a number of islands as we decided to focus on islands where more fossil data were available.
The historic and pre-historic data sets (i.e. including extinct species) were built using a range of literature sources (see Appendix S2 for details). The functional traits of extinct species were initially sourced through measurements made on specimens in various museums and literature searches. For 135 of the 158 extinct species, we were able to acquire at least one measurement from skin or skeleton (or both) specimens in museums, with body mass being estimated for the remaining species (see Appendix S2). All gaps were then imputed using Bayesian Hierarchical Probabilistic Matrix Factorization (Schrodt et al., 2015). We ran the imputation ten times, averaging the imputed values across the ten runs. As a sensitivity test, we re-ran the analyses using a randomly selected individual imputation run rather than averaging. Extinct species were also grafted onto our consensus phylogeny. Appendix S2 provides a detailed description of the extinct species data collection, trait inference and phylogeny grafting.
We then re-ran the power and linear model fitting for the historic and introduced species data sets, storing the power model z-values for the ISAR, IFDAR, IPDAR and the slopes (of the linear model) of the two ES-area relationships. We compared the values with those from our main analyses (i.e. current fauna excluding introduced species) using paired Wilcoxon signed-rank sum tests. We also re-fitted the models using the prehistoric data sets.
Unless otherwise stated, all analyses were undertaken in R (Version 4.2.0; R Core Team, 2019), and the analyses were run on a 500GB cluster using 51 cores (~2000 core-hours).

R E SU LT S
In total, we sourced 51 data sets (26 true island and 25 habitat island archipelagos), incorporating 1051 islands and 2111 species (1953 extant and 158 extinct species) (Appendix S1). The size of habitat islands ranged from 0.004 to 1592 km 2 , and true islands from 0.001 to 150,437 km 2 . A map of the locations of these data sets is provided as Figure S1 in Appendix S3. All best fit models, and the power and linear model parameters, for all five relationships across all data sets are provided in Table S3 in Appendix S3.

Q1 and Q2: ISAR, IFDAR and IPDAR model form
The non-asymptotic convex-upward Kobayashi, power and logarithmic models were always the three models with the highest mean AIC c weight values for the ISAR, IFDAR and IPDAR (but not always in the same order), across the 51 data sets (Figure 2). Inspecting the plots of model fits provided further evidence for the convexupward nature of most of the ISARs, IFDARs and IPDARs (e.g. Figure 3). In terms of the number of best fits (i.e. cases of lowest AIC c for a given IDAR and data set), the top model was always the power model, with the linear, logarithmic and Kobayashi models alternating in F I G U R E 2 Generally, the shape of IDARs consistently had a convex-upward nature with some variation in exact model shape between the ISAR, IFDAR and IPDAR. Regarding the ES-area relationships, the intercept model had the higher mean AIC c weight for both FD and PD, meaning a lack of relationship between FD.ES and PD.ES and area for many data sets. The bar charts show the mean model AIC c weights across all data sets in which a model fit converged, for the five IDARs. The total number of data sets is 51. Full model names can be found in Table S2. For the FD.ES-area and PD.ES-area relationships, the two models were fitted in semi-log space, for the other IDARs in untransformed space. second and third position ( Figure S2). The results were similar when looking at true and habitat islands separately ( Figures S3-S6 in Appendix S3).
The power model provided a reasonable approximation of the form of the three IDARs (mean R 2 across all data sets and the three IDARs = 0.62). In general, for a given data set, the z-value of the ISAR was larger than that of the IPDAR, which was slightly larger than that of the IFDAR, and these differences became more pronounced the steeper the ISAR was (Figure 4). Using a mixed-effect model with the diversity type as a fixed effect and the data set as a random effect revealed that the z-values significantly differed between the ISAR (mean z = 0.19), IFDAR (mean z = 0.14) and IPDAR (mean z = 0.16) (Type II Wald χ 2 test for the categorical fixed effect, χ 2 = 163.5, P < 0.001). This was also the F I G U R E 3 Some island systems exhibited positive FD.ES and PD.ES-area relationships, and others negative relationships. The top two rows show the IDARs of a data set of birds (number of species = 54) in 11 true islands in the Galápagos, generated using five diversity metrics: species richness (ISAR), functional diversity (IFDAR), phylogenetic diversity (IPDAR), and the FD (FD.ES-area) and PD (PD.ES-area) effect sizes. The bottom row shows the two ES-area relationships for a data set of birds (number of species = 101) in 77 true islands in the Aegean (Simaiakis et al., 2012). In the top row plots, the different coloured lines represent the fits of up to twenty competing models, and the thick black line represents a multi-model averaged curve generated using the AIC c weights of the individual model fits. In the middle and bottom row plots (left and middle), the dark green line is the fit of a standard linear model, while the light grey line is the fit of an intercept-only model. For the FD.ES-area and PD.ES-area relationships, the two models were fitted in semi-log 10 space, for the other IDARs in untransformed space. Increasing ES values from zero denote greater overdispersion, while decreasing values from zero denote greater clustering. The two bird photographs show example species from each archipelago: the middle right plot shows a lava gull (Larus fuliginosus), the rarest gull in the world and a species endemic to the archipelago, and the bottom right plot shows a Rüppell's warbler (Sylvia rueppelli), a species that breeds in Greece, Turkey and the Aegean Islands. Middle row photograph by Andy Morffew and under licence (https://creat iveco mmons.org/licen ses/by/2.0/); bottom row photo by Mick Sway and under licence (https://creat iveco mmons.org/licen ses/by-nd/2.0/). case when considering only true islands or only habitat island data sets. Figure 5 provides the results of the exploratory modelling of correlations between IDAR slopes and archipelago features. When considering all data sets, ISAR, IFDAR and IPDAR slopes were significantly positively correlated with Ni (number of islands), and significantly negatively correlated with (maximum) elevation and temperature. When considering only the true island data sets, there was still a significant negative correlation between elevation and the slopes of the three IDARs. There were also negative correlations with ArchArea, although for the ISAR this was not significant. When controlling for ISAR slope, there were no significant correlations, either for all data sets or just true island data sets.

Q3: Avifaunal community assembly: FD and PD effect sizes and their scaling relationships
The avifauna of most islands (87% for FD, and 79% for PD) exhibited random structure regarding FD.ES and PD.ES values, with a small proportion being characterized as significantly clustered (11% for FD and 19% for PD). Very few island avifaunas were significantly overdispersed (2% for both metrics). Mean ES values were −0.51 for FD and −0.89 for PD, indicating a slight tendency toward clustering ( Figure S7a). FD.ES and PD.ES significance results were equivalent for most islands, but there were notable exceptions ( Figure 6); for example, 126 of the islands had significantly negative PD.ES values, but non-significant FD.ES values ( Figure 6).
Across all data sets, the intercept-only model had the higher mean AIC c weight and provided the best fitting candidate model the most times, for both the FD.ES and PD.ES-area relationships (i.e. lowest AIC c in 34 and 37 out of 51 data sets, respectively). However, there were notable exceptions, with some FD.ES and PD.ES-area relationships exhibiting positive and negative linear relationships ( Figure 3). When looking at true and habitat islands separately ( Figures S3-S6 in Appendix S3), it was apparent that, for true islands, the relative performance of the linear model, regarding both ES-area relationships (but particularly PD.ES), improved.
Considering cases where the linear model provided the best fit, there were nine positive and eight negative relationships for the FD.ES-area, and eight and six respectively for the PD.ES-area relationship. The majority of significant linear cases were true island data sets (11 cases for both the FD.ES and PD.ES relationships) (see Appendix S3). The median slope of the linear model across all data sets was 0.02 (−0.03 and 0.03 for true and habitat island data sets, respectively) for the FD.ES-area and 0.15 (0.13 and 0.15) for the PD.ES-area relationship ( Figure S7b). Interestingly, when only focusing on the ten volcanic oceanic island data sets, the median linear slope values were higher: 0.35 and 0.55 for the FD.ES and PD.ES-area relationships, respectively (see Figure S8 in Appendix S3).
The slope values from the FD.ES and PD.ES-area relationships significantly differed according to a

F I G U R E 4
The difference between the ISAR z-value and the IPDAR and IFDAR z-values increases with increasing ISAR z-value. The figure shows the relationships between the z-values of three IDARs, plotted as a function of the ISAR z-value rank (higher rank = steeper ISAR): the ISAR (black lines and points), the IPDAR (blue lines and points) and the IFDAR (red lines and points). Different symbols are used for habitat (circles) and true island (triangles) systems. The z-values were generated from fitting the non-linear power model to the bird IDARs of 51 island data sets. mixed-effects model when considering all data sets together (χ 2 = 4.5, p = 0.03), but not true and habitat islands separately. Considering all data sets, there were significant positive correlations between the FD.ES and PD.ES-area relationship slopes and (maximum) elevation ( Figure 5). Considering only the true island data sets, there were significant positive correlations between both slopes and MeanDist, isolation and elevation, and a significant negative correlation between PD.ES-area slope and AreaScale ( Figure 5).

Q4: The effect of including extinct and introduced species on diversity scaling relationships
The power model z-value for the ISAR, IFDAR and IPDAR followed an interesting and relatively consistent pattern across the three data set types: historic fauna (A), current fauna excluding introduced species (B) and current fauna including introduced species (C) (Figure 7). For these three IDAR types, z decreased or remained roughly constant between A and B, and then generally increased between B and C. This pattern was stronger for certain data sets (e.g. Society Islands, Marianas) compared to others (Figure 7). The (paired) Wilcoxon signed-rank tests indicated that the differences between A and B were significant for the ISAR (p = 0.03), IFDAR (p = 0.02) and IPDAR (p = 0.02). The differences between B and C were also significant for all three IDAR types (p = 0.02, 0.02 and 0.04 for the ISAR, IFDAR and IPDAR, respectively), while the differences between A and C were non-significant. For the FD.ES-area and PD.ES-area relationship slopes, there were significant decreases in slopes between A and B (p = 0.01 and 0.04) and A and C (p = 0.04 and 0.01), but the differences between B and C were not significant (p > 0.05) (Figure 7).
Comparing models for the prehistoric and current avifaunas (excluding introduced and marine species) for five data sets, the z-values decreased or remained relatively constant for the ISAR, IFDAR and IPDARs, with the exception of Hawaii, for which z-values increased (Figure 8). For the FD.ES and PD.ES-area relationships, with two exceptions (Marianas for FD.ES and New Zealand for PD.ES) the slope of the relationships decreased between the two time periods (Figure 8).

Sensitivity analyses
The full results of all sensitivity analyses are presented in Appendices S4-S7. First, re-running the analyses using body-size corrected traits to construct the functional dendrogram resulted in very similar findings (Appendix S4). Second, re-running the analyses after subsetting the data sets to only include land birds also generated mostly similar results (Appendix S5). The main differences here related to the exploratory correlations (e.g. no significant associations involving the number of islands or isolation), and the introduced and extinct species analysis: while some data sets followed the same pattern as the main results, the general pattern was less clear and none of the ISAR, IFDAR or IPDAR F I G U R E 5 Some characteristics of archipelagos are correlated with IDAR slopes. The figures show Pearson's rank correlation heatmaps, with IDAR slope on the x-axis and various archipelago-level predictors on the y-axis. For the ISAR, IFDAR and IPDAR, slope was measured as the power model z-value. Correlations for IFDAR and IPDAR z-values were also undertaken using partial correlation using ISAR zvalues as a covariate. For the ES-area relationships, the slope was the slope of a linear model fitted in semi-log 10 space. Correlations were undertaken twice, once using all 51 data sets (a), and once using only the 26 true island data sets (b). Significant coefficient values (p < 0.05) are indicated using black circles. Cell colour indicates correlation strength. Grey cells indicate a correlation was not undertaken for that variable combination. Predictor acronyms are GA, Gamma; C, power model c-value; AA, ArchArea; AS, AreaScale; NI, number of islands; MD, MeanDist; IS, isolation; EL, elevation; and TP, temperature. Note that for the IFDAR and IPDAR correlations, Gamma was the total functional or phylogenetic richness.
paired Wilcoxon tests were significant, although this is perhaps expected given the smaller number of data sets involved (Appendix S5). Third, undertaking the model selection using residual assumption checks resulted in very similar results (Appendix S6). The power model passed the assumption checks for 41, 41 and 45 data sets for the ISAR, IFDAR and IPDAR, respectively. Using the z-values from the linear (log 10 -log 10 ) power model also generated similar results (Appendix S6). Fourth, using a randomly selected trait imputation run in combination with an individual Jetz et al. (2012) phylogeny resulted in very similar findings (Appendix S7).

Q1 and Q2: The form of island diversity-area relationships
In general, and in contrast to our prediction (Q1) that asymptotic models would provide a better relative fit to the IFDAR and IPDAR, the three island diversityarea relationships (IDARs; i.e. the island species-area relationship [ISAR], the island functional diversityarea relationship [IFDAR] and the island phylogenetic diversity-area relationship [IPDAR]) were all best modelled by non-asymptotic convex-upward models, although the linear model provided the best fit in certain cases (see also Triantis et al., 2012). Inspection of the model fit plots (Figure 3) also showed that the form of the three primary IDARs was generally convex-upward. This matches the recent findings of a study on habitat islands by Dias et al. (2020), but for a much larger number of data sets and broader range of island types.
As expected (Q2), we observed that, for a given data set, the ISAR was generally steeper than the IPDAR, which was in turn steeper than the IFDAR. These results indicate that as island area increases, more species are sampled from the archipelagic pool. These additional species initially add novel traits and phylogenetic branches to the island communities, but this process slows down with increasing richness as an increasing proportion of these species are functionally and to a slightly lesser extent phylogenetically, redundant (see also discussion in Karadimou et al., 2016;Schrader et al., 2021;Ferreira-Arruda et al., 2022).
Our exploratory correlation modelling indicated that IDAR slope was significantly negatively associated with elevation, when focusing on all data sets and true islands data sets separately. For true islands, this may seem counterintuitive as many of the relatively isolated (mostly oceanic) archipelagos (e.g. Hawaii, Cape Verde, New Zealand) have high elevation, and previous studies have theorized and shown that ISAR slope increases with isolation (Whittaker et al., 2017). This pattern could be specific to birds: due to their relatively high dispersal ability and the fact that many of the oceanic archipelago data sets are lacking very small islands, it is possible that many bird species are present on most islands, thus lowering IDAR slope. There could also be an effect of anthropogenic extinctions given that these archipelagos are also those that have likely experienced the most extinctions, and our results indicate that in many cases (when excluding introduced species) these extinctions have lowered IDAR slope. For the true island data set correlations, while few associations were significant, OLS regression models including all predictors explained a relatively large proportion of the variation in the slope of the ISAR, IFDAR and IPDAR (R 2 values: 0.72-0.76; adjusted R 2 values: 0.56-0.63). This matches the results of previous studies (e.g. Matthews, Rigal, et al., 2019;Triantis et al., 2012) and suggests that the lack of significant correlations here may be due to smaller sample sizes and thus a lack of power.

Q3: Island community assembly patterns and processes
Overall, we found that the majority of island avifaunas were classified as being randomly assembled in terms of functional diversity (FD) and phylogenetic diversity (PD), although a sizable minority were significantly clustered (11-19%). This could indicate that neutral dynamics predominate on most islands, as assumed within the core model of island biogeography (MacArthur & Wilson, 1967). It should be noted that our null model used the archipelago species list as the pool rather than a wider (mainland) species pool, as we were not focused on testing for the effects of mainland to island filters (Santos et al., 2016;Triantis et al., 2022). Nonetheless, if these results were viewed in isolation, it would be tempting to conclude that there were no patterns of interest beyond the 'null' observation that most islands had random functional and phylogenetic structure. However, analysis of the scaling of these assembly patterns (effect size [ES] values) reveals a more complex picture, at least in certain cases. For many data sets, the relationships between ES values (assembly patterns) and island area are indeed relatively flat. This indicates that, for these data sets, the convex-upward scaling of unstandardized FD and PD with area was primarily a result of increasing richness with island area, rather than changes in the dominant community assembly processes. However, there were numerous exceptions to this pattern, particularly regarding the FD.ES-area relationship, where the linear model provided a better fit for a third of the 51 data sets.
Based on previous work (e.g. Chen et al., 2020;Matthews et al., 2020), we had hypothesised that, owing to limited habitat availability, smaller true islands would be characterized by functional and phylogenetic clustering. In contrast, larger true islands, with a wider range of habitat types, were expected to display functional and phylogenetic neutrality or overdispersion (Carvajal-Endara et al., 2017;Matthews et al., 2020). Together, this would result in a positive linear relationship between richness-corrected FD and PD and island area (i.e. less clustering with increasing area). However, approaching half of the data sets where the linear model provided the best fit exhibited a negative relationship (i.e. more clustering with increasing area). This is the opposite of our theoretical prediction (Figure 3), but has been observed previously, such as for exotic plants in US National Parks  and mammals on oceanic islands (Si et al., 2022). One explanation for this pattern can be found in Diamond's (1975; for a review see Whittaker & Fernández-Palacios, 2007) work on assembly rules, which argues that very small islands can only support one bird species per guild (e.g. one fruit pigeon) due to limited niche space and increased interspecific competition. As island area increases, the number of species per guild or habitat type (which will be relatively functionally redundant) able to coexist on an island also increases. If the amount of niche space and number of guilds increases with area at a slower rate than for the number of species, this will increase the amount of functional and to a lesser extent phylogenetic redundancy on islands. Following this logic, smaller island assemblages would be expected to be overdispersed, with clustering increasing with island size, ultimately resulting in a negative relationship between richness-corrected FD and PD and island area (see also Si et al., 2017). Future research could test this theory by analysing the density of species per guild or habitat type for archipelagos that have negative ES-area relationships. The scaling of speciation rate with island area (see Whittaker & Fernández-Palacios, 2007) could theoretically also result in this pattern, at least for archipelagos where speciation is a source of new bird species. Specifically, if a small number of colonizers radiate on the larger islands into numerous closely related species without substantial trait disparification, this could result in increased clustering on the larger islands.
Interestingly, boxplots of the linear slope values across island types indicate that the average slope of both the FD.ES and PD.ES-area relationships was larger on oceanic islands (relative to other true islands and habitat islands; Figure S8), particularly for PD.ES, and, in our exploratory correlations restricted to true islands, isolation, elevation and MeanDist had positive correlations with the slope. Thus, it could be the case that our theoretical prediction is more applicable to large isolated oceanic island systems, than to other island types. However, it should be noted that, while the median slope was relatively high, for most of the oceanic data set ES-area relationships the best model was in fact the intercept-only model. This could partly be because most of these data sets have relatively few data points (a characteristic of many oceanic archipelagos), reducing the power of the test and increasing the effect of noise in the data (e.g. due to anthropogenic impacts). For habitat islands, there was a larger proportion of cases where the intercept-only model provided the best fit. This likely reflects the fact that they are often relatively noisy systems (Matthews, 2021) and can vary substantially in terms of various properties. For example, large forest fragments may contain a range of habitat types or be relatively homogenous (e.g. see Figure 1 in Willrich et al., 2019;dos Anjos et al., 2022).
It is important to note that there are several limitations associated with the community assembly framework . Each pair of coloured circles joined by a black line represents two different data sets for the same archipelago: pre-historic fauna including extinct species, and current fauna excluding introduced species. In this analysis, all marine species were removed from the data sets prior to model fitting and some islands were removed due to a lack of fossil data (i.e. the Modern z-values may differ from those in Figure 7). Note the different y-axis scales in the different plots. used here (see Münkemüller et al., 2020, for a review). These include (i) a focus on specific patterns ignores the reality that a given FD / PD pattern can be produced by multiple processes (Mayfield & Levine, 2010); and (ii) that, as mentioned above, defining the species pool as the archipelago species list can only underestimate filtering by excluding species unable to reach or persist on the archipelago (Carvajal-Endara et al., 2017;Si et al., 2022), but (iii) also underestimates competition by ignoring 'dark diversity' (i.e. species excluded from the archipelago due to past competition will not be present in the pool; Münkemüller et al., 2020).

Q4: The effect of anthropogenic introductions and extinctions on IDARs
Humans have introduced hundreds of species to islands . Our results supported our prediction (Q4) that the inclusion of introduced species would generally lead to steeper IDARs. These increases in z-values were relatively modest in absolute terms but were statistically significant. It is known that larger islands tend to have higher rates of anthropogenic colonization pressure and thus experience more introductions Blackburn et al., 2021), which will, all else being equal, have the effect of steepening the ISAR.
Particularly on certain oceanic island archipelagos, human colonization precipitated a wave of avian extinctions (Boyer, 2008;Boyer & Jetz, 2014;Matthews et al., 2022;Sayol et al., 2021;Triantis et al., 2022;Whittaker & Fernández-Palacios, 2007). However, while recent work has started to look at the effects of species extinctions on island diversity and FD (e.g. Boyer & Jetz, 2014;Matthews et al., 2022;Si et al., 2022;Sobral et al., 2016;Triantis et al., 2022), how these extinctions have affected IDARs has not been fully evaluated. Interestingly, and in contrast to our theoretical prediction (Q4), we found that the slopes of the ISAR, IFDAR and IPDAR significantly decreased between the historic and current period (excluding introduced species). This could indicate (i) that extinctions were more prevalent on the larger islands due to greater human impact, or (ii) that there is a bias resulting from greater knowledge of the historic fauna on larger islands. One other caveat is that Franklin & Steadman's (2008) conceptual model was based on all extinctions, but for most archipelagos we lack adequate data at the island-level for species that went extinct prior to 1500 CE. However, many archipelagos are known to have suffered numerous extinctions prior to 1500 CE (Hume, 2017;Sayol et al., 2021;Steadman, 2006). Our analysis of five data sets that did include pre-1500 extinct species (Figure 8) broadly confirmed the decrease in the ISAR, IFDAR and IPDAR slope between the pre-human colonization avifauna and the current avifauna excluding introduced species, although Hawaii was an exception (Figure 8).
Taking both the above findings together, it appears in several cases that introduced species have 're-calibrated' the IDARs such that the slopes are more like the historic period including extinct species; indeed, there were no significant differences in the z-values between the historic and current (with introduced species) assemblages. A similar pattern is observed with the power model R 2 values, with higher values observed for the historic (average R 2 of power model across ISAR, IFDAR and IPDAR = 0.66) and introduced data sets (0.63), with lower values for the current fauna without introduced species (0.53), for all three IDARs. Interestingly, Hawaii again provides an exception to this pattern, with the inclusion of introduced species lowering or not changing IDAR slopes. This could indicate that introductions to that archipelago (or at least the islands in the archipelago that comprise our data set) are occurring more independently of island area.
We also observed that extinctions resulted in a statistically significant decrease in slope for the FD.ES and PD.ES-area relationships (Figure 7), a pattern also apparent in the analysis of the five pre-1500 data sets ( Figure 8). One interpretation of this is that anthropogenic extinctions are leading to more random patterns of community assembly, or even greater clustering due to the selective extinction of certain types of species (e.g. large-bodied; Boyer, 2008;Boyer & Jetz, 2014;Hume, 2017;Matthews et al., 2022).
A notable caveat is that the prehistoric and historic data sets analysed here likely underestimate the true island composition at these time periods. First, there are known biases in the (sub)fossil record, such as largebodied species being more likely to leave material evidence than small-bodied species. Second, the fossil record is likely incomplete for almost all islands. Third, several studies present data on which islands extinct species occurred on, with fewer presenting data on the past distributions of extant species. Finally, the trait estimation and imputation, and phylogeny grafting, procedures obviously involve a certain degree of uncertainty. Appendix S2 provides a more detailed discussion of these issues.

Concluding remarks
Overall, we have shown that increasing richness with island size is the main driver of the IPDAR and IFDAR for most data sets, although there are numerous exceptions to this pattern. We also find that archipelagos with the steepest ISARs exhibit the biggest differences between ISARs and the IFDAR/IPDAR. These results indicate that, within a given archipelago, there is an increasing amount of functional and phylogenetic redundancy on larger islands. As a next step, it is necessary to test whether the patterns observed here are consistent across taxonomic groups, particularly those with lower dispersal ability compared to birds, given that ISAR slope has been shown to vary between taxa (Triantis et al., 2012). In addition, as more data on extinct island species distributions become available, it will be necessary to evaluate further how anthropogenic extinctions, in combination with introductions, have affected IDAR form and slope, and whether this re-calibration effect is a general pattern. This will ultimately improve our knowledge of the 'island biogeography of the Anthropocene'. All authors contributed to writing, editing and discussion of ideas.

AC K NO W L E DGE M E N T S
The comments of Anderson Saldanha Bueno, Jonathan Chase and three anonymous reviewers greatly improved the paper. The computations described in this paper were performed using the University of Birmingham's BlueBEAR HPC service, and the data compilation was supported by the University's GEES Research Support Fund. FS was supported by a Beatriu de Pinós postdoctoral fellowship (2020 BP 00067) from the Ministry of Research and Universities of the Government of Catalonia.

PE E R R EV I EW
The peer review history for this article is available at https://www.webof scien ce.com/api/gatew ay/wos/peerrevie w/10.1111/ele.14203.

OPE N R E SE A RC H BA DGE S
This article has earned Open Data and Open Materials badges. Data and materials are available at: doi. org/10.5281/zenodo.7632639.

DATA AVA I L A BI L I T Y STAT E M E N T
All code and data are available on GitHub (txm676/ DARs) and archived on Zenodo (doi.org/10.5281/ zenodo.7632639).