Identifying the Drivers of Spatial Taxonomic and Functional Beta-Diversity of British Breeding Birds

Spatial variation in community composition may be driven by a variety of processes, including environmental filtering and dispersal limitation. While work has been conducted on the relative importance of these processes on various taxa and at varying resolutions, tests using high-resolution empirical data across large spatial extents are sparse. Here, we use a dataset on the presence/absence of breeding bird species collected at the 10 km × 10 km scale across the whole of Britain. Pairwise spatial taxonomic and functional beta diversity, and the constituent components of each (turnover and nestedness/richness loss or gain), were calculated alongside two other measures of functional change (mean nearest taxon distance and mean pairwise distance). Predictor variables included climate and land use measures, as well as a measure of elevation, human influence, and habitat diversity. Generalized dissimilarity modeling was used to analyze the contribution of each predictor variable to variation in the different beta diversity metrics. Overall, we found that there was a moderate and unique proportion of the variance explained by geographical distance per se, which could highlight the role of dispersal limitation in community dissimilarity. Climate, land use, and human influence all also contributed to the observed patterns, but a large proportion of the explained variance in beta diversity was shared between these variables and geographical distance. However, both taxonomic nestedness and functional nestedness were uniquely predicted by a combination of land use, human influence, elevation, and climate variables, indicating a key role for environmental filtering. These findings may have important conservation implications in the face of a warming climate and future land use change.


INTRODUCTION
Biodiversity is currently facing a multitude of global-scale threats from human activity (Dirzo et al., 2014;McGill et al., 2015). As the human footprint on the natural world grows, it is becoming increasingly important to understand how these factors are impacting ecological communities in order to inform conservation efforts and make predictions about impacts under future scenarios (Newbold, 2018;Soininen et al., 2018). Analyzing spatial variation in species diversity is a powerful means of assessing the impact of different environmental factors on biodiversity as it provides us with information on what is currently limiting species ranges and occupation of sites. The analysis of spatial variation in biodiversity generally involves focusing on taxonomic changes between sites in the form of alpha (α) or, to a lesser extent, beta (β) diversity (Field et al., 2009;Calderón-Patrón et al., 2016;Soininen et al., 2018).
Comparing the alpha diversity of two communities separated in space provides a measure of the difference in the number of species between these sites but ignores species replacement/turnover (i.e., a species being extirpated from a site and another species colonizing) and can therefore mask biodiversity change (Gonzalez et al., 2016). In contrast, betadiversity provides a measure of community dissimilarity between sites (Whittaker, 1960;Koleff et al., 2003). Various metrics have been proposed to measure beta-diversity, which can be grouped into variance-based approaches (the focus of the present study) and diversity-partitioning approaches (Legendre and de Cáceres, 2013;Matthews et al., 2019). Recently, several variance-based metrics (e.g., the Sørensen index) have been partitioned into constituent components, such as species replacement/turnover (that is independent of richness differences) and species richness differences or 'nestedness' (Baselga, 2010). It has been argued that the study of these partitions provides insight into the drivers of compositional difference between sites (Baselga and Leprieur, 2015). Nestedness in this context is not 'true' nestedness [e.g., as measured by the nestedness metric based on overlap and decreasing fill (NODF)], but rather nestedness resultant dissimilarity that allows for the separation of dissimilarity due to turnover from that of nestedness (Baselga, 2012). For ease, we henceforth use the term 'nestedness' to describe nestedness resultant dissimilarity.
Standard beta-diversity metrics (herein termed 'taxonomic beta-diversity') assume all species are the same in terms of the role they play within an ecosystem (Sekercioglu, 2006), thereby ignoring the vital role that functional diversity plays in assemblage dynamics (Devictor et al., 2010;Eskildsen et al., 2015). Recently, several beta-diversity metrics have been expanded to incorporate functional information (Cardoso et al., 2014;Baselga and Leprieur, 2015) and can be used to shed light on the biotic/abiotic factors that cause variation in functional diversity between sites (Villéger et al., 2013;Cardoso et al., 2014). This evidence can be used to inform conservation activities such as protected area design and biological corridor selection (Socolar et al., 2016), and help protect ecosystem services (Şekercioǧlu et al., 2004;Cardinale et al., 2012;Galetti et al., 2013). For example, bird communities play an important role in providing services such as pollination, pest control, and carrion removal (Whelan et al., 2008;Wenny et al., 2011), which are critical services for humans and other taxa. Using a beta-diversity measure that incorporates species traits, and hence functionality, is therefore essential to gain a better understanding of biodiversity change and its consequences (Devictor et al., 2010;Jarzyna and Jetz, 2018;Tobias and Pigot, 2019;Carvalho et al., 2020).
Biological communities are predicted to vary spatially in the absence of anthropogenic influences (i.e., the natural pattern of distance decay in species similarity; Nekola and White, 1999) due to both dispersal limitation and niche filtering, among other factors (Lomolino et al., 2010). Dispersal limitation is hypothesized to impact the spatial variation in community composition by restricting the range of species through distance alone (Dambros et al., 2017), independent of environmental differences between the communities (Hubbell, 2001). Niche-filtering occurs when environmental gradients constrain communities to those species adapted and able to persist in local conditions (Weiher and Keddy, 1999;Cornwell et al., 2006). In addition, human-induced change (e.g., land-use change and climate change) are likely important drivers of both spatial taxonomic and functional beta-diversity (Devictor et al., 2007;Davey et al., 2012;Barnagaud et al., 2017).
The effect of dispersal on spatial beta-diversity can be assessed through the analysis of the geographic distance between sites, while the effect of environment can be tested by evaluating measures of habitat and land use types (hereafter 'land use') and climate (Luck et al., 2013;Wieczynski et al., 2019;Fluck et al., 2020). However, due to the spatial structuring of environmental gradients (i.e., a distance decay in environmental conditions), it is difficult to partition the unique effects of each (Leibold and Chase, 2017). Relating current land use and climate to spatial variation in community composition will also enable inferences to be made on how increases in the relative intensity of these drivers may impact spatial variation in the future (Barbet-Massin and Jetz, 2015).
The effect of land use on community composition is mostly a result of niche filtering, where species which are adapted to a specific land use type are unable to survive in contrasting land use types (Weiher and Keddy, 1999;Cornwell et al., 2006). Rapidly growing human populations (Tratalos et al., 2007) have facilitated considerable land use changes via increasing urbanization in some areas of the world (Seto et al., 2012), and conversion of natural land to agriculture, as well as an overall intensification of agricultural practices (Zabel et al., 2019). There is strong evidence that these practices have disrupted communities, leading to pools of generalist species in heavily modified habitats via the exclusion, through filtering, of species with narrower environmental requirements (i.e., specialists) (McKinney, 2006;Vellend et al., 2007;Clavel et al., 2011;Flohre et al., 2011;Barnagaud et al., 2017;Hagen et al., 2017). Thus, in anthropogenic landscapes (such as those that occur across much of the United Kingdom), turnover is predicted to be low across large spatial scales, and communities are predicted to become more nested, with high redundancy in functional diversity (Liang et al., 2019;Weideman et al., 2020).
Climate is also an important environmental filter and can drive high spatial beta-diversity between regions due to differences in energy availability (energy richness hypothesis; Hutchinson, 1959;Currie, 1991;Hurlbert and Haskell, 2003), variation in species tolerance (physiological tolerance hypothesis; Root, 1988), and variation in speciation rates (the speciation rates hypothesis; Currie et al., 2004;Hua and Wiens, 2013). While there is mixed support for these hypotheses, substantial evidence exists showing that differences in species composition between sites are often correlated with climatic gradients (Currie et al., 2004). Additional filtering impacts may occur where land use change interacts with climate (Auffret and Thomas, 2019).
Better understanding the role of climate in driving spatial beta-diversity is essential in order to accurately predict the effects of future climate change on community composition. For example, species populations have been found to expand or contract their ranges in response to changing climatic conditions (Fox et al., 2014;Batt et al., 2017). As well as range shifts, shifting phenologies across the trophic web can lead to disruption of communities through cascade effects due to altered species interactions (Bell et al., 2019). Other impacts of a changing climate, such as more frequent severe weather events, are also increasingly recognized as significant drivers of spatial betadiversity (Maxwell et al., 2019).
There are thus many potential drivers of spatial beta-diversity. However, few studies exist assessing the relative roles of these different drivers (e.g., land use, climate, human impacts) in terms of both taxonomic and functional beta-diversity, the aim of the present study. We use generalized dissimilarity models (GDMs) in combination with a dataset containing presence/absence data of British breeding birds collected at the 100 km 2 scale over the entirety of the British Isles. We aim to (1) test the effect of geographic distance and a range of environmental variables (e.g., land use type, climate) on the spatial taxonomic and functional beta-diversity of British breeding bird communities and (2) evaluate the role of human influence on spatial beta-diversity patterns. While we do not set out to test the niche-filtering and dispersal-limitation hypotheses directly, we interpret increasing dissimilarity due to geographical distance as an indication that dispersal limitation may play a role in the structuring of communities. In contrast, increasing dissimilarity due to climate, land use, or human influence would point to a role for niche filtering.

Species Composition Data
Data showing the summer (breeding) distributions of the British avifauna (Gillings et al., 2019) were collected during April-July over the period 2008-2011 (BA2010) by volunteers on behalf of the British Trust of Ornithology (BTO) and the Scottish Ornithologists' Club (SOC). Some fieldwork effort was permitted out of this field season, with specific instructions given on what evidence was permitted (see Gillings et al., 2019, for further information). FIGURE 1 | Map displaying the study location (Great Britain) and its location within Europe. The lower proportion of the island is gridded with the quadrats used to sample the avifauna (only a sample is shown here and the whole island was sampled).
The dataset summarizes the presence/absence of British bird species within 10 km × 10 km (100 km 2 ) quadrats covering the British Isles on a continuous grid (Figure 1). Only species designated as being "confirmed" or "probable" breeders (Gillings et al., 2019) were retained here. Vagrant and pelagic species were excluded, but we retained introduced breeding species for the analysis (McInerny et al., 2018). While some introduced species' occupied ranges may reflect in part their initial introduction sites, many of these species are now established, so their presence exerts an influence on community structure, resource use, and competition (Lennon et al., 2000). Species under threat from human persecution (particularly hunting or egg-collecting) were also removed from the analyses as data for such species were provided at larger spatial grains (i.e., larger than 100 km 2 ) or their locations were omitted entirely (Gillings et al., 2019). All quadrats with less than 50% land and all island regions that were considered disconnected from the mainland were removed. A total of 2257 100 km 2 quadrats remained with a species pool containing 169 species (Supplementary Table 1).

Trait Data
Continuous morphometric variables were measured from museum specimens or extracted from literature and used to characterize the functional diversity of each community (defined as all the species present in each quadrat). We selected eight morphological traits to represent the functional role of birds: two estimates of beak length (culmen from tip-to-skull and tip-to-nares), beak width, beak depth, tarsus length, wing chord length, tail length, and body mass, with evidence showing all of these traits provide useful information about avian dietary niche, locomotion and ecological function (Trisos et al., 2014;Tobias and Pigot, 2019;Pigot et al., 2020). Further information on measurements, including sampling per species and methods, are published separately (Pigot et al., 2020).
A principal components analysis (PCA) was performed using all eight traits, with the full eight axes extracted. All the axes were then standardized to a mean of 0 and a standard deviation of 1, producing a trait matrix (species x traits) with eight trait axes for each species. To test for the effect of raw trait variability, we also standardized the traits prior to running the PCA. The results using both the raw traits and the standardized traits in the PCAs were comparable, so we report only the results using the raw traits within the PCA here. All eight axes were included as it has been shown that all the axes provide useful information, with even minor axes capturing significant variation in traits (Pigot et al., 2020).

Climate Data
Monthly temperature and precipitation data were downloaded from the United Kingdom Met Office, which provides climate data interpolated from local weather stations onto a 1 km × 1 km grid across the United Kingdom (Hollis et al., 2019). Data were downloaded from the period 2000-2011. For the breeding season (defined as the start of May to the end of July, as many arriving migrants in April will not yet be breeding), key climate variables were selected a priori, and averages calculated. Precipitation (mm) was summed for each 100 km 2 quadrat over the breeding season for each year. The average temperature ( • C) was calculated as the daily average temperature across the quadrat and the breeding season. The range in temperature was also calculated as the average mean maximum daily temperature over the breeding season minus the average mean minimum daily temperature. The mean of each of these climatic variables was then calculated over the 2000 -2011 period to reduce the influence of yearly variation, leaving three measures of climate (Tavg mean , Prec mean , and Range mean ). Climate averages were also constructed over a more extended period  to test if birds were responding to longer-term climate variation. The majority of the GDMs fitted using the shorter period had a better fit, and thus we only report the results using 2000-2011 climate here.

Land Use Data
Data on land use were obtained from the EDINA environment digimap service for 2007 (Land Cover Map, 2007). These data provide land cover (23 land use classes) for the British Isles at a 25 m scale. From these data, the percentage cover for each land use within each 10 × 10 km quadrat was calculated. The woodland classes (coniferous woodland and broadleaved woodland) were grouped into one variable named 'forest, ' as were 'grasslands' (grouped from the improved grasslands and semi-natural grasslands categories), and 'urban areas' (grouped from the suburban and urban categories). Arable land was also included as a predictor variable. Shannon's diversity index was calculated for each quadrat as a measure of land-use heterogeneity (hereafter called Shan).

Human Influence Index (HII)
The Human Influence Index (HII) was used to assess the contribution of human impact on the variation in community composition (Wildlife Conservation Society-WCS, and Center for International Earth Science Information Network- CIESIN-Columbia University, 2005). The HII is derived from multiple data sources on population density, infrastructure (railroads, urban development, night-time lighting), and landcover ranging in date collected 1994-2005 (although in this version about half of the measures were collected around 2000 instead of 1995, as was the case for the first version). The measures are each weighted differently in the methodology and then standardized giving a measure of human impact ranging from zero (no human impact) to 100 (maximum human impact possible using the methodology). HII values were extracted from each of the 1 km 2 grid squares within each 100 km 2 quadrat. The average was taken over these values to obtain the mean HII within each quadrat. It is important to note here that there is some temporal disparity in the period the HII was developed over (see above) and the period the atlas was conducted (2008)(2009)(2010)(2011). However, even with this small disparity, HII should provide a robust indication of the impact human influences have on spatial variation in taxonomic and functional composition.

Elevation Data
Elevation data were obtained from the shuttle radar topography mission (SRTM; Jarvis et al., 2008). For each 100 km 2 quadrat, data were extracted using 400 equally spaced points. The mean (Mean elev ) and the standard deviation (SD elev ) were then calculated from these data as measures of elevation and variability in the elevation across the area.

Testing for Multicollinearity
Pearson's and Spearman's correlations were used to test for multicollinearity between the predictor variables. SD elev was removed due to the variable being strongly correlated with multiple other variables (Elevation, Tavg mean , and Prec mean ). The climatic variables were found to be collinear with one another, and with other variables (Supplementary Figures 1a,b). Therefore, the climatic variables were combined using a PCA. The PCA yielded three axes that explained all the variation of the original three climatic variables [hereafter; Climate 1 (81.90%), Climate 2 (13.54%), and Climate 3 (4.56%)]. All the axes were retained in the models to capture all the variability that could be explained by climate. Scatter plots and correlations between the PCA axes and the raw climate variable showed that Climate 1 was positively correlated with average temperature and negatively correlated with average precipitation (Supplementary Figure 2). Climate 2 and Climate 3 had less clear and more complex correlations with the original variables (Supplementary Figure 2).
After substituting the climate variables with the PCA axes, all variables had correlations < 0.70, with two exceptions: Elevation and Climate 1, and urban land use and HII ( Supplementary  Figures 1c,d). As a result, urban land use was removed from the analysis and the human influence index, which is a composite measure including urban land use, was retained. Both elevation and Climate 1 were retained because (1) both variables are known to be significant predictors of spatial variation in breeding avian communities, (2) the correlation was still below 0.8, and 3) GDM is known to be robust to multicollinearity to a certain degree (Glassman et al., 2017). A variance inflation factor test (VIF) was also performed, with all remaining variables having VIF values < 5 (Neter et al., 1983;Gareth et al., 2013).

Spatial Taxonomic and Functional Beta-Diversity
To assess taxonomic dissimilarity between the assemblages, pairwise taxonomic beta diversity was calculated for each 100 km 2 quadrat using the function beta.pair from the package 'betapart' (Baselga and Orme, 2012). This function computes the dissimilarity (here measured using Sørensen's dissimilarity index, βsor; Baselga, 2010;Koleff et al., 2003) between an assemblage and every other assemblage present in the dataset to create a pairwise dissimilarity matrix.
Using Sørensen's dissimilarity, total beta-diversity can then be partitioned into its two constituent components: dissimilarity due to turnover (BD TURN ) and nestedness resultant dissimilarity (BD NEST ), with BD TOTAL = BD TURN + BD NEST . Turnover is the proportion of dissimilarity due to species replacement between two assemblages, whereas nestedness is the proportion of the dissimilarity due to one assemblage being a nested subset of another assemblage through either species loss or gain (Baselga, 2010). It is important to note that, unlike BD TURN , BD TOTAL and BD NEST are not independent of species richness changes, as the measurements are dependent upon species richness gradients (Baselga and Leprieur, 2015).
A measure of functional beta-diversity was then calculated using Sørensen's dissimilarity index and Baselga's (2010) partitioning framework (Phylosor). For this approach, a global functional dendrogram was created containing all the United Kingdom breeding species, using a Euclidean trait distance matrix and the agglomerative hierarchical clustering method (UPGMA). This method produces a rooted tree with a constant weight assumption (i.e., where the distance between the root to all tips is equal), and this then describes the functional relationships between species (Petchey and Gaston, 2002). The phylo.sor function in the 'betapart' package (Baselga and Orme, 2012) was used to calculate functional dissimilarity based on the shared branch length of the functional dendrogram between each assemblage and every other assemblage (hereafter called FD TOTAL ). Although this method is usually used on phylogenies, here it is used on a functional dendrogram to give a functional measure analogous to taxonomic beta-diversity, allowing for straightforward comparison. In addition, using a convex hull approach (the standard Baselga functional metric) was not possible here due to the size of the dataset and the computational demands of such an approach. FD TOTAL was also partitioned into its constituent components of nestedness resultant dissimilarity (FD NEST ) and turnover (FD TURN ).
A Pearson's correlation was performed between the Euclidean distances (in the trait distance matrix) and the cophenetic distances (in the dendrogram) . The resultant correlation was high (Pearson's r = 0.97), indicating that the dendrogram provides an adequate measure of the functional distances between species.

MNTD (Mean Nearest Taxon Distance) and MPD (Mean Pairwise Distance)
As an alternative to Baselga's functional beta-diversity framework, mean nearest taxon distance (MNTD) and mean pairwise distance (MPD) were calculated. While MNTD is also sensitive to species richness differences, MPD is a measure that is mostly independent of species differences between sites (Miller et al., 2017). MNTD represents the mean distance (smallest nondiagonal value) between species in a community and is most sensitive to changes at the 'tips' of a dendrogram (Webb et al., 2002). MPD is a similar measure but is calculated as the mean between all non-diagonal elements between species within a community (Webb, 2000;Webb et al., 2008), and so it is more sensitive to changes at the roots of the functional dendrogram. Here, the beta-diversity versions of MNTD and MPD, that calculate the same measures but between assemblages are used. For ease, we refer to these as MNTD and MPD (Miller et al., 2017). MPD and MNTD were calculated using the comdist and comdistnt functions, respectively, in the R package 'picante' (Webb et al., 2008). MNTD and MPD were standardized by dividing each pairwise measure by the largest pairwise measure to produce dissimilarity bounded between 0 and 1. Standardizing MPD and MNTD in this way allowed the measures to be modeled using GDMs.

Modelling Variation in Spatial Beta-Diversity
As a first step, multidimensional scaling (MDS) was applied to each of the pairwise measures. The first axes from the MDS were taken and plotted. These were then assessed visually for any pattern in the dissimilarity/similarity between assemblages.
Generalized dissimilarity modeling (GDM) was then used to model functional and taxonomic beta-diversity. GDM is a statistical technique (an extension of matrix regression) that can be used for assessing the relationship between environmental gradients and variation in community composition (Ferrier et al., 2007). The modeling accommodates non-linearity that is present in ecological datasets over large extents (Ferrier, 2002;Ferrier et al., 2007) and can also incorporate geographical distance. This is vital to include, as dispersal limitation modulated by distance is known to be an important driver of community composition (Keil et al., 2012). GDM can deal with higher multicollinearity among predictor variables than many commonly used regression models (Glassman et al., 2017), and uses monotonic I-splines that constrain the coefficients of the regressions to be positive for a non-decreasing fit and non-positive for a non-increasing fit. The I-splines allow the evaluation of predictor effects on the dissimilarity metrics through the height and slope. The maximum height represents the total deviance explained by the predictor while holding all the other predictors constant, while the slope displays the rate of compositional change across the predictor's range (Fitzpatrick et al., 2013;Fitzpatrick and Keller, 2015). We applied a modeling framework using GDM, aimed at assessing the unique and shared roles of both geographic distance and environmental factors on our measures of beta diversity (Supplementary Figure 3).
A separate GDM was fitted with each of the taxonomic and functional beta-diversity metrics (i.e., BD TOTAL , BD TURN , BD NEST , FD TOTAL , FD TURN , FD NEST , MNTD, and MPD) as response variables using the 'gdm' package in R (Fitzpatrick et al., 2020). As a first step, matrix permutation was used to assess the model significance and variable importance using the gdm.varImp function. Due to the large amount of memory and time required to model the full site-pair combinations (N = 2,545,896) only a subset of the data could be used for the matrix permutation and variable importance process (as recommended for datasets with a large number of sites when calculating variable importance; Fitzpatrick et al., 2020). First, 60% of the sites were randomly removed, leaving 407,253 sitepairs. Another 60% of the site-pair combinations were then removed from the remaining site-pairs, leaving a total of 162,901 site-pairs for analysis. The removal of the site-pairs, after the initial site removal, removes site-pairs randomly but does not remove all site-pair combinations (Fitzpatrick et al., 2020). The predictor variables were not scaled, which allows assessment of the impacts each of the predictors has along actual environmental gradients (e.g., Fitzpatrick et al., 2013;Heino et al., 2019).
For the matrix permutation process, a GDM was first run using all of the predictor variables. The rows of the environmental data were then permuted, and a GDM model was fitted to those data. Significance was then evaluated by comparing the deviance of the model with unpermuted data to the model using permuted data. Variable importance was assessed by permuting each of the variables in turn while holding the other variables constant (unpermuted). Variable importance was then assessed as the difference in deviance explained using the permuted and unpermuted variable, with more important variables explaining a larger proportion of the deviance when unpermuted. The process was then repeated after dropping the least important variable (backward elimination) with variable importance and significance recalculated. All variable importance scores reported in the text were calculated at this point. The first model where all variables were significant (p < 0.05) was identified, and this model was then fit using all the sites (i.e., the full dataset) (final model). One hundred permutations were used for model significance testing and variable importance scores (Ferrier et al., 2007;Heino et al., 2019). Uncertainty in the I-splines was then evaluated using a bootstrapping approach (Shryock et al., 2015). A total of 50% of the sites were first removed randomly from the dataset. From these data, a GDM model was fit, and the I-spline coefficients extracted. For the bootstrapping, a further 80% of the sites were randomly removed, and a GDM model fit. The process was repeated 100 times. The I-splines were then plotted with error bands showing the standard deviation from the permutation process.
Variance shared between geographical distance and the environmental variables was calculated for each model in turn using the formula: where V s is the explained variance shared between the geographic predictor and the environmental variables, V full is the variance explained by the full model, V g is the variance explained by the geographic distance only model, and V e is the variance explained by the model containing only the environmental variables (Ray-Mukherjee et al., 2014).

Taxonomic and Functional Beta-Diversity of the British Avifauna
BD TOTAL was higher on average (0.313 ± 0.142) than FD TOTAL (0.285 ± 0.126). In both cases, overall beta-diversity was determined mostly by the turnover component (0.225 ± 0.121 (71.88% of the total) and 0.195 ± 0.108 (68.42% of the total), for BD TURN and FD TURN , respectively). Nestedness was responsible for a smaller proportion on average for both measures (Figures 2A,B). Average MPD (0.799 ± 0.045) was higher than MNTD (0.255 ± 0.125) (Figure 2C), as is to be expected.

Spatial Variation in Taxonomic and Functional Beta Diversity
Heat maps of the community metrics showed clear spatial patterns in the different beta diversity metrics. A northsouth divide was present for BD TOTAL ( Figure 3A) and FD TOTAL (Figure 3E), and for the turnover component of each ( Figures 3B,F, respectively). Alternatively, this could also be interpreted as a divide along the classic 'Tees-exe line' that roughly divides the uplands and lowlands of Britain (Prakash and Rumsey, 2018). The west of Wales was more similar to the north of England and parts of Scotland than the south of England for total beta diversity and turnover for both taxonomic and functional metrics, representing elevation changes between upland and lowland regions (Figures 3A,B,E,F). The eastern coastal regions in Scotland were more congruent with southern regions than with inland and west coast Scottish assemblages (Figures 3A,B,E,F). Nestedness patterns generally mirror the other metrics but the patterns were patchier (i.e., the divides between regions (north/south and Tees-exe line) was less delineated). The south-west of England and the west coast of Wales were a closer match with most of Scotland than they were with the majority of southern England (Figures 3C,G), again closely matching the 'Tees-exe line' (Prakash and Rumsey, 2018).
The pattern for MNTD was largely the same as that found for total taxonomic and functional beta-diversity and turnover ( Figure 3D). MPD showed a pattern similar to that of nestedness for both taxonomic and functional beta diversity ( Figure 3H). FIGURE 2 | Boxplots of pair-wise dissimilarity measures of community composition for British breeding birds. These are all on the same row (A) displays taxonomic spatial dissimilarity (turnover, nestedness resultant dissimilarity, and total), (B) is functional pairwise dissimilarity (turnover, nestedness resultant dissimilarity, and total), and (C) shows standardised mean pairwise distance (MPD) and standardised mean nearest taxon distance (MNTD). The horizontal line within the box represents the median, the box indicates the inter-quartile range (IQR), and the whiskers show data 1.5 times the IQR. Points highlight outliers.
FIGURE 3 | Heat maps of community dissimilarity in British breeding birds based on the first axis from a Principal co-ordinate analysis analysis (PCoA) for taxonomic and functional beta diversity, turnover, and nestedness resultant dissimilarity. Colors represent the ordering scores obtained from the PCoA, with areas displaying similar colors more similar and areas with differing colors less similar in terms of community composition. The first three maps on the first row (A-C) are for taxonomic beta-diversity (BD TOTAL , BD TURN , and BD NEST , respectively) and the first three maps on the second row (E-G) are for functional beta-diversity (FD TOTAL , FD TURN , and FD NEST , respectively). Mean nearest taxon distance (D) and mean pairwise distance (H) are the last maps on each row, respectively.

Modelling the Drivers of Spatial Beta-Diversity
Overall, the variance explained by all the final models (bar those for MPD, FD NEST , and BD NEST ) was high, with between 55.61% and 68.45% variance explained ( Table 1). The variance explained by the final GDM models for MPD, FD NEST , and BD NEST was lower (between 9.60% and 12.09%) ( Table 1). The deviance explained for the functional metrics was lower than for the taxonomic metrics (Table 1).
For all final models, the deviance explained by GDM models using only geographical distance as a predictor of dissimilarity was lower than the deviance explained by the models run using only the significant environmental variables (Table 1). However, there was overlap in the deviance explained by geographical distance and the environmental variables. Between 19.31% and 23.42% of the variance explained was shared within the MNTD, taxonomic, and functional beta-diversity models, excluding the FD NEST and BD NEST models ( Table 1). For the nestedness components, the geographical distance between sites explained a low percentage of deviance (Table 1). Geographic distance was non-significant in the MPD model.

Drivers of Taxonomic Beta-Diversity
The following results (and those in the 'Drivers of functional beta-diversity' section) relate to the final models (i.e., the models that have been simplified using the permutation approach described in the methods).
BD TOTAL was most impacted by geographical distance, Climate 1, Elevation, and HII (in order of decreasing importance, Table 1). BD TOTAL rose gradually with geographical distance and Also included are the results from models analyzing mean pairwise distance (MPD) and mean nearest taxon distance (MNTD) between sites. Variable importance scores were calculated using a subset of the data through permutation of each one of the predictors in turn, while holding all other predictors constant. The variable importance is then the mean difference in variation described by the model including the non-permuted variable and the permuted variable. Therefore, the higher the importance score the more important that variable is to explaining variation in the community dissimilarity metric. For individual variable descriptions, see the main text.
Elevation, with a sharper rise observed for Climate 1 initially, followed by a leveling out ( Figure 4A). HII had an initial effect on total beta-diversity, with a small but sharp increase observed over initial environmental dissimilarity, but then leveled off and remained relatively constant (Supplementary Figure 1). Arable, forest, and grass cover were also in the final BD TOTAL model (Supplementary Figure 4). Geographical distance, Climate 1, and arable land cover (in order of decreasing importance) had similar impacts on BD TURN as they did on BD TOTAL (Figure 4B). Elevation and forest cover were the two other most important variables in regard to this response variable, with turnover increasing sharply with forest cover initially before leveling off, and a near-linear increase with elevation (Supplementary Figure 5). HII, grass cover, and Climate 2 were also included in the final BD TURN model (Supplementary Figure 5).
Geographical distance was relatively unimportant for predicting BD NEST , although it was included in the final model. The best predictor (the predictors with the highest variable importance values; Table 1) for BD NEST was Climate 1, and its relationship was similar to that found for BD TOTAL (Figure 4C). Forest cover, HII, and elevation were also in the final model (Supplementary Figure 6).

Drivers of Functional Beta-Diversity
FD TOTAL had a similar relationship with geographic distance, Climate 1 and HII as was observed with BD TOTAL, and they were also the three most important predictors ( Figure 5A). Elevation had an almost linear relationship with FD TOTAL (Supplementary  Figure 7). Arable, grass, and forest cover, along with the other two climate axes (Climate 2 and Climate 3), were also included in the final FD TOTAL model (Supplementary Figure 7).
The main difference between the FD TURN and BD TURN final models was the inclusion of Climate 3 for FD TURN , and a slight difference in the ordering of the variable importance (Table 1, Figure 5B and Supplementary Figure 8). FD NEST , as with BD NEST , was also unaffected by geographical distance ( Table 1). FD NEST was impacted by Climate 1, forest cover, HII, climate 3, and geographical distance ( Figure 5C and Supplementary Figure 9).
MNTD was mainly impacted by Climate 1, geographical distance, HII, Elevation, and forest cover (Figure 6A and Supplementary Figure 10). Overall, the final MNTD model contained the same significant predictors as in the final models for BD TOTAL and FD TOTAL , and similar relationships were found with all the predictors, although the order was slightly different (Supplementary Figure 10). MPD was the only measure not predicted by geographical distance (Table 1). MPD was found to sharply increase with Elevation before leveling off around 100 m ( Figure 6B). HII was the next best predictor of MPD and its relationship was different to that found with the other measures, with a curvilinear increase observed midway through the gradient (Figure 6B). The relationship between Climate 1 and MPD was also different to that found between the predictor and the other response variables, with a threshold effect found near the tail end of the gradient (Figure 6B). FIGURE 4 | Plotted I-splines of the three most important variables (determined from the variable importance score) from generalized dissimilarity models analyzing the relationship between environmental and geographic gradients, and spatial community composition in British breeding birds. Plots on row (A) are for total Sorensen's beta-diversity (BD TOTAL ), (B) are for the turnover component (BD TURN ), and (C) are for the nestedness resultant dissimilarity (BD NEST ). Climate 1 is the first axis from a principal component analysis calculated from the average temperature, range of temperature, and precipitation over the months May to July, across the 2000-2011 period. Elevation is the average elevation across each 100 km 2 quadrat. Forest and Arable are the percentage of each land use within each quadrat. The human influence index (HII) is the average human influence across each quadrat. Geographic is the geographic distance between sites. Curves show the relationship between the gradients and community dissimilarity obtained using I-splines. The most important variables are on the left with decreasing variable importance to the right. Blue lines show the I-Spline correlations, with standard deviation (gray shaded area) calculated through bootstrapping (100 permutations) on a portion of the dataset. A rug plot on the x axis shows the spread of the data. Note the varying y-axis for each measure.
It should be noted that a large majority of the predictors were close in their variable importance scores across the models ( Table 1). The geographical distance between sites was the exception across the models, as it was the most important predictor by a large proportion in all but the two nestedness, and MPD/MNTD models. It is also important to note that although multicollinearity was assessed, geographic structuring of certain predictors may mean that the permutation of the variables may not have resulted in large importance values if geographic distance was retained in the model.

DISCUSSION
We found that variation in overall spatial functional and taxonomic beta-diversity (and their turnover components) of British breeding avian assemblages is driven by a combination of geographical distance per se and environmental gradients. The environment-only models explained more deviance than the geographic distance only models, but there was also a relatively large, shared variance explained component. Previous studies have suggested that this shared component be considered as indirect effects of climate (climate distance, Mazel et al., 2017;Qian et al., 2020). This may suggest that geographic distance plays a comparatively small role in predicting compositional differences between assemblages in contrast to environment. This aligns with previous studies that have shown that spatial variation in community composition, for a wide range of taxa, can be attributed to a combination of deterministic (including contemporary and historical), and stochastic factors (Steinitz et al., 2006;Melo et al., 2009;Dobrovolski et al., 2012;Vicente et al., 2014;Glassman FIGURE 5 | Plotted I-splines of the three most important variables (determined from the variable importance score) from generalized dissimilarity models analyzing the relationship between environmental and geographic gradients, and spatial community composition in British breeding birds. Plots on row (A) are for total functional beta-diversity (FD TOTAL ), (B) are for the turnover component (FD TURN ), and (C) are for the nestedness resultant dissimilarity (FD NEST ). Climate 1 is the first axis from a principal component analysis calculated from the average temperature, range of temperature, and precipitation over the months May to July, across the 2000-2011 period. Forest and Arable are the percentage of each land use within each quadrat. The human influence index (HII) is the average human influence across each quadrat. Geographic is the geographical distance between sites. Curves show the relationship between the gradients and community dissimilarity obtained using I-splines. The most important variables are on the left with decreasing variable importance to the right. Blue lines show the I-Spline correlations, with standard deviation (gray shaded area) calculated through bootstrapping (100 permutations) on a portion of the dataset. A rug plot on the x axis shows the spread of the data. Note the varying y-axis for each measure. Carvalho et al., 2020). However, we also found that the differences in species loss/gain between sites (as measured by nestedness resultant dissimilarity) were explained mostly by differences in climate, and to a lesser extent, land use. Overall variation explained by the nestedness models was, however, much lower than for the total beta-diversity and turnover models. The lower explanatory power of the nestedness models implies that other drivers not included here may be impacting dissimilarity between sites due to species loss and gains.

Geographical Distance
Both BD TOTAL and FD TOTAL were driven mostly by the turnover component, highlighting that compositional differences between communities in Britain are mainly a result of different species being replaced. This is consistent with other beta-diversity studies, across multiple taxa (see Soininen et al., 2018). A slow initial increase followed by a steeper rise in the slope of the I-splines of geographical distance in respect to turnover highlights a northsouth divide, or a divide between the 'Tees-exe' line, in community dissimilarity ( Figure 1B). The observed impact of distance could be due to dispersal limitation, geographical barriers, or historical factors (Nekola and White, 1999;Soininen et al., 2007;Dobrovolski et al., 2012;Barnagaud et al., 2017). Here, we expect it is a combination of all these factors, as well as an intertwining of distance with climate and land use. FIGURE 6 | Plotted I-splines of the three most important variables (determined from the variable importance score) from generalized dissimilarity models analyzing the relationship between environmental and geographic gradients, and spatial community composition in British breeding birds. Plots on row (A) are for standardised mean nearest taxon distance (MNTD) and (B) are for standardised mean pairwise distance (MPD). Climate 1 is the first axis from a principal component analysis calculated from the average temperature, range of temperature, and precipitation over the months May to July, across the 2000 -2011 period. Elevation is the average elevation across each 100 km 2 quadrat. The human influence index (HII) is the average human influence across each quadrat. Geographic is the geographical distance between sites. Curves show the relationship between the gradients and community dissimilarity obtained using I-splines. The most important variables are on the left with decreasing variable importance to the right. Blue lines show the I-Spline correlations, with standard deviation (gray shaded area) calculated through bootstrapping (100 permutations) on a portion of the dataset. A rug plot on the x axis shows the spread of the data. Note the varying y-axis for each measure.
FD TURN was lower than BD TURN , indicating that the species that are being turned over across communities share some functional traits. Petchey et al. (2007) found that many co-occurring species in British breeding birds were functionally similar, indicating low functional diversity across Britain as a whole and within individual communities.
The low functional alpha diversity observed by Petchey et al. (2007) may partly explain the lower functional turnover compared to taxonomic turnover observed here.
Nestedness made up a lower proportion of the total beta-diversity for both the taxonomic and functional metrics. Functional nestedness was mostly on par with taxonomic nestedness (Figure 2), highlighting that species lost/gained between assemblages were not functionally redundant (Petchey et al., 2007). The functional distinction of species lost/gained spatially between communities' highlights that there is a significant difference between northern/southern communities. The observed pattern is partially congruent with elevational peaks in Britain (the Tees-exe line; Prakash and Rumsey, 2018). The pattern is also supported by MPD, a measure independent of species richness ( Figure 3H). It is interesting to note that of the variables that were significant in the MPD model, elevation was the most important and had a threshold effect with any communities over ∼100 m being dissimilar to those at lower elevations ( Figure 6B). The results found here align with previous evidence that indicate the importance of elevational gradients as environmental filters (Sanders and Rahbek, 2012;Pigot et al., 2016).

Climate
Climate1 was included in the top three most important variables for all the final models ( Table 1). The impact of climate was similar across all final models (except for the MPD model), with an almost linear relationship observed with dissimilarity before a leveling off. As the correlations with the original climate variables showed a positive correlation with average temperature and a negative correlation with precipitation, this relates to a divide between warmer, drier, and wetter, cooler regions. The near linear relationship between partial ecological distance and Climate 1 for many of the metrics therefore shows a divide in terms of assemblage composition between these two types of regions. Climate was also the best predictor of nestedness by a large margin. This may point to climatic filtering of certain species, with a divide between north and south Britain. However, the overall pattern for both measures was patchy (Figures 2C,G). With global warming predicted to change the seasonality and severity of rainfall within Britain in the coming decades (Watts et al., 2015), these results indicate that drier and hotter summers in the future could lead to shifts in community composition.

Land Use and the Human Influence Index
While the overall impact of the different land-use types differed, the form of the relationship between the various land-use predictors and all the measures of beta-diversity were similar. After an initial sharp increase in community dissimilarity over the first small portion of the range of the environmental gradient, the dissimilarity leveled out. This impact over the initial environmental dissimilarity highlights the difference in community composition between areas with none of that land cover, and those with just a small percentage.
HII impacted spatial variation in communities in much the same way as the land-use predictors. After an initial increase in dissimilarity between communities with increasing HII, there was a leveling off. There are several reasons why this could be the case. For example, Tratalos et al. (2007) found that richness of all species initially increased faster with household densities (one of the measures included in the HII) than urban adapted species, but these then declined significantly after peaking at a very low density. For urban areas, captured within HII, this could highlight the initial homogenization impact of urbanization, with urban areas similar after this initial disturbance driven by an increase in generalists and loss of specialists (Davey et al., 2012). However, as HII is a composite measure, this cannot be conclusively confirmed from these results alone. The human influence captured in this measure may also be masked at this scale as it is likely that a remnant of suitable habitat with lower human disturbance exists within a quadrat, or near to a quadrat (Fattorini et al., 2016). As this study is also only considering presence/absence, it also makes no inferences about abundance, which is likely more sensitive to human influence, and which could be very different between quadrats (Tratalos et al., 2007).

Future Climate and Land-use Change
The importance of both climate and land-use variables points toward potential future disruption of community composition if these drivers increase in intensity as is expected (Seto et al., 2012;Watts et al., 2015). The latitudinal divide between northern and southern (or highland and lowland) community composition also indicates that future warming could see species extirpations/extinctions from the colder, northern regions of the United Kingdom (Tayleur et al., 2016). This would also likely see an increased similarity of the northern communities with southern communities as southern species extend their ranges northward (Hickling et al., 2006). However, species extending their ranges into the UK from Europe could obscure this impact, so studies approaching this question should focus on species identities. Indeed, this has already been observed in the study of bird abundance within England, with resident and short-distance migrants increasing in abundance through time potentially at the expense of long-distance migrants, habitat specialists and coldassociated species (Pearce-Higgins et al., 2015). With many rare species also dispersal limited, future changes in land-use and climate could potentially extirpate some of the few rare bird species Britain has (Baur, 2014).

Limitations
While we consider a range of predictors, we have not included measures of biotic interactions. Competitive interactions and predator presence/abundance can all have an impact on spatial beta-diversity (Wittwer et al., 2015;Koròan and Svitok, 2018). Abundance differences between sites were also not directly considered due to a lack of appropriate data. Given that a species' abundance can be an important determinant of that species' influence on ecosystem functioning (Winfree et al., 2015;Gaston et al., 2018), future studies aiming to analyze spatial variation in community and functional composition should attempt to analyze measures of population size.

CONCLUSION
Spatial variation in both the taxonomic and functional composition of United Kingdom breeding birds is driven mainly by species turnover, which can be explained through a combination of geographical distance per se and environmental gradients. The unique variance explained by distance alone could reflect an important role for dispersal limitation in driving these patterns, but more work is needed, as this variance component could also be due to a process not considered here. In contrast, species loss/gain, observed through nestedness between sites, was driven mainly by environmental factors. Future climate warming and land-use change could lead to an increase in the loss of species, particularly cold-adapted or rare and dispersal-limited, from communities, particularly in the north and in the uplands. With turnover driving these patterns, broad-ranging conservation efforts would be preferable to conservation focused on target areas (Si et al., 2016).
Future work should look for potential synergies between climate and land-use in order to assess if future increases in both could potentially have larger than expected impacts on biodiversity based on the individual effects of each in isolation (Brook et al., 2008;de Chazal and Rounsevell, 2009;Mantyka-pringle et al., 2012;Frishkoff et al., 2016). Future comparison between the results presented here and results of similar tests from areas in different climate regions or in less disturbed regions than the United Kingdom would also be informative. For example, in tropical systems it would be expected that, as the United Kingdom is a postperturbation system and rates of habitat loss will be higher in the tropics, land use will play a larger role in driving community composition dissimilarity between sites than found here (Hansen et al., 2013). Another potential future research direction could be to assess whether other taxonomic groups within the British Isles have similar patterns of functional or taxonomic beta-diversity, or whether the patterns observed here are bird-specific.

DATA AVAILABILITY STATEMENT
Publicly available datasets were analyzed in this study. These data can be found here: http://www.bto.org/datasets.

AUTHOR CONTRIBUTIONS
JW, JS, TP, and TJM conceived the original idea. JT collected the trait data. JW performed the analysis with input from JS, TP, and TJM. JW wrote the manuscript with significant contributions from TJM. All authors contributed to the writing and commented on the draft.

FUNDING
Central OA block grant processed through the University of Birmingham.

ACKNOWLEDGMENTS
We thank British Trust of Ornithology (BTO) volunteers for their hard work in collecting the data. We are also grateful to the Sir Stanley Stapley Trust for financial support (to JW). The computations described in this manuscript were performed using the University of Birmingham's BlueBEAR HPC service.