What North America’s skeleton crew of megafauna tells us about community disassembly

Associated Data

Supplementary Materials

Appendix S1

rspb20162116supp1.docx

(3.8M)

GUID: 9B464BF6-DEF5-4DE8-A581-A343DD2E77F9

Table S1

rspb20162116supp2.xls

(203K)

GUID: 9E4845AB-E793-4BEE-8DF0-53C7367DBA50

Table S3

rspb20162116supp3.xls

(45K)

GUID: BBE1DCD0-3CCB-403D-AFC2-8AD3D455A6A8

Allrange

rspb20162116supp4.csv

(20K)

GUID: 56B4F3C6-BF82-44D6-80F0-0D0AE65F7D88

Davis Convex Hull

rspb20162116supp5.txt

(3.1K)

GUID: C01486F5-DB9A-4056-B215-C54E3C87283F

Davis Functional Diversity Random Simulation

rspb20162116supp6.txt

(6.9K)

GUID: 017667F4-0DB6-4807-8B82-FC07ED66A588

Davis Functional Diversity Simple

rspb20162116supp7.txt

(3.0K)

GUID: 352DD617-AE69-41F3-98BA-4946CE9F9D8A

Davis Functional Diversity Simulation Mass Cutoffs

rspb20162116supp8.txt

(11K)

GUID: C2F3B990-91C8-45C4-AC4D-C0D9704A9220

Davis Functional Diversity

rspb20162116supp9.txt

(6.8K)

GUID: 765562B0-E9F8-40C0-9FFB-A32F3AA97E27

Davis Nearest Neighbor

rspb20162116supp10.txt

(5.9K)

GUID: F0A7F710-A2A8-46AC-A563-C5B48257CA2E

Davis Occurrence Matrix

rspb20162116supp11.txt

(1.0K)

GUID: 5ECDF4CF-66B9-477F-8AE0-4EFADFB956CA

Davis Plot Parameters

rspb20162116supp12.txt

(2.0K)

GUID: A4BF7AD4-00A4-4B3C-BB3D-19DF564A3DC5

Davis Timebinning

rspb20162116supp13.txt

(1.7K)

GUID: 23ACBE95-187E-4D61-9BA5-81A232D3B3CB

Davis Unique Branchlengths

rspb20162116supp14.txt

(4.5K)

GUID: 2E66941B-A50D-4FA8-B13A-F82053A397D6

Davis Functional Diversity Simulation

rspb20162116supp15.txt

(11K)

GUID: 1CA9C054-4817-4351-94B5-22612AC3C9CE

Davis Disassembly Main Code Open First

rspb20162116supp16.txt

(91K)

GUID: 3B8B28F6-30B6-4973-9B96-4C2A5D56E2A2

Small Extinct Mammal Functional Data

rspb20162116supp17.csv

(13K)

GUID: E6DB5BE3-75F4-4C66-87D5-A40C9AFCE818

EltonTraits1.0 North American Mammals

rspb20162116supp18.csv

(114K)

GUID: 718B519A-E53B-4F56-BDAC-859924B742A0

Data Availability Statement

All data used for this study are freely available in the electronic supplementary material and via the Dryad Digital Repository at http://dx.doi.org/10.5061/dryad.1bj88 [51].

Abstract

Functional trait diversity is increasingly used to model future changes in community structure despite a poor understanding of community disassembly’s effects on functional diversity. By tracking the functional diversity of the North American large mammal fauna through the End-Pleistocene megafaunal extinction and up to the present, I show that contrary to expectations, functionally unique species are no more likely to go extinct than functionally redundant species. This makes total functional richness loss no worse than expected given similar taxonomic richness declines. However, where current species sit in functional space relative to pre-anthropogenic baselines is not random and likely explains ecosystem functional changes better than total functional richness declines. Prehistoric extinctions have left many extant species functionally isolated and future extinctions will cause even more rapid drops in functional richness.

Keywords:

functional diversity, Rancholabrean, community disassembly, extinction, megafauna, Pleistocene

1. Background

With the potential, impending Sixth Mass Extinction [1], ecologists and conservationists are increasingly using biodiversity metrics like phylogenetic and functional diversity to triage taxa and predict the path and consequences of community disassembly [2–5]. Although rarely explicitly defined [6], functional diversity is usually thought to represent the distribution of organismal or specific traits in an environment that influences ecosystem function [7]. With taxonomic richness alone, each species’ ecological contribution is implicitly identical. Functional diversity predicts ecosystem function more accurately because it allows these contributions to vary [8–10].

Impressive, global scale models use functional diversity to predict the effects of community disassembly on large taxonomic groups given expected climate change [11,12] but ground truthing these models is difficult. Despite decades of fruitful research on community assembly, we still have a poor understanding of community disassembly and how it relates to functional diversity [13]. The rules of disassembly, it seems, are not simply those of assembly run in reverse [14–16]. There are few long-term studies examining disassembly [13] and most experiments must rely on unrealistic, highly controlled plant communities [10] or space for time substitutions of natural environments that are already denuded by human impacts [14].

Palaeontologists have a rich literature examining functional change through major extinctions and faunal shifts (e.g. [17,18]), but these studies usually have coarse, million year or longer, time bins and evolutionary findings that may not be immediately applicable to the short timescales modern conservation biologists investigate [19]. The frontier for understanding community function’s fate is examining the interaction between realistic disassembly scenarios and functional diversity over large spatial scales with (relatively) small temporal scales [5], something rarely achieved [20]. The latest Quaternary of North America provides an excellent natural experiment for examining massive community disassembly on a continental scale because a high proportion of large mammal species went extinct, probably due to the same stressors impacting vulnerable species today: rapid climate change and growing human populations [21]. Although the still debated causes of this extinction have received the most attention historically, researchers now recognize the important ecological consequences of megafaunal disassembly, which included major changes in vegetation cover, nutrient cycling, and possibly atmospheric composition [21,22].

To better understand disassembly’s role in this functional change, I took the 94 largest mammals that lived in North America during the past 50 000 years (see complete taxa list in the electronic supplementary material, table S1), and calculated how functional diversity changed through the Pleistocene–Holocene boundary, comparing functional traits between extant and extinct species (figures and ). At stake is whether or how response traits (a species’ sensitivity to extinction drivers) correlate with effect traits (the species’ effect on ecosystem function) [14]. There are few direct investigations of this relationship but because certain traits like large body size likely confer unique functional roles (e.g. ability to disperse larger seeds, higher trophic levels) and are known to increase extinction risk across a wide range of taxa, functional losses during disassembly should be greater than expected given the proportion of taxa going extinct [14,21]. If effect and response traits are not correlated, then functional diversity change during a real extinction should be indistinguishable from change given any equal drop in taxonomic richness, regardless of the species, something that is not expected by theory [5] but has been found in some habitat fragmentation studies [16]. Not only would this slow ecosystem collapse but it would also allow the use of taxonomic richness in place of much more difficult-to-collect functional data as on average, extinction order/contingency will not greatly alter the functional outcome. I show that even during a strongly trait selective disassembly like the Pleistocene megafaunal extinction, functional richness decline (the reduction in the volume of trait niche space occupied) may be no different than that of random species loss. However, focusing on total functional richness loss may be a poor way to understand community disassembly as it can obscure non-random trait shifts that have the potential to greatly affect ecosystem function.

Open in a separate window

Open in a separate window

2. Material and methods

This study considered the 94 largest mammal species occurring in North America over the last 50 000 years from the 3.9 kg tayra (Eira barbara) to the 8000 kg Columbian mammoth (Mammuthus columbi). This includes ‘native’ species as well as humans (Homo sapiens) and seven of the most common species introduced by humans like dogs (Canis familiaris) and cows (Bos taurus). For a thorough discussion of taxonomic revisions and complete list of taxa, see the electronic supplementary material, table S1 and appendix S1.

First and last appearance dates were taken from the literature. If no reliable radiocarbon or historical dates were available, dates from Faunmap were used [23]. For analysis, species’ temporal ranges were binned every thousand years into presence/absence bins. To understand how potential future extinctions compare to previous losses, two ‘future’ bins were also included. In the first, six species listed as threatened (endangered or vulnerable) by the International Union for the Conservation of Nature [24] go extinct. In the second, four species listed at near threatened go extinct. These bins arbitrarily start at 1000 and 2000 years in the future, respectively, for graphing purposes only (electronic supplementary material, figure S1).

For functional traits, I recorded species averages for body mass; percentage diet of invertebrates, vertebrates, browse and graze; locomotor ability in running, climbing, swimming and digging; and social group size—all from a variety of literature sources and data types. All the literature sources, notes and discussions of trait values and uncertainty are given in the electronic supplementary material, table S1. For detailed trait-level descriptions, see the electronic supplementary material, table S2.

Body mass and social group size were log base 10 transformed before analysis. All traits were standardized by z-transformation so that each trait had a mean of 0 and standard deviation of 1. Traits were given an arbitrary but typically used weighting of one-third to mass, one-third to diet and one-third to locomotor traits and sociality [11,20,25]. Varying these weightings by a factor of two produced qualitatively similar results.

Pairwise functional distances were used to construct a functional dendrogram (cophenetic correlation = 0.74) using UPGMA clustering ( a). Petchey & Gaston’s [26] modified dendrogram-based functional richness (FD) was calculated as the total length of all branches needed to connect the species in a given time bin. FD was standardized to the length of the complete dendrogram including all 94 species (total functional richness).

Pairwise functional distances were also used to create a principle component space ( b) for calculating Villéger et al.’s [27] functional richness (FRic), Rao’s quadratic entropy (RaoQ), and functional dispersion (FDis). Only the first four principle components (90% variance explained) were kept for functional analysis as they were deemed non-trivial by the Random Average Under Permutation method, which has favourable performance in datasets with correlated variables [28]. FRic [27] is the volume of the multidimensional convex hull enclosing all species in a time bin and was standardized by the volume of the convex hull enclosing all 94 species in the analysis (total functional richness). RaoQ [29] is the average functional distance between two randomly selected individuals in a bin with replacement. FDis [30] is the average functional distance to the centroid in each bin.

FD and FRic measure functional richness, how much functional niche space is filled, and both are highly correlated with species richness [27,31]. For this reason, two different null scenarios were used to examine functional diversity indices compared to expectations given species richness alone. For functional diversity through time, taxonomic richness per bin was identical to the real scenario but the species going extinct between bins were completely random. Introductions were identical to the real scenario with humans and dogs, and later humans and European domestic animals, including dogs, entering North America as they would in the real scenario ( b–e). Another null model created to look at best- and worst-case extinction scenarios included only ‘native’ species, those 85 taxa that occurred in the earliest bin (49 999–49 000 BP), and species went extinct randomly until only 10 remained (electronic supplementary material, figures S2 and S3). Both null scenarios were repeated 10 000 times and medians and 95% confidence intervals were compared to actual extinction scenario values.

Functional distinctiveness ( ; electronic supplementary material, figures S4–S6) measures how much of the functional dendrogram is attributable to each species and was calculated identically to fair proportions evolutionary distinctiveness by summing the lengths of all branches necessary to connect a species to the root of a tree but with each branch length divided by the number of species ultimately subtending it [2]. Functional uniqueness measures the amount of unique function that will be lost if a species goes extinct. It was measured both as the length of the pendant edge of the functional tree for that species (electronic supplementary material, figures S7–S9) [2], and the decrease in convex hull volume (FRic) if that species were removed (electronic supplementary material, figures S10–S12). Functional distinctiveness and uniqueness were calculated both for native species only ( ; electronic supplementary material, figures S4, S7 and S10) and for the whole dataset (electronic supplementary material, figures S5, S6, S8, S9, S11 and S12) and were standardized to the total functional richness of each bin.

Open in a separate window

The functional similarity between each bin’s actual and null model communities (electronic supplementary material, S13) was calculated with a tree based measure of the Sørensen Index [32] adapted for dendrogram measures of functional diversity [11].

To generate a metric of extinction risk of native species to compare to functional distinctiveness, hierarchical partitioning was used to determine predictive variables for extinction in a logistic regression. Mass, climbing ability and digging ability, which have previously been linked to extinction risk (e.g. [33,34]), were identified as potentially important traits for predicting extinction in this dataset. The model including mass and climbing ability had the lowest AICc score and so was used to generate a continuous prediction of extinction risk (electronic supplementary material, figures S2–S4). The effect of including smaller mammals in the dataset was simulated by repeating the analysis above but using successively smaller mass thresholds for species inclusion. All analyses were carried out in R v. 3.2.2 [35]. For more detailed descriptions of data and methods, see the electronic supplementary material, appendix S1 and R script in electronic supplementary material, appendix S2.

3. Results

Total functional richness for North American large mammals dropped substantially with the megafaunal extinction ( and b,c), but this drop was mostly within the expectations of random species loss. Actual FRic loss was completely within the bounds of the null model but FD dropped just below the 95% CI of expected values in the 9999–9000 BP bin after the main extinction pulse and then again after the extinction (or pseudo extinction, see the electronic supplementary material, appendix S1) of the ancient bison (Bison antiquus) from 3999–1000 BP. However, RaoQ and FDis, which measure functional divergence rather than functional richness like FD and FRic, both dropped below expected values ( d,e) signifying that species were packed more closely to each other and the centroid of trait space, respectively, than before the extinction. Extinct and extant species showed no major differences in functional distribution measures (effect traits) such as distance to centroid (electronic supplementary material, figure S14), functional distinctiveness ( a; electronic supplementary material, figures S5 and S6), functional uniqueness (electronic supplementary material, figures S7–S12) and distance to the nearest functional neighbour (electronic supplementary material, figure S15). The effect trait functional distinctiveness did not correlate with the response trait extinction risk whether risk was calculated by a model (electronic supplementary material, figure S4) or based on actual extinction dates ( b) causing real community disassembly to fall within the bounds of random extinction, well above results predicted by positive effect–response trait correlation (electronic supplementary material, figures S2 and S3). Simulations of functional space that included all small mammals showed a different pattern with actual extinctions falling outside the null model (see the electronic supplementary material, figures S27–S30, and appendix S1 for more detailed results).

Table 1.

bintaxonomic richnessrelative taxonomic richness (per cent)relative FD (per cent)relative FRic (per cent)Late Pleistocene (49 999–49 000 BP)94100100100Pre-Columbian (1999–1000 BP)46495762Modern (999–0 BP)52556265extinction of threatened species42455232Open in a separate window

4. Discussion

Despite the size selectivity of the megafaunal extinction (electronic supplementary material, figure S16), the amount of function lost was indistinguishable from (FRic, c) or very close to (FD, b) values predicted from random species extinctions because functional distinctiveness (effect trait) did not correlate with extinction risk (response trait; ; electronic supplementary material, figure S4). Classic megafauna like mammoths did have outsized contributions to functional diversity (functional distinctiveness = 1.06–2.30%) but many much smaller species at little risk of extinction like white-nosed coatis (Nasua narica) have unique trait combinations that make them even more functionally distinct (2.16%) than most megafauna. On average, extant and extinct species showed no difference in the amount of functional diversity they contributed to the North American megafauna.

This matches results from Pacific Islands, where despite trait selective extinctions after human arrival, functional richness loss in bird communities was no greater than expected given random species removal [20]. If these results are generalizable, community disassembly can be accurately modelled as a stochastic process, where contingency has little value and species have equal, albeit partially overlapping, functional worth. There are several reasons to interpret these results cautiously though.

The continental spatial scale of this study necessitated by the scant fossil record of some species is regrettably coarse and could create artificial communities by lumping together taxa that would never have interacted in real life, potentially inflating functional redundancy and moving the ‘real’ community functional richness loss away from a worst-case scenario where response and effect traits are positively correlated (electronic supplementary material, figures S2 and S3). However, Boyer & Jetz [20] found similar results with fossils despite using the much finer spatial grain of individual Pacific islands. Defining community boundaries is also not entirely straightforward [36], especially given that many of the highly mobile vertebrates in this dataset (e.g. cougars, Puma concolor) have ranges that cover large portions of North America [37] or had associations with other species during the latest Quaternary that are not analogous to recognized communities today [38]. Limited field experiments also reveal that larger spatial scales may actually capture factors important to community assembly like competition and habitat filtering that were previously thought to operate at finer spatial resolutions [39].

If measures of effect traits like functional distinctiveness are generated from multiple traits that are uncorrelated with traits like mass used to construct the predictive extinction model, there is a mathematical inevitability of low response–effect correlation as each functional trait added to the analysis will only increase noise (electronic supplementary material, figure S4). But even using actual extinction dates, which are methodologically independent of functional distinctiveness, as response traits, response–effect correlation is absent ( ). This can be a real signal of multiple, opposing niche-based processes cancelling each other out, which is the reason some are suggesting more nuanced analyses where assembly/disassembly models are tested against individual traits [25,39]. However, I would argue that a multiple trait approach better reflects our conceptualization of niche space as an N-dimensional hypervolume [40,41] and the amount of biodiversity that is necessary to conserve ecosystem function [42]. Extinction does not just remove one trait or effect but the sum of a species’ impact on ecosystem function.

Another reason there is no effect–response correlation is that large extinctions may not follow the same rules that guide community assembly or background extinction [43]. Two species should not be able to exist in the same niche space [44,45], especially during times of high stress, so it is reasonable to expect that species with close functional neighbours would be more extinction prone; but this does not bear out as measured by functional distinctiveness ( a), uniqueness (electronic supplementary material, figures S7–S12), or nearest-neighbour distance (electronic supplementary material, figure S15). Species with close functional neighbours may actually be less likely to go extinct (electronic supplementary material, figure S15). This might be because species are avoiding niche overlap by separating on another unmeasured axis like spatial range or activity period. For example, the Canadian lynx (Lynx canadensis), bobcat (Lynx rufus) and jaguarundi (Herpailurus yagouaroundi) are functionally indistinguishable by dendrogram measures ( a) but their modern ranges have minimal overlap, dividing up North America in clear latitudinal bands [37]. For the highly trait selective megafaunal disassembly, however, the proximity of a functional neighbour was probably less important than where in functional space both species were located. The modern fauna was created by chopping off an entire region of Pleistocene functional space ( b) but current spatial richness gradients seem mostly caused by how densely species are packed into a given niche space, not an expansion or contraction of that niche space [46]. This trait selectivity is also probably different to the mechanisms inferred in most previous disassembly studies, which relied on the inevitable richness–area relationship and gradual habitat fragmentation caused by factors like rising sea levels [16,47]. All things being equal, North America should have gained species during the last 50 000 years as land area lost to sea-level rise was more than offset by the retreat of the massive Cordilleran and Laurentide ice sheets.

The patterns found in this analysis might also be unique to megafauna. Although megafaunal effect and response traits were not correlated, simulations suggest that including smaller and smaller mammal species could eventually lead to a positive effect–response trait correlation (electronic supplementary material, figures S29 and S30). This is likely because, compared to megafauna, small mammals are more numerous and functionally redundant (electronic supplementary material, figure S30). They would lower community functional richness much less than megafauna would if they went extinct. But surprisingly, this expectation does not hold in Boyer & Jetz’s [20] analysis of Pacific island bird extinctions, despite their inclusion of small birds. Even though selective pressures were similar to the North American megafaunal extinction, they found no effect–response correlation when considering all known terrestrial birds [20] (a mass range of five orders of magnitude). Island birds may be a fundamentally different system than continental mammals causing different results but the limited functional data available for small mammals might also be insufficient to capture true functional diversity relationships (see extended discussion in the electronic supplementary material, appendix S1).

That functional richness loss for North American megafauna was no worse than expected might hearten some. It suggests that if disassembly follows similar patterns in the future, we will avoid a worst-case scenario where the most distinct species are lost first. But because functional richness loss was no greater than expected given species richness does not mean that functional loss was non-existent ( and ). And although the total volume of function lost was no different than null expectations, the region of functional space lost was not random. Without the size selectivity of the real extinction, random extinction scenarios would produce a modern North American megafauna about 70% similar to our current function (electronic supplementary material, figure S13). RaoQ and FDis were also lower than expected showing that the distribution of species within functional space was not random ( d,e). Although most trait changes are within the bounds of a random extinction, running, climbing and vertebrate meat consumption are currently higher and mass is much lower than expected because large herbivores that could not dig or climb were more likely to go extinct (electronic supplementary material, figures S17–S20). This non-random trait loss better explains the changes in ecosystem function seen after the Pleistocene extinction [22] than gross functional richness and contradicts suggestions that long-term disassembly may produce smaller communities with relatively similar functional group ratios to the initial pre-extinction state [16]. Functional richness, no matter how it is measured, is an abstract value that is heavily dependent on methodological choices and difficult to compare between studies [31,48]. Rather than focus on the total amount of functional space filled, it is more useful to study how communities will move through that functional space and how different trait distributions will change underlying ecosystem function.

5. Conclusion

During the disassembly of the large North American mammal fauna, species did not go extinct based on how much function they contributed to community trait space causing total functional richness to drop no more than expected given species loss. Because the extinction was trait biased though, which regions of functional space were lost were not random. If future disassembly follows the same course, species will be lost based on which functions they possess, not how much function they contribute to a community. The largest functional richness losses may be yet to come because functional redundancy is low for many taxa including vulnerable species like the giant anteater (Myrmecophaga tridactyla) that now has no functional equivalent (electronic supplementary material, figures S5–S12). Short-term experiments will likely under predict these losses because they overestimate redundancy [49], but long-term studies of the recent fossil record may provide some relief as they show that some species can evolve rapidly to fill in vacated niche space [50]. A long view that seeks to understand these past functional changes may be the best way to predict how future disassembly will impact community function [22].

Supplementary Material

Table S1:Click here to view.(203K, xls)

Supplementary Material

Table S3:Click here to view.(45K, xls)

Supplementary Material

Allrange:Click here to view.(20K, csv)

Supplementary Material

Davis Convex Hull:Click here to view.(3.1K, txt)

Supplementary Material

Davis Functional Diversity Random Simulation:Click here to view.(6.9K, txt)

Supplementary Material

Davis Functional Diversity Simple:Click here to view.(3.0K, txt)

Supplementary Material

Davis Functional Diversity Simulation Mass Cutoffs:Click here to view.(11K, txt)

Supplementary Material

Davis Functional Diversity:Click here to view.(6.8K, txt)

Supplementary Material

Davis Nearest Neighbor:Click here to view.(5.9K, txt)

Supplementary Material

Davis Occurrence Matrix:Click here to view.(1.0K, txt)

Supplementary Material

Davis Plot Parameters:Click here to view.(2.0K, txt)

Supplementary Material

Davis Timebinning:Click here to view.(1.7K, txt)

Supplementary Material

Davis Unique Branchlengths:Click here to view.(4.5K, txt)

Supplementary Material

Davis Functional Diversity Simulation:Click here to view.(11K, txt)

Supplementary Material

Davis Disassembly Main Code Open First:Click here to view.(91K, txt)

Supplementary Material

Small Extinct Mammal Functional Data:Click here to view.(13K, csv)

Supplementary Material

EltonTraits1.0 North American Mammals:Click here to view.(114K, csv)

Acknowledgements

I thank Marcel Cardillo and three anonymous reviewers for helpful comments on this manuscript.

Data accessibility

All data used for this study are freely available in the electronic supplementary material and via the Dryad Digital Repository at http://dx.doi.org/10.5061/dryad.1bj88 [51].

Competing interests

None declared.

Funding

Partial funding for this project came from the Yale Institute for Biospheric Studies, the Geological Society of America, The American Society of Mammalogists and a Smithsonian Institution Predoctoral Fellowship.