Non-technical Summary
Paleo-ecological niche modeling (paleoENM) is used to estimate the environmental conditions and distribution of extinct species, relying on fossil occurrence data and environmental climate simulations. While general circulation models (GCMs) are commonly used to recreate ancient climate conditions, the elevation models used in these simulations can influence the predictions of species’ past distributions. We examine two different paleo-elevation model-based GCMs (Getech and Scotese) and compare their effects on predicting the niche of a group of shelly marine invertebrates, the tower snails (Cerithoidea: Turritellidae) before and after the Cretaceous–Paleogene (K-Pg) mass extinction. The results show that the GCM using the Getech-based elevation model provided more accurate predictions of suitable habitats, while the Scotese-based model tended to overestimate the distribution. Both models suggested that the ecological niches of Turritellinae remained similar across the two periods. This study highlights how variations in elevation models can impact paleoENM outcomes and offers a new perspective on how marine invertebrates may have persisted through mass extinctions. It also provides important insights into the reliability of different models used in paleoecological research and lays the groundwork for future studies on fossil invertebrates.
Introduction
Here we test the hypothesis that the Cretaceous–Paleogene (K-Pg) mass extinction impacted the niche preferences of turritelline gastropods using fossil occurrence and general circulation model (GCM) data within an ecological niche modeling (ENM) framework. We project ENM results onto geography with the aim of detecting changes in the biogeographic distribution of Turritellinae between the Maastrichtian and Danian and identify the environmental niche variables that may have been associated with those distributional changes. Furthermore, we test the robustness of the paleoENMs developed with GCMs based on differing paleo-elevation models, which include differences in coastlines, to better understand the impact of model parameterization choices on model fit and realism.
ENM
ENM is a popular technique utilized by ecologists to estimate the abiotic niche dimensions and associated geographic distribution of a species by correlating observational (species’ occurrence) data with spatially explicit environmental data (Guisan and Zimmerman Reference Guisan and Zimmerman2000; Guisan and Thuiller Reference Guisan and Thuiller2005; Zurell et al. Reference Zurell, Franklin, König, Bouchet, Dormann, Elith, Fandos, Feng, Guillera‐Arroita and Guisan2020). Spatially explicit occurrence data (latitude and longitude) can be obtained from natural history collections, citizen science observations, or abundance data based on field sampling, whereas environmental data (raster files) can be derived from field stations or remote sensing satellites in the form of climate variables. Such modeling has become an important tool for answering ecological, evolutionary, and biogeographic questions and addressing issues in conservation biology and climate change (Guisan and Zimmerman Reference Guisan and Zimmerman2000; Guisan and Thuiller Reference Guisan and Thuiller2005; Pearman and Weber Reference Pearman and Weber2007; Peterson et al. Reference Peterson, Soberón, Pearson, Anderson, Martínez-Meyer, Nakamura and Araújo2011; Anderson Reference Anderson2012, Reference Anderson2013). Although ENMs allow us to potentially predict future species’ ranges constrained by anthropogenic effects on climate, the application of ENM to deep time (paleoENM) has been limited due to inherent biases within the fossil record and the challenges of reconstructing paleoenvironments (Myers et al. Reference Myers, Stigall and Lieberman2015). That said, it is also paleoENM that possesses the greatest potential to fill gaps in spatially patchy fossil records (Myers et al. Reference Myers, Stigall and Lieberman2015). Furthermore, paleoENM can also provide equally pertinent information about past changes in species’ distributions due to climate change at a variety of spatiotemporal scales that are not available using modern data.
PaleoENMs are hindered in accuracy and predictive power by two distinct factors: the completeness and abundance of focal fossil taxa and the generation of high-resolution paleoclimatic data (Myers et al. Reference Myers, Stigall and Lieberman2015). In recent years, the centralization of large biodiversity datasets (Global Biodiversity Information Facility, iNaturalist, Paleobiology Database) has allowed more accurate and large-scale modeling of species niches and their geographic distributions. However, accurate representation of fossil species’ ranges remains variable across different taxa and ultimately incomplete due to a litany of preservation, geologic, sampling, and taphonomic biases that result in uncertainty in estimating paleo-distributions of species (Myers et al. Reference Myers, Stigall and Lieberman2015). As such, the niches of fossil taxa are only robustly obtained if fossils are well-preserved, easily identifiable, and abundant (Waterson et al. Reference Waterson, Schmidt, Valdes, Holroyd, Nicholson, Farnsworth and Barrett2016). Despite these challenges, paleo-occurrence data offer the strength of preserving taxa within the environments in which they lived over millions of years, as opposed to data for modern-day taxa, which only provide a single timeframe of where the species lives presently.
ENMs and Paleoclimates
Modern ENMs typically rely on high-resolution (~1 km2 pixel size) global environmental layers (variables derived from temperature and precipitation measured at weather stations), which are readily available from online databases (e.g., Worldclim, Earthdata, MARSPEC); more recently, environmental layers have been able to reach back into the Pliocene epoch using stable isotope data (Brown et al. Reference Brown, Hill, Dolan, Carnaval and Haywood2018). In contrast, standard global climate layers within deep geologic time have not been readily available, necessitating paleobiologists to construct their own environmental layers using sedimentological and geochemical proxies (Stigall and Lieberman Reference Stigall and Lieberman2005; Maguire and Stigall Reference Maguire and Stigall2008; Dudei and Stigall Reference Dudei and Stigall2010; Stigall Reference Stigall2012; Brame and Stigall Reference Brame and Stigall2014; Purcell and Stigall Reference Purcell and Stigall2021). However, within the last 5 years, the incorporation of GCMs and state-of-the-art paleo-elevation data have allowed paleoclimatologists to simulate global climate variables through the Phanerozoic (Lunt et al. Reference Lunt, Farnsworth, Loptson, Foster, Markwick, O’Brien, Pancost, Robinson and Wrobel2016; Valdes et al. Reference Valdes, Armstrong, Badger, Bradshaw, Bragg, Crucifix, Davies-Barnard, Day, Farnsworth and Gordon2017, Reference Valdes, Scotese and Lunt2021; Farnsworth et al. Reference Farnsworth, Lunt, O’Brien, Foster, Inglis, Markwick, Pancost and Robinson2019). In brief, these models utilize a coupled ocean–atmosphere GCM in conjunction with paleogeographic digital elevation maps (DEMs), and information on past atmospheric carbon dioxide (pCO2) concentration and solar luminosity, to generate environmental raster data such as global maps of seasonal atmospheric and oceanic temperature and salinity. The HadCM3 family of GCMs has been used recently to estimate the distribution of theropods before and after the K-Pg mass extinction, as well as the evolution of palms, crocodilians, turtles, and birds (Waterson et al. Reference Waterson, Schmidt, Valdes, Holroyd, Nicholson, Farnsworth and Barrett2016; Chiarenza et al. Reference Chiarenza, Mannion, Lunt, Farnsworth, Jones, Kelland and Allison2019, Reference Chiarenza, Farnsworth, Mannion, Lunt, Valdes, Morgan and Allison2020, Reference Chiarenza, Waterson, Schmidt, Valdes, Yesson, Holroyd, Collinson, Farnsworth, Nicholson and Varela2023; Saupe et al. Reference Saupe, Farnsworth, Lunt, Sagoo, Pham and Field2019; Lim et al. Reference Lim, Huang, Farnsworth, Lunt, Baker, Morley, Kissling and Hoorn2022). Another potential rich source of well-preserved fossil taxa can be found within invertebrates, where in some cases (e.g., trilobites, bryozoans, brachiopods, and mollusks) the fossil record is so complete that large-scale patterns of turnover and evolutionary succession can be obtained (Lidgard et al. Reference Lidgard, McKinney and Taylor1993; Hopkins Reference Hopkins2013; Carlson Reference Carlson2016; Ferrari and Hautmann Reference Ferrari and Hautmann2022). Given that GCM-based paleoENMs of invertebrate fossils have generally only extended 3 Myr into the past (Saupe et al. Reference Saupe, Hendricks, Portell, Dowsett, Haywood, Hunter and Lieberman2014, Reference Saupe, Qiao, Hendricks, Portell, Hunter, Soberón and Lieberman2015), application of this technique to invertebrates in deep time is a promising avenue for examining patterns in ecology and evolution.
Although GCMs provide an invaluable opportunity to explore paleo-distributions of invertebrate fossil taxa, they are not perfect representations of paleoenvironments. For example, spatial and temporal resolution is low (i.e., in this model, 3.75° × 2.5° geographic pixel size and geologic stage-level time bins), and models often fail to predict interesting features of paleoclimate (e.g., the reduced latitudinal temperature gradient observed in the middle Cretaceous; Harper et al. Reference Harper, Farnsworth, Valdes, Markwick and Stockdale2022). Furthermore, one must be cautious to ensure that paleo-elevation and pCO2 data input into the GCM align with paleo-rotated localities of taxa from the present day back to the geologic past (Harper et al. Reference Harper, Farnsworth, Valdes, Markwick and Stockdale2022). This presents a conundrum of which paleo-elevation model to use in the context of fossil invertebrate distribution. Here in particular, we assess the sensitivity of GCMs to different paleo-elevation reconstructions and their impacts on paleoENM predictions.
Study System
We explore different GCM inputs to the analysis of niche change in a group of shelly marine invertebrates before and after the K-Pg mass extinction, specifically the gastropod subfamily Turritellinae (Cerithoidea: Turritellidae). Turritellid gastropods comprise some of the most abundant, geographically widespread, and recognizable marine macrofossil taxa throughout the Cretaceous and Cenozoic, with recent compilations suggesting approximately 150 valid living and 800 valid fossil turritellid species (Scott Reference Scott1970, Reference Scott1974; Squires Reference Squires1987, Reference Squires, Filewicz and Squires1988; Kollmann et al. Reference Kollmann, Decker and Lemone2003; Squires and Saul Reference Squires and Saul2007; Allmon Reference Allmon2011; Plotnick and Wagner Reference Plotnick and Wagner2018; Shin et al. Reference Shin, Allmon, Anderson, Kelly, Hiscock and Shin2020). Fossil turritellids have been recorded on every continent and are commonly the numerically dominant species in fossil assemblages and living communities (Plotnick and Wagner Reference Plotnick and Wagner2006, Reference Plotnick and Wagner2018; Allmon Reference Allmon2007, Reference Allmon2011; Allmon and Cohen Reference Allmon and Cohen2008). Turritellids are known to be diverse and abundant in both Maastrichtian (latest Cretaceous; 72.1–66.0 Ma) and Danian (earliest Paleocene; 66.0–61.6 Ma) strata but exhibit very high rates of nominal extinction at the boundary. Only three species of turritellids are reported to cross the K-Pg boundary, two in California (Turritella peninsularis adelaidana Merriam, Reference Merriam1941, Turritella webbi Saul, Reference Saul1983 [Saul Reference Saul1983]) and one in the Netherlands (Turritella plana Binkhorst, 1861 [Vellekoop et al. Reference Vellekoop, Van Tilborgh, Van Knippenberg, Jagt, Stassen, Goolaerts and Speijer2020]). Given that several deeply diverging lineages (e.g., members of the subfamilies Turritellinae and Pareorinae) are also known to occur both above and below the boundary, it is likely that several lineages survived the extinction.
Occurrence and abundance of turritelline species are hypothesized to be influenced by water temperatures, depth, salinity, sediment type (siliciclastic vs. carbonate), the presence of coral communities, and primary productivity (Allmon Reference Allmon1988). Data from turritelline-dominated assemblages suggest that turritelline gastropods have become less thermophilic and more associated with siliciclastic (as opposed to carbonate) environments since the Cretaceous (Allmon and Knight Reference Allmon and Knight1993), with the transition occurring not at the K-Pg boundary, but rather after the Eocene (Allmon Reference Allmon2007). This occurs despite the extremely high apparent taxonomic turnover at the K-Pg boundary.
Methods
Taxonomy of Turritellinae
Although turritellids are well known from the fossil record, the intrafamilial relationships among them remain an active area of taxonomic research. Turritellidae was divided into five subfamilies, Turritellinae, Protominae, Pareorinae, Turritellopsinae, and Vermiculariinae by Marwick (Reference Marwick1957), with Turritellinae including most nominal genera and subgenera. Subsequent analyses, informed by both molecular and morphological data, have found that the Vermiculariinae are well nested within the Turritellinae, which may also be the position of the Protominae (Anderson and Allmon Reference Anderson and Allmon2024). Many generic and subgeneric names exist in the literature, but very few are used consistently. Recent work supported by new molecular analyses (Lieberman et al. Reference Lieberman, Allmon and Eldredge1993; Sang et al. Reference Sang, Friend, Allmon and Anderson2019; Anderson Reference Anderson2018; Anderson and Allmon Reference Anderson and Allmon2024) has supported Marwick’s (Reference Marwick1957) prioritization of protoconch form, order of appearance of spiral ornamentation, and growth-line characters as indicative of taxonomic affinity, in the absence of significant morphological differentiation in the teleoconch (e.g., Vermicularia, Caviturritella) (Friend et al. Reference Friend, Anderson, Altier, Sang, Petsios, Portell and Allmon2023). However, the comparative lack of available well-preserved protoconch material leaves many species broadly assigned to “Turritella sensu lato” on the basis of plesiomorphic or convergent teleoconch form, making “Turritella” a true “wastepaper basket” genus (Allmon Reference Allmon1996, Reference Allmon2011; Plotnick and Wagner Reference Plotnick and Wagner2006; Allmon and Cohen Reference Allmon and Cohen2008; Hendricks et al. Reference Hendricks, Saupe, Myers, Hermsen and Allmon2014). Yet species assigned to Turritella s.l. (given their absence of teleoconch characters that assign them to other subfamilies) are all expected to be encompassed by a monophyletic Turritellinae. Furthermore, extant species within Turritellinae all have similar ecologies as shallow infaunal, suspension-feeding (and occasionally deposit-feeding) gastropods, with the overwhelming majority also having highly similar life histories (Allmon Reference Allmon1988, Reference Allmon2011; Anderson and Allmon Reference Anderson and Allmon2020). Therefore, we consider the Maastrichtian and Danian members of Turritellinae, a clade first appearing in the Jurassic (Das et al. Reference Das, Saha, Bardhan, Mallick and Allmon2018) and diversifying in the Cretaceous, to be broadly like the modern gastropod genus in terms of species richness, clade age, and ecological similarity among constituent species (Allmon Reference Allmon2011).
Unfortunately, confusion in the literature and large databases like the Paleobiology Database is not limited to the proper generic assignments of species. Large numbers of fossil occurrences of “Turritella sp.” do not represent members of monophyletic Turritellinae (or even Turritellidae) but may represent other high-spired gastropods mistakenly identified as “Turritella.” Furthermore, the name “Turritella" has been erroneously applied to entire rock layers such as “turritella agate” when such layers belong to the freshwater gastropod Elimia tenera Hall, Reference Hall1845 (Allmon and Knight Reference Allmon and Knight1993). Therefore, we chose to restrict our analyses to only those occurrences that had species-level identification.
Because of these taxonomic and temporal record uncertainties, we have decided to analyze the subfamily Turritellinae as a single taxon. We recognize that ENM techniques are best interpreted at the species level, however, taxonomic reality and spatiotemporal resolution supports a more robust analysis at the clade level (Hendricks et al. Reference Hendricks, Saupe, Myers, Hermsen and Allmon2014).
Occurrence Records
We acquired occurrence records of species belonging to Turritellinae within the Maastrichtian and Danian stages from the Paleobiology Database (PBDB; paleodb.org), Global Biodiversity Information Facility (GBIF; gbif.org, downloaded July 2022), additional occurrences from C.E.M. (Myers et al. Reference Myers, MacKenzie and Lieberman2013), as well as an exhaustive literature search for the U.S. Coastal Plain occurrences by W.D.A. and K.D.C. We selected genera belonging to Turritellinae based on taxonomic, morphological (Allmon Reference Allmon1988, Reference Allmon2011; Harzhauser and Landau Reference Harzhauser and Landau2019; Friend et al. Reference Friend, Anderson, Altier, Sang, Petsios, Portell and Allmon2023) and recent molecular phylogenetic evidence from Anderson (Reference Anderson2018) and Anderson and Allmon (Reference Anderson and Allmon2024) (Table 1). Phylogenetic evidence suggests a deep split between the majority of turritellid gastropod species, assigned to the subfamily Turritellinae, and the subfamily Pareorinae Marwick, Reference Marwick1957 which includes the genera Mesalia and Sigmesalia (Anderson Reference Anderson2018; Anderson and Allmon Reference Anderson and Allmon2024); we omitted both genera from the analysis (Table 1). Molecular evidence also suggests that two other traditional subfamilies, Vermiculariinae Dall, Reference Dall1913 and Protominae Marwick, Reference Marwick1957, are nested within the Turritellidae (Anderson and Allmon Reference Anderson and Allmon2024); however, neither of these have Maastrichtian or Danian fossil records, and therefore this inclusion does not affect the analyses. Because the genus Turritella s.s., is nested within Turritellinae, and fossils are commonly misidentified as “Turritella sp.,” we only selected taxa that possessed a species epithet (in the “accepted name” category within the PBDB) and that could be verified in the literature as being present within the selected time periods. Further occurrence filtering consisted of omitting specimens for which latitude and longitude could not be computed. We binned the time intervals to the Maastrichtian (72.1–66 Ma), which possessed 2032 occurrences, and the Danian (66–61.6 Ma), consisting of 1883 occurrences (See Supplementary Information).
Genera and species of Turritellidae included and omitted from ecological niche modeling (ENM) analysis within the Maastrichtian and Danian periods. The subfamily Turritellinae, considered by Marwick (Reference Marwick1957) and Allmon (Reference Allmon2011) as inclusive of the majority of extant members of Turritellidae, also properly encompasses several genera that have previously been assigned to other subfamilies (e.g., Zaria and Vermicularia), but nest within Turritellinae in molecular analyses (Lieberman et al. Reference Lieberman, Allmon and Eldredge1993; Anderson Reference Anderson2018; Sang et al. Reference Sang, Friend, Allmon and Anderson2019; Anderson and Allmon Reference Anderson and Allmon2024). We therefore list all genera and putative genera included in our concept of a monophyletic Turritellinae for the purpose of our analysis, regardless of their occurrence in the time period of interest. As many high-spired gastropods are mistakenly assigned to “Turritella sp.,” we included only occurrences identified with a species epithet (the “accepted name” category within the Paleobiology Database [PBDB]) and verified in the literature as a member of Turritellinae. Generic assignments presently employed by the PBDB for these species are reproduced here as downloaded for clarity, including species that may be represented in the PBDB in multiple ways. The genus Trobus was included, although it is uncertain whether it is most properly assigned to Turritellidae or Casiopidae (fide Bandel and Dockery Reference Bandel and Dockery2016)
Paleogeographic maps of the Phanerozoic predominantly follow the paleo-digital elevation model PALEOMAP, which estimates Earth’s past paleoceanography, the changing area of land, mountains, shallow seas, and deep oceans through time (Scotese Reference Scotese2016; Scotese and Wright Reference Scotese and Wright2018). Fossil databases such as the PBDB follow PALEOMAP for paleo-rotation of taxa to their estimated coordinates of deposition within the Phanerozoic. The GCM simulations of past climate from Valdes et al. (Reference Valdes, Scotese and Lunt2021) have utilized Scotese’s model of paleo-elevation. Paleo-elevation models such as those from the Getech Plc (Getech.com) have also been used in GCM simulations from Lunt et al. (Reference Lunt, Farnsworth, Loptson, Foster, Markwick, O’Brien, Pancost, Robinson and Wrobel2016) and Farnsworth et al. (Reference Farnsworth, Lunt, O’Brien, Foster, Inglis, Markwick, Pancost and Robinson2019); these provide alternative reconstructions of Earth’s past climate. However, differences in the paleo-coastlines of each model (particularly epicontinental sea boundaries) may result in the reconstruction of coastal invertebrate occurrences on land instead of in the ocean (Scotese and Wright Reference Scotese and Wright2018).
Using the rgplates package v. 0.3.2 (Kocsis et al. Reference Kocsis, Raja and Williams2023) in R v. 4.0 (R Core Team Reference Team2021), we paleo-rotated both datasets to their respective time periods using Scotese’s PALEOMAP (Scotese Reference Scotese2016) and Getech’s paleo-rotation model provided by Farnsworth et al. (Reference Farnsworth, Lunt, O’Brien, Foster, Inglis, Markwick, Pancost and Robinson2019). PALEOMAP consists of 120 unique paleo-digital elevation models (paleoDEMs) representing 3 Myr time slices roughly equating to different stratigraphic stages. We paleo-rotated the occurrences of the Maastrichtian to the 69 Ma time slice, and the Danian to the 66 Ma time slice. We repeated this process using Getech’s in-house model with the same time slices (Farnsworth et al. Reference Farnsworth, Lunt, O’Brien, Foster, Inglis, Markwick, Pancost and Robinson2019; see Supplementary Information for all paleorotated datasets).
Climate Data
We generated environmental layers using output from the HadCML3 climate model v. 4.5 (Valdes et al. Reference Valdes, Armstrong, Badger, Bradshaw, Bragg, Crucifix, Davies-Barnard, Day, Farnsworth and Gordon2017), with inputs from the solar luminosity, and both the PALEOMAP paleogeographic atlas (Scotese Reference Scotese2016; Scotese and Wright Reference Scotese and Wright2018) (hereafter referred to as the Scotese simulations) and from Getech Plc (Getech.com) (hereafter referred to as the Getech simulations) for the Maastrichtian and Danian periods. For both sets of simulations, we use the GCM “HadCM3LB-M2.1aD,” described in detail in Valdes et al. (Reference Valdes, Armstrong, Badger, Bradshaw, Bragg, Crucifix, Davies-Barnard, Day, Farnsworth and Gordon2017), for which surface resolution is 3.75° longitude × 2.75° latitude (grid box size of ~420 × 220 km at the equator, rescaling to ~200 × 280 km at 45° latitude).
The Scotese simulations are similar to the latest Maastrichtian (Map number 22) and Danian (Map number 20) simulations described in Valdes et al. (Reference Valdes, Scotese and Lunt2021). These are the simulations that prescribe a pCO2 concentration according to Foster et al. (Reference Foster, Royer and Lunt2017). Compared with the Valdes et al. (Reference Valdes, Scotese and Lunt2021) study, simulations were run for an additional 2000 years and with modified atmospheric and ocean physics by applying methods similar to those described in Sagoo et al. (Reference Sagoo, Valdes, Flecker and Gregoire2013), resulting in an improved representation of polar amplification in deep-time climates. In addition, they have islands defined correctly for the purposes of calculating the ocean barotropic stream function, according to Foreman (Reference Foreman2005).
The Getech simulations are identical to those described in Farnsworth et al. (Reference Farnsworth, Lunt, O’Brien, Foster, Inglis, Markwick, Pancost and Robinson2019). The pCO2 concentration was 1,120 ppmv for both the Maastrichtian and Danian Getech simulations, which is within the typical range of uncertainty in the Foster et al. (Reference Foster, Royer and Lunt2017) pCO2 data, and approximates the actual evolution of pCO2 through time with some uncertainty (see Foster et al. Reference Foster, Royer and Lunt2017: fig. 1).
Climate model output variables chosen within the Maastrichtian and Danian for both model sets were: annual mixed layer depth (meters), annual maximum and minimum (warmest and coldest) sea-surface temperature within a 3 month interval (season; Celsius), and monsoon seasonality index (difference in precipitation between the three driest and three wettest months of the year); monsoon seasonality is a proxy for preference of seasonal variability around the tropics. Furthermore, we acquired annual potential temperature (Celsius), and annual salinity (practical salinity unit [PSU]) averaged from 5 to 5800 m depths (for in-depth definitions of variables, see Valdes et al. Reference Valdes, Armstrong, Badger, Bradshaw, Bragg, Crucifix, Davies-Barnard, Day, Farnsworth and Gordon2017). Although modern turritellines are found predominantly between depths of 2 and 100 m (Allmon Reference Allmon2011), occurrences affected by paleo-rotation and placed at greater depths would be mistakenly removed if layers were restricted to shallow depths including epeiric seas found in North America and Eurasia. Finally, we incorporated bathymetry (meters), derived from the digital elevation maps (paleoDEMs) that underlie both the Getech and Scotese GCM simulations. These variables were chosen to encapsulate the abiotic restrictions of modern turritelline species (Allmon Reference Allmon1988, Reference Allmon1996, Reference Allmon2011; Allmon and Knight Reference Allmon and Knight1993; Allmon and Cohen Reference Allmon and Cohen2008). We rotated all variables from 0 to 360 longitudinal format to −180 to 180 longitudinal format to avoid cutoff at the international date line.
Sedimentology Data
We generated two distinct sedimentological, interpolated datasets derived from lithology and depositional environment proxies at fossil localities within the Maastrichtian and Danian. We downloaded all marine occurrences within the K-Pg boundary stages from the PBDB and assigned each collection to either carbonate or siliciclastic using the primary lithologic data associated with the occurrence, recorded as “lithology1” field in the collection record (see Hopkins [Reference Hopkins2014] and Foote [Reference Foote2006] for criteria). Within the PBDB, each lithological term possesses specific definitions so that the collection is consistently described in ways relative to all other collections in the database, allowing for consistent translation to the two lithological categories. We also assigned each collection to different depositional environments recorded from the “environment” field to determine preferred depositional settings for turritellines within the Maastrichtian and Danian. In total, we used 58 depositional categories. We supplemented the PBDB lithological and environmental datasets with localities from Paleo Reefs (https://www.paleo-reefs.pal.uni-erlangen.de/; Kiessling et al. Reference Kiessling, Flügel and Golonka1999), and the Sedimentary Geochemistry and Paleoenvironments Project (SGP; https://sgp.stanford.edu/; Farrell et al. Reference Farrell, Samawi, Anjanappa, Klykov, Adeboye, Agic, Ahm, Boag, Bowyer and Brocks2021). We classified all paleo-reef collections as being carbonate, and only selected collections from the SGP belonging to one of the two lithological categories. We categorized paleo-reef collections as “reef,” while the SGP only possessed “basinal,” “fluvial,” “inner shelf,” “lacustrine,” and “outer shelf” environmental categories, all of which are consistent with the PBDB “environment.” Definitions for lithologies and depositional environments can be found at https://paleobiodb.org/public/tips/lithtips.html and https://paleobiodb.org/public/tips/environtips.html, respectively. It should be noted that depositional environment categories as defined and used in the PBDB are not strictly independent of one another and are defined in a hierarchical fashion. Thus, our analysis provides an example of a very simplified use of this data. In the Supplementary Information, we provide a detailed description of how these categories can be condensed into more directly comparable categories and return to this topic in the “Discussion” (Supplementary Table S1).
We paleo-rotated the lithological and depositional environment datasets to the Getech and Scotese GCM paleo-rotation models within the Maastrichtian and Danian. We assigned each lithological and environmental category an integer (e.g., 1 = carbonate, 2 = siliciclastic) and performed a nearest-neighbor interpolation using the interpNear function in the R package terra (Hijmans et al. Reference Hijmans, Bivand, Forner, Ooms, Pebesma and Sumner2022) to generate categorical raster layers. When lithological or depositional categories vary within a single pixel, the nearest-neighbor interpolation assigns the pixel’s value based on Euclidean distance—selecting the occurrence closest to the pixel’s center. Nearest-neighbor is an efficient method for interpolating categorical data, as it preserves original categories, avoids blending (i.e., averaging), and is less computationally intensive (Johnson and Clarke Reference Johnson and Clarke2021). We resampled the categorical rasters to the same resolution of the GCM layers using the resample function in the raster package in R (Hijmans et al. Reference Hijmans, van Etten, Mattiuzzi, Sumner, Greenberg, Lamigueiro, Bevan, Racine and Shortridge2021). Climate model, lithological, and depositional environment output variables can be obtained in our Supplementary Information via Dryad digital repository number: https://doi.org/10.5061/dryad.w9ghx3g17.
ENM
Before modeling, we spatially thinned occurrences to the resolution of the raster layers to reduce sampling bias, artificial clustering, and subsequent spatial autocorrelation (Veloz Reference Veloz2009) using the spThin package (Aiello-Lammens et al. Reference Aiello-Lammens, Boria, Radosavljevic, Vilela and Anderson2015), which resulted in more even sample sizes of occurrence points across time periods (see Table 2). ENM modeling is more accurate when occurrences are thinned to the resolution of pixels within our raster layers, such that each occurrence gets a singular value of each of the climate and sedimentological values used (Araújo et al. Reference Araújo, Anderson, Barbosa, Beale, Dormann, Early and Garcia2019). Similarly, after thinning, only a single occurrence may exist in any single pixel, which is an important step to prevent model oversampling of the same pixels within the study extent (Aiello-Lammens et al. Reference Aiello-Lammens, Boria, Radosavljevic, Vilela and Anderson2015). Given that the collated occurrence data often contained already-assigned environmental conditions (e.g., lithology and depositional environment in PBDB-based records), we checked whether (1) multiple occurrences within the same pixel had the same environmental designations and (2) thinning occurrences caused discrepancies between the environmental values derived from our lithology and depositional environment layers and the original PBDB-based record designations. We found no discrepancies between the PBDB environmental designations linked to occurrences within a pixel and the interpolated pixel values from our lithology and depositional environment layers.
Sample size of occurrence localities before and after paleo-rotation spatial thinning and maximum entropy (MaxEnt) ecological niche modeling (ENM) parameterization and statistics for Getech and Scotese general circulation models (GCMs) within the Maastrichtian and Danian time periods. Parameters shown are feature class (features) and regularization multiplier (regularization). Statistics shown are mean validation area under the curve (AUCval), mean omission rates for 10-percentile training values (OR10), Akaike information criterion (AICc), and delta Akaike information criterion (delta.AICc), and number of nonzero coefficients (ncoef).
Turritellines possess a global distribution (Allmon Reference Allmon2011); thus we chose a study extent to encapsulate the entire preserved geographic range of Turritellinae in both the Maastrichtian and Danian time periods. Furthermore, to include potentially undersampled areas despite dispersal limitations of individual species (Peterson and Soberón Reference Peterson and Soberón2012), we defined a bounding box buffered by 8.12° (~100,000 m2), as this is the maximum estimated dispersal capability of modern-day turritellines (“M” training region as defined by Soberón and Nakamura [Reference Soberón and Nakamura2009]; Allmon Reference Allmon1996, Reference Allmon2011). From the total region encompassing all thinned occurrences enveloped by the bounding box, we randomly sampled background points for modeling (n = 888 for Scotese and n = 1842 for Getech study extents); background points differ between niche models due to the number of occurrence points differing between paleo-rotation models that produced the resultant bounding box. Background points are a random sample of the environmental conditions across our entire study area that help define what is suitable and unsuitable habitat in model output. Our models compare environmental conditions at our presence (occurrence) points with those of our background points to identify which conditions are more strongly associated with the presence of the taxon (Phillips and Dudík Reference Phillips and Dudík2008; Elith et al. Reference Elith, Phillips, Hastie, Dudík, Chee and Yates2011; Phillips et al. Reference Phillips, Anderson, Dudík, Schapire and Blair2017; Phillips Reference Phillips2021). We extracted environmental values from our background points to calculate correlations between variables using the vifcor function in the usdm package (Naimi Reference Naimi2017) and filtered out variables with correlation coefficients higher than 0.7. In ENM analyses, the balance of correlation threshold against the number of environmental variables retained for analysis is an area of active study, with coefficients ranging from 0.5 to 0.9 (Graham Reference Graham2003; Yan et al. Reference Yan, Feng, Zhao, Feng, Wu and Zhu2020). Because fewer variables create underfitted models, we chose a moderate correlation threshold (0.7), which acted as a natural break in which higher thresholds significantly reduced the number of environmental variables retained. We used the same set of variables for subsequent model building, retaining only those that remained after assessing correlations between GCM type and time period.
We used the machine learning algorithm MaxEnt v. 3.4.4. (Phillips et al. Reference Phillips, Anderson, Dudík, Schapire and Blair2017), which remains one of the top-performing algorithms for fitting presence-background ENMs (Valavi et al. Reference Valavi, Guillera‐Arroita, Lahoz‐Monfort and Elith2021). Even though a substantial number of occurrences were removed during thinning, sample size was large enough to produce robust results (Shcheglovitova and Anderson Reference Shcheglovitova and Anderson2013; Table 2). At small spatial and taxonomic scales, fewer than 10 occurrences can be sufficient to create highly accurate models (Hernandez et al. Reference Hernandez, Graham, Master and Albert2006; Pearson et al. Reference Pearson, Raxworthy, Nakamura and Peterson2007; van Proosdij et al. Reference van Proosdij, Sosef, Wieringa and Raes2016). A recent study by Qiao et al. (Reference Qiao, Peterson, Ji and Hu2017) suggested that incorporation of occurrences of related species increases model fidelity on undersampled taxa. Nonetheless, sample size should be considered carefully, because the number of occurrence points varies greatly depending on the study system, and too few paleoENM studies have been published to generate a consensus (see Myers et al. Reference Myers, Stigall and Lieberman2015: supplementary information). Within ENMs, partitioning data involves splitting the dataset into two subsets: the training dataset is the subset of data used to build the model, while the testing dataset is used to evaluate the performance of each model iteration (Phillips and Dudík Reference Phillips and Dudík2008; Elith et al. Reference Elith, Phillips, Hastie, Dudík, Chee and Yates2011; Phillips et al. Reference Phillips, Anderson, Dudík, Schapire and Blair2017; Phillips Reference Phillips2021). Partitioning reduces the risk of model overfitting and assesses its performance with unseen data. We used the “checkerboard2” strategy of spatial partitioning, which divides up our study extent into a grid of spatial cells whereby occurrence points are then alternatively assigned to different partitions (the default setting used here assigns four partitions).
All final models were fit to the full datasets (training + testing) (Phillips et al. Reference Phillips, Anderson, Dudík, Schapire and Blair2017). As the combination of two key complexity settings in MaxEnt models, feature classes and regularization multipliers, can strongly influence model outputs (Warren and Seifert Reference Warren and Seifert2011; Radosavljevic and Anderson Reference Radosavljevic and Anderson2014), we tuned model complexity to find optimal settings. For tuning, in order of increasing complexity, we chose the feature classes linear (L), quadratic (Q), and hinge (H), as well as regularization multipliers 1 through 5 (higher numbers penalize complexity more) (see Supplementary Information for a more detailed explanation of MaxEnt model building). In brief, feature classes determine the shape of the model fit, while regularization multipliers control how much complexity is penalized—this can result in predictor variable coefficients shrinking to 0 and thus dropping out of the model (Phillips and Dudík Reference Phillips and Dudík2008).
We assessed optimal models using sequential criteria that included threshold-dependent (omission rate) and threshold-independent (area under the curve [AUC]) performance metrics. We first filtered models that possessed a delta Akaike information criterion (with correction for small sample sizes, AICc) less than 2 to compare models between time periods and GCM type (Warren and Seifert Reference Warren and Seifert2011). We then chose models that possessed the smallest 10-percentile omission rate. If a series of models possessed near-identical 10-percentile omission rates, we chose the simplest model or the one with the fewest nonzero lambda values (model coefficients). For each time and GCM model, we documented variable importance and plotted marginal response curves to better understand the modeled relationships between the predictor variables and the data. We recorded the permutation importance metric output by MaxEnt, which is calculated by randomly permuting the values of all environmental variables but one, building a new model, then calculating the difference between each model’s training AUC and that of the empirical model (Phillips et al. Reference Phillips, Anderson, Dudík, Schapire and Blair2017). Marginal response curves are generated by constraining all predictor variables to their means except for one, then making model predictions along the full range of the focal variable associated with the training data. These curves show the modeled relationship of each variable individually with the occurrence data when all other variables are held constant and are affected by the complexity of the model settings (Phillips et al. Reference Phillips, Anderson, Dudík, Schapire and Blair2017).
We calculated a multivariate similarity surface (MESS; Elith et al. Reference Elith, Kearney and Phillips2010) to detect the degree of similarity in extracted empirical occurrence (presences) environmental values and extracted background environmental values for each time period and GCM type. Increasingly negative values in MESS scores suggest dissimilar environmental conditions (non-analogue environments resulting in model extrapolation), while positive values suggest more similar environmental conditions (interpolation). In the case of our data, we calculated MESS plots by extracting environmental values from occurrence points and projected them to their respective time period. However, to prevent over-extrapolation, we performed an “informed” MESS analysis by clamping each environmental variable a priori based on the completeness of response curve (e.g., extrapolation can be used if the response curve for an environmental variable reaches 0 or 1 for full modeled behavior). Code for our informed MESS analysis (as well as all other aspects of model building) is provided in our Supplementary Information, which was modified from the mess function in the R package dismo (Hijmans et al. Reference Hijmans, Phillips, Leathwick and Elith2017).
We made habitat suitability predictions for Turritellinae using our environmental predictor variables. We reclassified continuous prediction output to binary (0 = unsuitable; 1 = suitable), determined by the maximum sum of sensitivity and specificity (MaxSSS) threshold calculated from model evaluation metrics. MaxSSS remains one of the best threshold selection methods for presence/absence data (Liu et al. Reference Liu, Newell and White2016). We also compared our MaxSSS predictions to those calculated from the 10-percentile omission rate which is one of the strictest binary thresholds (Phillips et al. Reference Phillips, Anderson, Dudík, Schapire and Blair2017). Finally, we generated continuous predictions by transforming our raw MaxEnt predictions to a scale of 0 – 1 to approximate the probability that Turritellinae will be present at a particular location (“cloglog” transformation) (Phillips et al. Reference Phillips, Anderson, Dudík, Schapire and Blair2017).
Niche Overlap
We compared niche overlap of the Maastrichtian and Danian occurrences using an ordination framework by first reducing dimensionality within the datasets via a principal component analysis (PCA). Using the espace_pca function in the ade4 package v. 1.7 in R (Dray and Dufour Reference Dray and Dufour2007; Kass et al. Reference Kass, Vilela, Aiello‐Lammens, Muscarella, Merow and Anderson2018), we generated convex hulls of principal component space using the extracted environmental variables from the Turritellinae occurrences within both time periods and plotted with correlation loadings to infer the degree of influence particular environmental variables possess in the distribution of background points within niche space (see Supplementary Information for convex-hull reconstructions). The espace_occDens function in the package ecospat v. 3.2 (Di Cola et al. Reference Di Cola, Broennimann, Petitpierre, Breiner, d’Amen, Randin, Engler, Pottier, Pio and Dubuis2017) was used to estimate an occurrence density grid for both the environmental values at each occurrence point and background extent points using a kernel density estimation approach. Niche overlap between occurrence density grids of environmental values at each occurrence point and background points were compared using Schoener’s D (Schoener Reference Schoener1968) using the espace_nicheOv function.
Finally, using the ecospat.plot.overlap.test function, we conducted niche similarity and equivalency tests to assess the degree of niche differentiation of Turritellinae between time periods. For the similarity test, a null distribution was generated by randomly shifting the occurrences of Turritellinae from one time period (Maastrichtian) within the combined background extent of both time periods, while keeping the niche from the other time period (Danian) fixed. This process is permuted 1000 times to assess whether the observed niche overlap was greater than expected by chance. This process was then repeated by swapping time periods, shifting occurrences from the Danian while keeping the Maastrichtian niche fixed. A p-value < 0.05 indicates that the niches from both time periods are more similar than expected by chance. In contrast, the equivalency test pools the occurrences of Turritellinae from both time periods and randomly assigns them within their combined background extent to assess whether the niches are statistically indistinguishable; this process is also permuted 1000 times. A p-value < 0.05 suggests that the niches are significantly different from each other, indicating they are less equivalent than expected by chance.
Results
Optimal Model Settings
Optimal settings for the Getech and Scotese MaxEnt models differed in complexity between model type. Due to the results of the collinearity analysis, we used the following seven predictor variables: annual potential temperature, annual potential salinity, annual mixed layer depth, bathymetry, monsoon seasonality, depositional environment, and lithology. Most models were complex, with linear quadratic features (curvilinear with splines) and regularization multipliers of 4–5 (high complexity penalty). Overall, the Getech MaxEnt method showed higher model performance, possessing high mean validation AUC values, low 10% omission rates, low amounts of nonzero lambdas (coefficients), and low AICc values, suggesting low amounts of model overinflation (Table 2). The Scotese MaxEnt models performed poorly, with the Maastrichtian model having a low omission rate of 10% (0.12), but an unacceptably low AUC value (0.50), indicating performance no better than a random model (Table 2).
Suitability Predictions
Within our Getech MaxEnt models, predictions based on MaxSSS revealed suitable areas along the coasts of north and southeastern North America, north and western Africa, throughout the Caribbean and Neo-Tethys Ocean, northern India, and southeastern Australia. Furthermore, within our Maastrichtian model, we observed suitable areas within the Western Interior Seaway, the inland waterways of Eurasia, and off the coast of North America, which are either closed or absent within the Danian model (Fig. 1A,B).
Thresholded maximum entropy (MaxEnt) predictions for Turritellinae derived from our Getech general circulation models (GCMs) within the Maastrichtian (A) and Danian (B) and Scotese GCMs within the Maastrichtian (C) and Danian (D). Predictions were derived from the optimal model using the criterion of area under the curve (AUCtest) values being the highest and 10% omission rate being the lowest. Threshold predictions were calculated from the sum of sensitivity and specificity (MaxSSS) from our model evaluation. Yellow areas are predicted suitable; purple areas are predicted unsuitable; orange circles delineate occurrence data used to train the model.
In contrast, the Scotese MaxEnt models predicted an overly large extent of suitable habitat, encompassing nearly the entire globe (Fig. 1C,D). We could not compare differences between the time periods due to the large swaths of suitable area found within both time periods, with suitability within the Maastrichtian encompassing the entire extent, suggesting model overinflation (Fig. 1C,D).
Suitability predictions based on 10-percentile omission rate and cloglog transformation produced near-identical patterns for both GCM types and time periods (Supplementary Figs. S1, S2).
Permutation Importance and Response Curves
Optimal models retained 4–6 predictor variables after regularization and relied mainly on contributions from bathymetry, annual potential temperature, depositional environment, lithology, monsoon seasonality index, and annual mixed layer depth (Table 3). Response curves of climate variables and bathymetry revealed negative relationships with suitability (i.e., higher preferences to shallower depths and stronger seasonal variation). Depositional environments that possessed the highest suitability included shallow subtidal, transition zone, shoreface, deltaic, and marine substrates, consistent with modern turritelline ecology. Lithological data possessed minimal importance within the models (Table 3, Fig. 2, Supplementary Figs. S3, S4).
General circulation model (GCM) type, geologic stage, and environmental variable within each optimal maximum entropy (MaxEnt) model, permutation importance of each environmental variable, and the resultant response curve behavior expressed with that variable in relation to suitability or lithology/depositional environment. MaxEnt calculates permutation importance by changing the values of each environmental variable at random, then calculating the difference using the area under the curve (AUC) from the “training data.” The values of each environmental variable are randomly permuted within the training presence and background data, the resultant drop in training AUC is calculated, then normalized to percentages. Response curves show how each environmental variable individually affects the MaxEnt prediction in terms of increasing or decreasing suitability. The response curve behavior is dictated by model complexity (feature classes and regularization multipliers)
Environmental response curves and permutation importance of the Maastrichtian (top) and Danian (bottom) for our Getech general circulation model (GCM) maximum entropy (MaxEnt) models. We recorded the permutation importance metric output by MaxEnt, which is calculated by randomly permuting the values of all environmental variables but one, building a new model, then calculating the difference between each model’s training area under the curve (AUC) and that of the empirical model (Phillips et al. Reference Phillips, Anderson, Dudík, Schapire and Blair2017). Response curves show how each environmental variable individually affects the MaxEnt prediction in terms of increasing or decreasing suitability. Behavior of response curve is dictated by model complexity (feature classes and regularization multipliers). Dark blue represents probability densities of environmental values for our thinned occurrence (presence) dataset, while yellow represents the densities of our background points.
Niche Overlap and MESS Analysis
Ordination analysis revealed substantial overlap between both models, with environmental variables showing similar loadings in principal component space (Fig. 3, Table 4). The first principal component (PC 1) explained ~35% of the variance, while PC 2 accounted for ~24%. In both models, annual potential temperature, salinity, and monsoon seasonality index loaded strongly on PC 1, whereas annual mixed layer depth and bathymetry loaded strongly on PC 2 (Table 4).
Principal components analysis (PCA), correlation circle, and niche equivalency and similarity tests showing turritelline niche differences within the Maastrichtian and Danian between our Getech (top row) and Scotese (bottom row) general circulation model (GCM) variables. Solid contour lines in the PCA illustrate full range (100%) of climate space (“fundamental niche”), while dashed lines indicate 50% confidence intervals. Contour lines of the Maastrichtian are in orange; the Danian is blue. Shading shows the density of species’ occurrences per grid cell of the kernel density analysis (“realized niche”), and violet pixels show niche stability (climate conditions occupied in both time periods). Orange shading indicates climate conditions only occupied in the Maastrichtian, while blue indicates climate conditions only occupied in the Danian. Correlation circle indicates climactic variable loadings on the PCA space. Length and direction of arrows indicate influence and distribution of variables within PCA space, respectively. Histograms represent observed (red bar) and randomly simulated overlaps of pairwise niche equivalency and similarity. Pairwise niche similarity tests also showed observed overlaps exceeding 95% of simulations confirming high niche similarity. However, niche equivalency tests in both models indicated overlaps were lower than 95% of simulated values, supporting niche equivalency between the two time periods.
General circulation model (GCM) type, loadings, importance of variables, percent overlap, percent of unique niche overlap from both time periods, environmental similarity (Schoener’s D), and results of pairwise similarity, and equivalency tests from our ecospat analyses. D = Danian; M = Maastrichtian; PC 1/PC 2 = principal component 1/2.
In the Getech model, niche overlap between the Maastrichtian and Danian was high (Schoener’s D = 0.37), with pairwise similarity tests indicating observed overlaps exceeded 95% of simulated overlaps (p = 0.03), strongly supporting niche similarity (Fig. 3). The Scotese model showed similarly high overlap (Schoener’s D = 0.31), with Maastrichtian occurrence density grids being smaller and substantially overlapping with the Danian (Fig. 3). Pairwise niche similarity tests also showed observed overlaps exceeding 95% of simulations (p = 0.02), confirming high niche similarity (Fig. 3). Similarly, niche equivalency tests in both models indicated overlaps were lower than 95% of simulated values (p = 0.83 for Scotese, p = 0.61 for Getech), supporting niche equivalency between the two time periods.
Finally, the MESS analysis revealed similar environmental conditions within the Maastrichtian and the Danian among the Getech variables, with larger areas of environmental similarity within the Maastrichtian (Fig. 4A,B). These results were more apparent among the Scotese GCM variables, in which environmental similarity of niche space was drastically higher in the Maastrichtian than the Danian (Fig. 4C,D).
Multivariate environmental similarity surface (MESS) analysis described in Elith et al. (Reference Elith, Kearney and Phillips2010) and calculated using the dismo package in R using our Getech (top) and Scotese (bottom) variables for the Maastrichtian (left) and Danian (right) stages. MESS plots were calculated using environmental values extracted from occurrence points from each time period and projected to the same time period to infer the degree of environmental similarity, but clamped in accordance with the response curves of each resultant model (“informed” MESS). Negative scores (shown in red) indicate novel climate conditions (i.e., environmental values that fall outside the range of environmental variables in the training region). Positive scores (shown in blue) indicate similar climate conditions to the training dataset.
Discussion
Optimal Model Settings
Overall, the Getech MaxEnt models produced more robust results in terms of model performance (e.g., omission rate and AUC) and reasonable suitable habitat predictions, which were explained by similar combinations of linear and quadratic feature classes. The Getech-based models also showed high regularization for model complexity, resulting in similar model outputs regardless of geologic stage. In contrast, the Scotese-based MaxEnt models showed suboptimal performance, as evidenced by low AUC values and overinflated predictions. Differences in model performance may be partially due to the data lost during the paleo-rotation and occurrence-thinning stages of model building. The precise time for which a paleo-elevation model is reconstructed, along with any associated spatial averaging at the resolution of the model, coupled with the time averaging of fossil occurrence data across periods of different eustatic and local sea levels, can cause taxa from genuine marine paleoenvironments to appear to be falsely classified as terrestrial occurrences. Furthermore, low resolution of environmental layers can result in substantial loss of sample size after thinning. In this analysis, the Maastrichtian Getech MaxEnt model retained occurrences located at the mouth of the Western Interior Seaway (n = 8) that were removed by rotation and spatial thinning within the Scotese (n = 0) simulations (Fig. 1A,C). We obtained similar results when partitioning the data using the jackknife strategy of data partition, which is used for small sample sizes, as it allows for the largest training datasets possible for model building during cross-validation (Pearson et al. Reference Pearson, Raxworthy, Nakamura and Peterson2006; Shcheglovitova and Anderson Reference Shcheglovitova and Anderson2013).
Permutation Importance and Response Curves
For most models, bathymetry, annual mixed layer depth, depositional environment, and monsoon seasonality index possessed the highest importance for predicting the range of turritellines (Table 3). Previous work has suggested that the biogeographic distribution of turritellid gastropods was, and continues to be, strongly influenced by seasonal changes in temperature (Allmon Reference Allmon1988, Reference Allmon2007; Allmon and Knight Reference Allmon and Knight1993). However, our temperature variables showed minimal importance (Table 3, Fig. 2). Although ENMs involving depth-related variables are currently being explored (Owens Reference Owens2022), future research would require interpolation of more depth-related temperature variables, while accounting for geographic changes of past continents—a nontrivial task. However, incorporation of 3D ENMs may reveal more temperature influence in the distribution of turritellines.
Occurrences of Turritellinae exhibited a broad bathymetric range (~2–5800 m), with suitability being highest at ~5 m within the response curves (Fig. 2, Supplementary Fig. S3), however, values considered suitable ranged between 5 and ~3665 m within the threshold (MaxSSS) predictions. Depth preferences for fossil turritellines remain poorly constrained by the models, with living turritelline species generally inhabiting depths of 10–100 m (Allmon Reference Allmon2011). However, species have been found living in very shallow water <10 m or depths as great as 1200–1500 m (Allmon Reference Allmon1988). Whereas these results overlap with the known depth ranges of modern turritellines, we are skeptical of occurrences exceeding ranges of modern-day species. Instead, the coarseness of depth resolution here may account for the wide range of depth values among turritelline localities. This is further supported by thinned occurrences possessing high bathymetric estimates, but shallow-water depositional environmental categorizations (see Supplementary Tables S2–S5).
Annual mixed layer depth remained consistent across models, with the highest habitat suitability occurring at approximately 10 m in the response curves (Fig. 2, Supplementary Fig. S3) and ranging from 10 to 90 m in the thresholded predictions. Because the mixed layer is the ocean’s uppermost zone, where temperature and salinity are relatively uniform due to wind and wave action, its depth plays a crucial role in shaping the thermal conditions that marine communities experience (Graham et al. Reference Graham, Kinlan, Druehl, Garske and Banks2007). Additionally, annual mixed layer depth can change throughout the year, being influenced by regional climate and location (Kara et al. Reference Kara, Rochford and Hurlburt2003). Turritellines are likely to thrive within the observed mixed layer depths, as greater depths may lead to cooler conditions that fall outside their preferred temperature and salinity limits (Allmon Reference Allmon2011).
Suitable habitat for annual potential temperature peaked at 34°C in the response curves and ranged from 5°C to 34°C in the thresholded predictions. This temperature range is reasonable, as most occurrences were clustered within the tropics, where waters are naturally warmer. Turritellines are predominantly tropical to temperate species (Wilf et al. Reference Wilf, Johnson and Huber2003; Allmon Reference Allmon2011; Fig. 1).
The monsoon seasonality index possessed the highest suitability at −1.6, while values ranged from −1.6 to 1.5, within our thresholded predictions (Fig. 2, and Supplementary Figs. S3–S5). Monsoon seasonality is the difference between summer and winter precipitation (precipitation of 3 driest months − precipitation of 3 wettest months), suggesting that turritellines prefer stronger seasonal variation, with suitable habitats within areas of higher precipitation regardless of season. Modern turritellines are often abundant and diverse in localities with strong seasonal nutrient delivery, either from upwelling (which can be associated with changes in wind direction that also correspond to dry/wet seasonality; e.g., Anderson et al. Reference Anderson, Hendy, Johnson and Allmon2017) or episodes of strong fluvial input associated with nutrients from terrestrial runoff (associated with a wet season near river mouths; e.g., Scholz Reference Scholz2020; Scholz et al. Reference Scholz, Petersen and Anderson2024).
Depositional environment possessed the most informative data pertaining to habitat preferences of fossil turritellines. Lithology was relatively unimportant in determining suitable facies for turritellines, as both carbonate and siliciclastic depositional environments were supported, suggesting no significant preference (Fig. 2, Supplementary Figs. S3–S5). Throughout both models and time periods, suitable habitats of turritellines were located within siliciclastic environments such as deltas/fluvial systems, as well as carbonate environments such as shallow subtidal, and transition zone/lower shorefaces; all facies are located within or adjacent to marine settings where fossil and extant turritellines have previously been found. The Getech Danian MaxEnt model possessed freshwater settings such as fluvial/deltaic and lacustrine substrates as suitable; however, Shin et al. (Reference Shin, Allmon, Anderson, Kelly, Hiscock and Shin2020) observed that modern-day turritellines are found at “salinities as low as 10–15 psu,” which includes brackish, but not truly freshwater, settings. This is consistent with work by Scholz et al. (Reference Scholz, Petersen and Anderson2024), who found that most turritellines are fully marine, although they can be abundant in environments with riverine nutrient input and salinities that are well below fully marine conditions, at least seasonally. Consequently, given the lack of lacustrine or fully freshwater modern Turritellinae, these attributions are likely erroneous and may result from time and spatial averaging between occurrences of turritellids in deltaic environments and adjacent freshwater environments.
Suitability Predictions
We observed unreasonably large areas of suitable habitat within the Scotese MaxEnt predictions, suggesting those model sets are unreliable and potentially inaccurate in their predicted distribution of the subfamily (Fig. 1C,D), despite AICc values and the number of nonzero lambdas (parameters) suggesting minimal overinflation (Table 2). Explanations for overfitting within these models could be due to differences in Getech- or Scotese-based paleo-rotation, which could influence the distribution and potential removal of points discussed above.
Previous distribution studies of fossil taxa involved generating MaxEnt models using modern-day climate layers and occurrences, with fossil record predictions being generated by projecting modern-day ENM predictions onto paleoclimate layers (Waterson et al. Reference Waterson, Schmidt, Valdes, Holroyd, Nicholson, Farnsworth and Barrett2016; Jones et al. Reference Jones, Mannion, Farnsworth, Valdes, Kelland and Allison2019; Saupe et al. Reference Saupe, Farnsworth, Lunt, Sagoo, Pham and Field2019; Lim et al. Reference Lim, Huang, Farnsworth, Lunt, Baker, Morley, Kissling and Hoorn2022); the advantage of this method is the increase in sample size and resolution of predictions. Because the taxonomy of Turritellinae is not as resolved as those of other vertebrate, invertebrate, or plant taxa within the fossil record, restricting the occurrences to the Maastrichtian and Danian allowed for more fine-tuned predictions of the subfamily; constructing modern-day ENMs would likely encounter non-analogue environments when projected to the Maastrichtian and Danian, which MaxEnt would classify as unsuitable habitat.
Getech MaxEnt predictions revealed minor differences in suitable habitat between both time periods; however, we observed no major habitat shifts, expansion, or contraction, suggesting that Turritellinae possessed large enough niche breadths to find suitable habitat under conditions before and after the K-Pg mass extinction. Minor differences observable in Figure 1 could be due to the stochastic nature of the probability distribution that best fits our observed data.
Studies of the effects of volcanism and asteroid impact at the K-Pg boundary indicate solar dimming would have generated global cooling of roughly 20°C (Artemieva et al. Reference Artemieva, Morgan and Party2017). The results of the asteroid impact had devastating ecological cascading effects on much of the end-Cretaceous biota (e.g., Chiarenza et al. Reference Chiarenza, Mannion, Lunt, Farnsworth, Jones, Kelland and Allison2019, Reference Chiarenza, Farnsworth, Mannion, Lunt, Valdes, Morgan and Allison2020). However, this short-term climatic change does not appear to have affected the long-term environmental space occupation of Turritellinae, as predictions were the same before and after this dramatic climate shift.
Suitability predictions of turritellines across time periods overlap substantially with shallow-water tropical environments (restricted to the coastal regions of the tropics near northern Africa, India, South America Asia, and North America; Fig. 1), regardless of substrate, supporting the findings of previous research that this subfamily had not shifted to cooler waters or away from carbonate habitats substantially before the Neogene (Allmon Reference Allmon1992, Reference Allmon2007; Kiel and Wagreich Reference Kiel, Wagreich and Wagreich2002). Alternatively, our models may not have captured the sudden changes in suitability after the K-Pg boundary because the environmental variables are averaged across geologic stages that each span millions of years, whereas conditions associated with the extinction event itself was likely short-lived on evolutionary timescales. Niche shifts may have occurred on the timescale of the extinction event and its immediate aftermath, but the event did not redirect or truncate long-term occupation of environmental space for the Turritellinae.
An interesting observation is the presence of suitable habitat observed within the Western Interior Seaway in the Getech Maastrichtian MaxEnt model prediction, yet occurrences of turritellines are only found within the southernmost portion of the seaway (Fig. 1A). This may be the result of either undercollection/underreporting of fossil occurrences in these environments or due to aspects of these environments that are not parameterized by the models. Underreporting of gastropods in epicontinental seas may be the result of taphonomic biases, including more frequent occurrences of aragonitic taxa as molds (Dean et al. Reference Dean, Allison, Hampson and Hill2019), leading to collection or identification biases (Hollis Reference Hollis2008). That said, the relative importance of taphonomic biases and habitat heterogeneity within epeiric seas remains untested (Miller et al. Reference Miller, Brett and Parsons1988; Van Iten et al. Reference Van Iten, Lichtenwalter, Leme and Simões2006). Epeiric seas are predicted to have high habitat heterogeneity and be subject to differing environmental conditions than open ocean-facing environments (Peters Reference Peters2007; Miller and Foote Reference Miller and Foote2009; Miller et al. Reference Miller, Aberhan, Buick, Bulinski, Ferguson, Hendy and Kiessling2009; Anderson et al. Reference Anderson, Pisani, Miller and Peterson2011), including differences in circulation, nutrient supply, and more rapid response to eustatic sea-level change, all of which can impact biodiversity (Lagomarcino and Miller Reference Lagomarcino and Miller2012). Another interesting observation is the presence of suitable habitat within the northern pole for our Getech MaxEnt predictions (Fig. 1A,B). No occurrences for Turritellinae have been observed within this region, suggesting inaccurate model extrapolation to these conditions or suitable environments outside of the dispersal capability of the species. The former is more likely given MESS map results, wherein climatic similarity is negative at the poles for both time periods (Fig. 4A,B). Consequently, we interpret these regions as unsuitable habitat.
Niche Overlap and MESS Analysis
Our ecospat and MESS analyses provide the opportunity to explore the environmental space occupied by turritellines within a Hutchinsonian framework (Hutchinson Reference Hutchinson1957) to visualize the position and breadth of niches and to explore otherwise omitted variables as evident by the MaxEnt permutation analysis. Ordination analyses revealed near-identical niches of the Turritellinae between time periods and model types, suggesting niche stability across the K-Pg mass extinction within this clade. Although occurrence density grids of both models noted larger niche breadths for the Danian period, niche similarity and equivalency tests do not support statistically significant niche expansion. These findings are supported by the MESS analysis and broadly indicate that environmental conditions were shared between the Danian and Maastrichtian (within the Getech GCM) (Fig. 4A,B). We observe the same patterns of stability within both GCM types, suggesting paleo-rotation is not confounding the analyses. However, the number of occurrences used in these analyses could have affected these results.
Turritellids are considered highly diverse, “bloom taxa” in the Danian (Hansen Reference Hansen1988; Allmon Reference Allmon1996; Hansen and Josefson Reference Hansen and Josefson2004), and this might traditionally be hypothesized to associate with significant niche differentiation among species. However, recent work using simulated species to develop null models of diversification under various scenarios of niche change found that the highest diversity can be produced under simulations with perfect niche conservatism (i.e., no change in niche conditions with taxonomic evolution; Qiao et al. Reference Qiao, Peterson, Myers, Yang and Saupe2024). This suggests that the notable post-recovery diversity in turritellids does not necessarily require significant changes in niche space or breadth. It is possible that facultative deposit feeding contributed to the survival of turritellines at low local abundances, facilitating a bloom in fossil occurrences once marine productivity recovered. Therefore, the finding that niche occupation remained very similar across this event aligns with Paleogene diversification and evidence from high-density assemblages, which show that turritellines as a clade exhibit niche stability—particularly regarding sediment type and bathymetric affinity—until after the Paleogene (Allmon Reference Allmon2007), despite high rates of taxonomic extinction at the K-Pg boundary. Moreover, prior studies have identified a trend of phylogenetic niche conservatism (PNC; retention of ancestral niches of closely related species over evolutionary time) among marine mollusks, further supporting the notion of niche stability in this clade (Foote et al. Reference Foote, Ritterbush and Miller2016; Tomašových and Jablonski Reference Tomašových and Jablonski2017).
Future Considerations of GCMs in PaleoENMs
Throughout data collection and model building, several key methodological aspects required in-depth discussion and evaluation, which can be utilized and refined in future work of GCMs in the context of paleoENMs. The first aspect is the pixel size (i.e., spatial grain) of our environmental variables, specifically bathymetry. The coarseness of our DEMs may have resulted in the loss of information about the bathymetric preferences of Turritellinae, as the transition from coastal to deep-water conditions was represented by only a few pixels. An avenue of future research could be to generate interpolated layers of bathymetry using facies/sedimentology/ichnology/sequence stratigraphy and geological context from PBDB occurrences or direct measurements from existing literature or new fieldwork. The bias in this method is the uneven sampling of fossils, in which some areas may possess more fine-resolution bathymetric readings than other localities, which prompted us to use a global evenly sampled DEM underpinning our GCM layers. A compromise could be to compare global distributions of a fossil taxon using low-resolution environmental data from a GCM to a regional study of the same taxon using high-resolution PBDB-generated variables, including bathymetry, depositional environment, and lithology (e.g., Purcell and Stigall Reference Purcell and Stigall2021; Purcell et al. Reference Purcell, Scuderi and Myers2023).
Another aspect to consider when generating paleoENMs from GCMs is the number and nature of the categories utilized in our depositional environment layer. Utilizing the 58 categories from the PBDB may have created noise or bias if most pixels in our raster layers were of a single category, especially because most occurrences possessed a few select underdetermined categories (“marine indet.” and “carbonate indet.”; see Fig. 2, Supplementary Fig, S3, Supplementary Tables S1–S4). The PBDB categories are also hierarchical in the degree of environmental specificity they described and so are difficult to compare directly. To briefly assess the impact on the results, we condensed the set of categories into fewer, more directly comparable categories based on the PBDB definitions (e.g., deep subtidal indet., deep subtidal ramp, deep subtidal shelf = deep subtidal; see Supplementary Table S1), excluding the most underdetermined. Overall, the model outputs produced similar results, in which lithology, depositional environment, and depositional setting have minimal importance within model building, as well as similar results to our original analysis in terms of response curves and permutation importance (Supplementary Figs. S4, S5). Future studies should consider clustering or omitting depositional categories that do not apply to their study systems (e.g., removing coastal/deltaic/lacustrine categories if the study taxon lives in deep water) or removing depositionally ambiguous categories (e.g., carbonate indet.). Overall, geological proxy data should be generated based on the most up-to-date knowledge about the study system, to ensure biological realism in model-building.
An additional aspect to consider in GCM-based paleoENMs is the use of multiple lines of evidence to assess model fidelity. Incorporating supplemental analyses such as response curves, ordination analyses, and MESS maps provides clarity as to where the model is over- or under-extrapolating, which can be interpreted as either model error or accuracy. Response curves demonstrate whether the full environmental range of the species has been sampled for each environmental variable used in the model; MESS maps give a proxy for over- or under-fitting in model predictions, whereas ordination analyses show variable importance outside spatial context.
A final consideration of the analysis is the unknown structure of lower-level taxon niche dimensions within the larger clade. Hendricks et al. (Reference Hendricks, Saupe, Myers, Hermsen and Allmon2014) warn of the use of higher-level taxonomic groups as arbitrary units of evolution, as species-rich genera could have variable life-history strategies that would inflate geographic ranges, niche breadths, or taxon durations (i.e., where the structure of species’ niches may not be not overlapping within the composite generic niche). Greater attention must be used when studying genera or subfamilies if the attributes under study reflect biological properties and processes at the species-level (Hendricks et al. Reference Hendricks, Saupe, Myers, Hermsen and Allmon2014). We argue that the niche studies of Turritellinae do reflect biological realism at the subfamily level, due to the relatively unchanged ecology and limited variety of life-history strategies utilized by Turritellinae throughout their evolutionary history (Allmon Reference Allmon2011). Although turritellines possess global distribution, large-scale patterns of evolution have been noted at the genus- and species-level within the geologic past (Allmon Reference Allmon2007, Reference Allmon2011).
Studies of fossil taxa including birds, non-avian dinosaurs, and palms (Palmaceae) have also utilized family or subfamily groups and argue that changes within or among constituent species result in subsequent changes in species composition of higher groupings (Chiarenza et al. Reference Chiarenza, Mannion, Lunt, Farnsworth, Jones, Kelland and Allison2019; Saupe et al. Reference Saupe, Farnsworth, Lunt, Sagoo, Pham and Field2019; Lim et al. Reference Lim, Huang, Farnsworth, Lunt, Baker, Morley, Kissling and Hoorn2022). This issue is essentially a question of the prevalence and degree of PNC within clades (Wiens and Graham Reference Wiens and Graham2005; Kozak and Wiens Reference Kozak and Wiens2006, Reference Kozak and Wiens2010; Wiens et al. Reference Wiens, Ackerly, Allen, Anacker, Buckley, Cornell, Damschen, Davies, Grytnes and Harrison2010; Peterson Reference Peterson2011; Xu et al. Reference Xu, Dimitrov, Shrestha, Rahbek and Wang2019; Diniz-Filho Reference Diniz-Filho2023; Tuya et al. Reference Tuya, Martínez‐Pérez, Fueyo and Bosch2024). By assuming that higher-level taxa represent lower-level (species-level) phenomena, these studies assume PNC as a matter of course. However, in the case of our data, the similarity observed in post-recovery niche occupation in Turritellinae could be explained by: (1) survival of lineages that included a variety of environmental tolerances in the clade (non-PNC, but a wide range of species-level niches), (2) survival of lineages with broad environmental tolerances (non-PNC), and/or (3) the rapid radiation of surviving species into previously occupied niche space, whether these taxa began with many non-overlapping narrow or few broad niches (PNC). Limited earliest Paleocene fossiliferous marine deposits complicate the determination of how many lineages survived across the boundary without preservation of boundary-crossing species (T. Hansen et al. Reference Hansen, Farrand, Montgomery, Billman and Blechschmidt1987; Hansen Reference Hansen1988; Keller and Barrera Reference Keller, Barrera, Sharpton and Ward1990; H. Hansen et al. Reference Hansen, Rasmussen, Gwozdz and Walaszczyk1992; T. A. Hansen et al. Reference Hansen, Farrell and Upshaw1993). Careful phylogenetic analyses across the boundary will be necessary to discriminate among these possibilities. Performing paleoENM analyses at the genus- or species-level in a modern phylogenetic context with higher spatial and temporal resolution could determine whether similar post-extinction niche occupation is the result of multiple constituent lineages with differing surviving niche dimensions or whether turritellines all shared a similar broad niche. Nevertheless, the strong phylogenetic hypothesis of monophyly among species within Turritellinae provides first-order support for PNC (Anderson and Allmon Reference Anderson and Allmon2024), and further tests comparing niche similarity to phylogenetic distance, as well as examining niche dimensions in the context of tree topology of surviving species, would be an interesting next step in unraveling turritelline–environment interactions over evolutionary timescales.
Conclusions
The combination of atmosphere–ocean GCM simulations and presence-based ENMs provide the opportunity to capture species’ environmental associations within deep time (Waterson et al. Reference Waterson, Schmidt, Valdes, Holroyd, Nicholson, Farnsworth and Barrett2016; Valdes et al. Reference Valdes, Armstrong, Badger, Bradshaw, Bragg, Crucifix, Davies-Barnard, Day, Farnsworth and Gordon2017; Chiarenza et al. Reference Chiarenza, Mannion, Lunt, Farnsworth, Jones, Kelland and Allison2019, Reference Chiarenza, Farnsworth, Mannion, Lunt, Valdes, Morgan and Allison2020; Jones et al. Reference Jones, Mannion, Farnsworth, Valdes, Kelland and Allison2019; Saupe et al. Reference Saupe, Farnsworth, Lunt, Sagoo, Pham and Field2019). Environmental layers generated from sedimentological or depositional data through large biodiversity aggregators such as the PBDB allow paleobiologists to test hypotheses of biological meaning in the long-term persistence of extinct taxa (Foote Reference Foote2006; Hopkins Reference Hopkins2014; Hopkins et al. Reference Hopkins, Simpson and Kiessling2014; Sclafani et al. Reference Sclafani, Congreve and Patzkowsky2021; Bault et al. Reference Bault, Balseiro, Monnet and Crônier2022). However, GCM boundary conditions and the underlying paleo-rotation model of occurrences involved must be in congruence to derive meaning from habitat suitability for fossil lineages (Saupe et al. Reference Saupe, Farnsworth, Lunt, Sagoo, Pham and Field2019).
This study presents the first utilization of high-resolution paleoclimatic data on a well-sampled invertebrate fossil taxon to generate a paleoecological model of suitable habitat through time. We explored differences in suitable habitat predictions generated from ENMs using climatic layers extracted from differing paleo-elevation datasets. Our findings provide a framework for paleoENM whereby the incorporation of climatic data from GCMs can be utilized in conjunction with historically useful sedimentological and geochemical proxy data. We also agree with previous findings of MaxEnt model analyses whereby data quality greatly affects model output (Owens et al. Reference Owens, Campbell, Dornak, Saupe, Barve, Soberón, Ingenloff, Lira-Noriega, Hensz and Myers2013; Araújo et al. Reference Araújo, Anderson, Barbosa, Beale, Dormann, Early and Garcia2019). More broadly, our findings contribute to understanding how the environmental dynamics of mass extinction and its recovery impacted the abiotic niche and potentially shaped Cenozoic and modern niche occupation of an important group of marine mollusks. Although nearly all latest Cretaceous turritelline species did not survive to enter the earliest Paleogene fossil record, the K-Pg extinction was not associated with reduction or change in clade niche dimensions. Turritelline niche stability in the wake of a severe mass extinction supports macroecological and macroevolutionary studies indicating that high species-level extinction may be decoupled from ecological change (Edie et al. Reference Edie, Huang, Collins, Roy and Jablonski2018) and that niche breadth of a higher taxon may be decoupled from species-level diversity (Qiao et al. Reference Qiao, Peterson, Myers, Yang and Saupe2024).
Acknowledgments
Thank you to members of the Paleobiology Database community members who entered data used in this study, particularly M. Clapham, A. Hendy, K. Layou, M. Kosnik, J. Marcot, and J. Sessa, as well as W. Kiessling, A. Miller, B. Kröger, M. Foote, and M. Patzkowsky. We also thank Nick Freymueller for his guidance of interpolation of sedimentological datasets. B.M.A.’s research has been supported by the lab of Elizabeth Petsios (Baylor University), and this material is based upon work supported by the National Science Foundation under award no. 2225014 to J. Hendricks and W.D.A., supporting B.M.A. C.E.M. and W.D.A. acknowledge support from NSF SGP 1924807 and TCN 1601878. A.F. acknowledges Natural Environment Research Council of the UK (NE/V011405/1), Leverhulme Research Project Grant RPG-2019-365 and Chinese Academy of Sciences Visiting Professorship for Senior International Scientists (2021FSE0001). D.J.L. acknowledges NERC grant NE/X000222/1. We would like to thank reviewers A. Tomašových, A. Stigall, and J. Saulsbury as well as two anonymous reviewers.
Competing Interests
The authors declare no competing interests.
Data Availability Statement
Supplementary Data including figures, code, tables, and additional analyses are available from the Dryad Digital Repository: https://doi.org/10.5061/dryad.w9ghx3g17.
Open access
