Selective extinction against redundant species buffers functional diversity

(a) Diversification dynamics

We assessed the diversification dynamics of Caribbean molluscs by quantifying speciation and extinction rates, net diversification, and species richness during the last 23 Myr. To do so, we applied a Bayesian approach that accounts for sampling and preservation biases and dating uncertainties (see Material and methods and electronic supplementary material) to fossil occurrences identified at the species level. We found eight (95% credible interval, CI: 4–15) significant rate shifts in speciation rates and seven significant shifts (95% CI: 4–12) in extinction rates, which are punctuated across the evolutionary history of Caribbean molluscs (electronic supplementary material, figure S2). After a gradual decrease during the early Miocene, speciation rates were generally constant during the middle and late Miocene. Then, there was an increase during the early Pliocene, peaking around the early Pliocene/late Pliocene boundary (approx. 3.7 Ma), and sharply decreasing thereafter (figure 1a). Extinction rates were roughly constant and lower than speciation during the Miocene (figure 1b). During the Pliocene and Pleistocene (from approximately 5 Ma onwards), there was a progressive increase in extinction rates, reaching maximum values around the Pliocene/Pleistocene boundary, approximately 2.6 Ma (figure 1b). Extinction rates remained consistently high throughout the Pleistocene (figure 1b). Net diversification rates were positive throughout the Miocene and until the early Pliocene (figure 1c), followed by a negative phase from the late Pliocene to the Pleistocene, between 3.6 and 2.9 Ma (grey band in figure 1). Finally, there was an almost continuous increase in species richness during the Neogene, with maximum diversity reached in the late Pliocene and resulting in a 49% net loss from the late Pliocene to the Pleistocene (figure 1d).

• Download figure
• Open in new tab
• Download PowerPoint

Our results indicate that molluscan communities experienced major diversity changes during the last 23 Myr, which were characterized by two phases: (i) an increase in diversity during most of the Neogene (between 23 and 5 Ma, peaking around 3.7 Ma) and (ii) an extinction phase at end of the Neogene and throughout the Quaternary (from 5 Ma onwards, peaking around 2.6 Ma). During the first phase, diversity gradually increased over an approximately 18-million-year period of constantly low speciation and extinction rates (figure 1a,b). This phase culminated in the early Pliocene when speciation rates peaked, allowing species richness to reach maximum values (figure 1d). The diversity increase phase was coeval with the closing of the CAS [28–34], which gradually provided additional substrate, thereby increasing habitat complexity [38,39].

The phase of positive diversification was followed by a 5-million-year phase of high extinction (figure 1b) and negative net diversification (grey band; figure 1c), which resulted in a net loss of 49% of species. This phase roughly falls within the time range of elevated extinction rates previously reported during the formation of the Isthmus of Panama (figure 3d in [24]). However, our results contrast in terms of when the peak took place. While we found the extinction rates to be highest starting approximately 2.6 Ma, O'Dea et al. [24] found extinction to peak between 2 and 1 Ma. Contrasting results can be due to differences in geographical range covered (Panama versus Caribbean), taxonomic resolution (genus versus species), and analytical methods used to estimate extinction rates (number of last occurrences divided by the total number of taxa present per million years versus Bayesian inference that accounts for sampling and preservation biases; see Material and methods and electronic supplementary material). Of particular relevance is the difference in geographic range, as the use of the Caribbean area as a whole allowed us to sample greater habitat heterogeneity during the entire time series.

The time of the previously proposed extinction peak was interpreted as a lagged response to the environmental changes attributed to the closure of the CAS, between 2.25 and 3.45 Ma [24]. Based on our results and on the new evidence suggesting that the CAS constricted earlier, during the late Miocene (approx. 10 Ma) [28–34], we offer an alternative interpretation. We suggest that the high extinction rates of Caribbean molluscs at the end of the Neogene were likely a response to the environmental changes associated with the Northern Hemisphere glaciation, which started in the late Pliocene and persisted throughout the Pleistocene [40–43]. These climatic changes caused an abrupt transition to a new climatic state dominated by colder and more variable global temperatures [40–43]. Importantly, these climate changes produced large sea-level oscillations [40,42] and a significant reduction in the extent of sublittoral areas available [7], decreasing habitat for benthic communities in the Caribbean.

In sum, based on our results and recent evidence concerning the timing of the closure of the CAS, we interpret the positive diversification that took place during most of the Neogene as related to the additional substrate in the Caribbean Sea as the CAS was closing [38,39]. By contrast, the extinction of nearly half of Caribbean molluscs by the Quaternary was likely related to the unstable and reduced habitats that those species faced as a consequence of the Northern Hemisphere glaciation from the late Pliocene onwards [7,40,42].

(b) Functional diversity

In order to understand the ecological effects of the diversification of Caribbean molluscs, we quantified changes in functional diversity through time. We did so by performing trait-based analyses following the methods described in Mouillot et al. [1] and as applied to the fossil record in Pimiento et al. [7]. Accordingly, we assigned five ecological traits that characterize molluscan functional diversity (i.e. body size, life habit, locomotion, environment, and diet; see electronic supplementary material, table S1 and methods) to (i) quantify the number of functional entities (i.e. unique trait combinations; hereafter FEs) and (ii) construct a trait space per time bin (see Material and methods). We found a general increase from the Miocene to the early Pliocene in the number of FEs (electronic supplementary material, figure S3; 58% increase); the mean number of species per FE, hereafter functional redundancy (figure 1e and electronic supplementary material, figure S4; 68% increase); the over-representation of some FEs (see Material and methods), hereafter over-redundancy (electronic supplementary material, figure S3; 18% increase); and in the convex hull volume of the trait space (see Material and methods), hereafter functional richness (figure 1f; 31% increase). All of these metrics peak between the early Pliocene and late Pliocene (between 4 and 3 Ma). While functional redundancy sharply decreases afterwards (by 31%), the number of FEs, over-redundancy, and functional richness remain virtually unchanged (9%, 3%, and 3.5% decrease, respectively; figure 2e,f and electronic supplementary material, figure S2). We further compared these results against null expectations (see Material and methods) and found that all functional diversity values were no different than expected for most of the Neogene (electronic supplementary material, figure S3; electronic supplementary material, table S2). Conversely, from the late Miocene to the Quaternary, we found that: (i) functional redundancy values were higher than expected in the late Miocene, and mostly as expected thereafter; (ii) over-redundancy was mostly as expected, and (iii) the number of FEs and functional richness values were generally significantly higher than expected (electronic supplementary material, figure S3; and table S2).

• Download figure
• Open in new tab
• Download PowerPoint

Our analyses suggest that Caribbean molluscan assemblages experienced an overall increase in functional diversity during the Miocene and Pliocene. Such an increase of functional diversity was characterized by the appearance of new functions (as evidenced by the increased number of FEs), as opposed to the over-representation of pre-existing ones (as evidenced by a lower increase of functional over-redundancy during the same time period; electronic supplementary material, figure S3). The general increase in functional diversity during most of the Neogene was coupled with the growth of species richness that took place during the same time (figure 1). Indeed, our simulations indicate functional diversity values were as expected given the number of species present in the assemblage (electronic supplementary material, figure S3; and table S2). Nevertheless, from the late Miocene onwards, when extinction rates started to increase and species richness started to decline, the number of FEs and functional richness were largely maintained and therefore, decoupled from species richness (figure 1f; electronic supplementary material, figure S3). As such, at the time of the highest extinction rates and consequent species loss, FEs and functional richness are significantly higher than expected (electronic supplementary material, table S2).

Our results indicate that the extinction phase of Caribbean molluscs, which started in the late Miocene (approx. 5 Ma) and peaked approximately 2.6 Ma, had a disproportionally small effect on functional diversity (figure 1f; electronic supplementary material, figure S3; and table S2). These results are consistent with previous studies documenting negligible changes in global functional diversity in the face of mass extinction events [3,4,6,44]. The combined evidence suggests that the functional resilience of benthic assemblages to extinction is consistent across spatial scales (regional to global) and extinction events (moderate to large, i.e. ‘mass extinctions'). Further, the decoupling of species loss and functional diversity suggests that extinction was non-random.

(c) Extinction selectivity

We tested if extinction followed diversity-dependent dynamics (i.e. if there was extinction selectivity) by developing a novel multi-trait-dependent extinction model that we introduce here (see Material and methods). Briefly, we tested the effect of functional redundancy (i.e. number of species per functional entity) and ecological traits (i.e. body size, life habit, locomotion, environment, and diet; electronic supplementary material, table S1) on overall extinction rates. We found strong evidence that extinction did not affect all species equally (figure 2). Functional redundancy and body size had the highest effect on extinction rates (strong statistical support; log Bayes factor greater than 10; red line in figure 2) followed by environment (i.e. water depth; positive statistical support; log Bayes factor greater than 6; orange line in figure 2a). Accordingly, very high redundancy (species sharing their FEs with more than 32 species; figure 2b; also see Material and methods) had the most significant effect on extinctions (their entire range of relative effect on extinction was higher than expected; figure 2b). High redundancy (species sharing their FEs with 10 to 32 species; figure 2b) also had an important, albeit relatively weaker effect on extinctions. Similarly, very small and small size (less than 10 mm; figure 2d) were associated with higher extinction rates. Finally, shore-face environments (less than 5 m of depth; figure 2c) were less affected by extinction than all other environments (figure 2c,d). None of the other traits—such as diet, life habitat, and locomotion—had a significant effect on extinction rates (figure 2a).

Our multi-trait dependence models indicate that redundancy and size had a much stronger effect on extinctions than all other ecological traits, with small, highly redundant species being significantly more affected by extinctions than those that were large or functionally unique (figure 2). This result, combined with our finding that during the period of elevated extinction rates (from 5 Ma onwards; figure 1) Caribbean molluscs maintained significantly higher functional richness than expected by chance (electronic supplementary material, figure S3), suggests that selectivity against redundant species conferred functional resilience to extinction across assemblages. This is because the loss of redundant species likely allowed the spread of species among functional groups, preventing increases in over-redundancy, generally stabilizing the community, and ultimately enabling new FEs to persist. Indeed, extinction selectivity against rich FEs, or asymmetric functional selectivity, had been previously hypothesized as a possible explanation for the disproportional persistence of FEs, even during mass extinction events [6]. Similarly, the extinction susceptibility of highly redundant species suggests that competition may have played a major role in driving extinctions. This could be due to the fact that highly redundant species (i.e. those sharing traits and therefore, ecological functions) would have been competing more for resources than unique species, especially during times of large sea-level oscillations and habitat loss. Accordingly, we propose that biotic interactions could have mediated the response of species to the environmental changes that took place during times of high extinction rates.

We also found that body size had an important effect on extinctions (figure 2), with small-sized molluscs being more susceptible than medium- and large-sized species (figure 2d). The association between body size and extinction risk has been widely studied in the literature [e.g. 45–51]. In general, large size has been associated with historical extinctions and current risk (but see [46]), whereas large- and small-sized organisms have been equally at risk in the geological past [45]. In the fossil record of marine invertebrates, the effect of body size on extinctions has been found to vary across extinction events [45]. Our finding of selectivity against small-sized Caribbean molluscs reinforces our hypothesis of competitive interactions having an important role on extinction (figure 2) because smaller organisms tend to be more redundant [52,53] and therefore, face higher competition than larger organisms.

Consistent with our current knowledge on the habitat loss and sea-level oscillations that took place as Caribbean molluscs entered the extinction phase [7,20,40,42], our models found statistical support for an environmental effect on extinction rates (figure 2a). Within the neritic zone (i.e. above 200 m of depth, the zone in which all our samples lived; also see electronic supplementary material, dataset), species living in shore-face environments (less than 5 m of depth; also called intertidal; electronic supplementary material, table S1) had the lowest extinction risk. We suggest that these species were better equipped to tolerate the dramatic eustatic changes that took place from the late Pliocene onwards, as they would constantly face large sea-level changes during their lifetime [54] in comparison to species living in other environments (e.g. [55]). Although we found that environment, as defined here in terms of water depth, was an important factor in the extinction of Caribbean molluscs, our results show that redundancy and body size were the major determinants of extinction risk.

Diet, life style, and locomotion have been proposed as drivers of the extinction of molluscs at the end of the Neogene in the Caribbean, with infaunal suspension and predatory feeders being more affected than attached epifaunal deposit feeders [20,21,23,27]. These results led previous authors to hypothesize that changes in benthic diversity at the time were caused by a decrease in productivity (as indicated by an increase of suspension feeder abundance) and an increase in the development of seagrass meadows and coral reefs (as indicated by an increase in abundance of attached deposit feeders) [20,21,23,27]. By contrast, we found no effect of diet, life style, or locomotion on extinction rates (figure 2). Nevertheless, our finding that species living in environments between 0 and 5 m of depth had the lowest extinction risk (figure 2c) is consistent with the hypothesis of seagrasses enabling the survival of Caribbean molluscs. It is important to note that although statistically significant, the effect of environment (i.e. water depth) on extinction rates in Caribbean molluscs during the last 23 Myr shows no consistent trend (figure 2c) and should therefore be interpreted with caution. We further suggest the contrasting results presented here relative to previous studies could be due to the fact that we used different traits and trait resolution [20,21,23,27].

(d) Conclusions

Our study of Caribbean molluscs over the last 23 Ma reveals two distinct phases in terms of diversification and functional diversity: an approximately 18-million-year long phase of coupled taxonomic and functional diversity increase that resulted in maximum levels of functional redundancy; followed by an extinction phase that started in the late Miocene (approx. 5 Ma), peaked around 2.6 Ma, and resulted in negative net diversification (starting between 3.6 and 2.9 Ma) and the loss of 49% of species. Importantly, the extinction phase had negligible impacts on functional richness (figure 1). Further, our novel multi-trait dependence model identified selective pressure against small, functionally redundant species (figure 2). Based on our results we propose that (i) the response of species to environmental change at the end of the Neogene was determined by species competition within redundant FEs and that (ii) non-redundant species were generally released from competition, enabling the maintenance of functional richness in molluscan communities during times of extinction and thereby conferring ecological resilience.

(a) Dataset

We compiled a dataset consisting of 34 011 fossil occurrences, representing 960 genera and 4422 species of Caribbean molluscs ranging from the Neogene to the Pleistocene/Holocene boundary (23.03–0.01 Ma; electronic supplementary material, dataset). The data was collected at the species level from the published literature and unpublished field observations available from the Paleobiology Database (https://paleobiodb.orgl) and museum-deposited species available from iDigBio (https://www.idigbio.org; both downloaded 23/12/2018). While each data portal maintains a taxonomic backbone, the nomenclature of these independent datasets was standardized using the World Register of Marine Species taxon match service for extant species and genera. This eliminated many synonyms, misspelling, and consistently enforced taxonomic ranks. Extinct genera and species were accommodated within this taxonomic framework, with care taken to eliminate spelling variations. We defined the study area of the Caribbean following the marine ecoregions of the world [56] (as available at www.worldwildlife.org). Accordingly, our dataset only includes fossils from the ‘Tropical Northwest Atlantic' province (between 5.7° N and 35.7° S). To identify relevant records, we used the speciesgeocodeR package v. 1.8 in R [57]. In a few cases ([11]; figure 1), we assigned sites that are located outside any of the provinces and are terrestrial today, to the geographically closest province, using the ‘dist2Line' function of the geosphere package v. 1.5–10 [58]. The age range of each occurrence was further refined using published references to place them into the finest stratigraphic resolution possible. To do so, we followed the latest geological timescale [59] and the estimate of 2.58 Ma for the Pliocene/Pleistocene boundary [60].

(b) Speciation and extinction rates

We estimated speciation and extinction rates of Caribbean molluscs through a Bayesian analysis of all Neogene and Quaternary fossil occurrences identified to species level, while accounting for incomplete sampling due to fossilization and preservation processes and dating uncertainties. The stratigraphic ranges assigned to each occurrence were treated as dating uncertainty [61]. We resampled the age of each occurrence randomly from a uniform distribution spanning its minimum and maximum age and repeated this procedure 10 times to generate 10 resampled datasets. All analyses described below were repeated on the 10 replicates and their results combined in a joint posterior sample to incorporate dating uncertainty in the results [61]. Three main sets of parameters were simultaneously estimated using PyRate v. 2 [61,62]: (i) preservation rates quantifying the expected number of fossil occurrences per sampled lineage per time unit (1 Myr), (i) speciation and extinction times for each species, which usually extend beyond first and last fossil occurrences, and (iii) speciation and extinction rates for each species (expected number of speciation and extinction events per lineage per Myr) and their temporal variation. Details on how preservation was modelled and on how rate shifts were calculated can be found in the electronic supplementary material.

(c) Functional diversity analyses

We selected five traits (i.e. body size, life habit, locomotion, environment, and diet; electronic supplementary material, table S1) that characterize the ecological roles of Caribbean molluscs [63–65]. Traits were assigned to species using authoritative taxon-specific texts, online databases, and expert assessments based on both extant relatives and the fossil record (electronic supplementary material, dataset). There were limited differences in the number of species occurrences per trait category (electronic supplementary material, figure S6). Based on the trait assignments, we calculated a series of functional indices in order to evaluate changes in functional composition and structure across the Neogene and Quaternary (23.03–0.01 Ma). To make this comparison, we used the longevity of each species estimated in PyRate (i.e. the time of origination and extinction of each species; see above) to record the presence and absence of species in 23 time bins of 1 Myr each (from 23–0 Ma). Accordingly, all functional diversity analyses described below were made per time bin. Given that the age of each species at a given time bin was assigned based on their longevity, no sampling differences between time bins were detected. Nevertheless, in order to incorporate age uncertainties, we resampled the estimated longevity of each species 100 times and performed 100 iterations of all functional diversity calculations. For each time bin, we evaluated: (i) the functional composition of Caribbean molluscs based on FEs or unique trait combinations (which reflects the diversity of ecological roles) and (ii) their functional structure based on their functional space (which reflects the range of species' functional strategies).

(d) Multi-trait-dependent extinction model

To assess whether extinction rates varied among lineages as a function of one or multiple discrete traits, we developed a novel multi-trait-dependent extinction (MTE) model, which was included in the PyRate software. We express the likelihood of the estimated longevity of species i, here indicated with t_i, based on the birth-death model described by Silvestro et al. [61] and assuming a lineage-specific extinction rate μ_i. Thus, the likelihood of an extinct lineage i is computed as

where t_i is the lineage duration until extinction, while the likelihood of an extant lineage j is computed aswhere t_j is the lineage duration until the present, or until the most recent boundary of the time window of interest.

In the MTE model, we parameterized the extinction rate as a function of a mean extinction rate (μ₀) and Dirichlet-distributed multipliers m. The multipliers express the trait-dependent component of the model, and quantify the deviation of extinction from the mean associated with the lineage-specific trait state v. The resulting rate for a lineage i in a set of S species is

Under this parameterization, μ₀ equals the mean extinction rate across all lineages in a given time window, i.e.

Here, we extended the model to allow for multiple traits to have an effect on extinction rates, by making the lineage-specific extinction rates a function of the combined trait states:

where V is the number of traits included in the model.

We implemented the MTE model in a Bayesian framework to estimate all the free parameters, namely the mean extinction rate (μ₀) and the Dirichlet-distributed multipliers (m). We used the Markov chain Monte Carlo (MCMC) algorithm to obtain posterior samples of all parameters.

We used an exponential prior on the mean extinction rate with rate parameter equal to 1. Efficient sampling of Dirichlet-distributed variables is typically quite problematic; therefore, we used here a re-parameterization to (i) improve MCMC convergence and (ii) facilitate the use of hyper-priors on the shape parameters of the Dirichlet distribution (α₁, …, α_N). We followed the parameterization proposed by Gelman et al. [69] and expressed the trait-dependent multipliers as

We fixed y₁ = 0 (as only N − 1 values are necessary to define N Dirichlet-distributed values) and used normal prior distributions on the y parameters with equal means (i.e. the prior expectation is that extinction rates are constant across lineages and independent of the trait) and with precision parameter τ:

When multiple traits are analysed (V > 1) each trait is assigned an independent precision parameter τ. The precision parameter(s) τ were assigned a gamma hyper-prior:

with fixed shape parameter α = 1.5, and rate β. The rate β was itself considered as a free parameter and assigned a vague exponential hyper-prior:

This parameterization improves substantially the efficiency of the posterior sampling because it allows us to sample several (hyper)parameters directly from their conjugate posterior distribution. For a given set of y values and β, we sample the precision parameters τ directly from a gamma conjugate distribution:

Likewise, posterior values of the rate parameter β of the gamma prior on τ are sampled directly from its gamma distribution for a given set of precision parameters:where V is the length of the vector τ (i.e. number of traits analysed).

The inclusion in the model of traits that have no effect on the extinction rates may result in over-parameterization that can yield spurious results and false positives. We addressed this issue in two ways. First, our prior and hyper-prior settings (equations (3.3)–(3.5)) result in the Bayesian shrinkage of the y parameters around 0, therefore shrinking the multipliers m toward equal values (i.e. no effects of a trait on extinction). Specifically, as β shrinks to smaller values (as determined by the exponential hyper-prior), the precision parameters τ tend to become larger (i.e. smaller variance), therefore constraining the y parameters to values close to 0. Second, we implemented a Bayesian variable selection algorithm where each trait v is assigned an indicator I_v which can take a value of 0 or 1 (see Silvestro et al. [17] for a similar implementation). The indicators are used as multipliers to the y vectors so that when I_v = 0 all resulting m are equal (equation (3.2)). The indicators are treated as unknown parameters and estimated from the data through MCMC. We set the prior probability to P(I = 1) = 0.05 (therefore placing 0.95 probability on the absence of an effect). Under this prior a posterior frequency of I = 1 greater than 51.4% is considered as significant (2 log Bayes factor greater than 6 [69]).

We ran 10 replicated analyses using the resampled longevities of species inferred in the PyRate analyses of speciation, extinction, and preservation described above. For each analysis we ran 2 million MCMC iterations sampling every 1000 iterations. After combining all the output files, we computed the posterior frequency of each indicator to assess which traits were estimated to have a significant effect on extinction. We summarized any trait-dependent effects based on the sampled multipliers m. Because extinction is measured as a rate, i.e. a per-species value, the estimated relative effects on extinction is independent of the number of species in a state.

The results from the MTE model did not appear to be affected by species occupancy (a proxy of rarity; see electronic supplementary material, figure S8).