Diversification dynamics, species sorting, and changes in the functional diversity of marine benthic gastropods during the Pliocene-Quaternary at temperate western South America

Functional diversity based on species traits is a powerful tool to investigate how changes in species richness and composition affect ecosystem functioning. However, studies aimed at understanding changes in functional diversity over large temporal and spatial scales are still scant. Here we evaluate the combined effect of diversification and species sorting on functional diversity of fossil marine gastropods during the Pliocene-Quaternary transition in the Pacific coast of South America. We analyzed a total of 172 species in 29 Pliocene and 97 Quaternary sites. Each species was characterized according to six functional traits: body size, feeding type, mobility, attachment, life-habit, and larval mode. Functional diversity was estimated according to four indexes (functional richness, evenness, divergence and dispersion) based on functional traits measured. Extrapolated species richness showed a slight yet not significant decrease from the Pliocene to the Quaternary despite the fact that a large faunal turnover took place; furthermore, a large extinction of Pliocene species (61–76%) was followed by a high pulse of appearances (49–56%) during the Quaternary. Three out of four indices of functional diversity (evenness, divergence and dispersion) increased significantly towards the Quaternary which is more than expected under a random turnover of species. The increase in functional diversity is associated with a loss of large-sized carnivore forms, which tended to be replaced by small-sized grazers. Hence, this trait-selective species turnover, even in the absence of significant changes in species richness, likely had a large effect and has shaped the functional diversity of present-day assemblages.


Introduction
The current threats to global biodiversity have increased the need to understand the role of species diversity on ecosystem functioning [1][2][3][4]. While the importance of species richness on PLOS  ecosystem functioning has been demonstrated experimentally [3,5,6], its application across very large scales is impractical, limiting the generality of conclusions. A simple yet powerful approach to understand ecosystem functioning is through the concept of functional diversity (FD hereafter). FD measures "functional trait diversity, where functional traits are components of an organism's phenotype that influence ecosystem level processes" [7]. A large body of literature has focused on theoretical or methodological aspects of FD analyses [7][8][9][10][11].
Unlike spatial patterns of variation in FD that have been increasingly studied in multiple taxa and systems [12][13][14], patterns of temporal variation are much less understood or have been documented over shorter timescales at small spatial scales [15]. The fossil record provides an excellent arena for the study of temporal variation in FD at much larger timescales and over larger spatial scales [15,16]. Temporal changes in FD at evolutionary timescales can occur due to shifts in species richness (due to differential origination/extinction/dispersal rates) or species composition (due to species sorting). Species sorting (or species selection) refers to the preferential persistence of species according to the interaction between their intrinsic traits and the environment [17]. Therefore, non-random extinctions or originations may have important effects on FD, especially if the redundancy of the system is low [4,8,18]. Vast volumes of paleontological studies have documented the existence of species sorting in marine taxa [17,19]. The selectivity of other traits seems sensitive to extinction intensity; for instance, a long list of different factors-including local abundance, reproductive mode, body size, trophic strategy, and life habit-were not selective at genus level during the K/Pg mass extinction [20,21]. However, many of these factors have shown high selectivity during background extinction times or large extinctions at regional scales [22][23][24]. Whether the diversification trends and species sorting is translated into changes in FD remains much less studied [15,16].
The Pliocene-Quaternary mollusc turnover provides an interesting temporal frame to evaluate these ideas. This faunal turnover, characterized by enhanced extinction and origination rates, affected different regions around the globe with variable intensity [22,[25][26][27][28][29] and was associated with major climatic and oceanographic generalized re-organization (e.g. cooling and changes in marine productivity) [30]. The impact of these transformations is well illustrated along the temperate coast of western South America where for more than a century a large Plio-Quaternary marine faunal turnover has been described [31] which altered the standing species richness and biogeographic patterns of several taxa [22,32,33]. If the faunal turnover was also selective for the ecological and life-history traits of these species, as shown in bivalves [22] and vertebrates [33], then an evolutionary change in the FD should be expected -a possibility that remains untested. Here we test these ideas, using benthic gastropods, a diverse and well-documented group with an excellent fossil record in the region [32,34], as a study model. Our goals were: a) to quantify the magnitude of the Pliocene-Quaternary turnover in marine benthic gastropods along western South America, and b) to determine the impact of species sorting during the Pliocene-Quaternary turnover on the functional diversity of gastropod assemblages.

Database
We gathered a database documenting the occurrence of 172 gastropod species from 126 fossiliferous sites (29 Pliocene, 97 Quaternary) along the temperate Pacific coast of South America from northern Peru to southern Chile. Although the paleoenvironments have been little studied, the environments represent shallow-water habitats (i.e., coastal shelf) in wave-exposed and wave-protected areas [34][35][36]. A total of 1,553 occurrences were compiled (see S1 Table). The information was collected from our own field studies and complemented with an exhaustive literature search [32,. For each species, we gathered information about the body size (maximum shell length), locomotion (actively mobile, facultative mobile, slow-moving), attachment (attached, non-attached), life-habit (epifaunal, infaunal), feeding type (carnivore, grazer, suspension feeder), and larval mode (aplanktonic, planktonic feeding, planktonic non-feeding). Maximum shell length and larval mode were obtained from our own collections and observations and complemented with information from literature (see above and Supplemental information). Information for locomotion, attachment, life-habit and feeding type was gathered from Paleobiology Database (https://paleobiodb.org/), and its use has been validated by previous studies [59][60][61]. If information was not available for a given species, traits were inferred or generalized from species belonging to the same genus or family, given that ecological and life-history traits are often phylogenetically conserved [12,62]. The relationship among these ecological and life history traits is not orthogonal (i.e. not all the combinations of these factors are observed/possible) because some traits are correlated in part. For instance, carnivores are typically mobile forms and tend to be large-sized. However, this possible multicollinearity effect is properly addressed by functional diversity analyses (see below). The database of species occurrences and functional traits used in our analyses are provided in the supporting information (S1 and S2 Tables).

Analyses
We evaluated the intensity of the Pliocene-Quaternary turnover, quantifying the standing species richness and the intensity of extinction and appearance of new forms. To evaluate the standing species richness in Pliocene and Quaternary assemblages we used the non-parametric extrapolation Chao 2 index [63] based on species-sites presence-absence matrices. We also estimated the species richness at a local scale using a subset of sites for which abundance (number of individuals per species) was available, using the Chao 1 index [64]. Analyses were carried out in the library vegan in R [65]. Species extinction was estimated as the percentage of Pliocene species not found during the Quaternary; conversely, species appearance was estimated as the percentage of Quaternary species making a first appearance during the Quaternary. These estimates, however, assume that sampling effort was similar between Pliocene and Quaternary, which clearly is not the case. In order to account for this possible bias, we re-estimated the Chao 2 index, extinction, and appearance using similar sampling effort for Pliocene and Quaternary assemblages (n = 29 sites in 10,000 runs).
Temporal changes in the prevalence of the six functional traits were evaluated using t-tests and χ 2 tests. We estimated the importance of taxonomic selectivity for both extinction and apparition of taxa by comparing the observed proportion of families that went extinct and appeared with the expectations of a random null model which assumed that the faunal turnover was independent of the taxonomic identity of the species (10,000 runs).
All six functional traits were used to evaluate changes in the functional diversity and were characterized through four complementary metrics: a) functional richness (FRic), b) functional evenness (FEve), c) functional divergence (FDiv), and d) functional dispersion (FDis). These four metrics emphasize different aspects of the functional diversity of assemblages [15,66]. FRic, FEve and FDiv are based on principal coordinates analyses which are based on a dissimilarity matrix built on Gower's index [67]. FRic is measured as the convex hull volume of the multivariate trait space, where higher values indicate a larger functional space. FEve is measured as the minimum spanning tree linking all the species within the trait space, where high values of FEve imply a more regular filling of the functional space [67]. FDiv measures the species deviance from the centroid of the multivariate trait space, and it is higher if species are far from the centroid, implying the existence of more extreme trait values in the functional space [67]. FDis is the average distance to the centroid of the multivariate trait space, where higher values indicate a higher dispersion of species over the trait space [66]. Importantly, FD indices are not affected by the possible multi-collinearity among functional traits; indeed, the principal coordinate analysis transforms a number of correlated variables into a reduced subset of uncorrelated variables (PCO axes). We pooled species across all sites in Pliocene and Quaternary assemblages and estimated FD indices for each time bin. We estimated the change in FD according to each metric for: i) the Pliocene extinction (Pliocene FD-extinct species FD), ii) Quaternary appearance (appearing species FD-surviving species FD), and iii) between the Quaternary and Pliocene. We tested whether the change in each FD metric was different than expected under a null model of random species turnover by re-shuffling the species identity in the matrix of traits (10,000 runs). Analyses were carried out using the library FD [66] and scripts written in R [65].
The prevalence of different ecological traits was variable in time (Fig 2). There was a significant reduction in the body size of species from the Pliocene (median = 40 mm) to the Quaternary (median = 29 mm) (t-test on log-transformed data, P = 0.007, Fig 2A). The prevalence of feeding types also experienced changes, particularly a reduction in the proportion of carnivores (from 65% to 46%) and an increase in the herbivore grazers (from 25% to 41%) (Fig 2E). For other traits, no significant changes were detected (χ 2 test, P > 0.05, Fig 2).
Neither the extinction nor the apparition of species was taxonomically selective, as the proportion of extinct/appearing families was not different than expected by the null model (Fig 3). Patterns were, in general, similar within the most speciose families (Table 1). For instance, in Muricidae-the most speciose family during both Pliocene and Quaternary-the observed proportion of extinct species (0.64, 23 out of 36) was not different than expected by chance (95% CI: 0.41-0.74). However, the proportion of appearing Muricidae during the Quaternary (0.38, 8 out of 21) was lower than expected by the null model (95% CI: 0.48-0.79).
The four metrics of functional diversity, depicted by the two first axes of a principal coordinate analysis, showed an increase from the Pliocene to the Quaternary (Fig 4), revealing changes in functional evenness, divergence and disparity, which was significantly higher than expected under a random null model of a trait-independent diversification (P < 0.05, Table 2). The Pliocene extinction produced a significant increase in all indices, except in the case of functional richness ( Table 2). The Quaternary appearance of species also increased the functional divergence (P = 0.03), but the change for other indices were non-significant (Table 2).

Discussion
Our results show the FD of benthic gastropod assemblages along the western South America, indicating important changes during the Pliocene-Quaternary transition coupled with a large species turnover. Despite this species richness remaining relatively similar across time, the trait-selective replacement of a large number of taxa had a marked impact on the FD of gastropods. The existence of a large faunal turnover during the Pliocene-Quaternary transition at western South America has been hypothesized for more than a century [31,68]. Our results show the existence of a massive removal of species-confirming preliminary evaluations [51,53,56] -followed by the appearance of many new forms. These trends seem robust to differences in sampling intensity between the Pliocene and Quaternary and to the completeness of fossil inventories. Similarly, higher levels of biotic turnover have been estimated for bivalves and marine vertebrates [33]. The ultimate causes of the turnover are likely related to the onset of the modern conditions of the Humboldt Current and include several mutually non-exclusive factors, including a drastic decrease in sea temperature [69,70], increasing primary productivity [71], destruction of wave-sheltered coastal habitats due to tectonic subsidence [34,72], increasing anoxia due to increasing primary productivity [22], and habitat destruction and creation implied by the advance/retreat of glacial lobes in southern Chile during glacial/interglacial cycles [32]. Thus, the same factors promoting the disappearance of many forms also acted as drivers for the appearance of new species. The rapid appearance of new forms during the Quaternary has also been inferred based on phylogenetic analyses [73].
Evidence from the fossil record has shown drastic changes in FD through events of major biotic turnover [74], but patterns are hard to compare given the many methodological approaches used to assess patterns of functional diversity [75,76]. Different FD trends can be inferred according to different methodological approaches [16, 77,78]. Indeed, our analyses reveal that changes in FD can be inferred based on three out of four metrics. The functional richness index, sensitive to changes in species richness [67], increased slightly towards the Quaternary; however, the overall trait space occupied by species remains relatively similar, as revealed by the lack of significance in changes in functional richness. Indeed, no functional trait/group disappeared/appeared during the faunal replacement-a trend seen even during global mass extinction events [16]. However, the increase in functional evenness, divergence, and disparity indicates that the functional trait space is more homogenously filled and more dispersed towards the Quaternary. Berke et al. [12] have shown that functional evenness in marine bivalves increases towards higher latitudes which is attributed to an evolutionary accumulation of taxa due to biogeographic dynamics. The paleoenvironmental shift from subtropical to temperate conditions and the loss of many molluscs of tropical affinity during the late Neogene [79] are consistent with this pattern.     A strong species sorting during the faunal turnover was responsible for the increase in FD over time. Both extinction and appearance of forms were selective according to ecological and life-history traits, as commonly seen in the fossil record [17,19], though not driven by changes in the taxonomic composition of assemblages. Changes in FD can be attributed to the preferential loss of large-sized carnivores and the subsequent increase in the proportion of smallsized herbivore grazers and suspension-feeders. On the one hand, the loss of carnivory could be explained by two non-exclusive processes: i) the elevated losses during the Pliocene extinction, and ii) the reduced appearance of new forms during the Quaternary. Five carnivore families went extinct during the Pliocene, compared to only three families appearing in the Quaternary (Fig 5). In parallel, the appearance of new species within Muricidae, the most speciose family analyzed, was lower than expected by chance. On the other hand, the increasing importance in herbivory could be attributed to a null extinction of families, the diversification of surviving forms (e.g. limpets and turbid snails), and to the appearance of seven new families (Fig 5). Interestingly, the overall reduction in median body size did not occur within surviving families; i.e. this was remarkably stable across the Pliocene-Quaternary transition, and it is rather attributed to the appearance of small-sized herbivore families (Fig 5). Indeed, in agreement with previous studies [56,80] a few surviving families (e.g. Volutidae, Muricidae, and Fissurellidae] show an increase in median body size. However, this is not consistent with the so-called "Lilliput Effect", i.e. a miniaturization of species body size immediately after a large extinction event [81,82], because the decline in body size is largely due to newly originated species.
Previous studies have shown that the late Neogene extinction of bivalves and marine vertebrates along western South America were also trait-selective [22,33]. This combined evidence suggests that the deep environmental transformations that took place during the Pliocene-Quaternary transition [71,83] not only altered patterns of species diversity and composition but also might have brought a major shift in the functioning of marine ecosystems. Our results highlight the usefulness of FD indices to provide more detailed understanding of patterns of biotic turnover across evolutionary timescales.
Supporting information S1