Preferential dissolution of the biogenic carbonate polymorph aragonite promotes preservational bias in shelly marine faunas. While field studies have documented the impact of preferential aragonite dissolution on fossil molluscan diversity, its impact on regional and global biodiversity metrics is debated. Epicontinental seas are especially prone to conditions that both promote and inhibit preferential dissolution, which may result in spatially extensive zones with variable preservation. Here we present a multifaceted evaluation of aragonite dissolution within the Late Cretaceous Western Interior Seaway of North America. Occurrence data of mollusks from two time intervals (Cenomanian/Turonian boundary, early Campanian) are plotted on new high-resolution paleogeographies to assess aragonite preservation within the seaway. Fossil occurrences, diversity estimates, and sampling probabilities for calcitic and aragonitic fauna were compared in zones defined by depth and distance from the seaway margins. Apparent range sizes, which could be influenced by differential preservation potential of aragonite between separate localities, were also compared. Our results are consistent with exacerbated aragonite dissolution within specific depth zones for both time slices, with aragonitic bivalves additionally showing a statistically significant decrease in range size compared with calcitic fauna within carbonate-dominated Cenomanian–Turonian strata. However, we are unable to conclusively show that aragonite dissolution impacted diversity estimates. Therefore, while aragonite dissolution is likely to have affected the preservation of fauna in specific localities, time averaging and instantaneous preservation events preserve regional biodiversity. Our results suggest that the spatial expression of taphonomic biases should be an important consideration for paleontologists working on paleobiogeographic problems.
Introduction
While the fossil record provides a unique window into past life on Earth, it is well known that it is both pervasively and non-uniformly biased (Raup 1976; Koch 1978; Foote and Sepkoski 1999; Alroy et al. 2001; Allison and Bottjer 2011a). Geologic, taphonomic, and anthropogenic biases (such as the amount of available fossiliferous rock for sampling, variation in fossilization, and the degree to which the available rock record has been sampled) skew or remove information from the fossil record, leaving the remaining catalog of data uneven and incomplete. Although biomineralized remains have an increased preservation potential compared with soft-bodied tissues (Allison 1988; Briggs 2003), they are still influenced by various geologic and taphonomic processes (Kidwell and Jablonski 1983; Kidwell and Bosence 1991; Kidwell and Brenchley 1994; Best 2008; Hendy 2011). Shelly marine faunas are especially susceptible to misrepresentation due to preferential dissolution of biogenic carbonate polymorphs. It is well established that aragonite, a polymorph of CaCO3 found within the biomineralized shells of many invertebrates, dissolves more rapidly than the more common form of CaCO3, calcite, and at a higher pH (Canfield and Raiswell 1991; Tynan and Opdyke 2011). While both polymorphs can be destroyed by adverse conditions near the sediment–water interface (Best and Kidwell 2000) and the effects of dissolution can vary between fauna (due to microstructure surface area, morphology, and shell organic content; Walter and Morse 1984; Harper 2000; Kosnik et al. 2011), it is still the case that aragonitic shells are more likely to dissolve than calcitic remains (Brett and Baird 1986). As mineral composition of mollusks is usually conserved at the family level (Carter 1990), this has the potential to skew the record of molluscan diversity and trophic structure through time (Cherns and Wright 2000; Wright et al. 2003; Cherns et al. 2011) and negatively affect subsequent work that relies on the relative abundance and distribution of shelly marine fauna (Kidwell 2005). Cherns and Wright (2000) argued that early-stage dissolution could be substantial and referred to the phenomenon as the “missing mollusk” bias. Subsequent work at a multitude of temporal and spatial scales (Wright et al. 2003; Bush and Bambach 2004; Kidwell 2005; Crampton et al. 2006; Valentine et al. 2006; Foote et al. 2015; Jordan et al. 2015, Hsieh et al. 2019) has debated the magnitude of this bias; however, there is a broad agreement on the potential for dominantly aragonitic shells to suffer greater postmortem diagenetic destruction in the taphonomically active zone (TAZ) (Davies et al. 1989; Foote et al. 2015). While the effect of dissolution on the global macroevolutionary record of mollusks has been found to be limited, possibly due to the potential of aragonite to recrystallize to calcite (Kidwell 2005; Paul et al. 2008; Jordan et al. 2015), it is conceivable that local or regional conditions could severely impact perceived patterns of biodiversity in restricted areas (Bush and Bambach 2004). In a regional study of Cenozoic mollusks, Foote et al. (2015) found evidence to suggest that aragonite dissolution was both enhanced in carbonate sediments and insignificant within siliciclastic sediments, with similar preservation potential of aragonitic and calcitic fauna within the latter. They further emphasized the fact that scale is an important factor in determining the observable impacts of aragonite dissolution, which will strongly vary between local (potentially consisting of an individual bed), regional, and global studies. However, to date, research has focused on assessing the influence of aragonite bias on temporal trends of biodiversity and has ignored the potential for direct spatial expression.
Early-stage dissolution occurs within modern environments as a result of microbially mediated reactions increasing local acidity (Walter et al. 1993; Ku et al. 1999; Sanders 2003). Bacterially mediated decay of organic material within the upper sedimentary column occurs in a series of preferential redox reactions. By-products of these reactions, such as solid-phase sulfides from sulfate reduction and CO2 from aerobic oxidation, result in changes to local pore-water saturation of calcium carbonate (Canfield and Raiswell 1991; Ku et al. 1999). Additionally, oxidation of H2S above the oxycline increases acidity at that boundary; if this occurs at the sediment–water interface, it can adversely affect the preservation of shelly marine fauna (Ku et al. 1999). As such, dysoxic sedimentary environments might have a predisposition for dissolution of biogenic carbonate and enhance the effect of the “missing mollusk” bias (Jordan et al. 2015). Epicontinental seas, marine water bodies that form by the flooding of continental interiors, are especially prone to strong water-column stratification and sea-level variation and have a predisposition to seasonally anoxic or dysoxic conditions (Allison and Wells 2006; Peters 2009). As such, they have the potential to be more prone to both preferential aragonite loss and preservation than modern oceans. Cherns et al.'s (2011) model for taphonomic gradients of aragonite preservation along a shelf-to-basin transect can be readily applied to epicontinental sea settings (Fig. 1). If we assume the center of a seaway was stratified with at least a seasonally anoxic basin floor, we would expect enhanced dissolution to occur in the seaway margins, likely in the mid- to outer-shelf setting (Cherns et al. 2011). We would expect to see enhanced preservation in the anoxic basin center, as an aragonitic skeleton residing on the surface sediment in an anoxic water column would not be susceptible to dissolution from H2S oxidation (Jordan et al. 2015); however, we would not expect to see abundant benthos in such a setting because of bottom-water toxicity. It is apparent this could result in spatially expansive zones with conditions predisposed for heightened aragonite dissolution and preservation (Fig. 1); it is important to note that we do not expect all aragonitic fauna to be missing from any region of the seaway—merely that a lower relative proportion of aragonitic mollusks be found, due to a reduced probability of an individual site recording their occurrence. How these hypothesized basin-margin to basin-center zones could influence long-term patterns of mollusk distribution, preservation, and recovery remains to be examined. As epicontinental seas contain the majority of our Phanerozoic fossil record (Allison and Wells 2006), it is imperative that we understand systematic biases that may specifically affect these settings.
Figure1.
Diagram showing potential model of spatial aragonite bias within the Western Interior Seaway. Within the outer shelf, preferential dissolution of aragonitic fauna is common, which has the potential to be expressed spatially. Within the basin center, anoxia limits benthic organism development, but allows for preservation of aragonitic material. Modified after Cherns et al. (2011).
Here we present a spatial investigation of aragonite dissolution within the Late Cretaceous Western Interior Seaway (WIS) of North America, using sampling probability estimates and multiple logisitic regression to evaluate patterns of spatial distribution in preserved calcitic and aragonitic fauna. We address two key questions: (1) Does aragonite bias exhibit systematic spatial variation across the seaway, and (2) if so, does this influence perceived patterns of diversity?
Materials and Methods
Time Intervals and Paleogeography
The two stratigraphic intervals or time slices (Cenomanian–Turonian and early Campanian) were selected: (1) because of purported dysoxic conditions within their duration; and (2) due to their differences in environment, oceanography, and preserved lithology, allowing for comparison of taphonomic regimes. The first interval covers the Cenomanian/Turonian boundary, spanning the Dunveganoceras pondi Zone to the Collignoniceras woollgari Ammonite Zone (∼94.7–93 Ma) (Cobban et al. 2006). The second interval spans the early Campanian, from the Scaphites leei III Zone to the Baculites obtusus Ammonite Zone (∼83.5–80.58 Ma) (Cobban et al. 2006). The geologic context of stratigraphic intervals is detailed in Supplementary Information 1.
A global atlas of 1:20,000,000 scale paleogeographic maps, compiled by GETECH PLC, formed the basis for new regional-scale, high-resolution interpretations for the selected time intervals. The original paleogeographic maps (Markwick 2007) are underpinned by the GETECH plate model (v. 1), which is outlined further in Supplementary Information 1. High-resolution mapping involved synthesis of stratigraphic, sedimentologic, and paleontological information to produce 1:5,000,000 scale paleogeographies with suggested paleobathymetry. A full list of decisions on paleogeographic reconstructions and key references for each time interval are provided in Supplementary Information 1.
Landward-to-basinward arrangements of a priori binned zones for each time slice were based on average paleobathymetry (Fig. 2). Bathymetric reconstructions were divided into four bins, each of which covers a specific interpreted depth range: nearshore (<50 m), proximal offshore (50–100 m), distal offshore (100–150 m), and basin center (>150 m). These designations were based on the previously constructed paleobathymetry for the WIS produced by Sageman and Arthur (1994), but match the paleobathymetry in our new maps and represent a reasonably high resolution without being compromised by large changes in shoreline position within our chosen time slices.
Distance from paleoshoreline zones (Supplementary Fig. S1) were constructed based on 50 km intervals from the time-averaged paleoshoreline position until reaching the basin center, with number of occurrences, collections, and total outcrop area plotted per zone. These were generated by constructing a fishnet of points in ArcGIS (ESRI 2010) using the Fishnet tool, which were selected by the Select By Location tool with increasing distance in 50 km intervals from the paleoshoreline: the positions of the most-basinward selected points were used for the bin boundary. Results for depth zones are used in the main body of this paper; distance from paleoshoreline zones are available in Supplementary Information 1 and Supplementary Figures S2, S4, S6, and S7.
Fossil Data Set
A presence-only fossil occurrence data set of bivalve and ammonite taxa was produced for the selected stratigraphic intervals, collated from personally provided digitized collections from the U.S. Geological Survey (USGS) and Smithsonian Museum of Natural History (NMNH) and downloads from the Paleobiology Database (PBDB; http://paleobiodb.org) and iDigBio ( http://www.idigbio.org). Each occurrence includes taxonomic and geographic locality data, an associated collection with lithologic and geologic information, and modern latitudinal and longitudinal coordinates. Data were extensively screened for problematic records and to ensure taxonomic validation (see Supplementary Information 2 for the latter).
The resultant Cenomanian–Turonian data set contains 5867 occurrences from 2409 localities, with 207 genera, 1549 species, and 3886 specimens identifiable to the species level. The early Campanian data set comprises 2544 occurrences from 1186 localities, recording 163 genera, 1405 species, and 1405 specimens identifiable to the species level. Generic-level taxonomic diversity was used for all tests; species-level results can be found in Supplementary Information 1 and in Supplementary Figures S3–S6. Full information regarding downloads, sources, and screening of data can be found in Supplementary Information 1, and the full data set can be found in Supplementary Information 2.
Mineralogy
Bivalve shells are a composite of layered mineral crystallites, which are sheathed by a refractory organic matrix of fibrous protein (Taylor 1969). As these mineral layers can be composed of both calcite and aragonite, variation in overall mineral composition must be taken into account when assigning a predominant mineralogy to a specific bivalve genus. Different scoring mechanisms have been adopted by previous workers to address this issue. Kidwell (2005) used a five-point decimal scoring system from entirely aragonitic (1) to entirely calcitic (3), with three permutations of mineralogy between. Crampton et al. (2006) adopted a simple and effective system of counting organisms as calcitic if they contained a calcitic element that would allow them to be identified to the species level. We utilize a combination of these approaches—organisms were scored using the system of Kidwell (2005) to maintain the maximum amount of data, but simplified into binary categories afterward based on whether they contained sufficient calcitic parts to enhance preservation potential. Note that we have not included either the inner myostracal layer or periostracum in our assignments of mineralogy.
Figure2.
Paleogeographic zoned maps of the Western Interior Seaway used in this study. Depth-based zones are designated as nearshore, proximal offshore, distal offshore, and basin center (Fig. 1). A, Paleobathymetric map of the Cenomanian–Turonian; B, paleobathymetric map of the early Campanian.
Information on shell composition was predominantly gathered from a personally provided data set from S. Kidwell (Kidwell 2005), with additional information from studies by Taylor (Taylor 1969; Taylor and Layman 1972), Majewske (1974), Carter (1990), Schneider and Carter (2001), Lockwood (2003), Hollis (2008), and Ros Franch (2009) and many papers focused on single genera or families. For genera for which information regarding shell mineralogy was not available, composition was assigned based on the dominant mineralogy of the family, as composition is highly conservative both among species within a genera and among genera within a family (Taylor 1969). In total, 124 bivalve genera were assigned a mineralogy, of which 41 (∼33%) were achieved using familial relation (Supplementary Information 2).
Life Habits
Life habits of bivalves were assembled to allow additional interrogation and interpretation of environmental and sampling regimes. Life habits were separated into the following categories: relation to substrate, mobility, and diet. Data for each genus of bivalve were primarily gathered from the NMiTA Molluscan Life Habits Database (Todd 2017) and the PBDB, with further data collected from the wider literature (Supplementary Information 2).
Outcrop Area
Relevant rock outcrop area was plotted per zone to evaluate broader-scale bias influencing patterns of fossil distribution. Outcrop areas for the selected time slices were generated by combining statewide digitized geologic maps from publicly available USGS downloads and selecting shape files that matched formations found within those time slices. Some state surveys grouped relevant formations with other partially contemporaneous formations that spanned multiple stages: we chose to include these designations in order to present the maximum possible sampling extent in terms of outcrop area. Outcrop was projected in ArcGIS (ESRI 2010) using the USA Contiguous Albers Equal Area conic projection, to minimize distortion of distances. Outcrop areas per zone were created by using the Intersect tool in the Geoprocessing toolbar in ArcGIS, and area (km2) calculated using the Calculate Geometry function in the attribute table. Outcrop was split into depth zones by using the Intersect tool in ArcGIS (ESRI 2010). Outcrop area for each zone was calculated by summing the total area of all outcrop polygons within that zone. Collections per zone were counted by exporting occurrences selected in zones in the seaway as shapefiles, then using the arcgisbinding package to view and organize the data in R v. 3.0.2 (R Core Team 2017).
Dominant Lithology
Each collection was assigned a dominant lithology to allow for comparative testing. If these data were not available, a lithology was assigned from the dominant lithology of the formation, with reference to USGS formation records. Collections were assigned one of the following lithologies (primarily based on original USGS records): siliciclastic mudstone, siliciclastic siltstone, siliciclastic sandstone, conglomerate, ironstone, calcareous mudstone and siltstone, marl, calcarenite, limestone, and chalk.
Range Size
If the presence of preferentially destructive zones is affecting the spatial distribution of aragonitic fauna, we might expect to see overall smaller range sizes for aragonitic organisms compared with calcitic organisms (Fig. 3). As such, range-size estimates were produced for calcitic and aragonitic bivalves and compared to test whether aragonite bias influenced perceived range of aragonitic organisms. Note that ammonites were excluded from this test due to the difference in life habit between them and bivalve fauna: ammonites have a pelagic to nektono-benthic mode of life (Ritterbush et al. 2014), while bivalves are predominantly epifaunal and infaunal.
Geographic locality data for the selected fauna were visualized in ArcGIS (ESRI 2010). Faunal occurrences were paleo-rotated using the Getech Plate Model to match the paleogeography of the appropriate stages of the Late Cretaceous. This ensures that tectonic expansion and contraction of the North American Plate from the Mesozoic to the Recent has a negligible effect on propagating estimation error in range-size reconstructions. Fossil occurrences were projected into ArcGIS using the using the USA Contiguous Albers Equal Area conic projection. A 10 km buffer was additionally applied to each occurrence point to control for any error in paleogeographic or present position of fauna. ArcGIS (ESRI 2010) was then used to construct convex-hull polygons for each taxon, and the Spatial Analyst tools from this software calculated the area of each reconstructed polygon. We did not account for landforms within the ranges of any organisms, and thus ignored their area when calculating overall area of ranges. Several vertices for range-size polygons appeared on what is classified as land within our paleogeographies; due to rapid changes in shoreline position within the WIS, we decided to keep using these fauna for range-size estimations. Myers and Lieberman (2011) showed that relative range sizes for vertebrates in the WIS were not overly affected by resampling occurrence points—consequently, we have not carried out a similar test for this study. Comparisons between the ranges of aragonitic and calcitic fauna were carried out using the Wilcoxon-Mann-Whitney test with continuity correction (Brown and Rothery 1993). Geographic range data for all applicable taxa are provided in Supplementary Information 2.
Sampling Probability and Multiple Logistic Regression
To further observe differences between aragonitic and calcitic organisms throughout the seaway, we employed a modified version of the sampling probability method used by Foote et al. (2015; after Foote and Raup 1996). In this method, the sampling probability of a time bin was generated by compiling a list of all fauna with originations older than that bin and extinctions younger than that bin, and then dividing the total number of species found within the bin by that figure. This allows for a sampling probability to be estimated on a per-bin, per-group basis. Here we devised three variants on this method for application in the spatial realm. It should be made clear that the modified methods used in this work come with the caveat that in the spatial realm it is impossible to know whether a species was present in a precise location in the past: for instance, if zones A, B, and C are designated with increasing distance away from a paleoshoreline, it cannot be assumed that because an organism exists in zones A and C that it was ever present in zone B. Consequently, the probabilities generated from the methods described below are relative and cannot be taken as “true” probabilities. However, the methods used were designed to be as inclusive as possible and to deliver a strongly conservative estimate of true sampling probabilities between groups; these methods therefore provide a useful estimate on the relative likelihood of sampling aragonitic or calcitic fauna. Furthermore, sampling probabilities through time based on regional studies such as those used by Foote at al. (2015) rely on the assumption that groups were not genuinely absent from the study region at a particular time and that other geographic variables do not have an effect—as such, the use of these metrics to evaluate the distribution of fauna across the WIS is validated.
Figure3.
Diagram showing potential model of apparent range-size reduction due to spatially variable aragonite preservation. Assuming that calcitic and aragonitic species of bivalve were both living at four separate localities, but aragonitic dissolution strongly influenced one of those locations (A), the resulting convex hull for the aragonitic fauna drawn from surviving fossil occurrences would likely be smaller than that for the calcitic organisms (B).
Three methods were devised for dealing with the issue of unknown “correct” distribution of species across the seaway and to correct for differences in the number of collections between zones: method 1 finds two bins either side of the current bin and generates a list of the total number of possible species across those five bins; method 2 finds all formations that appear in the selected bin that contain specimens of the selected group (e.g., calcitic bivalves), and then finds the total number of species for that group from those formations; method 3 finds all formations in the current bin and the two adjacent bins that contain specimens of the selected group, and subsequently finds the total number of species from those formations. For all three methods, the total number of sampling opportunities per bin was generated by multiplying the number of potentially recoverable species by the number of collections to standardize for differences in collecting intensity. The low number of depth-based bins could potentially result in flattening the curve of sampling probability using method 3, and thus method 2 is employed in the main body of this paper for depth-based results.
Figure4.
Plots of generic-level M2 sampling probabilities for the Cenomanian–Turonian (A, C, E, G) and lower Campanian (B, D, F, H) time slices across depth zones, split into carbonate and siliciclastic sampling opportunities. All results are plotted with percentage of carbonate collections per depth zone. A, Cenomanian–Turonian generic-level sampling probability, plotted with percentage of carbonate collections per depth zone; B, lower Campanian generic-level sampling probability, plotted with percentage of carbonate collections per depth zone; C, Cenomanian–Turonian ammonite sampling probability; D, lower Campanian ammonite sampling probability; E, Cenomanian–Turonian aragonitic bivalve sampling probability; F, lower Campanian aragonitic bivalve sampling probability; G, Cenomanian–Turonian calcitic bivalve sampling probability; H, lower Campanian calcitic bivalve sampling probability.
To determine the primary controls on sampling probability between the two stages, we used multiple logistic regression, coding sampling opportunities as the response variable and mineralogy, lithology, life habits (mobility, relation to substrate, and feeding style), and depth zone as the predictor variables. Multiple logistic regression allows for the use of binomial nominal values by using the odds ratio, a measure of the relationship between the odds of an outcome, in this case sampled (1) or not sampled (0), along with multiple potentially explanatory ecological or physiographic variables. A full model is generated that incorporates all potential variables, and a null model is defined that includes none. Stepwise addition or deletion from the null or full models, respectively, and analysis in the change of likelihood and of respective Akaike information criterion (AIC) scores contribute to a final predictive model of explanatory variables and respective statistical significance.
Sampling opportunities were tabulated as the presence or absence of each recoverable genus per collection, per depth zone. Each sampling opportunity was assigned a lithology based on collection lithology, as well as all ecological attributes related to that genus. To test for multicollinearity between variables, correlation tests were run using Spearman's rank correlation using the Performance Analytics package in R. Explanatory variables that showed a strong (>0.7) statistically significant correlation were excluded from further analysis (Supplementary Information 1).
Interaction terms were also added to explore the possibility of multiple confounding factors and increased model complexity. These terms were restricted to a combination of lithology and mineralogy, so as to test for specific interactions between the two (e.g., whether preservation of aragonite was specifically enhanced within limestones). We also partitioned the data to be able to fully explore the influence of various contributing factors on sampling probability per depth zone and include all organisms in the data (ammonites were excluded from analyses involving life habits, as discussed previously). Both effect sizes of individual factors and AIC values of models are presented for statistically significant interactions. All methods were written and implemented using R.
Occurrences, Raw Diversity, and Shareholder Quorum Subsampling
To establish the potential influence of aragonite bias on diversity of shelly taxa, total occurrences of organisms were counted per zone using the Select By Location tool in ArcGIS (ESRI 2010) and were used to generate landward-to-basinward profiles of raw occurrences and raw and subsampled diversity estimates. Shareholder quorum subsampling (SQS; Alroy 2010), a method for standardizing taxonomic occurrence lists based on an estimate of coverage, was implemented in R using script provided by J. Alroy (personal communication 2015) for each faunal group. Calcitic and aragonitic groups were evaluated for statistically significant differences using the chi-squared test for nonrandom association (Brown and Rothery 1993). All statistical tests were implemented in R. Results pertaining to patterns within raw occurrences can be found within Supplementary Information 1 and Supplementary Figure S1.
Results
Sampling Probability
Cenomanian–Turonian.—For generic-level sampling probability (Fig. 4A), aragonitic bivalves and ammonites show a similar trend for the first three depth zones. After this, sampling probability drops to 0 for aragonitic bivalves (as none were recovered), while it increases to a peak for ammonites. Calcitic bivalves record a higher sampling probability than ammonites or aragonitic bivalves in all zones and show a basinward increase in sampling probability.
Campanian.—In the lower Campanian (Fig. 4B), ammonites have the highest sampling probabilities, showing a level trend across the seaway with a pronounced trough in the distal offshore. Aragonitic bivalves record a relative high sampling probability in the nearshore, followed by a sharp decline for both proximal and distal offshore zones and an increase toward the basin center. Calcitic fauna have a consistently higher sampling probability than aragonitic bivalves, but lower than ammonites; they also show a level trend across the seaway, experiencing a peak in the distal offshore.
Sampling Probability between Lithologies
Cenomanian–Turonian.—For the Cenomanian–Turonian (Fig. 4C,E,G), ammonites show the same trends and relatively little difference in absolute values between carbonate and siliciclastic sampling opportunities; the greatest difference appears in the basin center, where sampling probability is higher in carbonates. Aragonitic bivalves show a much larger difference, with siliciclastic opportunities scoring consistently higher than carbonate opportunities, even during the large decline within the proximal offshore. Calcitic bivalves show virtually no difference in sampling probability until the basin center, where sampling probability within carbonate sampling opportunities increases substantially.
Table 1.
Table for multiple logistic regression results for all bivalves within the Cenomanian–Turonian across the whole seaway, using the model with lowest Akaike information criterion (AIC) score. M, mineralogy; D, diet; L, lithology; Z, depth zone; DOS, distal offshore; NS, nearshore; POS, proximal offshore. *indicates statistical significance at p = < 0.05.
Campanian.—For the Campanian (Fig. 4D,F, H), siliciclastic opportunities of ammonites score higher than carbonate, except for within the nearshore. Aragonitic bivalves are not sampled within carbonate collections in either the nearshore, distal offshore, or basin center; their sampling probability curve is virtually entirely made by appearances in siliciclastic sampling opportunities. Calcitic bivalves show a decoupled trend between lithologies, with carbonate sampling opportunities showing higher on average sampling probabilities that increase toward the basin center, compared with the fairly low-scoring, level trend in siliciclastic.
Multiple Logistic Regression
Results of the logistic regressions are shown in Tables 1–8 and summarized in Figure 5. When interpreting these, note that calcitic mineralogy is compared with aragonitic mineralogy, so that positive regression coefficients indicate greater odds of sampling calcite. As lithology has multiple parameters, these were compared against the relative sampling probability of mudstone, which is used as a baseline. We are primarily interested in reporting effect sizes, which are gauged by the magnitude of regression coefficients.
Table 2.
Table for multiple logistic regression results for all bivalves within the Cenomanian–Turonian for each depth zone, using the models with lowest Akaike information criterion (AIC) scores. M, mineralogy; D, diet; L, lithology. * indicates statistical significance at p = <0.05.
AIC scores are utilized in choosing ideal model fit when comparing models with and without two-way interactive terms (a combination of effects between explanatory parameters: e.g., the relative odds of sampling calcitic fauna within a specific lithology), with lower scores indicating a better model fit. Only models with the lowest AIC scores are presented here, and we only report factors with statistically significant results (p < 0.05); full results can be found within Supplementary Information 2.
Table 3.
Table for multiple logistic regression results for all organisms (including ammonites) within the Cenomanian–Turonian across the whole seaway, using the model with lowest Akaike information criterion (AIC) score. M, mineralogy; L, lithology; Z, depth zone; DOS, distal offshore; NS, nearshore; POS, proximal offshore. * indicates statistical significance at p = < 0.05.
Cenomanian–Turonian.—Mineralogy, lithology, feeding style, and depth zone all influence the preservation potential of fauna in the seaway (Table 1); lower AIC scores when an interactive term is added suggest this provides a better model fit than when it is excluded. The odds of sampling calcitic fauna are shown to be 4.6 times (the exponential of the coefficient; 1.52) higher that of aragonitic fauna, with analysis of variance (ANOVA) results showing mineralogy contributing the most toward deviance from the null model. Limestone environments are shown to be detrimental to the sampling probability of fauna, whereas sandstones and siltstone enhance sampling probability. The positive interaction between mineralogy and limestone lithologies shows that aragonitic fauna have comparatively strongly reduced odds of being sampled within limestone environments. All depth zones are shown to have decreased sampling probability compared with the basin center, with nearshore and proximal offshore zones showing the worst sampling potential. Chemosymbiont deposit feeders are shown to have an increased preservation potential compared with other feeding styles.
Table 4.
Table for multiple logistic regression results for all organisms (including ammonites) within the Cenomanian–Turonian for each depth zone, using the models with lowest Akaike information criterion (AIC) scores. M, mineralogy; L, lithology. *indicates statistical significance at p = < 0.05.
Table 5.
Table for multiple logistic regression results for all bivalves within the lower Campanian across the whole seaway, using the model with lowest Akaike information criterion (AIC) score. M, mineralogy; D, diet; L, lithology; Z, depth zone. * indicates statistical significance at p = < 0.05.
We additionally partitioned the data into each depth zone to test for differences with increased bathymetry across the seaway (Table 2). The nearshore zone exhibits an increase in the odds of sampling calcitic fauna, although this effect is reduced compared with results across the whole seaway. Sandstones are also shown to exhibit increased sampling probability. The proximal offshore shows a significant increase in the odds of sampling calcitic bivalves relative to aragonitic bivalves (6.17 compared with 1.88 for the nearshore), as well as increased sampling probability in marl depositional environments and for chemosymbiotic deposit feeders. Limestone negatively impacts the sampling probability of bivalves; the positive interaction between calcite and limestone consequently suggests that this negative impact is related to the sampling probability of aragonitic bivalves. The distal offshore shows a similar pattern, although the relative odds of each are reduced compared with the proximal offshore. The basin center shows increased odds of sampling bivalves within calcarenite, calcareous shale, and marl environments, but no other statistically significant terms.
Table 6.
Table for multiple logistic regression results for all bivalves within the lower Campanian for each depth zone, using models with lowest Akaike information criterion (AIC) scores. M, mineralogy; D, diet; L, lithology; Z, depth zone. *indicates statistical significance at p = < 0.05.
We also assessed depth zones for the inclusion of all organisms (Table 3). When ammonites are included, the odds of sampling aragonitic fauna increase (calcitic bivalves show odds of 2.1 higher sampling probability). Sandstone shows reduced odds of sampling any fauna, the opposite of previous results. The interaction between mineralogy and lithology shows increased sampling probability of calcitic organisms within limestones, sandstones, calcarenites, and calcareous mudstones, suggesting this effect is predominantly produced by the addition of ammonite fauna.
Table 7.
Table for multiple logistic regression results for all organisms (including ammonites) within the lower Campanian across the whole seaway, using model with lowest Akaike information criterion (AIC) score. M, mineralogy; L, lithology; Z, depth zone; DOS, distal offshore; NS, nearshore; POS, proximal offshore. * indicates statistical significance at p = < 0.05.
When zones are assessed independently (Table 4), nearshore sampling probabilities are only controlled by mineralogy, although again with lower odds than reported elsewhere (1.56). In the proximal offshore, results show an increased sampling probability of calcitic fauna within sandstones and calcareous mudstones. The distal offshore also shows strong interactions between sampling probability of calcitic fauna and lithology, with strongly positive coefficients for sandstone, limestone, calcareous shale, and calcarenite two-way interactions. Overall, the sampling probability of calcite compared with aragonitic fauna is high, although reduced compared with the proximal offshore. Within the basin center, mineralogy is not listed as a statistically significant interactive term on its own, but calcitic fauna exhibit increased sampling probability for interactive terms with calcarenites, calcareous mudstones, limestones, and marls.
Campanian.—Models for all bivalves in the Campanian (Table 5) show comparatively few statistically significant contributors to sampling probability. By comparison with the Cenomanian, bivalve samples from the Campanian are only weakly influenced by mineralogy (showing odds of 2.16 increased likelihood of sampling calcitic organisms). Additionally, only sandstone and interactions between sandstone and limestone with calcitic organisms are shown to exert any other influence on sampling probability.
Table 8.
Table for multiple logistic regression results for all organisms (including ammonites) within the lower Campanian for each depth zone, using models with lowest Akaike information criterion (AIC) scores. M, mineralogy; L, lithology; NULL, null model.
Figure5.
Graph summarizing multiple logistic regression model results (T1–T8 = Tables 1–8). Final models are presented within columns, whereas factors are presented along rows for those models. Results for factors are presented as either comparatively positive or negative odds of sampling compared with reference factor, not statistically significant, or not included in the final model. The following factors are used as a baseline for comparison: mineralogy, aragonite; lithology, mudstone; lithology with mineralogy, aragonite:mudstone; zone, basin center (BC); diet, carnivore. Note that the magnitude of regression coefficients is not presented within this graph. NS, nearshore; POS, proximal offshore; DOS, distal offshore.
This trend continues when partitioning the bivalve data into depth zones (Table 6). The nearshore zone has no statistically significant individual factors contributing to sampling probability. The proximal offshore includes statistically significant effects due to mineralogy and lithology, particularly limestones and calcareous mudstones, where sampling probability is enhanced. Mineralogy, sandstone, and the interaction between mineralogy and sandstone are reported as statistically significant factors for the distal offshore; mineralogy has a relatively high positive coefficient (odds of 3.16 in favor of calcitic organisms). Sampling probability is enhanced in sandstones overall, but negatively influences the odds of recovering calcitic organisms: it therefore follows that aragonitic bivalves show particularly enhanced sampling within sandstones. Model results for the basin center suggest that only calcareous shale has a statistically significant positive impact on sampling probability.
When all organisms are assessed (Table 7), mineralogy and depth zone are the only contributors to the full model. Surprisingly, aragonitic organisms have a higher sampling probability than calcitic organisms using the full model, with mineralogy only contributing to a very small amount of deviance from the null ANOVA model. As this result is not observed when assessing bivalve fauna, it is likely that ammonite occurrences are principally contributing to this effect.
Depth zones were also evaluated for all organisms (Table 8). Only the proximal offshore supported a model other than the null, which reported mineralogy as a contributing factor; unusually, calcitic fauna are shown to have a reduced sampling probability compared with aragonitic fauna.
Range Size
Cenomanian–Turonian.—Box plots were generated on a log scale to show differences in mean range sizes between calcitic and aragonitic organisms (Fig. 6A). There is a visible difference in variability of range size between groupings, with calcitic fauna showing an average larger range than aragonitic fauna. The Wilcoxon-Mann-Whitney test also showed a statistically significant difference between the range sizes for the two groups (p-value = 0.00405), with a reported difference in median range size of 48,694 km2. As sample size varied between the groups, resampling measures were carried out to test the accuracy of these results. A randomized bootstrap with replacement calculating the difference between the means of range sizes was implemented 10,000 times in R (Fig. 6C). Our recorded difference in the mean was shown to have an associated p-value of 0.0172, showing statistical significance.
Campanian.—Box plots were generated to show differences in mean range sizes between early Campanian calcitic and aragonitic organisms (Fig. 6B). Calcitic bivalves show higher variability in mean range size than aragonitic bivalves. However, the Wilcoxon-Mann-Whitney test showed no statistically significant difference between the two groupings (p-value = 0.504) with a recorded difference in median range size of 13,540 km2, and a randomized bootstrap (Fig. 6D) with replacement recovered an associated p-value of 0.1527 (non-statistically significant).
Raw Diversity and SQS
Cenomanian–Turonian.—Within the Cenomanian–Turonian, broadly similar patterns of diversity occur in all groups (Fig. 7A,C)—peak diversity is within the distal offshore, with lowest values in the nearshore and basin center. Calcitic bivalves show proportionally enhanced diversity in the proximal offshore compared with the other faunal groups. These patterns closely align with the number of collections within each zone, but show limited similarity to zoned outcrop area.
Subsampled ammonite and calcitic bivalve diversities show a broadly similar pattern to their raw taxic diversity signals (Fig. 8A,E). The record of aragonitic bivalves (Fig. 8C) is too poor to resolve subsampled diversity for the basin center; however, a slight decline in subsampled generic richness exists in the proximal offshore.
Campanian.—Calcitic bivalves and ammonites exhibit a similar pattern in diversity (Fig. 7B,D) although the latter show an increase in the proximal offshore. Aragonitic bivalve diversity has a similar peak in the proximal offshore but declines toward the basin center. None of these trends show similarity to the distribution of collections or outcrop area throughout the seaway.
When subsampled (Fig. 8B,D,F), calcitic and aragonitic bivalves are most diverse within the proximal offshore, falling to relative lows within the distal offshore and basin center. Ammonites are most diverse in the nearshore, followed by a decline to a flat profile.
Figure6.
Range-size plots for the Cenomanian–Turonian and lower Campanian. A, Cenomanian–Turonian box plots of range size for both aragonitic bivalves and calcitic bivalves on log scale; B, lower Campanian box plots of range size for both aragonitic bivalves and calcitic bivalves on log scale; C, randomized bootstrap for Cenomanian–Turonian mean range sizes—recorded difference in the mean is shown to be statistically significant; D, randomized bootstrap for lower Campanian mean range sizes—recorded difference in the mean is not shown to be statistically significant.
Discussion
Sampling Probability and Multiple Logistic Regression
Our results from estimations of sampling probability and subsequent multiple logistic regression suggest that aragonite bias may be present within distinct depth zones of the seaway during the Cenomanian–Turonian. Mineralogy has a strong and statistically significant impact on sampling probability within the proximal and distal offshore bathymetric zones and shows the highest contribution to deviance from the null model. This is further supported by the fact that while all aragonitic taxa have lower sampling proportions overall, both aragonitic bivalves and ammonites disproportionally decrease in sampling probability within the proximal offshore compared with calcitic bivalves. Ammonites, while still showing reduced sampling probability compared with calcitic fauna, are more likely to be sampled than aragonitic bivalves; a potential explanation for this difference could be that aragonite dissolution acts differently upon ammonites compared with bivalves. Body sizes of ammonites and bivalves differ, with ammonites generally having larger forms (Jablonski 1996). This has been known to influence preservation potential and the extent of aragonite dissolution: Wright et al. (2003) showed that ammonites are affected less severely than aragonitic bivalves by early-stage aragonite dissolution, often exhibiting poor preservation rather than complete removal. Our results have the potential to be partially related to this effect.
Figure7.
Plots of generic-level diversity plots for the Cenomanian–Turonian and lower Campanian within depth zones, plotted with number of collections and outcrop area. A, Generic diversity and number of collections for the Cenomanian–Turonian; B, generic diversity and number of collections for the lower Campanian; C, generic diversity and outcrop area for the Cenomanian–Turonian; D, generic diversity and outcrop area for the lower Campanian.
Aragonitic bivalves have lower absolute sampling probabilities in carbonate environments than in siliciclastic environments, supporting the results of Foote et al. (2015). However, when examining the proximal offshore zone, we can see that sampling probability within siliciclastic lithologies falls dramatically. As this zone records the largest difference in odds of sampling between calcitic and aragonitic taxa, it can be argued that aragonite bias can influence fauna within siliciclastic deposits in epicontinental seas, in contradiction to Foote et al. (2015). The absolute sampling proportions of calcitic bivalves remain relatively consistent (at ∼2% of genera per collection) throughout the seaway until the basin center, where they increase dramatically within carbonates compared with siliciclastics. Foote et al. (2015) reported that calcitic organisms experienced higher sampling probabilities in carbonate-rich intervals, which is especially enhanced in limestones. As carbonates make up ∼93% of total sampling opportunities within this zone, our results align fairly closely with previous findings. While Foote et al. (2015) singled out lithology as an important factor for aragonite dissolution, they did not investigate whether differences in grain size significantly influenced results. Within this study, sandstone and siltstone are consistently shown to have better odds at preserving aragonitic fauna than mudstone. This is unsurprising, considering that coarser, oxidized sediments are likely to contain lower quantities of organic matter than finer sediments, and thus provide less material for the microbial decay that ultimately controls the dissolution of aragonite within the TAZ (Cherns et al. 2008). However, siltstone appears to have higher odds than sandstone, potentially a reflection of increased quality of preservation in lower-energy settings. It should be noted, however, that only a few models include both siltstone and sandstone and therefore allow for comparison of sampling probabilities of these two lithologies.
Figure8.
Plots of generic-level SQS results for depth zones in the Cenomanian–Turonian and lower Campanian, set at 0.4, 0.5, and 0.6 quora. A, SQS results for ammonites in the Cenomanian–Turonian; B, SQS results for ammonites in the lower Campanian; C, SQS results for aragonitic bivalves in the Cenomanian–Turonian; D, SQS results for aragonitic bivalves in the lower Campanian; E, SQS results for calcitic bivalves in the Cenomanian–Turonian; F, SQS results for calcitic bivalves in the lower Campanian.
Potential ecological signals can also be parsed from the results reported here. Within the Cenomanian–Turonian data set, odds of sampling chemosymbiont deposit feeders within the proximal offshore were higher than for sampling other bivalves, forming a statistically significant part of the final model and accounting for the second highest deviance from the null model. Chemosymbiosis in bivalves occurs in a range of environments to cope with life in sulfide-rich environments, typically at deep sea vents or in sediments at the oxic/anoxic interface (Cavanaugh 1994). Combined with evidence for poor sampling probability of aragonitic fauna in siliciclastic lithologies, this lends credence to the likelihood of fluctuating benthic oxygen conditions within the proximal offshore, ideal for preferential aragonite dissolution. More broadly, several previous works have suggested that aragonite bias strongly influences perceived trophic communities within molluscan fauna, favoring preservation of specific life habits (Cherns et al. 2008; Cherns and Wright 2009). Unfortunately, very few statistically significant life habit factors contribute to our final models (Fig. 5), and thus we cannot draw any conclusions regarding preservational shifts in trophic structure. In the basin center, ammonites are more likely to be sampled compared with other organisms. This confirms expectations of enhanced preservation within a predominantly anoxic water column, where dissolution and predation have reduced impact on the removal of fauna emplaced by pelagic fallout (Jordan et al. 2015).
Within the Campanian, there is a somewhat contradictory pattern. Multiple logistic regression results show that mineralogy only has a strong, statistically significant impact on relative sampling odds when assessing bivalves within the proximal and distal offshore bathymetric zones, with only the latter showing a strong deviation from the null model in ANOVA results. When ammonites are added, the odds of sampling aragonitic fauna are actually higher than for calcitic organisms within the proximal offshore, and all other zones show no statistically significant contributions from mineralogy. This is reinforced when one considers the absolute proportions of mineralogies sampled: ammonites exhibit the highest overall sampling probability among fauna. A potential cause of this contradiction is preferential sampling bias. Ease of collecting and human interest can result in skewed sampling effort and intensity, potentially inflating (Foote and Sepkoski 1999) or reducing (Lloyd and Friedman 2013) the published records of certain taxa, locations, and time periods above others. The WIS has long been known for its abundance and diversity of ammonite fauna, and consequently ammonites have been used for systematic biostratigraphic correlation since the 1930s (Stephenson and Reeside 1938). An intensive effort to collect ammonites for stratigraphic purposes was carried out by a selection of workers through the latter half of the twentieth century to the present day (Scott and Cobban 1959; Gill and Cobban 1973; Cobban and Hook 1984; Cobban et al. 2006; Merewether et al. 2011). Consequently, it is likely that records for biostratigraphically important organisms have been overinflated compared with other mollusks and between localities. By comparing previously existing collections and newly collected records for the upper Cenomanian Sciponoceras gracile Zone (now the Vascoceras diartianum Zone and the Euomphaloceras septemseriatum Zone; Cobban et al. 2006), Koch (1978) showed that ammonites were better studied and more commonly reported than bivalve fauna. Parts of these collections have made up the majority of the publicly available records of fossil occurrences within the WIS that are used in this study. As such, it is possible that ammonites are overrepresented in the early Campanian data set and are skewing perceived results. However, it is still possible to suggest that a suppressed expression of spatial aragonite bias occurs in the distal offshore, albeit at reduced levels in comparison to the Cenomanian–Turonian interval.
Range Size
Range-size results reported a difference between calcitic and aragonitic bivalves across the two time intervals studied, with aragonitic fauna showing a significantly smaller range size during the Cenomanian–Turonian but not the Campanian. This variation is also expressed spatially (Fig. 9). Within the Cenomanian–Turonian time slice, aragonitic geographic ranges (Fig. 9A) are generally restricted to the western and northern edges of the seaway in comparison to calcitic geographic ranges, which extend farther to the center of the basin, as well as to the east and south (Fig. 9C). This same pattern is slightly different in the early Campanian interval (Fig. 9B,D); while aragonitic fauna still show a limited range, the difference between both bivalve groups is less pronounced. This pattern also matches with the distribution of carbonate deposition within the WIS: the Cenomanian–Turonian interval experienced widespread carbonate sedimentation—in the form of the Greenhorn Limestone Formation—in the basin center (Miall et al. 2008), while deposition in the basin center transitioned from limestones of the Niobrara Formation to the siliciclastic mudstones of the Pierre Shale in the early Campanian (McGookey et al. 1972; Da Gama et al. 2014). As our results confirm that carbonate environments can exacerbate the effects of aragonite dissolution, it is possible that the differences between the Cenomanian–Turonian and the Campanian are partially driven by the enhanced effects of aragonite bias in carbonate-rich environments, resulting in a lowered sampling probability within carbonate-dominated localities.
Occurrence and Diversity Results
Overall, there is some evidence of aragonite dissolution influencing patterns of pure occurrences and taxonomic and subsampled diversity for aragonitic fauna, as previously hypothesized. In the Cenomanian–Turonian, aragonite bias is most pronounced within the proximal offshore bathymetric zone, with a lesser impact within the distal offshore zone. While all fauna show a close correlation to collection counts for depth zones, both aragonitic and calcitic fauna deviate from this correlation in the proximal offshore zone, recording lower raw occurrences and diversity. The same is broadly observed in the Campanian: maximum disparity of sampling probability between calcitic and aragonitic fauna is observed within the distal offshore zone, where aragonitic occurrences and raw taxic diversity show a noticeable decline and subsequent deviation from sampling proxies. Foote et al. (2015) reported similar results when comparing sampling-corrected results to ones that previously displayed the proportion of aragonitic taxa (Crampton et al. 2006) and concluded that similarities existed between sampling probabilities and relative proportions of aragonitic species.
Despite the potential relationships discussed earlier, we cannot report conclusive evidence for aragonite bias influencing the sampled diversity of molluscan fauna within the WIS. This aligns with other recent studies showing that despite evidence of widespread aragonite dissolution during early shallow diagenesis, perceived diversity is not largely affected by these processes (Behrensmeyer et al. 2005; Kidwell 2005; Crampton et al. 2006; Hsieh et al. 2019). Hence, we must additionally look at external influences that might capture, enhance, or control the distribution of aragonitic faunas that would otherwise be lost to preferential dissolution.
Known human influences have potentially contributed to the suppression of aragonite bias on a spatial scale. While the extent to which aragonite dissolution may have influenced our perceived record of molluscan diversity within the WIS is unclear, it is apparent that these records closely correlate with established sampling proxies. Results of Spearman rank correlation tests of occurrences and raw taxic diversity against sampling proxies for distance from paleoshoreline zones (Table 9) correlate strongly and significantly. It is clear that broader-scale sampling trends related to collector effort strongly influence the pattern of faunal distribution across the seaway, potentially overwriting the effects of aragonite dissolution.
Figure9.
Paleogeographic maps shown with range sizes of calcitic and aragonitic bivalves for both time slices. A, Aragonitic bivalve range sizes for the Cenomanian–Turonian; B, aragonitic bivalve range sizes for the lower Campanian; C, calcitic bivalve range sizes for the Cenomanian–Turonian; D, calcitic bivalve range sizes for the lower Campanian. P-N indicates direction of paleo-north.
While there have been many cases of preferential aragonite dissolution within local studies, aragonitic molluscan fauna are relatively well represented in the global fossil record (Harper 1998). This paradox suggests that processes must occur that capture records of molluscan fauna at a higher frequency than the frequency at which they are capable of being destroyed. Cherns et al. (2008, 2011) describe “taphonomic windows” as events in the fossil record that capture an unbiased view of aragonitic faunas that have escaped preferential dissolution, and these authors detail many examples that may have operated within the WIS. One such window that is prevalent within the WIS is concretions, sedimentary mineral masses of varying chemical composition that often form at shallow burial depths early in diagenesis when mineral cement precipitates locally during lithification (Berner 1968; McCoy et al. 2015). These have the potential to preserve three-dimensional fossilized remains, often in exquisite detail (Dean et al. 2015; Korn and Pagnac 2017). Concretions are also a characteristic mode of molluscan occurrences within the WIS, with fossil-bearing concretions found commonly throughout the seaway (Landman and Klofak 2012); as such, they could further contribute to a potential anthropogenic bias in that they provide easily spotted locations to find fauna in otherwise barren strata (such as the Pierre Shale), skewing collection intensity between localities with concretions and those without. However, only ∼3% of USGS collections were obtained by selective collecting (Koch 1980), and as USGS records make up ∼55% of our finished data set, this suggests that sampling intensity bias might be partially mitigated. Sediment accumulation rate could exert a large influence on the potential for preferential aragonite dissolution to affect spatial zones of the seafloor. If sediment accumulation rates were low, fauna would remain within the TAZ for an extended period of time, and thus are more likely to be removed through physical reworking, bioerosion, and enhanced dissolution (Cherns et al. 2011). In contrast, if sediment accumulation rates were high, fauna are likely to have been rapidly buried and thus have escaped into the sub-TAZ region, where vulnerable bioclasts are likely to be stabilized by shallow-burial diagenesis (Melim et al. 2002, 2004). Sediment accumulation rates within the WIS varied both longitudinally within a stratigraphic interval (with higher sediment accumulation rates toward the western paleoshoreline) and with increased bathymetry in a single location (Arthur and Sageman 2005): accounting for this potential influence is problematic, and the extent of its effects is ambiguous. The result of these factors is a potential suppression of the spatial influence of aragonite dissolution bias on recorded faunal diversity within the WIS.
Table 9.
Spearman's rank correlations between generic diversity of faunal groups and various sampling proxies for distance from paleoshoreline zones within the Cenomanian–Turonian and lower Campanian. Sig., significance; * indicates statistical significance at p = < 0.05.
Spatial Scale and Influence of Bias
The issue of scale is key to understanding the spatial impact of aragonite dissolution (Kosnik et al. 2011). Foote et al. (2015) recorded preferential aragonite bias within carbonate-rich environments on the regional spatial (∼106 km2) and stage-level temporal (1–10 Myr) scales. However, others (Behrensmeyer et al. 2005; Kidwell 2005; Kiessling et al. 2008; Kosnik et al. 2011) using global-scale data have reported negligible influence of shell mineralogy on temporal trends or frequency of occurrences. Foote et al. (2015) reported three key differences between previous studies and their work: higher taxonomic level of occurrences, larger time bins, and the use of global data. These factors were inferred to “even out” spatial and temporal variations in sampling, mitigating the influence and effect of locally variable biases inherent to the fossil record. Foote et al. (2015) further suggested that as their taxonomic and temporal scales were consistent with previously published work, an increase in spatial scale may prove the most influential factor in demoting the influence of aragonite dissolution.
This result can be easily translated into the spatial expression of aragonite bias by comparing its potential on alpha (within-site), beta (between-site), and gamma (global) diversity. At the alpha level, the impact of aragonite bias on a single species will be at its most severe, particularly within single-bed assemblages (Wright et al. 2003; Bush and Bambach 2004; Cherns et al. 2008, 2011). However, at gamma levels of diversity, the probability of not recording an individual drops substantially due to the number of possible localities to sample from, where various taphonomic windows may result in aragonite preservation. As such, an increased number of localities in a spatial setting are likely to partially obscure localized aragonite dissolution. As we recorded an impact on zoned sampling probabilities and range size of aragonitic fauna in the WIS, but could not conclusively prove an influence on total diversity estimates, our data support the suggestions of Foote et al. (2015) that spatial scale is a dominant factor on the severity of aragonite bias.
While preferential aragonite dissolution is unlikely to influence diversity on a global scale, this study has shown that it has the capacity to govern the sampling probability of a species in geographic space, and thus can influence the “variation” definition of beta diversity (Anderson et al. 2011). As the preferential dissolution of aragonite is a process that is exacerbated by certain environments (Foote et al. 2015), its influence will impact localities with different environmental conditions to differing extents—a species will be lost at one site and recorded at another. Our results confirm this, showing aragonite bias has an effect on observed diversity between locations, at least during times of widespread carbonate deposition.
Consequently, when looking at the spatial signal of aragonite dissolution as a whole, we can see a sliding scale of influence: strong, environmentally dependent impact on alpha diversity; a potentially large influence on beta diversity; and a negligible impact on gamma diversity. Bush et al. (2004) grouped biases affecting spatially organized biodiversity in similar alpha, beta, and gamma levels, with alpha biases influencing within-site diversity and beta and gamma arising from failure to sample all available habitats or environments within a region. While it was noted in this study that taphonomic effects were not included in this definition, this system can be modified in the light of our results. Aragonite bias, while operating at an alpha bias (local) level, evidently has the capacity to systematically influence estimates of beta diversity. As such, the influence of some taphonomic biases may be dependent on the spatial scale at which they are observed. This is an important consideration for studies of the spatial distribution of bias in the fossil record (Barnosky et al. 2005; Vilhena and Smith 2013; Benson et al. 2016; Close et al. 2017), and for paleobiogeographic studies in general.
Conclusions
A multifaceted approach shows that preferential aragonite dissolution is spatially variable and impacts the relative likelihood, absolute sampling probabilities, and range sizes of aragonitic organisms within the Cretaceous WIS of North America for a time interval that straddles the Cenomanian/Turonian boundary. A similar but reduced effect is additionally observed within an early Campanian time interval. A combination of depositional lithology (a limestone-dominated basin within the Cenomanian–Turonian compared with a more siliciclastic setting in the early Campanian) and an anoxic basin center are hypothesized as drivers for this effect.
Carbonate environments enhance the effects of aragonite dissolution and the preservation of calcitic organisms, as has been previously demonstrated. However, in contrast to previous studies, siliciclastic environments are also shown to be influenced by preferential aragonite dissolution.
While similarities are observed between faunal distribution and absolute sampling probabilities, we cannot conclusively say that aragonite dissolution has influenced perceived diversity of mollusks within the WIS. Taphonomic windows act to preserve records of organisms that would otherwise be lost. Other anthropogenic and geologic biases appear to have a more obvious effect on the molluscan record, and likely mask the influence of aragonite dissolution.
While aragonite bias can be thought of as an alpha bias, results show it could have a systematic and severe impact on beta diversity. This suggests that taphonomic biases can act differently at different scales in the spatial realm.
Acknowledgments
We would like to thank A. Quallington and Getech Group PLC for permission to use and edit the paleogeographic reconstructions in this study. We are especially thankful to S. Kidwell for providing molluscan mineralogy information, J. Alroy for providing the SQS function in R, and C. McKinney and K. Hollis for providing valuable museum data sets. We thank W. Kiessling, A. Bush, and one anonymous reviewer for comments that improved an earlier version of this work. We would also like to thank P. Wright, M. Sephton, S. Maidment, and A. Chiarenza for fruitful discussions regarding the article, as well as contributors to PhyloPic.org and the PBDB. C.D.D. was supported by a National Environmental Research Council (UK) studentship (1394514). This is Paleobiology Database official publication number 347.
