Research ArticleECOLOGY

Erosion of global functional diversity across the tree of life

See allHide authors and affiliations

Science Advances  26 Mar 2021:
Vol. 7, no. 13, eabf2675
DOI: 10.1126/sciadv.abf2675

Abstract

Although one-quarter of plant and vertebrate species are threatened with extinction, little is known about the potential effect of extinctions on the global diversity of ecological strategies. Using trait and phylogenetic information for more than 75,000 species of vascular plants, mammals, birds, reptiles, amphibians, and freshwater fish, we characterized the global functional spectra of each of these groups. Mapping extinction risk within these spectra showed that larger species with slower pace of life are universally threatened. Simulated extinction scenarios exposed extensive internal reorganizations in the global functional spectra, which were larger than expected by chance for all groups, and particularly severe for mammals and amphibians. Considering the disproportionate importance of the largest species for ecological processes, our results emphasize the importance of actions to prevent the extinction of the megabiota.

INTRODUCTION

Earth is entering a sixth mass extinction period triggered by human activities (1) in which nearly 1 million species are estimated to be at risk of extinction (2). However, counts of threatened species do not fully reflect the ecological and evolutionary impacts of extinctions because species’ responses to environmental changes and their contributions to ecosystem functioning depend on their functional traits (3, 4). Consequently, extinctions of species with unique traits are likely to have more marked consequences than extinctions of species with redundant traits (5, 6). Yet, little is known about the impacts of the ongoing mass extinction on the global functional diversity for most organisms. Improving our understanding of the factors modulating species’ extinction risk is critical for conservation (79). These factors include functional traits—morphological, physiological, phenological, or behavioral features that govern the functional role of organisms and the effects the environment has on them (10, 11).

Trait variation among species is remarkable, encompassing differences of several orders of magnitude. For example, mammalian body size ranges from less than 2 g for shrews to more than 100 tons for whales, and plant seed mass ranges from less than 1 μg for some orchids to more than 15 kg for coconut. Despite all this variation, species’ ecological strategies resulting from trait combinations are constrained by physiological limits set by evolutionary history and trade-offs in resource allocation (12, 13). Accordingly, recent mappings of the global trait spectra of plants (14), and birds and mammals (15) have described them as two-dimensional surfaces, where occupation of the trait space is restricted compared to null expectations. There are still no characterizations of the global trait spectra of other groups of vertebrates, which precludes understanding how extinction risk is distributed within these functional spaces. So far, the few characterizations of the impact of species extinctions have only considered reductions in the amount of functional space occupied by each taxonomic group (1517). However, functional redundancy among species [i.e., adjacent species in the functional space (18, 19)] is widespread. This means that it is likely that many of the functional consequences of extinctions do not only affect the overall volume and boundaries of the functional spectra but deeply reorganize their internal structure. These changes can be better examined with probabilistic approaches that consider shifts in the density of occupation of the functional space (18, 20), hence fully accounting for the potential effect of functional redundancy.

Here, using traits for more than 75,000 species of vascular plants, mammals, birds, reptiles, amphibians, and freshwater fish, we describe the impacts of potential species extinction scenarios on functional diversity across the multicellular tree of life. By mapping extinction risk in the functional space of each group, we found that extinction risk is not randomly distributed but localized in certain areas of the functional space occupied by species with large size, slow pace of life, or low fecundity. We show that extinctions will lead to a denser functional aggregation [functional homogenization (16, 21)] of species at the global scale, along with increased vulnerability of large portions of the functional space. Potential extinctions will therefore cause marked erosion and rearrangement of ecological strategies across the tree of life.

RESULTS

Uneven redundancy in the global functional spectra

We used functional trait information from 39,260 species of vascular plants, 4953 mammals, 9802 birds, 6567 reptiles, 6776 amphibians, and 10,705 freshwater fish from different published databases (17, 2224). For each of these groups, we chose a set of fundamental functional traits associated with different key aspects of their ecology (table S1). Because the functional trait databases for none of the groups were complete, we performed a phylogenetically informed trait imputation procedure (see Materials and Methods). We first characterized the global functional space occupied by each group by means of principal components analyses (PCA) based on the compiled traits (14, 15). We explored the performance of the imputation procedure by artificially removing trait values and then imputing traits from species with complete trait information in each group and examining the difference between the original position of the species in the corresponding PCA and the position after imputations. Our results show that the imputation procedure performed well in retrieving the positions of species in the functional space for all groups, even when a high proportion of trait information was missing [table S2; see (25, 26) for references about imputation performance measurements].

For five groups, the resulting functional spaces consisted of two main dimensions, which captured 68 to 83% of the total functional trait variation. In the case of freshwater fish, the main functional space extended over four dimensions (table S3). In all cases, the first principal component was linked to traits related to the size of the organisms, such as plant height, body length, or body mass (Fig. 1). The second principal component was frequently related to traits linked to reproduction, such as the frequency and amount of offspring produced. The estimated functional spaces when considering only species with complete trait information were always very similar to the ones that considered imputed traits, both in terms of the loadings of the individual traits in the different components and in terms of the concordance in the position of the species in both spaces (Procrustes tests with 1000 permutations for all groups: P = 0.001).

Fig. 1 Global functional spectra.

Global functional spectra of plants (A), mammals (B), birds (C), reptiles (D), amphibians (E), and freshwater fish (F). Probabilistic species distributions in the spaces defined by the two first principal components (PC1 and PC2) of PCA (details in table S3) considering different functional traits for each group (see table S1 for definitions of each functional trait). Arrows indicate the direction and weighting of each trait in the PCA. The color gradient (red, yellow, and white) depicts different density of species in the defined space (red areas are more densely populated). Arrows show the loadings of the considered traits in the resulting PCA. Thick contour lines indicate the 0.5 (hotspots, see main text) and 0.99 quantiles, and thinner ones indicate quantiles 0.6, 0.7, 0.8, and 0.9. Silhouettes were downloaded from PhyloPic (www.phylopic.org). The legends within each panel show the amount of functional space (measured in SD units, SD2, except for freshwater fish, in SD4) occupied by the hotspots (0.5 quantile) and the 0.99 quantile distribution. sla, specific leaf area; ln, nitrogen content per unit leaf mass; sm, seed mass; la, leaf area; ph, plant height; ssd, specific stem density; ls, litter size; svl, snout-vent length; bm, adult body mass; long, longevity; wea, weaning length; fmat, time to reach female maturity; gest, gestation length; em, egg mass; fa, fledging age; inc, incubation time; am, age at maturity; bs, body size; os, offspring size; es, eye size; mp, mouth position; pp, pectoral position; ms, mouth size; elo, body elongation; cs, caudal peduncle throttling; ps, pectoral fin size; wid, lateral body shape.

To account for functional redundancy between organisms, we used a trait probability density (TPD) approach, which allows mapping the functional spectra of organisms as probabilistic surfaces (18, 20, 27). Instead of simply characterizing the boundaries of the spectra, the TPD approach reflects the abundance of species with similar suites of traits (14, 18). This method effectively represents the functional spectra as a landscape with “peaks” that reflect areas with high density of species and “valleys” where the density is lower and allows to test changes in aspects of trait space occupation other than volume. Specifically, we examined some properties of these spectra, such as the degree of aggregation of species in particular areas of the functional space [functional hotspots; (14)] and the positioning of these areas within the occupied space. We estimated the TPD of each group considering the dimensionality of the corresponding functional space (i.e., two dimensions for all groups except freshwater fish, with four dimensions).

Our analyses reveal high unevenness in the patterns of occupation of the functional space by species in all groups. On one hand, all groups presented a high degree of lumpiness and hence of functional redundancy (4, 20). As a result, the amount of functional space occupied was always much smaller than in equivalent multivariate normal distributions (fig. S1). In particular, the hotspots of the different groups [i.e., the smallest portion of functional space including 50% of the species (14)] consistently had a small extent. For example, in the case of birds, half of the species occupied only 9.8% of the total spectrum, whereas amphibians were less aggregated (18.2% of the total spectrum). The plants and mammals’ spectra displayed multiple distinct hotspots. For plants, one of the hotspots is occupied by grasses and herbs and the other by angiosperm tree species in agreement with (14) (Fig. 1A). Mammals showed two main hotspots separated along the second PCA axis, one corresponding to species that reproduce often and produce many offspring with relatively short gestation times, whereas the other included species with smaller reproductive outputs and longer gestation times. The third mammal hotspot, much smaller, included primarily primate species (suborder Simiiformes) with long life spans and late weaning times (Fig. 1B). Birds, reptiles, amphibians, and freshwater fish displayed a single hotspot, with skewed distributions of species in their size-related axes for all taxa but fish (even after log transformation before PCA). These spectra revealed large proportions of small species and much fewer species of larger size (Fig. 1, C to F), in accordance with previous observations (15, 28). Therefore, the hotspots tended to be close to the boundaries of the occupied space, as revealed by the high values of functional divergence (fig. S1). In contrast, the freshwater fish hotspot was centered within the functional space, revealing that median-sized species with generalist morphology represent the core of the freshwater fish fauna.

On the other hand, large proportions of the spectra for each group displayed extremely low or no redundancy. In this sense, the amount of the total space that was occupied by a single species was comparable to the size of the hotspots, ranging between 9.5% (mammals) and 16.7% (freshwater fish). Further, around one-third of the space (average 32.9%; minimum 26.5% for mammals and maximum 39% for freshwater fish) was occupied by five or fewer species (fig. S2). These areas of the functional space and the species occupying them could be considered to be of critical conservation importance, because losing them would imply the complete disappearance of their functional strategies from Earth.

Mapping extinction risk

In light of the previous results, we used International Union for Conservation of Nature (IUCN) categories to map extinction risk in the functional trait space of each group by means of generalized additive models [GAMs; (29)], which allowed us to discover substantial differences in extinction risk between functional strategies. The position of a given species along the multivariate trait space informs about its conservation status in all the taxonomic groups examined. Threatened species showed distinct occupation patterns of the different functional spaces, evidenced by the lower similarity between the global spectra and the spectra of threatened species than between the global spectra and the spectra of nonthreatened species in all taxa (table S4 and fig. S3). Accordingly, mapping of threat probabilities revealed substantial increases in the risk of being threatened in areas of the functional space occupied by species with larger sizes and slower or less copious reproductive outputs (Fig. 2 and fig. S4). In particular, for plants, the probability of being threatened was up to three times higher for woody species than for herbaceous ones (Fig. 2A). Mammal species with long weaning and gestation periods, with relatively large sizes (right end of the functional spectra), showed up to eight times higher threat risk than smaller species toward the left side of the spectrum (Figs. 1B and 2B). The pattern for birds was similar, with species with later fledging ages, longer incubation times, and larger sizes having up to six times higher threat risk than smaller species with faster breeding times (Figs. 1C and 2C). For two groups, freshwater fish and reptiles, in addition to large species, threat risk also increased toward species with smaller size (Fig. 2, D and F).

Fig. 2 Risk of extinction in the functional space.

Extinction risk in the functional spaces of plants (A), mammals (B), birds (C), reptiles (D), amphibians (E), and freshwater fish (F). Probability of species being classified as threatened (see Materials and Methods) according to GAMs (with binomial distribution) using the position of species in the functional space as predictors. Yellow tones indicate lower risk of extinction, whereas purple tones indicate high risk of extinction. Legend shows χ2 and P values of the GAM corresponding to each taxonomic group. For each group, the red contour lines indicate the average threat probability (proportion of species classified as threatened in the group). The gray line indicates the 0.99 quantile of the spectra of each group considering only species whose threat status is known.

We then used different simulated scenarios of extinction based on the IUCN species’ extinction risk assessments to define a continuous range of potential future global functional spectra after different numbers of threatened species have gone extinct. We compared these spectra with alternative scenarios where the same number of species goes extinct randomly. Matching the current functional spectra with the spectra after possible extinctions allowed us to reveal the degree of “erosion” that could be potentially experienced by the global functional diversity, exposing functional areas that will become particularly vulnerable to further extinctions.

Our extinction simulations revealed that the high degree of aggregation observed in the examined taxonomic groups buffers the effects of extinctions on the total functional space occupied. Accordingly, projected losses of functional space (functional richness) after extinctions were substantially lower than the proportion of lost species, ranging between 0.3% (reptiles) and 5.1% (freshwater fish; Fig. 3, A to F) in the central scenario. The effect of extinctions on functional richness followed linear trajectories between the most optimistic scenario considered (50% fewer extinctions than in the central scenario) and the central scenario. However, nonlinear responses appeared for some groups between the central and the most pessimistic (50% more extinctions than in the central scenario) extinction scenarios. These nonlinear responses reflected contrasting organizations of the most vulnerable species among groups, with reptiles’ functional space occupation starting an abrupt decline in the pessimistic scenarios (Fig. 3D) and plants and mammals becoming less sensitive to extinctions near the most pessimistic scenario (Fig. 3, A and B). Despite the apparently small effects of extinctions, the proportion of lost space after extinctions was higher than expected under a random species extinction hypothesis for plants, mammals, birds, and freshwater fish, revealing that the diversities of functional strategies of these groups are particularly vulnerable to extinction.

Fig. 3 Effects of simulated extinction scenarios on different functional diversity indicators of the different taxonomic groups.

Values of functional richness (top row), functional redundancy (central row), and overlap with the current TPD functions (bottom row) of the TPD functions estimated after simulated extinctions (current values of each index are denoted by the value 1 in all plots). The x axis of all figures correspond to different extinction scenarios (total of 101 scenarios with 100 repetitions each), ranging from 50% fewer extinctions than under the central scenario (based on IUCN 100-year extinction probabilities) to 50% more extinctions than under the central scenarios (numbers of extinct species per group in these three cases are shown in the bottom row). Lines represent the fits of GAM (mean ± 2 SE) in which smooth functions of the number of species gone extinct were fitted for both realistic extinctions (based on IUCN categories; orange line) and random extinctions (blue lines). In all cases, the AIC (Akaike information criterion) values of the model including two smooths were much smaller (difference > 10 AIC) than a model including a single smooth for the number of species extinct; hence, the two lines are always represented. sp., species.

The resistance of functional richness to extinctions was brought by the high functional redundancy of species, which makes redundancy a key factor to understand the response to extinctions in the global spectra. Changes in redundancy after extinctions reflected to a large extent the differences among groups in the positions of threatened species within the functional space. For example, threatened plant species tend to be within the woody-plant hotspot (Fig. 2A) and hence to be functionally redundant so that simulated extinctions reduced average redundancy more than expected by chance (Fig. 3G). By contrast, threatened mammals, birds, reptiles, or freshwater fish tend to be closer to the boundaries of the trait distributions (and hence to be functionally unique) so that their extinction leads to smaller reductions in redundancy than expected by chance (Fig. 3, H to J and L).

Nonetheless, projected extinctions will erode the functional spectra in ways that go beyond the amount of functional space being completely lost or the redundancy of the remaining species. This was revealed by the changes in overlap between the present-day TPD and the TPD after extinctions. In all groups, reductions in overlap with extinction (reflecting overall changes of the functional organization of the group at a global scale) were much larger than expected following random extinctions, and this difference tended to increase with the number of predicted extinctions (Fig. 3, M to R). To better understand these changes, we compared the shifts in the patterns of occupation of the functional space of all groups before (present-day) and after extinctions. These shifts were particularly notable for mammals, amphibians, and freshwater fish. In the case of mammals, most of the projected functional diversity erosion would take place close to the boundaries of the functional spectra. Thus, further extinctions would increase the risk of completely losing parts of the spectra corresponding to mammal species with high longevity, late sexual maturity, and long gestation and weaning periods (Fig. 4B). In particular, the hotspot formed by primate species is projected to largely decrease because of the high proportion of threatened species in that group (Fig. 4B and fig. S5B). The shift in the functional spectra of amphibians after extinctions, characterized by a decrease in the relative proportion of species with small reproductive outputs, was also remarkable (Fig. 4E and fig. S5E). For freshwater fish, functional changes were mainly clustered toward particular boundaries of the functional space, corresponding, for instance, to the large-sized species (Fig. 4F and fig. S5F), which are, as for other taxonomic groups, long-lived species.

Fig. 4 Shifts in the functional spectra after extinctions.

Shifts in functional spectra after simulated extinctions for plants (A), mammals (B), birds (C), reptiles (D), amphibians (E), and freshwater fish (F). Differences (expressed in quantile changes) between the functional spectra of species assessed by IUCN before and after simulated extinctions based on IUCN category–based probabilities of extinction. For each group, the plot represents the average of the TPD functions of the 100 simulations performed for the central scenario of extinctions (see Materials and Methods). Brown tones reflect areas in which the TPD quantiles are lower after extinctions (i.e., those traits become relatively less frequent at the global scale), and blue tones reflect areas in which the TPD quantiles increase (i.e., those traits become relatively more frequent at the global scale). Black areas show the parts of the functional space that would disappear from the global spectra after extinctions (trait combinations that go completely extinct). The legend for each group shows the average absolute change in quantiles across the spectra and the range of these changes. Figure S5 shows similar plots for the two extreme extinction scenarios considered (−50 and +50% extinctions).

In general, the higher risk of extinction of large, long-lived and slow-reproducing species will result in decreased functional redundancy for these species at the global scale. In particular, our simulations projected generalized increases across groups in the proportion of the total functional space that is occupied by a single species. These increases were particularly marked for mammals (17% increase in the central scenario and 39% increase in the most pessimistic scenario), amphibians (23 and 35% increase, respectively), and freshwater fish (20 and 31% increase; fig. S2). This result shows that, although functional redundancy is able to buffer the effect of extinctions on the amount of globally occupied functional space, the ongoing erosion of the functional spectra of the different taxonomic groups will result in overall higher susceptibility to future extinctions.

DISCUSSION

Using the most species diverse trait databases collected to date, we show that, in most groups, realized functional strategies are constrained to a single plane in which species are clumped around a few strategies that are rather prevalent (14, 15). Contemporary functional spectra are probably quite different from those in the relatively recent past because of extinctions in the last millennia (30, 31). For example, since the Pleistocene, larger animal species have gone extinct at much higher rates than smaller species (3133). This trend continues nowadays, boosted by the rise of human population and technologies that disfavors large and slow-pace-of-life species (8, 31, 34, 35). In addition, the smallest animal species are also among the most threatened, as already reported for fish and reptiles, because small animals often have low dispersal abilities and remain endemic from restricted areas, making them vulnerable to extinction (9, 36, 37). Our models show that large-sized, slow-paced, and slow-reproducing species are much more likely to go extinct across all groups and provide some support for the notion that small species are also more threatened in the case of reptiles and freshwater fish. Mammals and amphibians emerge as the groups most affected by these changes, experiencing shifts in the functional space toward faster-living and highly reproductive strategies, respectively. These results are in line with the notion that larger organisms are more sensitive to global change, a trend affecting plants, land, and aquatic animals that will probably be exacerbated in the future (30, 38).

We show that estimated extinctions will cause noticeable reorganizations in the functional spectra of most groups, making them more vulnerable to further extinctions in the future. Although the impacts of extinctions on multivariate functional diversity is increasingly being studied, research focus has so far been on quantifying losses in the amount of functional space occupied (1517). However, functional diversity encompasses all aspects of trait variation among organisms, of which total occupation (functional richness) is only one component (18, 39). Here, we show that species extinction impacts on functional richness alone appear to be much milder than the intense erosion observed when the whole spectra is considered.

We show that understanding the risks of losing unique functional strategies requires considering the interplay between the organization of species within the functional space (6) and the risk of extinction of species (35, 40). For example, we found that the proportion of functional space occupied by a single species in mammals is much smaller than in birds (fig. S2). However, the disproportionally higher extinction risk for mammals close to the boundaries of the functional space results in a much higher predicted reduction of the global functional space of mammals in the near future. A similar trend toward higher isolation of species in the functional space since Late Pleistocene has been observed using fossils of large North American mammals (41). Indicators that combine information about the functional uniqueness of species (i.e., which species occupy low-redundancy areas of the functional space) and their risk of extinction, such as the FUSE (functionally unique, specialized, and endangered) metric (35) appear as a major opportunity for ranking species in terms of their conservation priority to safeguard functional diversity.

Our findings are relevant for biodiversity management and conservation policies. For example, the planetary boundary (PB) framework (42) aims to link global change, biodiversity, and sustainability by assessing the degree to which human perturbations will destabilize the Earth system. In the PB framework, the current status of the functional role of biosphere integrity has been found to be hard to quantify because of the lack of an adequate control variable indicating how far the system is from suffering serious impacts and the lack of general metrics to assess functional diversity at regional or global scales (42, 43). To date, most analyses of functional diversity conservation have focused on the loss of total functional space occupied by different groups (15, 16, 35). However, here, we show that changes in functional richness are strongly buffered by functional redundancy among species. As a more comprehensive way to reflect the full impact of species extinctions on global functional diversity, we propose the use of the overlap between the TPD functions before and after human impacts (e.g., between before industrial era and in the present time or between present time and after extinctions, as done here) as a suitable control variable for functional diversity. For example, we show that the effect of predicted extinctions in the next 100 years is higher than expected from random extinctions in all studied groups, reaching values of more than 5% reductions in overlap in the case of amphibians (Fig. 3Q) under the central scenario of extinctions. In addition to being bounded between 0 (complete shift in the functional space occupation patterns between the compared times) and 100% (no change), overlap can be estimated at any spatial scale (18) so that the indicator can be used to estimate the integrity of individual biomes (43). In addition, TPD functions can incorporate changes in the population abundances of species, hence providing extreme flexibility to study any driver of biodiversity change in addition to species extinctions (44, 45).

Several improvements will be required to fully incorporate functional diversity assessments into conservation policies. First, the functional spaces that we considered here are by no means complete, being limited by factors such as trait availability. For example, the plant functional space does not incorporate any belowground traits although roots can represent up to 70% of total plant biomass (46) and are crucial for species coexistence, ecosystem processes, and symbiosis with other organisms (47). Despite recent increases in the availability of root traits data (48), the proportion of species with such information is still very limited. In the case of vertebrates, there are very large differences among groups in terms of which traits (and of what kind) are available (49). While for freshwater fish we used a set of morphological traits that reflects the ecological roles of fish species (50), in the case of terrestrial vertebrates, we aimed to use comparable sets of life history traits, which led to the definition of two-dimensional spaces with a size-related and a reproduction-related axis for each group. The future incorporation of other traits, such as diet or habitat information (15), as they become more widely available across groups will surely strength the link between changes in the occupation of trait spaces and changes in ecosystem functioning.

It is important to note that the extinction effects in the global spectra are likely even greater in local communities. For example, although extinction risk is much higher for tree species, the global spectrum of plants is not likely to experience marked shifts because of high functional redundancy. However, in local communities with fewer species, functioning could be markedly altered even by the extinction of a low proportion of species and amplified by cascades of secondary extinctions (51, 52). These impacts are likely to be particularly important if the lost species are the largest ones, given their disproportionate importance for ecological processes (38). Understanding how the effects of extinctions downscale from the global scale into smaller scales (continents, realms, regions, and local communities) emerges as an important priority to understand the impacts of extinction on ecosystem functioning. In this sense, groups such as prokaryotes, protists, insects, or fungi, which have great importance for ecosystem functioning, are absent from our assessment. Both trait and conservation information are lacking for a very large proportion of the species in these groups, but we hope that future initiatives will eventually lead to their inclusion in global assessments of functional diversity.

MATERIALS AND METHODS

Data collection and processing

Functional traits and phylogenies. We collected published information on functional traits for all the studied groups of organisms from different sources (see table S1 for detailed descriptions of each trait)

Vascular plants

We used six traits previously shown to capture global spectrum of plant form and function (14): plant height (ph, m), specific stem density (ssd, g/m3), leaf area (la, mm2), specific leaf area (sla, mm2/mg), nitrogen content per unit mass (ln, mg/g), and seed mass (sm, mg). We used publicly available data for these traits from the latest version of the TRY Plant Trait Database [version 5.0, www.try-db.org/TryWeb/Home.php, accessed April 2019 (53)]. Together, our dataset included over 955,000 trait measurements for 44,431 vascular plant taxa. In the analysis, each taxon was represented by an average trait value (excluding outliers with >3 SD). To account for within-species variation, the averages for each species-trait combination were calculated first within individuals (if multiple measurements were taken from a single individual), then within datasets (if multiple individuals were measured in the same location), and last, within species (if multiple individuals were measured in various locations).

Plant height data included 179,263 measurements of adult plant vegetative height for 15,008 taxa. In most datasets, this was represented as observed height or average of measurements. In some cases, plant height was represented as the maximum observation (12,942 records). Specific stem density (ssd) data included 28,723 measurements for 8840 taxa. As this trait is usually measured for woody species, we estimated ssd for herbaceous plants using leaf dry mass content information (127,067 measurements for 5764 taxa), following the procedures described in (14). Leaf area data included 119,172 measurements for 13,928 taxa. Different datasets in TRY reported various measurements of leaf area (e.g., leaflet or leaf, petiole included or excluded). To maximize our data coverage, we included both leaflet and leaf measurements and preferred measurements including petiole (if both data types, petiole included or excluded, were reported for the same individual). Specific leaf area data included 223,126 measurements for 10,674 taxa. Similarly to leaf area data, we preferred measurements that included petiole. Data for N content per unit leaf mass included 92,850 measurements for 10,530 taxa. Data for seed mass included 185,182 measurements for 25,394 taxa.

Mammals, birds, and reptiles

We used the Amniote database (24) including data for 4953 species of mammals, 9802 species of birds, and 6567 species of reptiles. The database includes data for 29 life history traits, but information is very incomplete for many of them (see table S1). Hence, for each group, we selected subsets of traits with sufficient information (at least 1000 species). In the case of mammals, we selected eight traits: litter size (ls, number of offspring), number of litters per year (ly), adult body mass (bm, g), longevity (long, years), gestation length (gest, days), weaning length (wea, days), time to reach female maturity (fmat, days), and snout-vent length (svl, cm). We also selected a total of eight traits for birds: clutch size (ls, number of eggs), number of clutches per year (ly), adult body mass (bm, g), incubation time (inc, days), longevity (long, years), fledging age (fa, days), egg mass (em, g), and snout-vent length (svl, cm). Last, we selected six traits for reptiles: clutch size (ls, number of eggs), number of clutches per year (ly), adult body mass (bm, g), incubation time (inc, days), longevity (long, years), and snout-vent length (svl, cm).

Amphibians

We used the AmphiBIO database (23) to get data for 6776 species of amphibians (see table S1). Within this dataset, we selected a total of four traits with clear correspondence with the traits selected for the other groups of terrestrial vertebrates: age at maturity (am, years); body size, measured in Anura as snout-vent length and in Gymnophiona and Caudata as total length (bs, mm); maximum litter size (ls, number of individuals); and offspring size (os, mm).

Freshwater fish

We used the last updated version of the most comprehensive database to get morphological traits, available for 10,705 species of strictly freshwater fish (17, 54). Morphological measures and conventions are detailed in previous studies (17, 55). All morphological traits have been measured on side view picture, using one specimen per species and applying conventional rules for unusual morphologies (e.g., species without tail and flatfishes) as defined in previous studies [see details in (17, 55)]. Briefly, this database encompasses 11 traits describing size and shape of body parts involved in food acquisition and locomotion (50). Fish body shape and weight were described through the size using the standard length (svl) and body mass (bm) taken directly from FishBase (56), body elongation (elo, ratio between body length and body depth), and body lateral shape (bls, ratio between head depth and body depth). The other traits describing the position and the size of each part of the fish were eye size (es) and position (ep), mouth size (ms) and position (mp), pectoral fin size (ps) and position (pp), and caudal peduncle throttling (cs).

Phylogenies and imputation of missing trait values

We obtained published phylogenies for each of the considered groups (5763). Species that were not present in the phylogeny were added to the root of their genus, using the “add.species.to.genus” from the R package phytools (64). Species for which we did not have evolutionary information were removed from the trait databases before the missing trait imputation procedure (see below).

Because none of the functional trait databases assembled were complete, we completed this information by performing a trait imputation procedure for each group using the missForest R package (16, 65). Before the imputation process, all traits were log10–transformed, centered, and scaled. The missForest package uses random forest techniques to impute trait data, which allows to include phylogenetic information in the trait imputation process, which is known to improve the estimations of missing values (26). We included the evolutionary relationships between species in the imputation process by including the first 10 phylogenetic eigenvectors in the matrix to be imputed, as recommended in (26). The final numbers of species with functional trait information used in each group are shown in table S4.

The accuracy of trait imputation procedures depends on factors such as the proportion of missing observations and the potential biases in trait measurements. Although some studies have analyzed the performance of imputation methods in the context of retrieving individual trait values (25, 26), the goal of our imputations was rather to characterize the position of species in the corresponding trait space established using PCA based on the individual traits. The facts that (i) many traits are strongly correlated and (ii) evolutionarily closely related species tend to be close in the functional space effectively mean that the position of species in the PCA should be easier to estimate than the individual trait values. Following Penone et al. (26), we explored the performance of the imputation procedure by artificially removing trait values from a subset (10%) of the species with complete information from each group. For each of the species with complete information that was selected for artificial removal of trait information, we selected one random species with incomplete information and superimposed its pattern of missing values. This way, we attained a pattern of missing values that was consistent with the one in the original dataset. We then combined the three subset of species (90% of species with complete trait information, 10% of species with complete trait information plus artificial missing values, and all species with noncomplete trait information) into a single dataset in which we performed the phylogenetically informed imputation procedure as described above. Then, we used the imputed traits to project species onto the functional space based on the full dataset (see the “Construction of the global spectra” section below). Last, we compared the position of the species for which we had artificially removed traits in the different dimensions of the trait space without removing data (real position of the species) and after artificial removal and imputation of trait information. We estimated the performance of the imputation by means of the normalized root mean square error (NRMSE), which expresses the average distance between real and imputed positions of species as a proportion of the range of values of species in the corresponding dimension (26). We repeated this examination 100 times for each taxonomic and attained an estimation of the average NRMSE value across repetitions.

Conservation status of species. We collected conservation status of species from the IUCN Red List (66) (retrieved 25 September 2019) using the R package ‘rredlist’(67). We reclassified the IUCN categories as “threatened” (including the “extinct in the wild”, “critically endangered,” “endangered,” and “vulnerable” categories) or “nonthreatened” (including the “near-threatened” and “least concern” categories).

Taxonomic standardization. Taxonomies from all the used sources (trait databases, phylogenies, and IUCN Red List), were standardized using the R packages taxize (for animals) (68) and Taxonstand (for plants) (69). In the case of animals, all names were resolved against the Global Biodiversity Information Facility (GBIF) Backbone Taxonomy (70), whereas in the case of plants, we used The Plant List (71).

Construction of the global spectra

We identified the main axes of functional trait variation by performing PCA on the log-transformed and scaled functional traits of each group. We used Horn’s parallel analysis in the R package paran (72) to determine the number of axes retained in these PCA; we will refer to these reduced spaces as functional spaces from now on. We checked the reliability of the functional spaces obtained with imputed functional trait values by comparing them with the spaces that were based only on species with complete functional information. We performed this comparison by estimating the correlation between distance matrices of the species that were common to the two spaces (space with imputed data and space with complete species only) through a Procrustes of each taxonomic group (14), using the ‘procuste.rtest’ function from the R package ade4 (73). To assess the significance of the correlation, permutation tests (9999 randomizations) based on Monte Carlo simulations were generated. All the Procrustes tests were highly significant (P = 0.0001 in all cases; see table S3), indicating a strong correspondence between the complete and imputed functional spaces; consequently, we used the PCA based on imputed trait data in the rest of analyses.

We estimated the probabilistic distribution of the species within the functional spaces by performing multivariate kernel density estimations with the ‘TPD’ and ‘ks’ R packages (20, 7476). The kernel for each species was a multivariate normal distribution centered in the coordinates of the species in the life history or functional space and bandwidth chosen using unconstrained bandwidth selectors from the ‘Hpi’ function in the ‘ks’ package (77). The aggregated kernels for all species in a group result into the TPD function (20, 27, 78) of that group in the corresponding space. Although TPD functions are continuous, to perform operations with them, it is more practical to divide the functional space into a D-dimensional grid composed of many equal sized cells (we divided the two-dimensional spaces in 40,000 cells, 200 per dimension, and the four-dimensional space in 810,000 cells, 30 per dimension; previous sensitivity analyses showed that 30 divisions per dimension in the four-dimensional case provide results that are virtually indistinguishable from those attained using larger numbers of divisions but require much less memory and computation time). Then, the value of the TPD function is estimated for each cell. The value of the TPD function in a given point of the space reflects the density of species in that particular area of the space (i.e., species with similar functional traits). For each of these spaces, we represented graphically the global TPD and the contours containing 50, 60, 70, 80, 90, and 99% of the total probability.

We compared the distribution of species within the different functional spaces with a null model considering that species are distributed following a multivariate normal distribution (14). For this, for each taxonomic group, we drew 199 samples of 1000 simulated species from multivariate normal distributions with the same mean and covariance matrix as the observed spectra. For each of these samples, we estimated a TPD function and measured functional richness [amount of space occupied by the spectra (18, 20)] at the 99 and 50% quantile thresholds and functional divergence [which indicates the degree to which the density of species in the functional trait space is distributed toward the extremes of the spectra; (20, 38)]. Then, we drew 199 samples of 1000 species from the observed global species pool of each taxonomic group and performed similar analyses. We compared the estimations of functional richness (at 50 and 99% quantile thresholds) and functional divergence of the observed and simulated data by means of two-tailed t tests (after checking that the data were normally distributed by means of Shapiro-Wilk tests).

Effects of extinctions on global functional diversity

We estimated TPD functions in the different functional spaces for the sets of species classified as threatened or nonthreatened, as well as the TPD functions for all the species assessed by IUCN (including both threatened and nonthreatened), using the same procedure described above. We then estimated the similarity between these TPD functions and the global spectra (the global TPD function considering also species not assessed by IUCN) as the overlap between the TPD functions (18, 78, 79). Compared to methods that consider exclusively the boundaries of the distributions [e.g., hypervolumes or convex hulls (15, 16, 80)], TPD-based probabilistic overlap considers also the differences in density within those boundaries. This approach provides a more complete idea of what the differences between the functional spectra are, particularly in cases where functional redundancy is high (6, 78). Given that a high proportion of the considered species might be clumped in particular areas of the considered space (14, 15), this methodological aspect can be particularly useful to detect differences in the occupation of functional spaces between groups of species with different conservation status. Estimating the similarity between the different groups allowed us to examine (i) whether there is any bias regarding which species have been assessed by IUCN (overlap between the global distribution and the IUCN TPD functions) and (ii) whether nonthreatened and threatened species occupy the considered space in different ways (overlap between the nonthreatened and the threatened TPD functions). In all vertebrate taxa, the species whose conservation status has been assessed by IUCN formed a relatively random subset of the functional spectra of the different groups (high similarity between the assessed species and the global spectra; table S4). Plants constituted an exception to this pattern, reflecting the bias toward trees in the IUCN Red List (fig. S3) (81), as well as the high proportion of species for which conservation status is not known (in our dataset, roughly 80% of the plant species with trait measurements have not been assessed or are classified as data deficient).

After examining the overlap of TPD functions, we mapped the conservation status of species within the functional spaces. For this, we considered only species assessed by IUCN. The relationship between conservation status (1: threatened and 0: nonthreatened), and the position in the corresponding space (PCA axes) was analyzed using a tensor product smoother-based GAM (29) with a binomial response [R package ‘mgcv’ (82)]. We then mapped the predictions of the models (including the 95% confidence intervals of the means) to visually examine how different combinations of functional traits affect the probability of species being threatened.

Then, we explored the potential effects of extinctions on the distribution of species in the functional spaces. For this, in each space, we estimated a TPD function considering all the species assessed by IUCN, which represents the present-day functional spectrum of the corresponding group. We then estimated the erosion of the global functional diversity of each group following different simulated scenarios of future extinctions. Following Cooke et al. (15), we created a “central” scenario in which we assigned probabilities of extinction over the next 100 years to species based on their IUCN categories: 0.999 for critically endangered species, 0.667 for endangered species (EN), 0.1 for vulnerable species, 0.01 for near-threatened species, and 0.0001 for least concern and data deficient. Using these probabilities, we estimated, for each group, the expected number of species that are expected to go extinct (e.g., our plant dataset included 396 species classified as EN, which corresponds with 264 projected extinctions). We then randomly selected from each IUCN category as many species as its expected number of extinctions and estimated a TPD function representing the distribution of species within the functional spaces after extinctions. We repeated this process 100 times to incorporate uncertainty about the species that might go extinct. In each repetition, we estimated functional richness, functional redundancy—which indicates how many species on average are present in each cell of the functional space (1820)—and the overlap between the simulated TPD function and the present-day TPD function. We performed the same estimations using 100 repetitions of a null model in which the same number of species were lost at random from the total set of IUCN-assessed species. This analytical strategy allowed us to ascertain whether losing threatened species affects the functional spectra of the different groups more or less than expected by chance in terms of total space lost (functional richness), vulnerability to further extinctions in the future (reductions in redundancy), and overall similarity with the present-day distribution of traits (overlap).

In addition, we created a series of alternative scenarios representing either fewer extinctions due to increased efforts in conservation, restoration, and climate action [corresponding to up to 50% reduction in the total number of projected extinctions (83, 84)] or increased extinctions due to lack of these actions (corresponding to up to 50% increase in the total number of projected extinctions). Between these two extremes, we simulated 101 of these scenarios using intervals of 1% in changes in projected extinctions and performed the same calculations as described for the central scenario (i.e., 100 repetitions using extinctions based on IUCN categories and 100 repetitions using random extinctions). We analyzed the changes in functional diversity across the considered scenarios by fitting, for each group and variable (functional richness, functional redundancy, and overlap with present-day TPD function) using GAM. In these models, we created smooth functions of the number of species extinct (from −50 to +50% extinctions) for each type of extinctions (based on IUCN categories or random). To examine whether there were differences between the IUCN-based and random extinctions, we also fitted a model including a single smooth for the number of extinct species and compared it with the previous model using the Akaike information criterion (AIC).

While changes in overlap between TPD functions after extinctions and present-day TPD functions give us a numerical estimation of the changes in global occupation of the functional space, they do not show which combinations of traits are likely to experience the largest impacts. To visualize the erosion of functional space after extinctions, for each group and extinction scenario, we averaged the 100 TPD functions of both the IUCN category–based extinctions and the random extinctions. Averaged TPD functions share the same mathematical properties of TPD functions, namely, they are probability density functions so that they integrate to one across the functional space, and the same operations can be applied to them. We then expressed the probabilities of each TPD in terms of quantiles to ease interpretability of the results. We represented the impact of simulated and random extinctions by subtracting, in each cell, the quantile value of the TPD function after extinctions (realistic or random) from the quantile value of present-day TPD function. Negative values in this index indicate a decrease in the relative abundance of the trait values corresponding to that cell and vice versa. To quantify the sensibility of the whole spectra to extinctions, we estimated for each cell the absolute value of this quantile difference and averaged these values across cells (mean absolute quantile change). While mean absolute quantile changes reflect sensitivity across the whole spectra, lost space quantifies the amount of the functional spectra that goes globally extinct. Last, we explored the expected relative change in the different areas of the functional space of the corresponding group by mapping quantile changes in this space.

SUPPLEMENTARY MATERIALS

Supplementary material for this article is available at http://advances.sciencemag.org/cgi/content/full/7/13/eabf2675/DC1

https://creativecommons.org/licenses/by-nc/4.0/

This is an open-access article distributed under the terms of the Creative Commons Attribution-NonCommercial license, which permits use, distribution, and reproduction in any medium, so long as the resultant use is not for commercial advantage and provided the original work is properly cited.

REFERENCES AND NOTES

Acknowledgments: The study has been supported by the TRY initiative on plant traits (www.try-db.org). The TRY initiative and database is hosted, developed, and maintained by J. Kattge and G.Boenisch (Max Planck Institute for Biogeochemistry, Jena, Germany). TRY is currently supported by Future Earth/bioDISCOVERY and the German Centre for Integrative Biodiversity Research (iDiv) Halle-Jena-Leipzig. Funding: C.P.C., A.T., M.P., and R.T. were supported by the Estonian Ministry of Education and Research (PSG293, IUT20-29, PRG609, and PSG505) and the European Regional Development Fund (Centre of Excellence EcolChange). M.V.-V. and R.G.-M. were supported by the Dora Plus Fellowship Programme (University of Tartu). S.B. was supported by “Investissement d’Avenir” grants (CEBA, ANR-10-LABX-0025 and TULIP, ANR-10-LABX-41). R.S.-G. was supported by NERC IRF NE/M018458/1. F.d.B. was supported by the Plan Nacional de I+D+i (project PGC2018-099027-B-I00). Author contributions: C.P.C., A.T., R.T., M.V.-V., and R.G.-M. conceived the study. C.P.C., A.T., and R.T. integrated and cleaned data. S.B. and A.T. contributed data. C.P.C. and A.T. designed and performed analyses. C.P.C. wrote the first draft of the manuscript, and all authors contributed to interpret and analyze results and contributed to revisions. Competing interests: The authors declare that they have no competing interests. Data and materials availability: The data and R scripts used to make all analyses are available in Figshare (https://figshare.com/s/134db1124b4890d9e67b). All data needed to evaluate the conclusions in the paper are present in the paper and/or the Supplementary Materials. Additional data related to this paper may be requested from the authors.

Stay Connected to Science Advances

Navigate This Article