Distance decay 2.0 - A global synthesis of taxonomic and functional turnover in ecological communities.

Aim: Understanding the variation in community composition and species abundances (i.e., β-diversity) is at the heart of community ecology. A common approach to examine β-diversity is to evaluate directional variation in community composition by measuring the decay in the similarity among pairs of communities along spatial or environmental distance. We provide the first global synthesis of taxonomic and functional distance decay along spatial and environmental distance by analysing 148 datasets comprising different types of organisms and environments. Location: Global. Time period: 1990 to present. Major taxa studied: From diatoms to mammals. Method: We measured the strength of the decay using ranked Mantel tests (Mantel r) and the rate of distance decay as the slope of an exponential fit using generalized linear models. We used null models to test whether functional similarity decays faster or slower than expected given the taxonomic decay along the spatial and environmental distance. We also unveiled the factors driving the rate of decay across the datasets, including latitude, spatial extent, realm and organismal features.
Results: Taxonomic distance decay was stronger than functional distance decay along both spatial and environmental distance. Functional distance decay was random given the taxonomic distance decay. The rate of taxonomic and functional spatial distance decay was fastest in the datasets from mid-latitudes. Overall, datasets covering larger spatial extents showed a lower rate of decay along spatial distance but a higher rate of decay along environmental distance. Marine ecosystems had the slowest rate of decay along environmental distances. Main conclusions: In general, taxonomic distance decay is a useful tool for biogeographical research because it reflects dispersal-related factors in addition to species responses to climatic and environmental variables. Moreover, functional distance decay might be a cost-effective option for investigating community changes in heterogeneous environments.

INTRODUCTION

Biodiversity on Earth is shrinking (IPBES, 2019). Understanding its distribution is therefore paramount to inform conservation efforts and to evaluate the links between biodiversity, ecosystem functioning, ecosystem services and human well‐being (Cardinale et al., 2012). The variation in the occurrence and abundance of species in space and time (i.e., β‐diversity) is at the heart of community ecology and biogeography because it provides a direct link between local (α) and regional (γ) diversity (Mori et al., 2018). Moreover, β‐diversity has become an essential currency in spatial (Kraft et al., 2011) and temporal (Blowes et al., 2019) comparisons of biodiversity patterns and their underlying drivers. The β‐diversity is also informative in the context of biodiversity conservation and practical management decisions in rapidly changing environments (Gossner et al., 2016).

A common approach to examine spatial β‐diversity is to consider directional turnover in community composition with distance (i.e., distance decay) (Anderson et al., 2011; Nekola & White, 1999). The similarity among the pairs of biological communities typically decreases (“decays”) with increasing spatial or environmental distance (Nekola & White, 1999). This pattern stems mainly from dispersal limitation [related to physical barriers and dispersal ability (Hubbell, 2001)] and species‐specific responses to spatially structured environmental variation [related to environmental filters and evolutionary processes (Cottenie, 2005)] and is well documented in observational (Astorga et al., 2012) and theoretical studies (Morlon et al., 2008) and in meta‐analyses (Soininen et al., 2007). Although the patterns and drivers of taxonomic β‐diversity are relatively well studied in the biogeographical literature, whether the same patterns occur for functional β‐diversity is much less understood (Villéger et al., 2012).

Understanding functional diversity relies on trait‐based approaches, which are built on the idea that the environment selects species based on their ecological requirements and that functional traits capture these requirements better than species identity (McGill et al., 2006). Thus, a trait‐based approach should reflect the functional response of biotic communities to environmental gradients better than an approach based only on the taxonomic identities of species and should predict better how biotic communities respond to environmental changes (Mouillot et al., 2013). Functional diversity has been investigated widely at the α‐diversity level (Buisson et al., 2013; Villéger et al., 2008), but our understanding of functional β‐diversity is much more limited and fragmented (Heino & Tolonen, 2017; Penone et al., 2016; Villéger et al., 2013). Comparing the patterns of functional and taxonomic β‐diversity across different biotic groups, ecosystems and geographical contexts has the potential to contribute greatly to a mechanistic understanding of the drivers behind the spatial variation in ecosystem functionality and shed further light on how environmental change might affect ecological communities.

Several ecological processes can be inferred from the correlation between taxonomic and functional similarity. For example, a coupling of taxonomic and functional distance decay might indicate that species from the regional pool have equal probabilities of reaching all sites, but local communities are assembled based on local habitat constraints on organisms present at each site (Sokol et al., 2011). This generates functional clustering (i.e., trait variability is smaller than expected given the taxonomic composition) at the site level, but overdispersion (i.e., trait variability is larger than expected given the taxonomic composition) at the regional level. This phenomenon has been observed, for example, for specific leaf area of tree communities along an elevational gradient (Swenson et al., 2011). However, high taxonomic β‐diversity does not always mean high functional β‐diversity (Leibold & McPeek, 2006; Mouillot et al., 2013) (Figure 1a). In fact, a functional decay stronger than expected given the taxonomic decay might occur if the species in the two communities are functionally more divergent than expected given the species pool. In contrast, a weaker functional decay than expected given the taxonomic decay might occur if local communities under same habitat constraints are subsamples of multiple regional pools with different species composition. Therefore, the most pressing question is whether functional similarity decays typically faster or slower along environmental gradients than expected given the taxonomic decay, as suggested by some earlier studies (Carvalho et al., 2020; Sokol et al., 2011; Swenson et al., 2011).

FIGURE 1

(a) Taxonomic and functional distance decay. Two scenarios of distance decay of taxonomic and functional similarities along spatial and environmental distance. In scenario 1 (for simplicity, we consider here replacement only), the replacement occurs among species that have different traits (i.e., colours), which leads to both taxonomic and functional distance decay. In scenario 2, the replacement occurs among species that have similar traits, which leads to zero functional distance decay measured by the slope. (b) Master hypothesis: spatial distance decay is stronger for taxonomic similarities than for functional similarities, whereas environmental distance decay is stronger for functional similarities. (c) Specific hypotheses (higher values indicate steeper slopes) across datasets. For latitude, spatial distance decay is flatter in the datasets from higher latitude and, more notably, for taxonomic similarities than for functional similarities. Environmental distance decay is steeper in datasets from higher latitude for functional similarities, whereas it does not vary notably with latitude for taxonomic similarities. For spatial extent, both taxonomic and functional spatial distance decay are flatter in the datasets covering a larger spatial extent, whereas environmental distance decay is steeper in datasets covering a larger extent. For realm, marine ecosystems show flatter spatial and environmental distance decay than terrestrial and freshwater systems. Abbreviations: FRE = freshwater systems; MAR = marine systems; TER = terrestrial systems

Hypotheses

Following the first comprehensive distance decay meta‐analysis (Soininen, McDonald, et al., 2007), our understanding of community turnover along spatial and environmental gradients has increased notably. Here, based on existing ecological literature and theory, and as an initial step towards synthesizing knowledge, we tested four hypotheses concerning the differences between taxonomic and functional distance decay along the spatial and environmental distance. The master hypothesis (H1a) is that the distance decay along spatial gradients is stronger for taxonomic similarity than for functional similarity. This is because spatial factors relate more to taxonomic than functional composition owing to dispersal processes, dispersal history and speciation (Soininen et al., 2016). Such a hypothesis should be valid when functional traits do not comprise dispersal‐related traits. In contrast (H1b), distance decay along environmental gradients is stronger for functional similarity than for taxonomic similarity because functional composition should respond more strongly to environmental variation (Meynard et al., 2011; Soininen et al., 2016) (Figure 1b).

Latitudinal gradients

We also generalize the effects of major geographical and environmental factors in the three hypotheses that are tested across the datasets. For example, latitudinal effect has been recognized as a relevant factor in meta‐analyses (Soininen et al., 2007) and case studies (Qian et al., 2009), and these studies suggest that β‐diversity should decrease with increasing latitude (Figure 1c). This is indicated by the faster latitudinal decline in γ‐diversity than in α‐diversity (Hillebrand, 2004; Soininen, 2010) and the decrease in slopes of the species–area relationships (proxy for turnover) with latitude (Drakare et al., 2005). Moreover, Rapoport's rule (Stevens, 1989) postulates that species range sizes are larger at high latitudes, leading to lower β‐diversity (but see Rohde, 1996). Therefore, we hypothesize (H2a) that the rate of taxonomic distance decay along spatial gradients is generally slower in the datasets that originate from higher latitudes. In contrast, functional distance decay along a spatial gradient might be faster in the datasets from higher latitudes because large‐scale environmental heterogeneity tends to increase towards the poles (Soininen & Hillebrand, 2007; Soininen, McDonald, et al., 2007; Terborgh, 1973). Thus, environmental filtering becomes stronger with increasing latitude (Jarzyna et al., 2021; Lamanna et al., 2014), leading to functionally clustered communities locally that become increasingly overdispersed along the regional environmental gradient. This would result in a faster rate of functional distance decay along environmental gradients at higher latitudes (H2b). An alternative hypothesis is that extreme climatic conditions at high latitudes decrease functional diversity because abiotic filtering limits the number of possible ecological strategies found in a biotic community (Cornwell & Ackerly, 2009), resulting in a relatively slow rate of functional distance decay.

Spatial extent

Distance decay is also likely to be affected by the spatial extent of a given study (Nekola & McGill, 2014; Steinbauer et al., 2012). It has been shown that distance decay has a power‐law shape at spatial extents that do not exceed regional species pools and an exponential shape when extent encompasses multiple species pools (Nekola & McGill, 2014). This suggests that the slope of the relationship becomes flatter with increasing spatial extent (Soininen, McDonald, et al., 2007), mainly because regional species diversity is limited with a certain upper boundary (Triantis et al., 2011). Furthermore, environmental heterogeneity affects the diversity of species (Pianka, 1966) and functional traits at a regional level (Questad & Foster, 2008), but such effects are likely to be scale dependent (Gazol et al., 2013; Laanisto et al., 2012). To summarize, we hypothesize (H3a) that the rate of distance decay along spatial gradients is generally slower in the datasets covering larger spatial extent. In contrast (H3b), the rate of distance decay along environmental gradients is generally faster when the spatial extent is larger, especially for functional similarities.

Realms

We also expect the patterns of distance decay to vary among realms. In general, marine ecosystems are environmentally more homogeneous than terrestrial or freshwater ecosystems, at least in the open ocean (Clarke, 1992), and typically show weaker dispersal barriers than terrestrial or freshwater ecosystems (Bierne et al., 2003; Cornell & Harrison, 2014). Therefore, we hypothesize (H4) that the datasets from marine ecosystems generally have slower rates of taxonomic and functional distance decay than the other ecosystems.

Here, we tested these hypotheses using datasets that cover a wide range of biotic groups, from unicellular diatoms to vascular plants, fungi, invertebrates, fish, birds, amphibians and mammals, and that originate from marine, terrestrial and freshwater ecosystems spanning broad latitudinal gradients (Figure 2). To account for major biological differences in biotic groups, we also investigated whether distance decay varied among different‐sized taxa or among taxa with different dispersal modes (Bie et al., 2012; Jenkins et al., 2007). By using such a comprehensive, multi‐realm and multi‐taxon dataset, we explore patterns at a more general level compared with case studies that have examined both taxonomic and functional β‐diversity but considered only a single or few biotic groups.

FIGURE 2

Study design highlighting (a) a map of the study sites coloured according to the realms (FRE = freshwater; MAR = marine; TER = terrestrial); (b) the number of datasets for major biotic groups; and (c) the distribution of the datasets with respect to spatial extent, number of study sites, functional γ‐diversity (log10 hypervolume SD3), taxonomic γ‐diversity (number of species), number of environmental variables and latitude

MATERIAL AND METHODS

Data collection

We gathered our data by directly contacting data owners or using the existing data sources, such as GrassPlot (Biurrun et al., 2021; Dengler et al., 2018), sPlot (Bruelheide et al., 2019) and CESTES (Jeliazkov et al., 2020). We included datasets that provided raw data of species abundances, functional traits, environmental variables and spatial coordinates of the study sites. A few datasets (n = 6) provided only species occurrence rather than abundance information (Supporting Information Appendix S1). The traits included in the datasets were chosen by data owners from a suite of traits that should respond well to environmental variation. For plant datasets compiled from the sPlot database, trait information was commonly derived from the TRY database (Kattge et al., 2011). Regarding the CESTES database, we compiled 48 datasets, specifically from: fish communities (Brind’Amour et al., 2011; Chong‐Seng et al., 2012; Cleary et al., 2016; Villéger et al., 2012), terrestrial vascular plants (Bagaria et al., 2012; Chmura et al., 2016; Eallonardo et al., 2013; Frenette‐Dussault et al., 2011; Fried et al., 2012; Jamil et al., 2013; van Klink et al., 2017; Meffert & Dziock, 2013; Pakeman, 2011; Raevel et al., 2012; Ribera et al., 2001; Robroek et al., 2017), aquatic macroinvertebrates (Díaz et al., 2007; Gallardo et al., 2009; Jeliazkov et al., 2014), terrestrial arthropods (Barbaro & van Halder, 2009; Barbosa et al., 2015; Bartonova et al., 2016; Dziock et al., 2011; Gibb et al., 2015; Gonҫalves‐Souza et al., 2014; van Klink et al., 2017; Krasnov et al., 2015; Lowe et al., 2018), birds (Barbaro et al., 2012, 2017; Barbaro & van Halder, 2009; Charbonnier et al., 2016; Cleary & Renema, 2007; Meffert & Dziock, 2013), bats (Charbonnier et al., 2016; Farneda et al., 2015), bryophytes (Robroek et al., 2017), butterflies (Bartonova et al., 2016; Robinson et al., 2014), corals (Rachello‐Dolmen & Cleary, 2007) and foraminifera (Cleary & Renema, 2007). We included only datasets with ≥10 sites, two environmental variables and three traits or trait categories. In some cases, more than one dataset, representing different taxonomic groups with different responses to the environment and dispersal abilities (e.g., stream macroinvertebrates and diatoms) was collected in the same study area. In total, 148 datasets representing 17 major biotic groups from terrestrial (n = 87), freshwater (n = 41) and marine (n = 21) environments were assembled, amounting to >17,000 study sites around the globe (Figure 2). Of the 148 datasets, 118 were published in peer‐reviewed journals (Appendix S1). Taxa were identified mostly to species or morphospecies; in a few cases, we used data at the genus level if existing taxonomic knowledge did not allow individual species to be distinguished. Within biotic groups, traits were generally the same or at least covered similar functional roles (Appendix S1). Finally, each dataset included (1) a sites‐by‐species abundance matrix; (2) a species‐by‐traits table; (3) a sites‐by‐spatial coordinates table; and (4) a sites‐by‐environmental variables table (Figure 3a). Detailed information about collected datasets can be found in Appendix S1 and information about data curation in Appendix S2. All the data curation and further analyses were performed in the software R v.4.0.2 using the appropriate R packages. We will consistently refer to the functions used and their respective packages from here on.

FIGURE 3

The analytical framework described in a stepwise manner: (a–c) hierarchical description of the methods performed at dataset level, including the estimation of similarities and distance in addition to the distance decay models of each dataset; and (d) description of the tests performed after the compilation of the metrics from all datasets. (a) The four objects used in the analyses: a species‐by‐traits table, a sites‐by‐species matrix, a sites‐by‐coordinates table and a sites‐by‐environment table. (b) The calculation of taxonomic and functional similarities and of spatial and environmental distance. In the first example, only species identities are considered, and because sites k and k do not share any species, community similarity (blue) equals zero. In the second example, the functional traits of species are considered, and community similarity (orange) is higher than zero. The third example shows how spatial distance was calculated as the geographical distance between pairs of sites using spatial coordinates. The fourth example illustrates how sites far from each other may show similar environmental conditions and therefore small environmental distance. Environmental distance was calculated as the Euclidean distance between pairs of sites considering the standardized environmental variables. (c) Illustration of the metrics extracted to study the distance decay across datasets. The strength of distance decay was measured from Mantel tests using Spearman correlations (Mantel r), and the rate of decay was measured as the slopes of generalized linear models following a quasibinomial family with a log link. The models were built separately for each response variable (taxonomic or functional similarity) and explanatory variables (spatial or environmental distance), totalling four Mantel r values and four slopes. Also, the data of marine fish from the Mediterranean Sea are shown as an example in which the distance decay of similarity along environmental distance is stronger (higher Mantel r) for functional similarity than for taxonomic similarity, irrespective of the rate of decay (slope). (d) Description of the analyses used to test the hypotheses and which metrics were considered for each analysis. The strength (Mantel r) of decay was used to test hypothesis H1, and the rate of decay (slope) was used to hypotheses H2–H4

Community similarity

When estimating community similarities, we used both occurrence and abundance data, because occurrences are informative about the drivers and patterns of communities along geographical gradients, whereas abundances inform patterns along environmental gradients well (Declerck et al., 2011). The similarity in species composition between two communities (hereafter, taxonomic similarity) was estimated with the function beta in the package BAT v.2.7.0 (Cardoso et al., 2020). For the similarity in trait composition (hereafter, functional similarity), we first represented the ecological niche spaces of species as n‐dimensional hypervolumes (Hutchinson, 1957) sensu Blonder (2018). We used probabilistic hypervolumes based on the Gaussian kernel density estimator instead of other widely used methods (e.g., convex hulls; Villéger et al., 2013) because it allows the inclusion of abundance information in the estimation of functional space and it gives different probabilities to more or less populated areas of the functional space, being less influenced by outliers (Mammola & Cardoso, 2020; Mammola et al., 2021). Community similarity (S) ranges between zero and one and is commonly calculated for pairs of communities as the sum of the unique features of each community divided by the sum of the shared features between communities and the unique features of each community. We estimated the taxonomic and functional similarities between communities j and k using the Sørensen similarity index as: where a is the sum of the shared features between the communities j and k, b is the sum of the unique features of community j, and c is the unique features of community k. When estimating taxonomic similarities, each species is a feature that may or may not be shared by two communities (Figure 3b). For functional similarities, each feature comprises the area of the hypervolume that is either shared between two communities or unique to each community (Figure 3b).

To construct the hypervolume of each community, we first calculated the pairwise difference between species using the Gower similarity index. We used Gower similarity as modified by Bello et al. (2021) because it gives balanced weights for categorical, continuous, dummy and fuzzy‐coded variables. The modified gower distance was calculated with the function gawdis within the package gawdis v.0.1.3 (Bello et al., 2021). Then, we ran a principal coordinates analysis (PCoA) on this distance matrix and summarized the trait data in three PCoA axes, following a trade‐off between computation time and information quality. Using the three PCoA axes, we constructed the hypervolume using the function kernel.build from BAT. Finally, we estimated the amount of overlap between two hypervolumes, in addition to the unique area of each community, using the function kernel.beta of BAT that builds on hypervolume_set of hypervolume v.2.0.1 (Blonder & Harris, 2019). For functional similarities, if two communities do not share any species, taxonomic similarity would be, by definition, lower than functional similarity if any continuous trait (e.g., body size; Figure 3b) is included. Details of the calculation of similarities using the Sørensen‐based indices for occurrence and abundance data can be found in the (Appendix S2). The main results are given for occurrence data in the main text, whereas abundance‐based results can be found in the (Appendix S3).

Spatial and environmental distance

For each dataset, spatial distance was calculated as the geographical distance (in kilometres) between the pairs of sites using the function earth.dist of the package fossil v.0.4.0 (Vavrek, 2011) (Figure 3b). Environmental distance was calculated as the Euclidean distance between all the pairs of sites, considering the standardized environmental variables within each dataset (Figure 3b), with the function vegdist of the package vegan v.2.5–6 (Oksanen et al., 2019). From the original environmental variables in each dataset, we only kept the continuous variables with <5% missing values (for data curation, see Appendix S2). Given that the datasets contained different numbers and types of environmental variables, the values of environmental distance were context dependent and not very informative for comparison across datasets. We therefore assumed that the environmental gradient increased with spatial extent and rescaled the actual environmental distance to range between zero and one in each dataset by dividing actual values by the maximum environmental distance of the dataset.

Distance decay of similarity

To comply with the assumption of nonlinearity in the distance decays, the strength of the distance decay relationship was assessed by ranked Mantel tests (using a Spearman correlation, i.e., Mantel r). The rate of the decay was modelled as the slope of generalized linear models (GLMs) following a quasi‐binomial family with log‐link (Millar et al., 2011), representing a negative exponential curve between the community similarity and distance (Figure 3). One of the main assumptions of the distance decay is that the slope of the relationship should be negative (Nekola & McGill, 2014), and positive slopes suggest either periodicity in the environmental gradient or a mismatch between the communities and the measured environmental variables (Nekola & White, 1999). Therefore, whenever datasets showed positive distance decay slopes, these were transformed to zero values. In total, five datasets had positive slopes for taxonomic similarities, whereas 11 datasets had positive slopes for functional similarities.

Statistical analysis

We tested our master hypothesis using two different approaches. Firstly, we investigated whether taxonomic or functional distance decay is stronger along spatial and environmental distance (H1) by performing Student’s paired t tests to compare Mantel r drawn from taxonomic and functional similarity for each dataset, considering both spatial and environmental distance (Figure 3d). Secondly, we generated a null distribution of functional similarity values by randomizing the names of the species across the trait table 999 times. At each iteration, the functional similarities were calculated and regressed against spatial and environmental distance. The slopes of these relationships were used to obtain a null distribution of slopes under the assumption that the rate of functional decay is random given the rate of taxonomic decay. Deviations from null distribution were tested using standardized effect sizes (SES; Gotelli & Graves, 1996); SES values >1.96 indicate that functional similarity decays faster than expected given taxonomic decay, whereas SES values <−1.96 indicate that functional similarity decays slower than expected given taxonomic decay (Swenson et al., 2011).

We also investigated the ecological and geographical factors driving the rate of the distance decay across datasets. Each dataset was characterized with respect to: (1) latitude, recorded as the absolute mean value of all the sites of the dataset; (2) spatial extent, expressed as the largest pairwise distance (in kilometres) between study sites; (3) realm, classified into freshwater, marine and terrestrial environments; (4) body size information drawn from literature (Hillebrand, 2004; Peters, 1983), estimated as the mean log10‐transformed fresh weight (in grams) of the biotic group included in the dataset; (5) dispersal mode, classified as active and passive modes and organisms dispersed by seeds; (6) taxonomic γ‐diversity, expressed as the total number of species in the dataset; (7) functional γ‐diversity, measured as the total volume of the union of the n‐dimensional hypervolumes estimated within the dataset (Mammola & Cardoso, 2020); (8) total number of study sites in the dataset; and (9) the number of environmental variables in the dataset. For body sizes, we note that although the size range within the biotic group can be large (up to five orders of magnitude), it is small compared with the overall variation obtained across organism groups (12 orders of magnitude). For more details on body size approximations, see the papers by Hillebrand (2004) and Drakare et al. (2005). The taxonomic γ‐diversity was included to study whether there is a typical positive relationship between γ‐diversity (taxonomic and functional) and β‐diversity (Kraft et al., 2011; Lamanna et al., 2014). Functional γ‐diversity was estimated using the function hypervolume_set of the package hypervolume v.2.0.12 (Blonder & Harris, 2019). Hypervolumes are expressed in units of SDs to the power of the number of trait dimensions used (i.e., three). The number of study sites and the number of environmental variables for each dataset were included to explore their potential effect on distance decay.

Finally, we used boosted regression trees (BRTs) to test the effects of latitude (H2), spatial extent (H3) and realm (H4) on the rate of taxonomic and functional distance decay along spatial and environmental distance across the datasets. The BRT parameters were selected to amplify the deviance explained by the model. We tested interaction depths of 0.1, 0.01 and 0.001, and learning rates between two and five. The best models were the ones with a learning rate of five and interaction depth of 0.001. Given that the datasets in this study have not always followed the same sampling methodology and show different functional traits and environmental variables, we fitted the BRT models following a Laplace distribution of the errors to reduce the absolute error loss from the variation among datasets. We included taxonomic and functional γ‐diversity, number of sites and number of environmental variables in the dataset as predictor variables for controlling the heterogeneity across datasets within the model (Figure 3d). The BRT models were fitted using the function gbm.step of the package dismo v.1.1–4 (Hijmans et al., 2017). Given that some predictors could show cross‐information (e.g., marine ecosystems have smaller organisms than terrestrial systems in our datasets), we tested whether there were significant interactions between predictors in BRTs using the functions gbm.interaction. For understanding how gbm.interaction works and for more in‐depth details about BRTs, we refer to the (Appendix S2 and references therein). Additionally, we performed sensitivity analysis to ensure that the patterns detected were not an artefact of sample size. In the sensitivity analysis, models were refitted using subsampled data with 90, 70 and 50% of the full observations.

We partitioned the similarities into replacement and richness difference components following the methodology described in the (Appendix S2). Replacement gives the variation resulting from the substitution of species (taxonomic replacement) or functional traits (functional replacement), whereas richness differences account for the variation resulting from net differences induced by the loss/gain of species or traits (Carvalho et al., 2012). We showed the full results of the distance decay using occurrence‐based total similarities (Equation 1), but we also used abundance‐based similarities and showed their main findings in the main text, with further details in the (Appendix S3). To keep the narrative concise, in the main text we show only the results of the partitioned components using occurrence data. All figures were built using the packages from the tidyverse suite (Wickham et al., 2019).

RESULTS

Strength of the distance decay

The distance decays showed a wide range of shapes, from very steep decays to almost flat relationships (Figure 4). The average Mantel r using occurrence data along spatial distance for taxonomic similarities was .254 (SD ±0.197) and .115 (SD ±0.143) for functional similarities. Spatial distance decays of taxonomic similarities were significantly stronger than the distance decays of functional similarities when considering both occurrence (Figure 4a; t = 13.124, p < .001, d.f. = 146) and abundance data (Appendix S3; Figure S3.3), supporting H1a that spatial distance decay is stronger for taxonomic than functional similarities (Figure 4a). In 31 datasets, the spatial distance decay of functional similarities was stronger than taxonomic similarities.

FIGURE 4

The distance decay along (a) spatial distance and (b) environmental distance. The light blue lines show the distance decay of taxonomic similarity, and the orange lines show the distance decay of functional similarity. The first and second columns show the rate (slope) of the taxonomic and functional distance decay, respectively; the third column shows the strength (Mantel r) of the distance decay of taxonomic and functional similarities; and the fourth column shows the standardized effect sizes of the slopes of each dataset

However, our results did not support H1b, because the distance decay for taxonomic similarities (mean Mantel r = .272, SD ±0.176) was also, on average, stronger than for functional similarities (mean Mantel r = .150, SD ±0.144) along environmental distance, considering both occurrence (Figure 4b; t = 10.342, p < .001, d.f. = 146) and abundance data (Appendix S3; Figure S3.3). Note, however, that 32 of 148 datasets had stronger distance decay of functional similarities than taxonomic similarities along environmental gradients.

Rate of the distance decay

The mean slope of the spatial distance decay was −0.009 (SD ±0.028) for taxonomic similarities and −0.004 (SD ±0.020) for functional similarities (Figure 4a). Null models showed that only 13 datasets had a functional distance decay significantly different from expected given taxonomic decay, with three datasets having faster slopes and 10 having slower slopes (Figure 4a). For environmental distance, the mean slope of the distance decay was −1.077 (SD ±1.066) for taxonomic similarities and −0.355 (SD ±0.031) for functional similarities (Figure 4b). Null models along environmental distance showed that only 11 datasets had a functional distance decay significantly different from expected given the taxonomic decay, from which only one was slower (Figure 4b).

Regarding the biotic groups, terrestrial plants had the steepest slopes along spatial distance for both taxonomic and functional similarities (Figure 5). Along environmental distance, corals had the steepest slopes along spatial distance, whereas foraminifera had the steepest slopes along environmental distance (Figure 5). Similar patterns were found for abundance‐based similarities, except for the biotic groups, where aquatic plants had the steepest slopes along spatial distance (Appendix S3; Figure S3.4).

FIGURE 5

The average rate of decay (slopes) of biotic groups using occurrence data along spatial and environmental distance. The vertical dashed lines highlight the zero rate (absence of decay), and the horizontal lines indicate the standard deviation of the mean. The blue circles show the rate of decay of taxonomic similarities, and the orange circles show the rate of decay of functional similarities. Large error bars are attributable to low sample size (i.e., a low number of datasets for a given taxon)

Across datasets, BRT explained 35.4% of the deviance of the slopes of the spatial distance decay for taxonomic similarities and 38.3% for functional similarities using occurrence data. For the distance decay along environmental distance, BRT explained only 14.4% of the deviance of the slopes of the decay of taxonomic similarities and 6.6% for functional similarities. Spatial extent, latitude and γ‐diversity contributed the most to the variation in slopes along either spatial or environmental distance using both occurrence‐ and abundance‐based similarities (Figures 6a and 7a; Appendix S3; Figures S3.5 and S3.6). We found interactions for the BRT only on taxonomic similarities along environmental gradients. These were between: (1) taxonomic γ‐diversity and body size (interaction size = 5.09); and (2) spatial extent and latitude (interaction size = 3.79). These are probably related to the higher number of species included in the datasets containing small biotic groups and the fact that larger datasets had a mean latitude near the equator. There was no evident effect of sample size on the main patterns according to the sensitivity analysis (Appendix S2; Figures S2.1 and S2.2).

FIGURE 6

Relative effects (expressed as percentages) of geographical factors on the rate of decay along (a) spatial distance decay and (b) environmental distance decay of the total component of taxonomic (TAX, light blue) and functional (FUN, orange) similarities using occurrence data across datasets. Partial dependence plots show the effects of a predictor variable on the response variable after accounting for the average effects of all other variables in the model. Positive values indicate an increase in the rate of decay (steeper slopes) compared with the mean rate, whereas negative values indicate a decrease in the rate of decay (flatter slopes) compared with the mean rate. Semi‐transparent lines represent the actual predicted effects; continuous lines represent LOESS fits to predicted values from boosted regression trees (BRTs). We show here only the variables related to the specific hypotheses [i.e., latitude, spatial extent and realms (FRE = freshwater; MAR = marine; TER = terrestrial)]

FIGURE 7

Relative effects (expressed as a percentage) of organismal variables and dataset features on the rate of decay along (a) spatial distance and (b) environmental distance, considering the total component of taxonomic (light blue lines) and functional (orange lines) similarities using occurrence data across datasets. Partial dependence plots show the effects of a predictor variable on the response variable after accounting for the average effects of all other variables in the model. Positive values indicate an increase in the rate of decay (steeper slopes) compared with the mean rate, whereas negative values indicate a decrease in the rate of decay (flatter slopes) compared with the mean rate. Semi‐transparent lines represent the actual predicted effects; continuous lines represent LOESS fits to predicted values from boosted regression trees (BRTs). We show here the organismal variables and the variables related to the dataset features

Latitudinal patterns

The slopes of spatial distance decay of both taxonomic and functional similarities were steepest in datasets centred at c. 35–45°, providing only partial support for H2a that distance decay was flatter at high latitudes (Figure 6a). Note that spatial distance decay decreased sharply towards the poles. The slopes of environmental distance decay were flattest in the datasets at c. 50° (Figure 6b). However, note that functional distance decay increased towards the poles, providing partial support to hypothesis H2b. Similar patterns were found for abundance‐based similarities (Appendix S3; Figure S3.5).

Spatial extent

The distance decay of taxonomic and functional similarities was flatter in the datasets that covered a larger spatial extent for both occurrence (Figure 6a) and abundance data (Appendix S3; Figure S3.5a), supporting hypothesis H3a that distance decay becomes flatter with increasing spatial extent. For environmental distance, distance decay was steeper in the datasets that covered larger spatial extents only for taxonomic similarities, but functional distance decay did not vary with extent. Thus, our results agreed, in part, with H3b that distance decay would become steeper with larger spatial extent.

Realms

Marine ecosystems typically had flatter slopes in comparison to freshwater or terrestrial ecosystems, thus agreeing with H4 (Figure 6). However, the importance of the realms in BRTs was low overall. A similar pattern emerged for abundance‐based similarities (Appendix S3; Figure S3.5).

Organismal variables and dataset features

Organisms relying on seed dispersal had steeper slopes along spatial and environmental distance than other dispersal types, but the overall importance of dispersal mode was low (Figure 7b). The slopes of both spatial and environmental distance decays were steeper for larger‐bodied organisms in taxonomic and functional similarities (Figure 7a,b). Taxonomic γ‐diversity had a U‐shaped relationship with slopes for distance decay along spatial and environmental distance (Figure 7b). Slopes of distance decay had an overall flattening trend towards higher functional γ‐diversity for both spatial and environmental distance (Figure 7a,b). Generally, taxonomic slopes were steeper in the datasets where the number of study sites was higher (Figure 7a), whereas the opposite was true for functional slopes. The slopes were flatter when datasets contained only a few environmental variables (Figure 7b).

Replacement and richness differences

The slopes of taxonomic and functional replacement along spatial distance decreased rapidly in the datasets above 35° (Appendix S4; Figure S4.7a). Along environmental distance, the taxonomic replacement increased towards higher latitudes, whereas the functional replacement had a U‐shaped pattern, with a decrease from low to mid‐latitudes (c. 45°) and a sharp increase towards the poles (Appendix S4; Figure S4.7b). For the richness differences component, the slopes of taxonomic similarities were steepest in the datasets at c. 45° for the spatial distance decay, whereas the slopes of functional similarities became notably steeper with latitude (Appendix S4; Figure S4.8a). For environmental distance, slopes became flatter from low to high latitudes up to c. 50° for taxonomic similarities, whereas for functional similarities, slopes did not vary along latitude (Appendix S4; Figure S4.8b). Both replacement and richness differences showed flatter spatial slopes with increasing spatial extent (Appendix S4; Figure S4.7 and S4.8). In contrast, environmental slopes became steeper with spatial extent for taxonomic replacement and flattened for functional replacement (Appendix S4; Figure S4.7b), whereas the functional slopes showed an opposite pattern (Appendix S4; Figure S4.8b). Furthermore, marine ecosystems showed the flattest slopes for taxonomic replacement along environmental gradients (Appendix S4; Figure S4.7b), whereas terrestrial ecosystems had the flattest slopes for richness differences (Appendix S4; Figure S4.8b). Details about the organismal variables and dataset features can be found in the (Appendix S4; Figures S4.9 and S4.10).

DISCUSSION

Community ecology, macroecology and biogeography have lacked a comprehensive evaluation of functional β‐diversity across different taxa and ecosystems globally. Earlier studies suggest that functional β‐diversity better reflects environmental variability compared with taxonomic β‐diversity and that focusing on functional β‐diversity might help us, for example, to gain an understanding of how humans impact ecosystems by modifying the local environment (Meynard et al., 2011; Sokol et al., 2011; Spasojevic et al., 2014; Weinstein et al., 2014). This is because functional traits should reflect best the ecological requirements of species. Using a comparative analysis across biotic groups, ecosystem types and realms, we show here that taxonomic distance decay is generally stronger along spatial gradients than functional distance decay and that the decay of functional similarities along environmental gradients is typically not stronger than the decay of taxonomic similarities, unlike previous suggestions.

The strength of the distance decay of taxonomic and functional similarities

The stronger signal of taxonomic than functional distance decay along space provides empirical evidence that taxonomic distance decay is a robust approach for ecological and biogeographical studies, supporting H1a. Compositional differences effectively summarize dispersal‐related factors in addition to species responses to climatic and other spatially structured environmental variables. However, spatial distance decay of functional similarities might not reflect the geographical differences in biotic communities well. This is likely to stem from the different roles played by deterministic and stochastic factors shaping community composition, because it has been shown that dispersal limitation or species pool effects should be more important for taxonomic than for functional composition (Soininen et al., 2016). Some morphological or morphometric traits are informative when exploring geographical patterns in functional composition (Soininen et al., 2016); for example, seed mass and wood density explained the variability of tree communities along broad spatial gradients better than species identity alone (Siefert et al., 2013). The type of dispersal is also an important trait to include when assessing community‐level patterns along spatial gradients (Bie et al., 2012). Unsurprisingly, we found that the datasets with larger spatial extent and species pool size were more likely to have a stronger distance decay of taxonomic than functional similarities (Appendix S5; Table S5.1). We argue that when the study extent is large it might cross several species pools, and it is more likely that species occurrences are also affected by historical and dispersal‐related factors and not only by environmental preferences (Soininen et al., 2016). Furthermore, we found that passive dispersers and datasets with higher functional γ‐diversity were less likely to have stronger decay of taxonomic similarities than functional similarities. In our datasets, passive dispersers were microorganisms and, therefore, efficient dispersers with a good ability to reaching sites with suitable environmental conditions (Bie et al., 2012; Fontaneto, 2019). Thus, for these taxa, functional distance decay should be more informative than taxonomic distance decay.

Regarding environmental gradients, functional distance decay was also weaker than taxonomic distance decay. Although counterintuitive, this relationship has been found previously (Heino & Tolonen, 2017; da Silva & Hernández, 2015; Teittinen & Virta, 2021). For example, Teittinen and Virta (2021) observed stronger distance decay of taxonomic than functional similarities along environmental gradients, which they attributed to the greater number of species than functional traits in their data. Also, Heino and Tolonen (2017) found similar results for macroinvertebrate communities of boreal lakes and related it to the trait resolution, which could probably be improved by the addition of several other physiological traits relevant for the organisms in question. Here, additional analysis showed that increasing spatial extent, species pool and the number of environmental variables significantly increased the probability of a dataset having stronger decay of taxonomic similarities compared with functional similarities (Appendix S5; Table S5.2). In fact, the ratio between taxonomic and functional decay depends on whether the species replaced from one community to another are a random subsample of functionally redundant species from the regional pool (Swenson et al., 2011). Also, species pool size and functional redundancy typically exhibit a positive correlation (Cannicci et al., 2021; Mouillot et al., 2014), which, in turn, should increase the functional similarities between sites (Jarzyna & Jetz, 2018). We suggest that within a large species pool, the functional redundancy of species increases, given the limited set of trait combinations and/or available niches. Therefore, smaller species pools are more likely to have functionally unique species and lower functional similarities than larger pools. In the case of large pools, we found that taxonomic decay was often stronger than functional decay. Furthermore, because species pool size increases with study extent (Drakare et al., 2005; Palmer & White, 1994; Triantis et al., 2011), the datasets with larger extents had slower functional distance decay even along environmental gradients, and taxonomic composition turned out to be the best descriptor of distance decay patterns. Another possible reasoning is that filtering on a given trait might filter other traits concomitantly, and if the focal trait is not included in the analyses, a mismatch between functional composition and the environment is expected. On the contrary, a dataset might comprise traits not affected by the environment, which tends to increase the functional similarity among sites. Therefore, because functional diversity patterns depend strongly on the traits measured (Zhu et al., 2017), the choice of traits should be planned carefully.

The effect of latitude on the rate of distance decay

In addition to our master hypothesis, we investigated whether the rate of distance decay showed consistent variation across realms, along geographical gradients and among major taxonomic groups. We did not find slower rates of decay along spatial distances in the datasets at higher latitudes, but we found a unimodal relationship with the highest decay rate at c. 30°. Similar results have been found earlier in terrestrial vertebrates when considering only the turnover component of β‐diversity (Castro‐Insua et al., 2016) and for the total β‐diversity of marine phytoplankton (Martin et al., 2021). It is noteworthy that our latitudinal patterns were related mainly to the replacement component for both taxonomic and functional decay (Appendix S4; Figure S4.7). Regarding environmental gradients, we found opposing patterns compared with spatial gradients, with the flattest rates of decay in the datasets near 50° and a notable increase from 60° towards the poles. A hump‐shaped relationship between functional diversity and latitude has also been found previously for aquatic macroinvertebrates (Múrria et al., 2020), also with the minimum at c. 50°. Múrria et al. (2020) were studying patterns in functional dispersion, whereas we found here that the breakpoint was related mainly to the differences in richness for taxonomic similarities and replacement for functional similarities.

Traditionally, latitudinal patterns of biodiversity have been explained by Rapoport's rule, positing that there is an increase in species range size towards high latitudes (Stevens, 1989), hence lower taxonomic replacement. However, the breakpoints found in our data suggest that some additional factors might have generated the patterns. For example, landscape fragmentation might increase β‐diversity (Jamoneau et al., 2012), especially at mid‐latitudes that showed the highest levels of human impact in this study (Halpern et al., 2008; Venter et al., 2016). Also, it has been suggested previously that the distance decay along spatial distances is stronger at mid‐latitudes than at the poles because northern communities result from postglacial recolonization processes, flattening distance decay relationships (Gómez‐Rodríguez & Baselga, 2018). Although inferring processes from observational data is difficult (Cadotte & Tucker, 2017), we would like to speculate on some possible mechanisms generating our breakpoint patterns. Strong seasonality, resource scarcity and climatic stress should select only the highly specialized taxa and modify the functional space towards the poles (Lamanna et al., 2014). Therefore, it is plausible that the climatic stress leads to an increase in richness differences in communities towards the poles, as observed in vertebrates elsewhere (Castro‐Insua et al., 2016). Moreover, as environmental heterogeneity increases towards the poles, and functional clustering is expected to be stronger at higher latitudes (Jarzyna et al., 2021; but see Kruk et al., 2017), we suggest that strong environmental filtering in datasets at higher latitudes (above 50°) selects for the species with different trait combinations between sites, thereby increasing the rate of functional decay. The latitudinal decrease in the rate of abundance‐based functional distance decay (Appendix S3; Figure S3.3) is further evidence of an optimal utilization of the functional space, as has been observed earlier exclusively for marine organisms (Edie et al., 2018). However, these potential explanations should be tested further.

The effect of spatial extent on the rate of distance decay

The rate of spatial distance decay was slower in the datasets covering a larger spatial extent, suggesting that regional species pools are limited and that new species are not found constantly at the same frequency when extent is larger. Lower decay rates in larger study areas could also result from repeated patterns in environmental variation; that is, environmental patchiness or natural periodicity in the environment (Nekola & White, 1999). In agreement with our hypothesis, we also found that the rate of taxonomic decay along environmental distance was higher in the datasets covering a larger spatial extent. These findings indicate that spatial distance decay is more affected by species pool effects and dispersal processes than environmental distance decay, possibly because the latter reflects more strongly the level of local deterministic environmental filtering processes. Similar evidence has accumulated from case studies conducted in various ecosystems (Meynard et al., 2011; Sokol et al., 2011; Weinstein et al., 2014; Zagmajster et al., 2014). The finding that the rate of distance decay along environmental distance was higher in the datasets covering larger extents indicates the stronger environmental filtering for larger study areas. We also note that, in our BRT models, the extent, latitude and γ‐diversity had by far the largest relative importance, suggesting that their interplay shapes distance decay to a great extent.

The effect of realm on the distance decay

We found evidence for a lower rate of distance decay in marine versus terrestrial or freshwater ecosystems. Overall, this finding agrees with an earlier meta‐review on β‐diversity (Soininen, McDonald, et al., 2007), suggesting that large‐scale diversity patterns are generally weaker in marine ecosystems (Bierne et al., 2003). Given that connectivity, energy flows, dispersal modes, body size structure and trophic dynamics differ substantially between dry and wet ecosystems (Shurin et al., 2006), it is vital to investigate possible differences in turnover among the realms more closely.

Organism size did seem to affect taxonomic or functional distance decay along spatial and environmental gradients, because the slopes typically increased with organism body size. This might be because β‐diversity should be low among the small microbial taxa with efficient passive dispersal (Soininen, McDonald, et al., 2007). The rationale behind this idea is that efficient dispersal homogenizes communities among sites, resulting in lower β‐diversity (Mouquet & Loreau, 2003). Body size is also a key driver of the biological complexity of organisms (Heim et al., 2017), and it might be that smaller organisms show a much more limited set of trait combinations than macroorganisms, leading to a lower functional redundancy among larger species. Furthermore, our knowledge about the taxonomy and functional traits of organisms is typically size dependent. For example, the identification of larger species is much easier than that of microorganisms, which also applies to the identification and measurement of soft functional traits (Hodgson et al., 1999; Martínez et al., 2021). Therefore, the values of functional β‐diversity of small organisms might typically be underestimated.

Patterns in environmental distance decay were relatively congruent with spatial distance decay regarding dispersal mode, suggesting that taxa dispersing passively do not seem to track environmental gradients more efficiently compared with less dispersive taxa. It might also be that small‐sized taxa were filtered along some unmeasured spatially structured environmental gradients, and the pattern was thus detected as spatial turnover even if caused by some underlying unmeasured environmental factors. Forthcoming studies would benefit greatly from disentangling the signal of unmeasured environmental variables from true dispersal limitation (Stegen et al., 2013).

Study design

There are also some possibly influential aspects in our study design that should be discussed. Although the study is global in its extent, the availability of datasets was not evenly distributed geographically. This is a well‐known problem in biodiversity research (Titley et al., 2017) that calls for complementary studies to verify that these trends hold true in poorly sampled regions. Also, we relied on the suite of traits and environmental variables included in the original datasets, hence the collection of traits and environmental variables used differed somewhat among datasets even for the same focal taxonomic groups. Although traits covered mostly the same functional roles of the species, the variation in traits and environmental variables across datasets increases the uncertainty on how environmental variables filter the functional structure of communities in different contexts and how strong the community–environment relationships might be. An alignment of key traits and environmental variables is therefore desirable but requires a suite of sister studies following the same protocol, which is, unfortunately, not yet available. Moreover, the fact that some of the biotic groups (e.g., bryophytes, corals, foraminifera) were underrepresented in our analysis, with only a few datasets included (Figure 2), or the total lack of some taxa (e.g., bacteria, and aquatic and terrestrial mammals), makes it more difficult to generalize distance decay across certain taxa.

Concluding remarks

We believe our analysis is an important step towards a more comprehensive understanding of patterns and drivers of functional β‐diversity, particularly in comparison to the patterns and drivers of taxonomic β‐diversity that have so far attracted much more research interest. Here, we found that functional distance decay is scale dependent and a product of large‐scale geographical factors (latitude) and taxonomic and functional γ‐diversity but is also driven by the biology of organisms to some degree. In general, taxonomic distance decay is a useful tool for many aspects of biogeographical research because it reflects dispersal‐related factors in addition to species responses to climatic and other spatially structured environmental variables. However, functional distance decay might be a cost‐effective option for investigating how species respond to the environment, especially for microorganisms (e.g., microalgae), which are typically difficult to identify to the species level. Overall, the present findings and data shed light into the congruence between the functional and taxonomic diversity patterns globally and provide useful new information to the field of functional biogeography.

CONFLICT OF INTEREST

The authors declare no conflicts of interest.

AUTHOR CONTRIBUTIONS

Caio Graco‐Roza and Janne Soininen contributed equally to the original idea, data analysis and the writing of the first draft. Jani Heino and Otso Ovaskainen advised on the main idea, analysis and commented on the first draft. Francesco Maria Sabatini coordinated the data compilation from the “sPlot” database and commented on the first draft. Martin Gossner coordinated the compilation of the data from the “Biodiversity Exploratories” project and commented on the first draft. All other authors shown in alphabetical order contributed data and commented on the draft.

Appendix S1

Click here for additional data file.

Appendix S2

Click here for additional data file.

Appendix S3

Click here for additional data file.

Appendix S4

Click here for additional data file.

Appendix S5

Click here for additional data file.