Research ArticlePALEONTOLOGY

Little lasting impact of the Paleocene-Eocene Thermal Maximum on shallow marine molluscan faunas

See allHide authors and affiliations

Science Advances  05 Sep 2018:
Vol. 4, no. 9, eaat5528
DOI: 10.1126/sciadv.aat5528


Global warming, acidification, and oxygen stress at the Paleocene-Eocene Thermal Maximum (PETM) are associated with severe extinction in the deep sea and major biogeographic and ecologic changes in planktonic and terrestrial ecosystems, yet impacts on shallow marine macrofaunas are obscured by the incompleteness of shelf sections. We analyze mollusk assemblages bracketing (but not including) the PETM and find few notable lasting impacts on diversity, turnover, functional ecology, body size, or life history of important clades. Infaunal and chemosymbiotic taxa become more common, and body size and abundance drop in one clade, consistent with hypoxia-driven selection, but within-clade changes are not generalizable across taxa. While an unrecorded transient response is still possible, the long-term evolutionary impact is minimal. Adaptation to already-warm conditions and slow release of CO2 relative to the time scale of ocean mixing likely buffered the impact of PETM climate change on shelf faunas.


The Paleocene-Eocene Thermal Maximum (PETM) [~56 million years (Ma) ago] was an episode of global warmth brought about by the rapid [<5 thousand years (ka) (1)] release of up to 10,000 gigatons (2) of carbon into the atmosphere from a sedimentary (3) and/or volcanic (2) source. Global marine temperatures rose by 5° to 8°C (4), and acidification in the mixed layer (2, 5) and deep ocean (6) was severe and often coupled with hypoxia both at depth (7) and on continental margins (8). The PETM likely underestimates the expected impact of ongoing combustion of fossil fuels (2, 9) but nonetheless remains the closest analog to the present offered by the geological record. As such, the interval has received intense scrutiny for insights on how the present-day climate system and biosphere might react to projected anthropogenic increases in atmospheric CO2. The evolutionary, ecological, and biogeographic responses of terrestrial faunas and floras, planktonic ecosystems, and the deep-sea benthos have all been documented in some detail (10), yet very little is known about the impact of the PETM on marine shelf macrofaunal assemblages (11). The dearth of known macrofossiliferous shallow-marine sections, compounded by the low probability of preserving such a short-duration event in complex depositional settings rife with small-scale hiatuses (12), continues to hinder interpretations.

The U.S. Gulf Coastal Plain (GCP) contains a thick, highly fossiliferous, and much-studied section spanning the Paleogene and, as such, provides one of the best overall records of marine mollusk assemblage change in shallow-shelf environments through this interval in the world (13). The PETM has been located in the section based on dinoflagellate biostratigraphy and a diagnostic carbon isotope excursion (CIE) (14), and while marine macrofossils are not known from the CIE—and hence the interval of maximum warmth—itself, rich shell beds closely bracket the interval. Despite an inability to see the immediate short-term response of mollusks to PETM CO2 release, a comprehensive comparison of faunas before and after can reveal the long-term evolutionary impacts of this transient environmental perturbation. In the same way that other biotic perturbations leave a lasting impact on the history of a group, as evidenced by delayed recovery of richness, peaks in first and last appearances, lasting shifts in ecologic structure and composition, or persistent changes in body size, we should expect to see the fingerprint of the PETM on the evolutionary and ecological history of mollusks, even while an assessment of its immediate impact remains elusive. We therefore bring a range of tools to bear on the rich molluscan fossil record of the GCP to explore whether the PETM had a significant lasting effect on the richness, turnover, or ecological structure of shelf faunas or the body size, growth rate, or life span of component taxa.


We derive the results discussed below from the analysis of taxonomic occurrence and abundance data sets with taxa coded by ecological guild, and body size and growth data for species in two common molluscan clades. The geologic setting, details of each data set, and analytical approaches are described in Materials and Methods, and data and additional results are provided in the Supplementary Materials.

Taxonomic diversity and turnover

We analyze a synoptic data set of species-level taxonomic occurrences for bivalves, gastropods, and scaphopods from five discrete shell beds spanning the Paleocene-Eocene transition in the central GCP (fig. S1) and an abundance data set comprised of the subset of those collections with counts of individuals (data file S1). Raw richness data suggest a pulse of new taxa after the PETM (Fig. 1A). However, as Dockery (13) recognized early on, differences in sampling intensity might severely bias the perceptions of turnover in the GCP section. Resampling all units to equivalent numbers of occurrences shows this pulse to be largely an artifact of heavy sampling within the earliest Eocene Bashi Marl (BM) and supports only minor changes in richness through the interval (Fig. 1A). Sampling standardization and correction for differences in interval durations show that the proportion of taxa making their first or last appearance at the PETM is lower than that at the other horizon boundaries within the Wilcox Group (Fig. 1, B and C). The short duration of intervening time yields higher turnover rates between the Greggs Landing Marl (GLM) and Bells Landing Marl (BLM), both part of the same upper Paleocene formation, than elsewhere in the section. The PETM transition therefore does not stand out as unusual in the GCP with respect to changes in standardized richness or turnover of molluscan species.

Fig. 1 Richness and turnover patterns derived from resampling pooled taxon occurrences compared to raw values in shell beds bracketing the PETM (dashed line).

The number of collections in each unit is provided beneath the unit name. (A) Standardized mean species richness for each shell bed ±1 SD for the five sampled horizons through time. Raw richness is depicted in pale gray bars and approximates that presented by Dockery (13, 42). Note difference in scales. (B) Proportion of taxa showing first appearances in a given unit per lineage million years (Lma), standardized for sampling intensity and normalized by duration of intervening time between each shell bed, in comparison to raw, nonnormalized proportions (gray). (C) Same as in (B), for last appearances. OTB, Ostrea thirsae beds; UH, upper Hatchetigbee Formation. Time scale is from Gradstein (60, 61).

Differences in sample-standardized richness can arise solely from differences in abundance structure, so we use sample coverage analysis (15) to compare units at a common level of sampling completeness (for example, 0.8 in Fig. 2A; see Materials and Methods). Analysis of pooled occurrence data reveals a relative richness relationship among horizons consistent with that recovered by sample standardization (Fig. 2A) and similarly supports somewhat higher diversity in the BM. However, richness curves for a subset of the best-sampled collections with abundance data from each unit exhibit overlapping variation, suggesting no significant difference in the richness or abundance structure of individual assemblages from the five horizons (Fig. 2B). Higher richness in the pooled BM synoptic data is likely due to the greater geographic spread of these samples (3.6° of longitude) in comparison to other units (2.0° to 2.9°) and, hence, greater incorporation of subtle facies or biogeographic variations.

Fig. 2 Comparisons of species richness versus sample coverage within and among sampled horizons, with 95% confidence intervals.

Only collections with 10 or more species present were included. Richness values should be compared at the same level of coverage. (A) Pooled occurrence data within each horizon from the synoptic data set. (B) Abundance data for the four best-sampled collections within each of the five sampling horizons.

The composition of well-sampled faunas of the BLM and BM, immediately straddling the PETM, reinforces the similarity exhibited by diversity data for individual collections. Turritelline, strombid, and volutid gastropods, as well as venerid, venericardiine, and corbulid bivalves, dominate both assemblages. Ten of the 25 most abundant species in each unit are the same or congeners, and nearly all have family-level counterparts in the top 25 taxa of the other (table S1).

Ecological composition

Abundance data are further used to evaluate changes in the proportional representation of ecological guilds through time, as per Bush et al. (16) and modified by Sessa et al. (17). Each species is assigned motility, trophic, and tiering traits based on living congeners or family-level descriptions from the Paleobiology Database (PaleoDB; and the Neogene Marine Biota of Tropical America (NMiTA; database, with some adjustments, yielding 28 possible unique ecological combinations (data file S2). We find little support for change in proportional representation of guild categories across the boundary. Variation in ecological guild structure between pairs of collections in different horizons is little greater than that between collections within horizons, including those that straddle the PETM, as represented by the low ANOSIM (analysis of similarities) R value and large overlap in dissimilarity ranks (Fig. 3A). While the associated P value is significant, ANOSIM is biased toward finding structure where none is present (18). There is some variation among collections (for example, fig. S2), but the mean relative abundance of ecological life modes within horizons shows no trend through time, nor is there substantial change across the PETM (Fig. 3, B to D), with the exception of an increase in the relative abundance of infaunal taxa at the expense of epifaunal and semi-infaunal taxa (Fig. 3D and fig. S2) and an increase in the abundance of chemosymbiotic bivalves in the Eocene. As is typical of Cenozoic fossil assemblages in general (19) and the GCP in particular (17, 19, 20), mobile, suspension-feeding bivalves and carnivorous gastropods remain important components of the ecosystem throughout. Therefore, despite some taxonomic turnover at the PETM, ecological structure of the faunas remains largely unchanged. A similar analysis counting ecologically coded occurrences in the synoptic data set rather than individuals in the abundance data set likewise shows no significant change through time (fig. S3). Ordination of ecologically coded abundance collections using nonmetric multidimensional scaling (NMDS) illustrates the significant overlap among collections from different horizons (fig. S4).

Fig. 3 Ecological structure, using collections with abundance data and taxa coded by ecological guild.

(A) ANOSIM results comparing ranks of dissimilarities of collection pairs within each of the five horizons against those between horizons (shaded in gray). Width of spindles reflects the number of pairwise comparisons in each, and constriction marks the mean rank for each. (B to D) Proportional representation of individuals coded by (B) trophic strategy, (C) motility, and (D) tiering among pooled collections from each sampled horizon. Error bars represent 95% confidence regions obtained from the 2.5 and 97.5% percentiles of the proportion posterior distributions and are associated with the most represented ecological guild in each plot. PETM is indicated by heavy dashed line. Fac. mobile, facultatively mobile.

Clade-specific patterns in body size and life history

Venericardiine bivalves and turritelline gastropods are dominant groups with broadly similar ecologies (shallow–to–semi-infaunal, slow-moving, and suspension feeding) and, hence, might be expected to respond in similar ways to any perceived PETM environmental perturbation. Both groups show moderate richness changes from the Paleocene into the Eocene, but not always in the same direction. Venericardiines in the GCP increase from three to seven species (only one of which crosses the boundary), while turritellines drop from four to two with no shared species. Turritellines exhibit a clear size decrease at the PETM with the loss of several large taxa (Fig. 4A, top), including the iconic carinate species Kapalmerella mortoni postmortoni; Eocene species in the GCP never again reach the impressive proportions of their Paleocene predecessors (21). Venericardiine bivalves, on the other hand, with effectively the same trophic strategy, show no overall size trend through the window of study (Fig. 4A, bottom). Size frequency distributions for three of five additional common mollusk genera show small but significant decreases in body size across the PETM, although Eocene forms all overlap the size range of their Paleocene relatives (fig. S5). The absence of a consistent pattern across taxa suggests a lack of strong generalizable selection for body size during the PETM. If size selection did occur, as in deep-sea benthic foraminifera (22) and terrestrial mammals (23), its effects were transient and reversed before deposition of the BM.

Fig. 4 Body size and growth data in turritelline gastropods and venericardiine bivalves on either side of the PETM.

(A) Means ±1 SD and maximum body size for species of turritelline gastropods (top) and venericardiine bivalves (bottom) in each of the five sampled horizons; dashed line denotes the position of the PETM. From left to right, turritelline species are Turritella praecincta (GLM), Turritella multilira, K. mortoni postmortoni, Turritella praecincta, Turritella eurynome (BLM), Haustator gilberti (BM), and H. gilberti (upper Hatchetigbee Formation, UH); venericardiine species are Venericor aposmithii, V. nanaplata, V. pilsbryi (Ostrea thirsae beds OTB), Claibornicardia alticostata, Venericor aposmithii, V. nanaplata, V. pilsbryi (GLM), V. aposmithii (BLM), V. horatiana, V. bashiplata, V. turneri, Venericardia gulielmi (BM), Venericor hatcheplata, and Venericor turneri (UH). (B) Serially sampled δ18O data plotted by distance along the growth axis, revealing annual (seasonal) cycles from two turritelline gastropods in the Paleocene BLM compared to two from the Eocene BM. VPDB, Vienna Pee Dee Belemnite. (C) Same for venericardiine bivalves using both δ18O and δ13C data to help resolve ambiguous annual cycles; our interpretation is indicated by pale vertical bars. Distance axes are the same scale for Paleocene and Eocene shells in each panel, allowing comparison of growth rates (shaded gray bars).

Stable isotope sclerochronologic data reveal unambiguous changes in growth rate from the Paleocene to the Eocene in both groups, but in opposite directions. Turritellines on both sides of the PETM lived for 1.5 to 2 years, but latest Paleocene shells attained sizes nearly three times larger than those from the Eocene, indicating much faster growth before the PETM (Fig. 4B). In contrast, early Eocene Venericor bashiplata grew faster early in life than did late Paleocene Venericor aposmithii, becoming roughly twice as large in its first 5 years (Fig. 4C). The number of potentially resolvable cycles constrained by multiple data points in the first 80 mm or so of shell growth is markedly greater for the Paleocene shell, making the pattern of faster Eocene growth hold even given significant uncertainty in identifying seasonal cycles, a consequence of low-amplitude temperature variation compounded by relatively slow growth in both taxa.


Taken as a whole, our results indicate that the long-term impact of the PETM on these shallow-water benthic communities was unremarkable. Unlike the deep-sea benthos (22), molluscan shelf associations on the GCP suffered little in the way of lasting biodiversity loss, taxonomic turnover, or persistent ecological restructuring relative to changes at earlier or later formation boundaries. Turritellines show a shift toward smaller body size, a prediction associated with both warming (24) and acidification (25), but the magnitude of change seen in this group is not matched by other important groups. Growth rates of the most common turritelline and venericardiine species do show pronounced change across the PETM, but in opposite directions. This lack of substantial or generalizable pattern in faunas that straddle the PETM suggests that any potential selection pressure imparted by environmental changes during the PETM must have been weak, taxon-specific, and/or short-lived, and ultimately inconsequential to overall molluscan evolutionary history.

One potentially selective pattern present in our data set is the statistically significant shift to infauna-dominated assemblages after the PETM, a trend documented in proportional abundance (Fig. 3 and fig. S2), but not in richness (fig. S3). This change could have been driven by selection during the PETM such that earliest Eocene faunas exhibit the relict consequences of that ecological reorganization. During the PETM, sediments of the U.S. Atlantic Coastal Plain (ACP) (8, 26) and GCP (14) both suggest at least seasonal hypoxia and eutrophication along continental margins, an increase in nutrients delivered with runoff, and remobilization of phosphate in sediments. Epifaunal organisms generally tend to be more sensitive to hypoxia than infauna (27), and infaunal bivalves are better suited to conditions of increased primary production (28), consistent with observed patterns. A trend toward increasing richness and abundance of infaunal suspension feeders and decreasing richness of epifaunal suspension feeders characterizes the Cretaceous and Paleogene section in the GCP overall (17, 20), but a marked drop in the abundance (but not richness) of epifauna coupled with an increase in abundance (but not richness) of infauna suggests a departure from the longer-term evolutionary trends, perhaps reflecting a more immediate response to PETM ecological forcing via hypoxia that is carried over into earliest Eocene assemblages. Turritellines dominate the semi-infaunal guild, and living species are known to decrease in biomass and abundance in association with hypoxia (29), consistent with their observed lower abundance and smaller body size in post-PETM assemblages. We are likewise tempted to ascribe significance to the proportional increase in both richness and abundance of chemosymbiotic bivalves following the PETM, as a chemosymbiotic ecology should be favored in hypoxic sediments with more active sulfate reduction (30). Kosnik (20) cautions that there is a good deal of heterogeneity in both the richness and abundance of chemosymbiotic taxa among formations in the Cretaceous-Paleogene GCP. Despite this, the appearance in our data set of seven species in five different genera from two families, together comprising a small but significant percentage of individuals [4.3 ± 0.003% (1 SD)], in the earliest Eocene GCP is perhaps notable when late Paleocene units each contain only one species at lower abundance. Like the shift toward infauna, this might also be a relic of selection favoring hypoxia-tolerant taxa during the PETM itself. Note that while chemosymbiotic taxa are also infaunal, they comprise only a small proportion of occurrences (6.8%) and individuals (6.5%) within the earliest Eocene infauna, and so do not themselves control the trend toward infaunalization.

Acknowledging the potential for some lingering effects of PETM hypoxia-driven selection in the earliest Eocene faunas, by and large the mollusk faunas of the GCP offer little sign of impact from the marked but short-lived environmental changes that transpired during the PETM globally. Work on the foraminiferal assemblage of the ACP suggests a similarly muted response in the long term, with almost no taxonomic turnover on the shelf (26). When viewed at a comparable temporal resolution (assemblages from before and after, but not during, the CIE), our mollusk record from the GCP and the foraminiferal record from the ACP shelf (26) look very much like the long-studied plant fossil record from the western United States—each shows no major lasting change. It was only after floras (or foraminiferans on the ACP) were discovered within the CIE in sections that fortuitously captured the event that the full impact of the PETM became apparent (31). A similarly significant yet transient biogeographic, ecophenotypic, and/or ecologic response is still possible and even probable in marine shelf mollusks as well, but the stratigraphic record of the GCP as yet does not allow us to test this hypothesis.

The absence of a notable P-E (Paleocene-Eocene) boundary peak in first or last appearances of marine macrofauna in a global Phanerozoic database (32) suggests that the GCP pattern is not atypical of marine shelf assemblages more broadly. If so, why is it that shelf macrofauna were evidently less affected in the long term by PETM events? Several factors could be involved, including historical “preconditioning” of faunas to warm temperatures and the comparatively slow rate of CO2 release that limited the acidification of surface waters. With respect to the former, organisms at Earth’s surface were already adapted to the warm conditions of the Paleocene, and with low equator-to-pole thermal gradients, a good deal of the planet likely experienced occasional summer temperatures that approached PETM-like warmth. Benthic mollusks, with multiyear life spans and essentially no ability to migrate, would have had to endure these episodes, fortuitously preparing them to emerge unscathed from the PETM itself. The only other shelf macrofaunas explicitly studied across the PETM are those from the tropics, where sea surface temperatures may have exceeded 36°C (33)—reef corals were replaced by larger benthic foraminiferans (11) and reef fishes experienced a major evolutionary turnover (34), having a lasting impact on the history of both groups. While tropical planktonic taxa, with the capacity for rapid biogeographic range shifts, were able to recover rapidly despite short-term precipitous drops in abundance and diversity (33), at least some macrofauna evidently suffered more lasting effects. Whether tropical mollusks experienced a similar perturbation is unknown, but our data show that GCP subtropical assemblages, at least, escaped relatively intact.

In addition to warming, ocean acidification is an anticipated consequence of the addition of CO2 to the atmosphere (35), and one with potentially serious implications for marine ecosystems (36). Carbonate dissolution is a prominent feature of deep-ocean PETM sediments (6), but carbon cycle modeling suggests that CO2 release was sufficiently slow that transfer to the large deep-sea reservoir via overturning circulation limited buildup in the atmosphere and surface ocean (2). Recent work has provided proxy evidence for PETM sea surface acidification [0.3 to 0.4 pH units or approximately 100% more acidic (2, 5)], but a biological response has yet to be demonstrated, even in the most sensitive of carbonate skeletal elements (8, 37). PETM mollusks would have been most vulnerable as larvae (36), but a comprehensive metadata analysis (36) of modern mollusks found minimal to no impact on juveniles or adults experiencing a drop of <0.4 pH units, at least on a time scale of years or less. Hence, the PETM drop in mixed-layer pH was evidently insufficient to hinder calcification or leave a lasting imprint on patterns of survivorship or body size in these skeletal invertebrates. Buffering of coastal waters by more alkaline runoff sourced from enhanced continental weathering in the aftermath of the PETM may have additionally helped to spare shelf ecosystems (38), but this process would not be instantaneous even on shelves, providing thousands of years for selection to be effective had forcing been important.

What, if any, implications do these results hold for the present and future response of shallow marine biota to ongoing global change? The rate of anthropogenic CO2 release to the atmosphere today is an order of magnitude greater than that during the PETM (2). Given the amount of carbon released to date and even the conservative projections for the future, the degree of warming, acidification, and hypoxia of the surface ocean associated with the PETM is likely to be a “best-case” scenario for what to expect in our geologically near future (9). The markedly faster rate of environmental change today combined with the fact that modern organisms are instead adapted to a cooler climate highlight the differences between the PETM and the present day. The minimal response on long (geological) time scales to a perturbation that took place over several thousands of years does not imply that modern shelf ecosystems are not at risk from the climate change playing out over the next few hundred years.


Experimental design

Our goals overall were to identify and sample shell beds that closely bracket the PETM, to quantitatively compare aspects of their diversity and ecology, and to assess taxon-specific properties of body size and life history in two prominent clades. This multipronged and interdisciplinary approach should bring to light any significant biological change potentially associated with PETM environmental change and place it into the context of “background” change in the section spanning this interval. Below, we detail the stratigraphic context of the shell beds sampled and step through the numerical analyses we used to examine taxonomic richness and turnover; ecological composition based on relative representation of guilds and taxon similarity; and body size, growth rate, and life span of turritelline gastropods and venericardiine bivalves straddling the boundary. All data and R code for faunal analysis are provided in the Supplementary Materials.

Geologic setting

The Paleogene section in the GCP consists of a thick sequence of mostly unconsolidated, mixed carbonate and siliciclastic deposits from marginal and fully marine, shelf environments interspersed with fluvio-deltaic sediments (39, 40). Marine mollusk faunas have been studied for well over a century and enjoy a consistent and comprehensive taxonomy (41). Paleontologists have long recognized clusters of first and last appearances in these deposits, the largest of which are associated with major regional regressions and form the basis for GCP lithostratigraphic group boundaries (42, 43). The Paleocene-Eocene boundary lies not at one of these turnovers, but within the Wilcox Group (Fig. 1), already suggesting that PETM taxonomic change may not be as significant as elsewhere in the section (13).

Sluijs et al. (14) have placed the PETM between the abundantly fossiliferous shell beds of the BLM of the Tuscahoma Formation below and the BM of the overlying Hatchetigbee Formation above; intervening sediments are sparsely to nonfossiliferous with regard to mollusks. They argue that the PETM CIE (<150 ka) (44), encompassing the most extreme environmental conditions, has been removed to varying degrees by subsequent erosion at the base of the BM across the GCP, consistent with expectation in dynamic sedimentary environments (12). Therefore, the PETM in the GCP is closely bracketed by diverse mollusk assemblages separated by no more than 1 million years (45, 46), but there are no mollusks yet known from the CIE itself. Nevertheless, while the transient, immediate response to warming is unavailable for study, any lasting impact of the PETM on these faunas should be readily apparent.

Taxonomic diversity and turnover

We tabulated species-level taxonomic occurrence and abundance data on gastropods, bivalves, and scaphopods from collections in five discrete shell beds in the central GCP (Alabama and Mississippi; fig. S1), three from the latest Paleocene and two from the earliest Eocene. All were deposited in closely similar shallow-shelf environments (17, 47, 48), are generally unlithified and so minimize bias against small and delicate taxa (49), and routinely preserve original aragonite [for example, (50, 51)]. Given biostratigraphic age constraints, shell beds must be separated by about 1 Ma or less. Note that while some authors consider the Eocene BM and upper Hatchetigbee Formation (UH) to be time-equivalent (52), their respective species of planicostate Venericardia are distinct and virtually nonoverlapping, and we followed Mancini and Tew (45) in treating them as successive units.

We analyzed a synoptic data set of 117 existing and new collections with species-level taxonomic occurrence data (492 taxa) and an abundance data set comprised of 64 of those collections (294 taxa) that contain counts of individuals (data file S1). The synoptic data encompass 1967 taxonomic occurrences that derive largely (but not exclusively) from Palmer and Brann (41) together with references below for abundance data. Taxa in this reference were originally presented with a list of localities at which each species was found; for our work, species were parsed into unique collections (in this case, localities), each with its own species list. There is no way to know how many times a particular locality referenced by Palmer and Brann was visited and (re)sampled, so each unique locality has only a single associated list of taxonomic occurrences. These collections differ somewhat from those associated with abundance data in that each of the latter represents a unique collection event, with multiple visits to, and collections from, the same locality possible. Although generic assignments for some taxa (for example, “Pleurotoma,” “Corbula,” and “Tellina”) have not been updated, the species treatments all stem from Palmer and Brann (41) and hence are internally consistent and will not affect taxon counting. We followed Dockery (13, 42) in counting all species, subspecies, varieties, and unnamed forms recognized by Palmer and Brann (41), as did the authors of the abundance data sets below. This practice increases richness and turnover rates over a species-only data set but does not affect the overall patterns revealed in the analyses described below. Abundance data include 11,986 specimens representing 1181 occurrences and were derived largely (but not exclusively) from Sessa et al. (17) and Toulmin (47) based on counts of sieved bulk samples. While collection methods were unclear for Toulmin (47) and Palmer and Brann (41), Sessa et al. (17) found general agreement among diversity and ecology data from multiple authors across a set of the same localities. The minimum number of individuals in a collection is generally given by the number of identifiable fragments containing an apex (for gastropods) or containing an umbo divided by two (for bivalves). All data are archived in the PaleoDB with collections tied to this paper (reference #34008), as well as their original published references, and are available in the Supplementary Materials. Note that the taxonomy for some species has yet to be completed in the PaleoDB, and hence, some older genus names were retained in data file S1. These have been updated to the degree possible in the text, figures, and other supplementary materials.

Shifts in taxonomic richness and turnover were first evaluated using resampling techniques on the synoptic data set to standardize for differences in sampling intensity. Taxonomic occurrences across all collections within each shell bed were pooled and subsampled, with replacement, to the number of occurrences associated with the least sampled unit (UH, 167 occurrences). For each set of random draws, we tabulated sampled richness and the proportion of taxa that make their first and last appearances in each horizon. The processes was repeated 1000 times, yielding a mean, sample-standardized richness for each unit along with turnover among them that can then be compared to each other and to patterns in the raw data. Standardized richness values were not normalized for the duration of intervening time because shell beds are taphonomically similar, of comparable stratigraphic thickness, and thin relative to their enclosing units, and so were thought to encompass approximately similar amounts of time. Standardized turnover statistics were normalized for the estimated duration of time separating shell beds, yielding the rate of first (or last) appearances per lineage million years for each interval. It should be noted that first and last appearances in our data set reflect an as-yet unknown combination of biogeographic range shifts and in situ origination and extinction, both facilitated by the significant changes in sea level recorded in the section (42).

Trends in standardized richness could still be affected by differences in abundance structure across units (53). To explore this possibility, the richness and abundance structure of faunas in each unit were further assessed via sample coverage analysis (15) on occurrence (incidence) data and abundance data. This approach acknowledges that differences in abundance distributions make conventional rarefaction—comparing the richness of assemblages at the same level of sampling—less likely to recover true differences in richness, and argues instead for comparing samples at the same level of completeness, or coverage (15, 53). Samples are equally complete at sample sizes for each that essentially make the probability of finding a new species with the addition of one new individual (or occurrence) equal. We rarefied to equal levels of sample coverage (R code S1), rather than sample size, to complement among-horizon richness patterns recovered through traditional resampling (described above). Occurrence frequencies for each taxon in pooled, well-sampled collections (arbitrarily defined as those containing ≥10 species) were treated as taxon abundances in a single reference sample for each horizon, and horizons were then compared using coverage-based rarefaction on incidence data. In addition to pooled incidence data, we compared richness values for the four largest (number of individuals) samples from each of the five horizons to illustrate differences in alpha diversity (collection-level richness) within and among horizons using coverage-based rarefaction and extrapolation (15). Taxonomic compositional change across the PETM was also illustrated through a tabulation of shared taxa in the uppermost Paleocene BLM and the lowermost Eocene BM (table S1).

Ecological composition

Proportional representation of guilds within and among horizons was compared using ANOSIM with a Bray-Curtis distance matrix (R code S2). To illustrate guild breakdown within the motility, trophic, and tiering categories (data file S2), individuals of taxa coded by guild were pooled within horizons, and the SD of the most common guild’s proportion for a horizon was calculated using a nested, mixed-effect model (fig. S2) (54), where collections were treated as random effects (R code S3).

Within a horizon, the mean proportion of an ecological category isEmbedded Imagewhere there are m collections in the horizon and ni taxa with the ecological model of interest in collection i. Each βi consists of a population fixed effect and random effectEmbedded ImageThe proportion of interest for the horizon is Embedded Image. Proportions for each collection were estimated using best linear unbiased predictors (BLUPs), Embedded Image. Data were fitted using the lme4 package in the R statistical package. Owing to the possibility that estimates of p could be close to 0 or 1, Bayesian confidence regions were estimated by simulating the posterior distribution of p and selecting the 2.5 and 97.5% quantiles. These computations were performed using the R package arm. Collections with individuals coded by life mode were also ordinated using NMDS so as to visualize overlap in ecological composition.

Clade-specific patterns in body size and life history

Clade-level analysis of body size and life history focused on turritelline gastropods (family Turritellidae) and venericardiine bivalves (family Carditidae), two taxa ubiquitous in the GCP throughout the Paleogene (21, 55). Body size of all available taxa in both groups was quantified through time, and we report both mean and maximum size measures for each (data file S3). For turritellines, body size was given as specimen lengths measured from collections at the Paleontological Research Institution (PRI) and our own field collections, corrected for the frequently missing apex (56). For venericardiines, we used centroid size based on nine homologous landmarks digitized in the lateral (internal) orientation from collections at the PRI, the National Museum of Natural History, the Mississippi Department of Environmental Quality, and our own field collections. In addition, to explore changes in GCP shell size more generally across the PETM, we measured the maximum dimension of specimens of the gastropod genera Athleta, Calyptraphorus, Natica, and Pleurotoma and the bivalve genus Tellina in the PRI collections on either side of the PETM (data file S3) and compared their size distributions using a t test of mean size for Paleocene versus Eocene taxa. Multiple species were measured for each genus (see data file S3).

Size is a function of growth rate and life span, and changes in these two important life history variables give rise to, and can be masked by, body size trends. To better understand the factors determining body size in turritellines and venericardiines, growth rates and life spans were evaluated through stable oxygen and carbon isotope analysis of shell material serially sampled along ontogenetic trajectories. We sampled two turritellines and one large venericardiine each from the BLM and BM, the units straddling the PETM, at high spatial resolution (data file S4). Taxa analyzed are closely related but were probably not direct ancestor-descendant pairs (21, 57). They are, however, the dominant species of their family in each unit. Turritellines were sampled around the spiral growth axis, and venericardiines were sampled from the umbo to commissure. Only one to two individuals from each formation were sampled, a consequence of the need for a large number of samples within an individual ontogenetic trajectory to constrain seasonal cycles. While a greater number of individuals would ideally be preferred, similar data from both turritellines (58) and venericardiines (59) in other settings demonstrate general consistency in size-age relationships among individuals of a species in a given unit.

Annual cycles were identified on the basis of regular deviations from the mean in oxygen and carbon isotope values along the growth trajectory. If growth is sufficiently fast and seasonal variation is large, the number of samples within any year is sufficient to recognize the expected annual sinusoid. However, slow or seasonally truncated growth combined with the muted range of seasonal temperature variation in low-latitude settings can make delineation of years difficult. While turritellines are fast-growing and so posed little problem, age determination in these venericardiine bivalves was especially difficult due to the limited resolution afforded by isotope data for both the aforementioned reasons. To minimize the problem, we sampled on the external surface of shells instead of on cross sections where the need for sufficient material for analysis limits the temporal resolution possible. This strategy, however, eliminated the possibility of identifying growth bands in cross section, which can help ascertain age if they can be shown to be annual. Ultimately, annual cycles in venericardiines were cautiously identified using variation in both oxygen and carbon isotope values, particularly by looking for runs of multiple samples that constrain cycles, but we concede a good deal of uncertainty in these assignments.

Samples were analyzed at the University of Michigan’s Stable Isotope Laboratory or at the University of Kansas Keck Paleoenvironmental and Environmental Stable Isotope Laboratory using a MAT 251 or 253 mass spectrometer coupled to a Kiel automated carbonate preparation system. All values are reported relative to the Vienna Pee Dee Belemnite standard, and precision of values was better than 0.05 ‰.


Supplementary material for this article is available at

Data file S1. Taxonomic data set.

Data file S2. Ecological classifications.

Data file S3. Body size data set.

Data file S4. Stable isotope data set.

R code S1. Chao

R code S2. Analysis of similarities

R code S3. Collection random effects

Fig. S1. Locality and stratigraphic information.

Fig. S2. Variation in the proportion of individuals classified as “epifaunal” or “shallow infaunal” plus “deep infaunal” in collections of the five horizons.

Fig. S3. Ecological structure using synoptic data and taxa coded by ecological guild.

Fig. S4. Ordination of collections by NMDS of abundance data, with taxa coded by life mode, using Bray-Curtis dissimilarity.

Fig. S5. Size frequency distributions for five taxa across the PETM.

Table S1. Twenty-five most abundant mollusk species from collections in the BLM (Paleocene, n = 23) and the BM (Eocene, n = 67) for a comparison of taxonomic similarity across the Paleocene-Eocene boundary.

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.


Acknowledgments: We thank D. Dockery from the Mississippi Department of Environmental Quality for providing and helping with additional bulk sample material, for assisting in the field on multiple occasions, and for his long-standing support of research on GCP paleontology. A number of students helped with sample processing, size measurements, and data entry, including E. Judd, M. Kosloski, U. Smith, P. Stoddard, L. Eccles, E. Fenlon, R. Lee, I. Sobalvarro, K. Ohman, D. Reed, K. McClure, R./T. Lee, and T. Schlossnagle. B. Gollands constructed the data entry portal for assembling collection lists from the Palmer and Brann monograph. Thanks to L. Ward (Virginia Museum of Natural History) for help in the field. E. Thomas, C. Jaramillo, J. Jackson, two anonymous reviewers, and members of the PaleoX seminar group at Syracuse University read the manuscript and provided much helpful feedback. This represents PaleoDB paper #314. Funding: Research was supported by collaborative NSF Earth Sciences awards to L.C.I. (0719645), W.D.A. (0719642), and R.L. (0718745). Author contributions: L.C.I., W.D.A., and R.L. conceived the project, secured funding with J.A.S., conducted field work, and oversaw student contributions. J.A.S. and C.P. processed and tabulated the bulk sample data. J.A.S. and L.C.I. tabulated monographic works into occurrence-based collections and archived the data with the PaleoDB. J.A.S. and C.P. extracted and tabulated all the faunal data from the PaleoDB. C.P. performed guild analyses and sampled venericardiines for stable isotopes. R.L. oversaw the body size and phylogeny work on venericardiines. W.D.A. oversaw the body size and phylogeny work on turritellines. L.C.I. oversaw life history analyses using stable isotopes. J.C.H. wrote and implemented the R code and performed all statistical analyses, except those associated with body size. L.C.I. and C.P. wrote the manuscript and drafted the figures. All authors read the manuscript and provided edits. Competing interests: The authors declare that they have no competing interests. Data and materials availability: New and existing taxonomic data for these analyses reside in the PaleoDB (; collections used are linked to this paper (PaleoDB reference #34008) and were derived largely from the works of K. W. Palmer and D. C. Brann, L. D. Toulmin, and J.A.S. All data and R codes are also provided in full in the Supplementary Material files associated with this paper or may be requested from the authors.

Stay Connected to Science Advances

Navigate This Article