Introduction
Global environmental changes impact ecosystems at global to regional scales, but biotic responses often vary based on local conditions (e.g., Leckie et al. Reference Leckie, Yureitch, West, Finkelstein, Schmidt, Dean and Arthur1998; Seddon et al. Reference Seddon, Mackay, Baker, Birks, Breman, Buck and Ellis2014; Danise et al. Reference Danise, Clémence, Price, Murphy, Gómez and Twitchett2019; Piazza et al. Reference Piazza, Ullmann and Aberhan2020). Generalizing faunal responses to a given event among geographically widespread localities is challenging, because emergent faunal changes can be disparate, even if driven by a common factor (Seddon et al. Reference Seddon, Mackay, Baker, Birks, Breman, Buck and Ellis2014). Basin-wide changes in sea level, however, can produce predictable shifts in biofacies that are correlatable despite differences in water depth and spatial-temporal dynamics among localities (Holland et al. Reference Holland, Miller, Meyer and Dattilo2001; Miller et al. Reference Miller, Holland, Meyer and Dattilo2001; Webber Reference Webber2004). In these studies, multivariate ordination of faunal censuses across multiple stratigraphic records spanning ~65–100 km elucidated a primary faunal gradient related to water-depth preferences, thus demonstrating consistent stratigraphic patterns of faunal variation across a continental interior basin in response to large-scale environmental change. Similarly, basin-wide changes in water mass properties can drive contemporaneous shifts in faunal communities across broad spatial scales. For example, quantitative analyses combining marine sedimentary core records across ~25° latitude permitted the interpretation of spatial patterns in benthic foraminiferal faunal changes during dysoxic events in terms of underlying oceanographic drivers (Sharon and Belanger Reference Sharon and Belanger2022). Here we co-analyze benthic foraminiferal assemblages from Western Interior Seaway (WIS) localities ranging from ~1800 km to ~150 km apart (Fig. 1) to test for spatial patterns in biotic responses to biogeochemical perturbations associated with Oceanic Anoxic Event 2 (OAE2) and use those responses to infer the underlying oceanographic drivers of the environmental change.
Figure 1. Location of the 10 study localities in the late Cenomanian U.S. Western Interior Seaway (WIS); shorelines based on Blakey (2016). Northern localities denoted by triangles: Bull Creek, WY (BC), Black Gap, SD (BG), and Hot Springs, SD (HS). Central localities denoted by circles: Horsetooth Reservoir, CO (HR), Rock Canyon, CO (RC), Hartland-Bridge Creek, KS (HBC), Bunker Hill, KS (BH), and Cuba, KS (CK). Southern localities denoted by squares: Carthage, NM (CNM), and Lozier Canyon, TX (LC).
OAE2 is associated with biotic turnover in the marine realm (e.g., Elder Reference Elder1989; Leckie et al. Reference Leckie, Yureitch, West, Finkelstein, Schmidt, Dean and Arthur1998; Friedrich et al. Reference Friedrich, Erbacher and Mutterlose2006; Corbett and Watkins Reference Corbett and Watkins2013; Coccioni et al. Reference Coccioni, Sideri, Frontalini and Montanari2016; Elderbak and Leckie Reference Elderbak and Leckie2016; van Helmond et al. Reference van Helmond, Sluijs, Papadomanolaki, Plint, Gröcke, Pearce and Eldrett2016; Boudinot et al. Reference Boudinot, Dildar, Leckie, Parker, Jones, Sageman, Bralower and Sepúlveda2020) and is generally recognized by a positive δ13C excursion in the latest Cenomanian (e.g., Scholle and Arthur Reference Scholle and Arthur1980; Pratt and Threlkeld Reference Pratt, Threlkeld, Scott and Glass1984; Erbacher et al. Reference Erbacher, Oliver, Wilson, Birch and Mutterlose2005; Caron et al. Reference Caron, Dall'Agnolo, Accarie, Barrera, Kauffman, Amédro and Robaszynski2006; Sageman et al. Reference Sageman, Meyers and Arthur2006; Elrick et al. Reference Elrick, Molina-Garza, Duncan and Snow2009; Jenkyns Reference Jenkyns2010; Laurin et al. Reference Laurin, Barclay, Sageman, Dawson, Pagani, Schmitz, Eaton, McInerney and McElwain2019; Jones et al. Reference Jones, Sageman, Selby, Jicha, Singer and Titus2021). OAE2 was accompanied by dynamic environmental changes in the WIS that could drive basin-wide ecological change but could also impart spatial heterogeneity in those changes. For example, microfossil studies reveal changes in the position of thermal fronts and upwelling and downwelling zones through the OAE2 time interval (e.g., Fisher and Hay Reference Fisher, Hay, Johnson and Barrera1999; Fisher Reference Fisher2003; Elderbak and Leckie Reference Elderbak and Leckie2016; Eldrett et al. Reference Eldrett, Dodsworth, Bergman, Wright and Minisini2017; Lowery et al. Reference Lowery, Leckie, Bryant, Elderbak, Parker, Polyak, Schmidt, Snoeyenbos-West and Sterzinar2018; Bryant et al. Reference Bryant, Leckie, Bralower, Jones and Sageman2021). Further, during OAE2, oxygenation and organic matter quality were influenced by local and regional changes such as increased terrestrial runoff and iron fertilization due to increased volcanic activity (e.g., Adams et al. Reference Adams, Hurtgen and Sageman2010; Barclay et al. Reference Barclay, McElwain and Sageman2010; Monteiro et al. Reference Monteiro, Pancost, Ridgwell and Donnadieu2012; Raven et al. Reference Raven, Fike, Bradley, Gomes, Owens and Webb2019; Boudinot et al. Reference Boudinot, Dildar, Leckie, Parker, Jones, Sageman, Bralower and Sepúlveda2020). Benthic foraminiferal assemblages track changes in dominant water masses, fluctuating oxygenation, and the amount and quality of organic matter exported to the seafloor (e.g., Bernhard Reference Bernhard1986; Gooday Reference Gooday1993; Jorissen et al. Reference Jorissen, Stigter and Widmark1995). Thus, we use benthic foraminiferal assemblages to understand ecosystem responses to OAE2 and to test for a consistent basin-wide response to this global perturbation.
In the WIS, OAE2 is also associated with increased density and diversity of benthic foraminifera relative to planktonic foraminifera; this interval is termed the Benthonic Zone (BZ) (Eicher and Worstell Reference Eicher and Worstell1970; Leckie Reference Leckie, Pratt, Kauffman and Zelt1985; Leckie et al. Reference Leckie, Yureitch, West, Finkelstein, Schmidt, Dean and Arthur1998; Friedrich et al. Reference Friedrich, Erbacher and Mutterlose2006). Due to its stratigraphic association with OAE2, the BZ is often used as an informal biostratigraphic marker (e.g., Eicher and Worstell Reference Eicher and Worstell1970; Elderbak et al. Reference Elderbak, Leckie and Tibert2014; Lowery et al. Reference Lowery, Corbett, Leckie, Watkins, Miceli Romero and Pradmudito2014; Bryant et al. Reference Bryant, Leckie, Bralower, Jones and Sageman2021). However, its use is potentially problematic, because the current definition of the BZ is informed by the density of benthic foraminifera in sediments within individual sections, which could be controlled by taphonomic processes rather than basin-wide ecological responses. Further, single-section studies use section-specific definitions of the BZ instead of employing a universal definition. Site-specific definitions are based either on an increase in the abundances of benthic foraminifera compared with planktonic foraminifera or on the occurrence of species that were found in previously defined BZ intervals (e.g., Eicher and Worstell Reference Eicher and Worstell1970; Frush and Eicher Reference Frush, Eicher and Caldwell1975; Leckie Reference Leckie, Pratt, Kauffman and Zelt1985; Leckie et al. Reference Leckie, Schmidt, Finkelstein, Yuretich, Nations and Eaton1991; West et al. Reference West, Leckie, Schmidt, Dean and Arthur1998; Friedrich et al. Reference Friedrich, Erbacher and Mutterlose2006; Denne et al. Reference Denne, Hinote, Breyer, Kosanke, Lees, Engelhardt-Moore, Spaw and Tur2014; Elderbak et al. Reference Elderbak, Leckie and Tibert2014; Lowery et al. Reference Lowery, Corbett, Leckie, Watkins, Miceli Romero and Pradmudito2014; Parker Reference Parker2016; Bryant et al. Reference Bryant, Leckie, Bralower, Jones and Sageman2021). Therefore, we use our faunal analyses to test whether the BZ samples, as defined by previous authors, have a spatially consistent, recognizable assemblage composition that supports its use as a biostratigraphic marker.
Geological Setting
The WIS was a shallow (<200 m) epicontinental sea that occupied an asymmetrical foreland basin in western North America through the Late Cretaceous (Kauffman Reference Kauffman and Westermann1984; Caldwell and Kauffman Reference Caldwell and Kauffman1993). The deepest axis of the basin was flanked by a forebulge to the west and a hinge zone to the east, which created a gradient from mud-dominated facies nearshore to carbonate-dominated lithofacies in basinward localities (Elder et al. Reference Elder, Gustason and Sageman1994; Sageman and Arthur Reference Sageman, Arthur, Caputo and Peterson1994). Sedimentary archives from the northern part of the seaway are recorded in strata from U.S. states such as Wyoming, North Dakota, and South Dakota, where boreal-dominated water masses have been detected (Fisher Reference Fisher2003; Polyak Reference Polyak2003; Lockshin et al. Reference Lockshin, Yacobucci, Gorsevski and Gregory2017; Lowery et al. Reference Lowery, Leckie, Bryant, Elderbak, Parker, Polyak, Schmidt, Snoeyenbos-West and Sterzinar2018). In the southern WIS (New Mexico and Texas), the water masses are, at times, influenced by tropical waters from farther south (Eicher and Diner Reference Eicher, Diner, Pratt, Kauffman and Zelt1985; Lowery et al. Reference Lowery, Corbett, Leckie, Watkins, Miceli Romero and Pradmudito2014, Reference Lowery, Leckie, Bryant, Elderbak, Parker, Polyak, Schmidt, Snoeyenbos-West and Sterzinar2018; Fig. 1). The basin was shallower here because of the Comanche Platform (Donovan et al. Reference Donovan, Staerker, Pramudito, Li, Corbett, Lowery, Romero and Gardner2012), a carbonate platform and bathymetric high that behaved as a sill. Combined with sea-level fluctuations, this paleogeographic feature resulted in changing relative contributions of southern and northern waters into the WIS as the sill was eventually breached (Fisher Reference Fisher2003; Lowery et al. Reference Lowery, Leckie, Bryant, Elderbak, Parker, Polyak, Schmidt, Snoeyenbos-West and Sterzinar2018; Bryant et al. Reference Bryant, Leckie, Bralower, Jones and Sageman2021). In the central region, represented herein by strata from Colorado and Kansas, there is evidence for water mass mixing, stratification, and a thermal front (Leckie et al. Reference Leckie, Yureitch, West, Finkelstein, Schmidt, Dean and Arthur1998; Fisher Reference Fisher2003; Elderbak et al. Reference Elderbak, Leckie and Tibert2014; Lockshin et al. Reference Lockshin, Yacobucci, Gorsevski and Gregory2017).
Despite the regional environmental variation in the WIS, correlation among sections in the seaway is facilitated by chemostratigraphy and lithostratigraphy. The positive carbon isotope excursion (CIE) that defines OAE2 is well defined in outcrop and core at Pueblo, Colorado, the Global Boundary Stratotype Section and Point (GSSP) for the Cenomanian/Turonian boundary (Pratt and Threlkeld Reference Pratt, Threlkeld, Scott and Glass1984; Caron et al. Reference Caron, Dall'Agnolo, Accarie, Barrera, Kauffman, Amédro and Robaszynski2006; Sageman et al. Reference Sageman, Meyers and Arthur2006). Marker Bentonites A–D and limestone beds of Elder (Reference Elder, Pratt, Kauffman and Zelt1985) and Cobban and Scott (Reference Cobban and Scott1972) are observed in association with the CIE at Pueblo and other WIS sites (e.g., Bowman and Bralower Reference Bowman and Bralower2005; Caron et al. Reference Caron, Dall'Agnolo, Accarie, Barrera, Kauffman, Amédro and Robaszynski2006; Cobban et al. Reference Cobban, Walaszczyk, Obradovich and McKinney2006; Sageman et al. Reference Sageman, Meyers and Arthur2006; Jones et al. Reference Jones, Sageman, Oakes, Parker, Leckie, Bralower, Sepúlveda and Fortiz2019; Bryant et al. Reference Bryant, Leckie, Bralower, Jones and Sageman2021), permitting correlations using either the CIE or lithology across the WIS. At Pueblo, the onset of OAE2 is stratigraphically consistent with Bed 63 or LS1 (base of the Bridge Creek Limestone), whereas the initial decrease of carbon isotope values is consistent with Bed 86 (base of the Turonian), and values return to pre-OAE2 values at the stratigraphic position of Bentonite C. Previous workers have applied this marker bed schema to identify correlative beds at the local and regional scales to constrain the timing of biotic events (e.g., Elder et al. Reference Elder, Gustason and Sageman1994; Leckie et al. Reference Leckie, Yureitch, West, Finkelstein, Schmidt, Dean and Arthur1998; Bowman and Bralower Reference Bowman and Bralower2005; Lowery et al. Reference Lowery, Corbett, Leckie, Watkins, Miceli Romero and Pradmudito2014; Jones et al. Reference Jones, Sageman, Selby, Jicha, Singer and Titus2021).
Materials and Methods
Data Selection
Benthic foraminiferal assemblages from 10 stratigraphic sections across the U.S. WIS were taxonomically standardized for co-analysis. Samples from all localities were exclusively collected from shales and mudstones, because indurated limestones prevent easy liberation of foram tests. For this study, we only included stratigraphic sections that have species-level foram records, span the Cenomanian/Turonian boundary, and include samples assigned to the BZ by the original authors (Supplementary Table 1). Samples from Bull Creek, Crook County, Wyoming (BC), Black Gap, Pennington County, South Dakota (BG), and Hot Spring, Fall River County, South Dakota (HS), represent the northern region of the WIS; samples from Horsetooth Reservoir, Larimer County, Colorado (HR), Rock Canyon, Pueblo County, Colorado (RC), Hartland-Bridge Creek, Hamilton County, Kansas (HBC), Bunker Hill, Russell County, Kansas (BH), and Cuba, Republic County, Kansas (CK), represent the central region; samples from Lozier Canyon, Terrell County, Texas (LC), and Carthage, Socorro County, New Mexico (CNM) represent the southern region (Fig. 1). Each region includes samples from before, during, and after the reported BZ range, because only BG lacks samples from after the interval, and just HR and BH lack samples from before (Supplementary Table 1).
For three of the data sets (CNM, LC, CK), the proportional abundances of taxa were calculated for each sample based on total benthic foraminifera reported for the sample (Elderbak et al. Reference Elderbak, Leckie and Tibert2014; Lowery et al. Reference Lowery, Corbett, Leckie, Watkins, Miceli Romero and Pradmudito2014; Bryant et al. Reference Bryant, Leckie, Bralower, Jones and Sageman2021). Among the other seven localities (BC, BG, HS, HR, RC, HBC, BH), relative abundances for some samples were reported in three bins: abundant, common, and rare (Eicher and Worstell Reference Eicher and Worstell1970; Supplementary Table 2). In all localities, except for HR, where all samples were reported as binned abundances, the authors reported proportional abundances for all BZ samples. Given that we know the number of rare, common, and abundant taxa in each sample and that all abundances in a sample must sum to 100%, we can use the totality of the data to solve a system of equations that provide a range of values for each abundance bin that applies to all samples. Thus, to make the data interoperable, we assign estimated proportional abundances of species in each abundance bin as follows: rare ≤ 1% of the assemblage, common = 2–19%, and abundant > 19%. The maximum number of common species in a given sample was five (Eicher and Worstell Reference Eicher and Worstell1970); thus, when a single species was the only common species reported in a given sample, it was assigned 19%. When multiple species were reported, each species was assigned a percentage by dividing 19% by the number of common species, but no species deemed common was assigned less than 2%. Similarly, rare species were all accommodated within no more than 2% of the assemblage. After accounting for rare and common species, the remaining ~78% was assigned to a single abundant species or distributed among up to three abundant species reported in a sample. For one sample, five abundant species were reported; in this case, common was assigned to 2% and each abundant taxon was 19.2% after accommodating rare taxa. Using this method, we ensured that no abundant species was assigned a lower proportional abundance than any common species in the full data set. Although imputing the binned data in this way brings the data closer to the interval nature of the data inherent in the underlying abundances, it would not be appropriate to analyze these transformed data using parametric statistics. However, we are interested here in the dissimilarity in assemblage composition among samples. Given our abundance estimates are constrained by the number of taxa assigned to each category and the necessity of summing of 100%, we are confident that the dominance structure of the assemblages is sufficiently replicated for the ordination analyses described in “Multivariate Analysis.”
Morphotype Assignments
For some analyses, we grouped congenerics and assigned each genus to one of 10 morphotypes to reduce the effect of taxonomic uncertainty (Table 1). WIS benthic foraminiferal assemblages usually feature unbroken, whole tests with minimal compactional deformation. But poor preservation of some tests complicates species-level identification, because matrix often obscures features like sutures and wall texture that are important for species-level identification. Morphotype assignments allow us to analyze our data using a taxon-free metric that accommodates the integration of the data herein, which were collected over more than 50 years and, thus, vary in their application of species concepts. We defined morphotypes based on overall test shape and general chamber arrangement to focus analysis on the ecological structure of foraminiferal communities; although genera are often differentiated by characteristics including aperture types, wall structure, and ornamentation features, test shapes are closely related to the sedimentary environment the foraminifer occupies (Bernhard Reference Bernhard1986; Corliss and Chen Reference Corliss and Chen1988; Murray et al. Reference Murray, Alve and Jones2011; Table 1). We used morphotype definitions from Patarroyo and Martinez (Reference Patarroyo and Martinez2021), Koutsoukos and Hart (Reference Koutsoukos and Hart1990), and Corliss and Chen (Reference Corliss and Chen1988) to create a list of consensus morphotypes: branching, spherical, milioline, concave-convex, discoidal, plano-convex, biconvex, flat-tapered, cylindrical-tubular, and cylindrical-tapered (Table 1); and we assigned genera to morphotypes based on plates from Eicher and Worstell (Reference Eicher and Worstell1970) and Tappan and Loeblich (Reference Tappan and Loeblich1988). We chose not to differentiate between agglutinated taxa and calcareous taxa, although they are often assigned separate morphotypes (e.g., Reolid et al. Reference Reolid, Rodríguez-Tovar, Nagy and Olóriz2008), to avoid the influence of taphonomic biases (i.e., dissolution of calcareous specimens and breakage of agglutinated specimens), which can affect the ratio of agglutinated to calcareous benthic foraminifera. In the WIS, differences in the relative abundance of agglutinated and calcareous taxa are often attributed to differences in temperature, salinity, and carbonate content (Leckie et al. Reference Leckie, Yureitch, West, Finkelstein, Schmidt, Dean and Arthur1998; Lowery et al. Reference Lowery, Leckie, Bryant, Elderbak, Parker, Polyak, Schmidt, Snoeyenbos-West and Sterzinar2018), however, in other studies, increased abundances of agglutinated taxa are interpreted as a taphonomic effect from increased carbonate dissolution (Murray and Alve Reference Murray and Alve1999; Matsumoto et al. Reference Matsumoto, Kuroda, Coccioni, Frontalini, Sakai, Ogawa and Ohkouchi2020).
Table 1. Morphotype assignments. For complete morphotype assignments see Supplementary Data Appendix. Example genera are well represented in the assemblages (> 2% of a given sample).
Multivariate Analysis
Ordination methods are commonly used to reduce the dimensionality of a multivariate dataset and summarize major faunal gradients in paleoecological data within stratigraphic sections and across geographic regions (Bush and Daley Reference Bush and Daley2008; Patzkowsky and Holland Reference Patzkowsky and Holland2012). Simultaneous ordination of multiple stratigraphic records further allows recognition of basin-wide patterns and correlation among sedimentary sections where faunal changes are driven by a common factor (Holland et al. Reference Holland, Miller, Meyer and Dattilo2001; Miller et al. Reference Miller, Holland, Meyer and Dattilo2001; Webber Reference Webber2004; Hendy Reference Hendy2013; Belanger and Garcia Reference Belanger and Garcia2014; Sharon and Belanger Reference Sharon and Belanger2022). However, many ordination studies are restricted to single stratigraphic sections, because faunal differences among localities often prohibit data integration. Previous studies were able to compare among stratigraphic sections by maintaining consistent taxonomic identifications in each record, some of which were necessarily low resolution (i.e., classes rather than species or genera) or based upon morphological forms (Miller et al. Reference Miller, Holland, Meyer and Dattilo2001; Webber Reference Webber2004). Seven of the datasets used herein are derived from a single reference (Eicher and Worstell Reference Eicher and Worstell1970) and thus are internally taxonomically standardized. The remaining three datasets (CNM, LC, CK) were created within a single laboratory group (Elderbak et al. Reference Elderbak, Leckie and Tibert2014; Lowery et al. Reference Lowery, Corbett, Leckie, Watkins, Miceli Romero and Pradmudito2014; Bryant et al. Reference Bryant, Leckie, Bralower, Jones and Sageman2021) using species concepts from Eicher and Worstell (Reference Eicher and Worstell1970). Thus, we are confident in our ability to quantitatively compare these 10 datasets at the species level. Extension of these methods to future sites that are not taxonomically standardized would require a taxon-free approach; therefore, we also analyze these data at the morphotype level, based on the assignment procedure described earlier (Table 1), and compare stratigraphic patterns in morphotype composition with the species-level results.
We use principal coordinates analysis (abbreviated here as “PCO”) of the Bray-Curtis dissimilarities between samples based on the proportional abundance of species or morphotypes to summarize major gradients in faunal composition. This method is most appropriate for this dataset, because the Bray-Curtis dissimilarity is unaffected by joint absences between samples, which are common when comparing faunal assemblages from disparate localities (Field et al. Reference Field, Clarke and Warwick1982). R v. 4.1.2 was used for all analyses described below. We calculate dissimilarities and PCO scores using vegdist and cmdscale, respectively, from the vegan package in R (Oksanen et al. Reference Oksanen, Blanchet, Friendly, Kindt, Legendre, McGlinn and Minchin2017). To determine whether the BZ is compositionally distinct, we group samples by position within or outside the BZ as defined by previous authors (Eicher and Worstell Reference Eicher and Worstell1970; Elderbak et al. Reference Elderbak, Leckie and Tibert2014; Lowery et al. Reference Lowery, Corbett, Leckie, Watkins, Miceli Romero and Pradmudito2014; Bryant et al. Reference Bryant, Leckie, Bralower, Jones and Sageman2021). Finally, we group samples by region (northern, central, or southern) to test whether there is a geographic control on assemblage composition and on the trajectory of assemblage change over time. We examine the clustering of groups in the ordinations and perform Mann-Whitney tests on the PCO scores among groups using the function wilcox.test in R (R Core Team 2021). We also evaluate the similarity between BZ and non-BZ samples within localities and the similarity among BZ samples from all localities using Mann-Whitney tests. To compare species- and morphotype-level PCO results, we perform a Mantel test on the Euclidean distances among samples in PCO space using the function mantel in the vegan package of R (Oksanen et al. Reference Oksanen, Blanchet, Friendly, Kindt, Legendre, McGlinn and Minchin2017). To examine potential taphonomic influence on faunal composition, we examine the position of calcareous and agglutinated species in PCO space and use Spearman rank order correlation to test the association between the percentage of calcareous specimens and PCO axis scores using cor.test in R (R Core Team 2021).
Stratigraphic Correlations
We correlate among the localities using established chemostratigraphic and lithostratigraphic indicators (Fig. 2) to compare the stratigraphic expression of PCO scores among localities. Carbon isotope records measured on organic matter exist for the northern (Bryant Reference Bryant2021; L. J. Robinson, K. S. George, C. P. Fox, J. E. A. Marshall, I. C. Harding, P. R. Bown, J. R. Lively, S. Marroquín, R. M. Leckie, S. Dameron, D. R. Gröcke, N. M. Papadomanolaki, N. A. G. M van Helmond, and J. H. Whiteside personal communication), central (Bowman and Bralower Reference Bowman and Bralower2005; Sageman et al. Reference Sageman, Meyers and Arthur2006), and southern (Lowery et al. Reference Lowery, Corbett, Leckie, Watkins, Miceli Romero and Pradmudito2014; Bryant et al. Reference Bryant, Leckie, Bralower, Jones and Sageman2021) regions that allow us to define pre-OAE2, OAE2, and post-OAE2 stratigraphic intervals. All of these records capture the 2.0–2.5‰ positive CIE that defines OAE2. Pre-excursion values range from −27‰ to −25‰, excursion values range from −24‰ to −23‰, and post-excursion values range from −26‰ to −25‰ (Bowman and Bralower Reference Bowman and Bralower2005; Sageman et al. Reference Sageman, Meyers and Arthur2006; Lowery et al. Reference Lowery, Corbett, Leckie, Watkins, Miceli Romero and Pradmudito2014; Bryant Reference Bryant2021; Bryant et al. Reference Bryant, Leckie, Bralower, Jones and Sageman2021). At localities where there is no existing carbon isotope record (Supplementary Table 1), we rely on prominent WIS marker beds for correlation between sections (e.g., Hattin Reference Hattin1971, Reference Hattin, Pratt, Kauffman and Zelt1985; Elder et al. Reference Elder, Gustason and Sageman1994; Jones et al. Reference Jones, Sageman, Oakes, Parker, Leckie, Bralower, Sepúlveda and Fortiz2019). The pre-OAE2 interval generally corresponds to the Hartland Shale and its equivalents, whereas the OAE2 interval is generally stratigraphically consistent with the first limestone of the Bridge Creek Limestone Member and equivalent strata. The end-OAE2 interval occurs stratigraphically between Elder's (Reference Elder1989) Limestone 9 and Bentonite C (Fig. 2). Our OAE2 intervals and identification of Limestone 1 (Bed 63) and its equivalents are consistent with the stratigraphic position of the last occurrence of Rotalipora cushmani, which previous authors have used to correlate WIS sections.
Figure 2. Chemostratigraphic and lithostratigraphic correlations from Billings, MT (Bryant Reference Bryant2021; L. J. Robinson, K. S. George, C. P. Fox, J. E. A. Marshall, I. C. Harding, P. R. Bown, J. R. Lively, S. Marroquín, R. M. Leckie, S. Dameron, D. R. Gröcke, N. M. Papadomanolaki, N. A. G. M van Helmond, and J. H. Whiteside personal communication), Pueblo, CO (Sageman et al. Reference Sageman, Meyers and Arthur2006), Carthage, NM (Bryant et al. Reference Bryant, Leckie, Bralower, Jones and Sageman2021), and Lozier Canyon, TX (Lowery et al. Reference Lowery, Corbett, Leckie, Watkins, Miceli Romero and Pradmudito2014). The Oceanic Anoxic Event 2 (OAE2) interval as defined by the positive carbon isotope excursion is shaded in gray. All carbon isotope records are of organic carbon. Marker beds of Elder (Reference Elder, Gustason and Sageman1994). Capital letters represent marker Bentonites; “a”, “b”, and “c” denote the carbon isotope curve inflections from Pratt and Threkeld (Reference Pratt1984), and small letters in parentheses correspond to the scheme by Caron et al. (Reference Caron, Dall'Agnolo, Accarie, Barrera, Kauffman, Amédro and Robaszynski2006). GSSP, Global Boundary Stratotype Section and Point.
Results
Species and Morphotype Assemblages
Across the full data set, we find 87 species in 55 genera, which we placed in 10 morphotypes (Table 1). Of the 87 species, 30 species in the data set are agglutinated. Agglutinated individuals are most abundant in the northern region and absent at CK, LC, and CNM (Supplementary Fig. 1). Northern localities have more agglutinated species (>8) in their assemblages, while HS, RC, and HBC have fewer agglutinated species (<4).
We focus on four ecologically meaningful taxa, Neobulimina albertensis, Tappanina laciniosa, Buliminella fabilis, and Gavelinella dakotensis (Table 2), which are abundant or common at most sites and within many samples. Notably, N. albertensis is absent only at BC. Neobulimina albertensis is rare before OAE2 and abundant to common during OAE2, except at LC, where it is abundant before OAE2. Tappanina laciniosa increases in abundance during the OAE2 interval at CNM, RC, and CK; increases after the OAE2 interval in the northern region; and is absent at HR and LC. Buliminella fabilis largely follows a pattern similar to T. laciniosa, including being absent at HR, but it also increases in abundance at HBC and BH during the OAE2 interval (Supplementary Fig. 1). Gavelinella dakotensis increases in relative abundance at southern and central localities during the OAE2 interval, while it increases in relative abundance at the northern localities after the OAE2 interval.
Table 2. Paleoecological attributes of indicator species. PCO, principal coordinates analysis.
In terms of morphotypes, the OAE2 interval is characterized by cylindrical-tapered and concave-convex forms, except at BC in the northern region, where cylindrical-tubular forms dominate (Supplementary Fig 2). After OAE2, cylindrical-tapered and flat-tapered forms appear at the northern localities. Among the central localities, cylindrical-tapered morphotypes dominate during the OAE2 interval, and the assemblages all feature cylindrical-tubular and flat-tapered forms. After OAE2, concave-convex forms dominate at HR, and other central sites are dominated by cylindrical-tapered and concave-convex forms (Supplementary Fig 2). Flat-tapered and cylindrical-tubular morphotypes appear in post-OAE2 assemblages at HR and RC. During OAE2, assemblages are dominated by cylindrical-tapered forms in the southern regions and flat-tapered forms first occur in the record at CNM during this interval. After OAE2, assemblages are dominated by cylindrical-tapered morphotypes with 25–40% of the assemblage being biconvex forms between 45 and 50 m at LC (Supplementary Fig. 2). Cylindrical-tubular and flat-tapered forms only occur in the post-OAE2 assemblages at LC.
Gradients in Assemblage Composition
We use PCO analyses to examine gradients in assemblage composition across all sites on a common scale that incorporates the full benthic foraminiferal fauna. PCO axis 1 summarizes 31.55% of the variance in species composition, PCO axis 2 summarizes 14.89%, PCO axis 3 summarizes 10.15%; all other axes summarize less than 4.72% of the assemblage variance (Fig. 3C). Given that 46.44% of the species variance is summarized by PCO axes 1 and 2, we focus our interpretations on those axes. PCO axes 3 and 4 summarize 10.15% and <5% of the variance, respectively (Supplementary Table 2), so we do not interpret those axes. Species do not group in PCO axis 1 or 2 space according to test wall material (Fig. 3C). Neobulimina albertensis has the most positive PCO axis 1 score (Fig. 3C, Table 2). Gavelinella dakotensis and T. laciniosa have more negative PCO axis 1 scores (Fig. 3C, Table 2). Gavelinella dakotensis has the most negative PCO axis 2 scores. Neobulimina albertensis has a slightly positive PCO axis 2 score, whereas T. laciniosa and B. fabilis have high positive PCO axis 2 scores (Fig. 3C, Table 2). In the species-level ordination, agglutinated and calcareous taxa are distributed throughout PCO space (Fig 3, Supplementary Table 5). The relative abundance of calcareous species is only weakly related to PCO axis 1 and 2 scores for both species and morphotype analyses (Spearman's rho = −0.4 to 0.5; Supplementary Table 3).
Figure 3. Principal coordinates analysis (PCO) scores of samples, species, and morphotypes. A, PCO scores of samples by species analysis, where filled points represent Benthonic Zone (BZ) samples. B, PCO scores of samples by morphotype analysis, where filled points represent BZ samples. C, PCO scores of species; indicator taxa are identified with filled down-facing triangles, B.f = Buliminella fabilis; T.l = Tappanina lacinosa; N.a = Neobulimina albertensis; G.d = Gavelinella dakotensis. D, PCO scores of morphotypes.
In the morphotype ordination, PCO axis 1 summarizes 31.39% of the variance, PCO axis 2 summarizes 14.60%, and PCO axes 3 and 4 summarize less than 7% of the variance (Supplementary Table 2), similar to the species-level analysis (Fig. 3D). The positions of samples within the species-level and morphotype-based ordination spaces are positively correlated (Mantel r: 0.6175, p = 0.001). The cylindrical-tapered morphotype has the only positive PCO axis 1 score in the morphotype ordination (Fig. 3, Table 1); this morphotype typically lives infaunally and is associated with dysoxic conditions (Koutsoukos and Hart Reference Koutsoukos and Hart1990; Jorissen et al. Reference Jorissen, Stigter and Widmark1995; Table 1). In contrast, cylindrical-tubular morphotypes have the most negative PCO axis 1 scores, and spherical forms, like Lagena spp., Globulina lacrima, and Saccammina alexanderi (Bernhard Reference Bernhard1986; Koutsoukos and Hart Reference Koutsoukos and Hart1990; Jorissen et al. Reference Jorissen, Stigter and Widmark1995), also have low PCO axis 1 scores in the morphotype analysis. Branching, milioline, planoconvex, and biconvex morphotypes also have low PCO axis 1 scores (Table 1). On PCO axis 2, the concave-convex morphotype has the most negative score, and the discoidal morphotype has the most positive score (Fig. 3D, Table 1). The cylindrical-tapered morphotype has more negative scores, while cylindrical-tubular forms have more positive scores on PCO axis 2 (Fig. 3D, Table 1). Cylindrical-tubular species such as N. albertensis (Bernhard Reference Bernhard1986; Koutsoukos and Hart Reference Koutsoukos and Hart1990) have positive PCO axis 2 scores.
Regional Patterns in Sample Scores through Time
Samples from each region have both negative and positive scores on PCO axis 1 and 2, such that regions overlap in species (Fig. 3A) and morphotype ordination space (Fig. 3B); however, there are regional differences in PCO scores through time (Figs. 4, 5). Stratigraphic patterns in PCO scores for both the species and morphotype analyses are largely consistent among regions; therefore, we describe the temporal patterns for the species-level analysis and discuss the morphotype analysis only where results differ.
Figure 4. Stratigraphic plots of principal coordinates analysis (PCO) axis 1 scores by species (solid line) and morphotype (dashed) for all localities: Bull Creek (BC), Black Gap (BG), and Hot Springs (HS), Horsetooth Reservoir (HR), Rock Canyon (RC), Hartland-Bridge Creek (HBC), Bunker Hill (BH), Cuba, Kansas (CK), Carthage, NM (CNM), and Lozier Canyon, Texas (LC). Refer to Fig. 1 caption for locality abbreviations. Oceanic Anoxic Event 2 (OAE2) interval is shaded; filled points represent Benthonic Zone (BZ) samples as defined by previous authors.
Figure 5. Stratigraphic plots of principal coordinates analysis (PCO) axis 2 scores by morphotype (solid line) and species (dashed) for all localities (see Fig. 1 caption for locality abbreviations). Oceanic Anoxic Event 2 (OAE2) interval is shaded; filled points represent Benthonic Zone (BZ) samples as defined by previous authors.
PCO axis 1 scores from the species-level analyses range from −0.544 to 0.332, giving a gradient length of 0.876. PCO axis 1 scores are high and stable before and through OAE2 at the southernmost site, LC (Fig. 4). PCO axis 1 scores are similarly stable, but lower, at the northernmost locality, BC (Fig. 4). Scores at BG, HS, and HR increase by >50% of the gradient length through the OAE2 interval, whereas RC, CK, BH, HBC, and CNM values increase by between 70% and 96% of the gradient through the interval (Fig. 4). Post-OAE2 samples predominantly include samples with high PCO axis 1 scores, except at BC, HR, and RC, where some samples have low scores (Fig. 4). Temporal patterns in scores differ between morphotypes and species analyses at BC, such that morphotype-level scores increase across nearly the full gradient length of 0.984 after OAE2, but this shift is not recorded in the species-level analysis (Fig. 4).
PCO axis 2 scores from the species-level analysis range from −0.584 to 0.326 for a total gradient length of 0.910. PCO axis 2 scores are stable through the OAE2 interval at LC, BC, and BG (Fig. 5). Sample scores within the OAE2 interval at HBC decrease by ~90% of the gradient length, while all other localities decrease by ~50–70% of the gradient (Fig. 5). Post-OAE2, scores decrease across the full gradient length at BC, but sample scores are mostly stable around a score of ~0.1 at other localities such as BG, HS, and CNM (Fig. 5). HR, CK, BH, and HBC have samples with both low and high PCO axis 2 scores (Fig. 5). Temporal patterns in scores differ between morphotypes and species analyses at RC, such that morphotype-level scores increase by >50% of the gradient (gradient length for morphotype level = 1.023) after OAE2, but this shift is not recorded in the species-level analysis (Fig. 5).
Distinctiveness of BZ Assemblages
BZ and non-BZ samples overlap in ordination space (Fig. 3), but there are regional differences in assemblage compositions (Fig. 6). PCO axis 1 scores for BZ samples from the northern and central regions are indistinguishable from each other in both the species and morphotype analyses, whereas PCO axis 1 scores for BZ samples from the southern region are significantly higher than the northern and central region scores (Fig. 6). PCO axis 2 scores for BZ samples are indistinguishable among regions in the species-level analysis; however, the central region has significantly higher PCO axis 2 scores for BZ samples in the morphotype-level analysis (Fig. 6D). Non-BZ samples are also distinct at the regional level (Supplementary Fig. 3). Within localities, BZ and non-BZ samples are significantly different in both PCO axis 1 and 2 scores at only two sites; however, the distinction is supported in only the species-level analysis for each site (Supplementary Table 4). The BZ samples, as defined by previous authors, occur during the OAE2 interval at central and southern localities but occur after the OAE2 interval among northern localities (Supplementary Table 1). BZ samples from BG, HS, HR, RC, CNM, and LC include samples with higher PCO axis 1 scores (Fig. 4). BZ samples from BC, BH, CK, and HBC have low PCO axis 1 scores. In contrast, BZ samples from all localities feature samples with low PCO axis 2 scores during all or part of the interval (Fig. 5).
Figure 6. Principal coordinates analysis (PCO) scores of Benthonic Zone samples from the northern, central, and southern regions. A, PCO axis 1 scores by species level. B, PCO axis 1 scores by morphotype level. C, PCO axis 2 scores by species level. D, PCO axis 2 scores by morphotype level. For each box, the bold line is the median PCO axis score, and the box encapsulates the interquartile range (IQR). Whiskers extend 1.5 times the IQR, and open points are PCO scores outside that range. Notches in boxes extend 1.57 times the IQR divided by the square root of the number of samples (n) to approximate a 95% confidence interval; notches that do not overlap suggest the medians are significantly different. U-values for significant pairwise Mann-Whitney tests are presented; bar under U-value extends between the included regions. *p < 0.05; ***p < 0.001.
Discussion
Assemblage Compositions Indicate Variation Controlled by Oxygenation and Organic Matter Quality
The position of species and morphotypes along PCO axis 1 suggests that faunal variation among the assemblages herein is primarily driven by oxygenation. For discussion, we use a relative oxygen scale in which “dysoxic” means the lowest oxygen conditions with fauna still present; we use “suboxic” to refer to conditions between 1 and 20 μmol O2/L, and we use “oxic” to refer to conditions with >20 μmol O2/L (Helly and Levin Reference Helly and Levin2004; Middelburg and Levin Reference Middelburg and Levin2009). Neobulimina albertensis has the highest PCO axis 1 scores among species, and we interpret its high abundance to indicate the lowest oxygen conditions recorded in our samples (Fig. 3C, Table 2). This species is associated with dysoxic conditions and with water-column stratification in the WIS (Leckie et al. Reference Leckie, Yureitch, West, Finkelstein, Schmidt, Dean and Arthur1998; West et al. Reference West, Leckie, Schmidt, Dean and Arthur1998; Lowery et al. Reference Lowery, Leckie, Bryant, Elderbak, Parker, Polyak, Schmidt, Snoeyenbos-West and Sterzinar2018; Boudinot et al. Reference Boudinot, Dildar, Leckie, Parker, Jones, Sageman, Bralower and Sepúlveda2020; Bryant et al. Reference Bryant, Leckie, Bralower, Jones and Sageman2021). On the western shore of the WIS, when samples have benthic foraminiferal assemblages dominated by N. albertensis, the sample typically also has nannofossil assemblages that suggest episodic productivity and biomarkers for euxinia (Parker Reference Parker2016; Fortiz Reference Fortiz2017; Boudinot et al. Reference Boudinot, Dildar, Leckie, Parker, Jones, Sageman, Bralower and Sepúlveda2020). Similarly, when N. albertensis dominates at CNM, geochemical records reconstruct elevated C/N ratios that are likely driven by dysoxic conditions (Bryant et al. Reference Bryant, Leckie, Bralower, Jones and Sageman2021). Furthermore, modern morphological analogues that share features with N. albertensis, such as cylindrical-tapered morphology, pronounced and serially arranged chambers, and a high, hooded aperture, are able to thrive in anoxic conditions through specialized physiology and relationships with nitrate-reducing bacteria (Risgaard-Petersen et al. Reference Risgaard-Petersen, Langezaal, Ingvardsen, Schmid, Jetten, M, Camp and Derksen2006; Bernhard et al. Reference Bernhard, Casciotti, McIlvin, Beaudoin, Visscher and Edgcomb2012; Suokhrie et al. Reference Suokhrie, Saraswat and Nigam2020). Thus, we are confident that samples with high PCO axis 1 scores in both the species-level and morphotype-level analyses were deposited in dysoxic conditions.
Species such as Tappanina laciniosa and all species belonging to the branching, spherical, milioline, and plano-convex morphogroups have negative PCO axis 1 scores (Table 1 and 2). Tappaniniids, including T. laciniosa, have been previously associated with the onset of better oxygenated conditions in Late Cretaceous strata (Friedrich and Erbacher Reference Friedrich and Erbacher2006; Friedrich et al. Reference Friedrich, Erbacher and Mutterlose2006) and with increased seafloor oxygenation in the WIS (Eicher and Diner Reference Eicher, Diner, Pratt, Kauffman and Zelt1985; Leckie et al. Reference Leckie, Yureitch, West, Finkelstein, Schmidt, Dean and Arthur1998). Further, epifaunal species, which are associated with oxic to suboxic conditions, also have low PCO axis 1 scores (Table 1). Thus, we interpret samples with lower PCO axis 1 scores in both the species-level and morphotype analyses as representing suboxic conditions.
Species and morphotypes are arrayed on PCO axis 2 consistent with a gradient in organic matter quality, suggesting the changes in carbon export impart a secondary control of assemblage composition in this dataset. Gavelinella dakotensis and T. laciniosa are typically associated with enhanced productivity and organic matter export; however, these species differ in the type of organic matter they prefer (e.g., Kuhnt and Wiedmann Reference Kuhnt, Wiedmann, Huc and Schneidermann1995; Holbourn and Kuhnt Reference Holbourn and Kuhnt2002; Gustafsson et al. Reference Gustafsson, Holbourn and Kuhnt2003). Gavelinella dakotensis feeds at the sediment–water interface and has a preference for labile organic matter associated with episodic productivity (Bernhard Reference Bernhard1986; Koutsoukos and Hart Reference Koutsoukos and Hart1990; Gooday Reference Gooday1993; Thomas and Gooday Reference Thomas and Gooday1996; Leckie et al. Reference Leckie, Yureitch, West, Finkelstein, Schmidt, Dean and Arthur1998; West et al. Reference West, Leckie, Schmidt, Dean and Arthur1998), whereas T. laciniosa, which lives deeper in the sediment, has a preference for more refractory organic matter (e.g., Bernhard Reference Bernhard1986; Koutsoukos and Hart Reference Koutsoukos and Hart1990; Jorissen et al. Reference Jorissen, Stigter and Widmark1995; Gustafsson et al. Reference Gustafsson, Holbourn and Kuhnt2003). During OAE2 in the WIS, increased relative abundance of G. dakotensis is stratigraphically consistent with a relative peak of δ13C values in the latest Cenomanian in the WIS and elsewhere (Elderbak and Leckie Reference Elderbak and Leckie2016; Boudinot et al. Reference Boudinot, Dildar, Leckie, Parker, Jones, Sageman, Bralower and Sepúlveda2020; Bryant et al. Reference Bryant, Leckie, Bralower, Jones and Sageman2021), suggesting a change in productivity. Further, elevated C/N ratios at CNM concurrent with an increase in carbonate content suggest a switch to pulsed productivity when G. dakotensis is abundant (Bryant et al. Reference Bryant, Leckie, Bralower, Jones and Sageman2021), and elevated biomarkers for euxinia at Big Water, Kane County, Utah, USA, are associated with high productivity and concurrent with abundant G. dakotensis (Boudinot et al. Reference Boudinot, Dildar, Leckie, Parker, Jones, Sageman, Bralower and Sepúlveda2020). In Colorado, Utah, New Mexico, and Texas, samples with increased relative abundance of G. dakotensis are also stratigraphically associated with planktonic foraminiferal and calcareous nannofossil assemblages indicative of high-fertility surface waters (Leckie et al. Reference Leckie, Yureitch, West, Finkelstein, Schmidt, Dean and Arthur1998; West et al. Reference West, Leckie, Schmidt, Dean and Arthur1998; Lowery et al. Reference Lowery, Corbett, Leckie, Watkins, Miceli Romero and Pradmudito2014; Elderbak and Leckie Reference Elderbak and Leckie2016; Parker Reference Parker2016; Bryant et al. Reference Bryant, Leckie, Bralower, Jones and Sageman2021). Thus, both geochemical and independent microfossil evidence suggest high abundances of G. dakotensis in the records we use here reflect times of high availability of labile organic matter and support our interpretation that low PCO axis 2 scores indicate episodic productivity.
In the modern, episodic productivity at oceanic fronts is effective at exporting large quantities of labile organic matter to the seafloor, in contrast to steady productivity away from fronts (Stukel et al. Reference Stukel, Aluwihare, Barbeau, Chekalyuk, Goericke, Miller and Ohman2017). Thus, the relative abundances of these species indicate changes in the style of WIS productivity or proximity to oceanic fronts. Our analyses show that G. dakotensis has a low PCO axis 2 score, while T. laciniosa has a high score (Fig. 3C, Table 2), and changes in PCO axis 2 therefore reflect relative contributions of episodic and steady carbon export. In addition, N. albertensis, which is also thought to tolerate low levels of labile organic matter (Lipps Reference Lipps, Tevesz and McCall1983; Bernhard Reference Bernhard1986; Koutsoukos and Hart Reference Koutsoukos and Hart1990), also has high PCO axis 2 scores (Fig. 3C, Table 2), further supporting that samples with high PCO axis 2 values were deposited during times of sustained, rather than episodic, productivity.
Given that this ordination analysis considers the full benthic foraminiferal assemblage and not just known indicator species, we can also infer the ecological preferences of the species that covary with the indicator species. Buliminella fabilis covaries with T. laciniosa in ordination space (Fig. 3C, Table 2), suggesting that they have similar ecological preferences, including a tolerance for suboxic conditions and refractory organic matter. This interpretation is consistent with that for other buliminiids, which are usually cylindrical and tapered and associated with high fluxes of low-quality organic matter (Koutsoukos et al. Reference Koutsoukos, Leary and Hart1990; Coccioni and Galeotti Reference Coccioni and Galeotti1993; Gooday Reference Gooday1993; Ortiz et al. Reference Ortiz, Alegret, Payros, Orue-Etxebarria, Apellaniz and Molina2011; Orabi and Khalil Reference Orabi and Khalil2014; Deprez et al. Reference Deprez, Jehle, Bornemann and Speijer2017; Arreguín-Rodriguez et al. Reference Arreguín-Rodríguez, Barnet, Leng, Littler, Kroon, Schmidt, Thomas and Alegret2021). Thus, B. fabilis can be used as an additional indicator of environments with refractory organic matter and suboxic conditions.
The faunal gradients captured in the species-level analysis were largely replicated in the morphotype-level analysis, despite the loss of taxonomic information. This replication is evidenced by the positive multivariate correlation between the PCO scores from each analysis and the consistency in species and morphotype position in PCO space described earlier. This consistency is largely because the dominant species in the WIS assemblages we analyze have distinct morphotypes. However, there are incongruities between the positions of species in the species and morphotype ordination spaces. When species of the same morphotype are well distributed in ordination space, the dominant species drives the PCO score in the morphotype-level analysis. For example, all other convex-concave species get grouped with G. dakotensis, which has a low PCO axis 2 score (Fig. 3C) and is common at most localities (Supplementary Fig. 1). Given concave-convex species with more positive or more negative PCO scores could be combined with G. dakotensis at the morphotype level, PCO axis 2 scores could diverge between the species and morphotype analyses depending on the species composition of the sample. This may explain the disagreement in PCO axis 2 scores at HR and RC after OAE2 (Fig. 5). In another instance, B. fabilis is grouped with N. albertensis because they both have a cylindrical-tapered morphotype; however, B. fabilis has a more negative PCO axis 1 score than N. albertensis in the species-level analysis. This suggests that B. fabilis and N. albertensis are ecologically different despite having the same overall morphotype and may explain brief disagreement in PCO axis 1 scores after OAE2 at BC (Fig. 4). Buliminella fabilis is smaller, less elongate, and more ovoid than N. albertensis; more elongate foraminifera are often interpreted as preferring deeper infaunal environments (Bernhard Reference Bernhard1986; Gooday Reference Gooday1993; Jorissen et al. Reference Jorissen, Stigter and Widmark1995).
Refining the morphotype approach to include ecologically relevant morphological features, such as size or aspect ratio, could further improve the ability of morphotype analyses to replicate species-level efforts. Despite this minor dissimilarity between the analyses, PCO axes from both analyses reconstruct the same environmental gradients, and thus our interpretations based on morphotypes are largely equivalent to those based on species (Figs. 4, 5). Morphotype analyses disregard test-type information; thus, changes in the relative contribution of agglutinated and calcareous taxa do not contribute to that ordination. The consistency between species-level and morphotype analyses therefore suggests that the faunal patterns we observe are not dominated by taphonomic processes, which differentially affect calcareous and agglutinated foraminifera. PCO sample scores are also only weakly related to the percentage of calcareous specimens (Supplementary Table 3), further indicating taphonomic processes do not dominate the patterns we observe. Faunal patterns are also not strongly affected by differences in data type given samples with ordinal and continuous data are dispersed in PCO space and localities with continuous data, ordinal data, and a mix of both data types reconstruct similar temporal changes in faunal composition.
Changes in Oxygen and Productivity through OAE2 are Driven by Water Mass Interactions
By placing benthic foraminiferal assemblages from all stratigraphic sections in a common ordination, we can quantitatively compare changes in oxygenation and organic matter quality among regions and through OAE2. These analyses also separate oxygenation and organic matter quality such that we can use the PCO scores to independently examine temporal patterns of changes in oxygenation and organic matter quality in response to OAE2, unlike previous index-taxon based assessments (Friedrich et al. Reference Friedrich, Erbacher and Mutterlose2006; Elderbak et al. Reference Elderbak, Leckie and Tibert2014; Lowery et al. Reference Lowery, Corbett, Leckie, Watkins, Miceli Romero and Pradmudito2014; Elderbak and Leckie Reference Elderbak and Leckie2016).
At the localities where we have samples from the pre-OAE2 interval, we find a south-to-north gradient of increasing oxygenation evidenced by PCO axis 1 scores that are highest at LC in the south, intermediate at central localities and CNM, and lowest among northern localities (Fig. 4). This reconstruction of paleo-oxygen is consistent with paleoceanographic models that show that the southern part of the seaway was influenced before the other parts of the basin by dysoxic waters from the proto–Gulf of Mexico to the south (Arthur and Sageman Reference Arthur, Sageman and Harris2005; Denne et al. Reference Denne, Hinote, Breyer, Kosanke, Lees, Engelhardt-Moore, Spaw and Tur2014; Lowery et al. Reference Lowery, Corbett, Leckie, Watkins, Miceli Romero and Pradmudito2014; Bryant et al. Reference Bryant, Leckie, Bralower, Jones and Sageman2021). During the pre-OAE2 interval, the Greenhorn Formation and Hartland Shale in the central region of the WIS experienced prolonged dysoxia due to stratification (e.g., Hattin Reference Hattin1971; Pratt Reference Pratt1984; Arthur et al. Reference Arthur, Dean, Schlanger, Sundquist and Broecker1985; Leckie Reference Leckie, Pratt, Kauffman and Zelt1985; Denne et al. Reference Denne, Hinote, Breyer, Kosanke, Lees, Engelhardt-Moore, Spaw and Tur2014). Other studies show that some northern and west-central parts of the seaway were strongly influenced by a boreal-sourced, better-oxygenated water mass with little to no contributions from southern waters due to a strong thermal front (Eicher and Diner Reference Eicher, Diner, Pratt, Kauffman and Zelt1985; Leckie et al. Reference Leckie, Yureitch, West, Finkelstein, Schmidt, Dean and Arthur1998; Fisher and Arthur Reference Fisher and Arthur2002; Fisher Reference Fisher2003; Lowery et al. Reference Lowery, Leckie, Bryant, Elderbak, Parker, Polyak, Schmidt, Snoeyenbos-West and Sterzinar2018). Before OAE2, CNM is the only locality with positive PCO axis 2 scores that accord well with paleoceanographic models that reconstruct environments prone to steady productivity there (Bryant et al. Reference Bryant, Leckie, Bralower, Jones and Sageman2021). Before OAE2, there are small decreases in PCO axis 2 values at LC that reflect times of episodic productivity that agree with previous models of paleoproductivity conditions in the southernmost WIS (Denne et al. Reference Denne, Hinote, Breyer, Kosanke, Lees, Engelhardt-Moore, Spaw and Tur2014; Lowery et al. Reference Lowery, Corbett, Leckie, Watkins, Miceli Romero and Pradmudito2014). Thus, the benthic foraminiferal faunas examined herein reflect the expected spatial organization of WIS water masses before OAE2.
Entering OAE2, PCO axis 1 scores increase in the southern, central, and northern regions, except at BG, suggesting a decrease in oxygenation. This finding is in contrast to previous studies that suggest a basin-wide ventilation event at this interval (e.g., Eicher and Worstell Reference Eicher and Worstell1970; Leckie Reference Leckie, Pratt, Kauffman and Zelt1985; Eicher and Diner Reference Eicher, Diner, Pratt, Kauffman and Zelt1985; Leckie et al. Reference Leckie, Yureitch, West, Finkelstein, Schmidt, Dean and Arthur1998). The change in oxygenation between pre-OAE2 and OAE2 intervals is more subdued at LC than at CNM and the central sites (Fig. 4), indicating that conditions were more stable in the southern WIS. This is consistent with the presence of the southern water mass at LC before OAE2. Coordinated large decreases in PCO axis 1 scores at BG, HS, the central localities, and CNM represent the northward movement of this southern, low-oxygen, water mass into the seaway at the onset of OAE2 (Arthur and Sageman Reference Arthur, Sageman and Harris2005; Lowery et al. Reference Lowery, Leckie, Bryant, Elderbak, Parker, Polyak, Schmidt, Snoeyenbos-West and Sterzinar2018; Bryant et al. Reference Bryant, Leckie, Bralower, Jones and Sageman2021). Previous studies have found substantial evidence in chalk and limestone units for a well-oxygenated seaway during OAE2 based on increased bioturbation and less total organic carbon (e.g., Hattin Reference Hattin1971; Pratt and Threlkeld Reference Pratt, Threlkeld, Scott and Glass1984; Denne et al. Reference Denne, Hinote, Breyer, Kosanke, Lees, Engelhardt-Moore, Spaw and Tur2014) and have attributed the ventilation to the late Cenomanian transgressive pulse that disrupted pervasive stratification. Although our mudstone-dominated records may reconstruct lower-oxygen conditions than chalk/limestone records, we are confident in the spatial patterns in relative oxygenation we reconstruct because of the lithological consistency of samples among the stratigraphic sections. Further, LC has the highest contribution of chalks but the lowest oxygen faunas, counter to the expectation that environments conducive to chalk formation were generally better oxygenated.
We do not find evidence of basin-wide oxygenation associated with the late Cenomanian transgressive pulse as hypothesized by others (e.g., Eicher and Worstell Reference Eicher and Worstell1970; Leckie Reference Leckie, Pratt, Kauffman and Zelt1985; Lowery et al. Reference Lowery, Leckie, Bryant, Elderbak, Parker, Polyak, Schmidt, Snoeyenbos-West and Sterzinar2018), but instead argue that downwelling, created by impingement of the southern water mass, drove localized ventilation near the front. Because this front was created by two distinct water masses, it could have been prone to water mass mixing through caballing (e.g., Hay et al. Reference Hay, Eicher, Diner, Caldwell and Kauffman1993; Fisher Reference Fisher2003; Elderbak and Leckie Reference Elderbak and Leckie2016; Lowery et al. Reference Lowery, Leckie, Bryant, Elderbak, Parker, Polyak, Schmidt, Snoeyenbos-West and Sterzinar2018; Bryant et al. Reference Bryant, Leckie, Bralower, Jones and Sageman2021), which would have created local downwelling in the central region and led to intermittent better-oxygenated conditions despite the impingement of the low-oxygen water mass. The effect of the southern water mass is subdued at northward sites like BG, HS, and HR, where it is followed by subsequent oxygenation, and no change in oxygenation is detected at northernmost BC (Fig. 4). These patterns imply that during OAE2, a thermal front in the northern region isolated the better oxygenated, fresher, boreal water mass in the northern localities (Frush and Eicher Reference Frush, Eicher and Caldwell1975; Leckie et al. Reference Leckie, Yureitch, West, Finkelstein, Schmidt, Dean and Arthur1998; Fisher Reference Fisher2003; Lowery et al. Reference Lowery, Leckie, Bryant, Elderbak, Parker, Polyak, Schmidt, Snoeyenbos-West and Sterzinar2018). The strong thermal front prevented this pattern from dominating in the northward sites and entirely prevented a change in oxygenation at BC, demonstrating the strength of the front during OAE2 time.
The thermal front weakened after OAE2 and permitted the progression of low-oxygen conditions northward. At BC, the species-level analysis does not show a change in oxygenation, primarily due to the absence of N. albertensis, which is the primary driver of high PCO axis 1 scores. The morphotype-level analyses combined N. albertensis and B. fabilis into a single variable, thus allowing increasing abundances of B. fabilis in post-OAE2 samples from BC to drive higher PCO axis 1 scores and thus contribute to the appearance of lower-oxygen conditions. The morphotype analysis may inappropriately combine B. fabilis, which we interpret as preferring suboxic environments like T. laciniosa, with N. albertensis, a taxon associated with dysoxic conditions. In contrast, both the species and morphotype analyses suggest oxygenation declined at BG and HS after OAE2, because it was driven by increases in abundance by N. albertensis (Supplementary Fig. 1). The persistence of better-oxygenated conditions at BC but not BG and HS would suggest that water masses remained disparate between these localities in the WIS even after OAE2. Weakening of the thermal front after OAE2 has been invoked by previous studies examining single stratigraphic sections (e.g., Fisher Reference Fisher2003; Arthur and Sageman Reference Arthur, Sageman and Harris2005; Elderbak and Leckie Reference Elderbak and Leckie2016; Lowery et al. Reference Lowery, Leckie, Bryant, Elderbak, Parker, Polyak, Schmidt, Snoeyenbos-West and Sterzinar2018; Bryant et al. Reference Bryant, Leckie, Bralower, Jones and Sageman2021), and our results support this interpretation up to BG and HS but suggest it does not extend to BC.
Productivity style changed in response to the thermal front in the central part of the WIS during OAE2. PCO axis 2 scores are stable through OAE2 at the northernmost (BC) and southernmost (LC) localities (Fig. 5), indicating little change in productivity style, which can be attributed to their distance from the front. We attribute the shifts to lower PCO axis 2 scores at localities such as RC to increased contributions from episodic productivity in the central region due to its proximity to the thermal front, where oceanic mixing can increase surface water fertility (Fisher and Hay Reference Fisher, Hay, Johnson and Barrera1999; Fisher and Arthur Reference Fisher and Arthur2002) and export of organic matter (Stukel et al Reference Stukel, Aluwihare, Barbeau, Chekalyuk, Goericke, Miller and Ohman2017). At the central localities, the deoxygenation at the onset of OAE2 predates the onset of episodic productivity, which suggests that the increase in pulsed productivity was likely not the sole driver of deoxygenation. CNM receives more terrestrial refractory organic matter, because the southwestern WIS had greater influx terrestrial material and fresher waters (Leckie et al. Reference Leckie, Yureitch, West, Finkelstein, Schmidt, Dean and Arthur1998; West et al. Reference West, Leckie, Schmidt, Dean and Arthur1998; Charbonnier et al. Reference Charbonnier, Desmares and Bardin2020; Bryant et al. Reference Bryant, Leckie, Bralower, Jones and Sageman2021); this difference in organic matter quality is reflected in the higher PCO axis 2 scores at CNM than at the central localities, where marine productivity dominates organic matter production. These intersite differences suggest foraminiferal assemblages are recording local changes in organic matter instead of changes in the global carbon cycle. Local changes in oxygenation and organic matter quality are controlled by the relative contributions of boreal and southern water masses and the subsequent, localized conditions at the thermal front.
Samples assigned to the BZ interval by previous workers have different faunal compositions and occur at different times with respect to OAE2 (Figs. 4, 5). In our analysis, we focus on the primary components of faunal composition that are shared across localities and find that BZ samples do not share a specific, recognizable, assemblage of benthic foraminifera as demonstrated by comparing PCO scores (Fig. 6). Further, BZ samples are indistinguishable from non-BZ samples across localities (Supplementary Table 6), even among localities where BZ samples have ratio data and non-BZ samples have ordinal data; this data difference would only serve to make BZ and non-BZ samples appear distinct, counter to our results. We do find, however, that within stratigraphic sections, faunas from BZ assemblages have a greater abundance of faunal elements associated with labile organic matter than non-BZ samples (Fig. 5), consistent with a response to episodic productivity near the water mass front. In general, previous workers interpret the elevated density of foraminifera in BZ samples as representing better oxygenation conditions, which they attribute to a brief seafloor ventilation associated with the late Cenomanian transgressive pulse (Eicher and Worstell Reference Eicher and Worstell1970; Frush and Eicher Reference Frush, Eicher and Caldwell1975; Eicher and Diner Reference Eicher, Diner, Pratt, Kauffman and Zelt1985; Leckie Reference Leckie, Pratt, Kauffman and Zelt1985; Leckie et al. Reference Leckie, Yureitch, West, Finkelstein, Schmidt, Dean and Arthur1998). However, low-oxygen settings can result in a higher density of benthic foraminifera, because predation upon foraminifera by metazoans and competition for organic matter decreases (e.g., Levin Reference Levin, Gibson and Atkinson2003; Gooday et al. Reference Gooday, Levin, Aranda da Silva, Bett, Cowie, Dissard and Gage2009; Enge et al. Reference Enge, Wukovits, Wanek, Watzka, Witte, Hunter and Heinz2016); thus, ventilation is not needed to explain the higher density of foraminifera. Instead, our results suggest that episodic productivity and increased carbon export likely drove the higher foraminiferal densities that define most BZ samples, as previously suggested for faunas at CNM (Bryant et al. Reference Bryant, Leckie, Bralower, Jones and Sageman2021). This shift to more episodic productivity occurred at different times in the WIS (Fig. 5), as it was controlled by the south-to-north migration of the front. Thus, the BZ is not a reliable temporal marker within this basin when defined as either a consistent, recognizable, faunal composition or as increased faunal density. Instead, we argue the BZ, as presently defined, is a spatially diachronous ecological response to water mass changes in the WIS.
Conclusions
The integration of benthic foraminiferal assemblage data from multiple localities across the northern, central, and southern regions of the WIS revealed that oceanographic responses to OAE2 are spatiotemporally dynamic across the WIS due to differences in local paleoenvironments driven by water mass movement. We show that regional and local mechanisms impart a signal on the response to global events, which is important in the context of understanding the impact of biogeochemical and hydrodynamic perturbations on the marine realm. This analysis allows us to separate gradients of oxygen and organic matter quality to better constrain local paleoenvironments and to reconstruct temporal trends in oxygenation and productivity style in the context of OAE2. Temporal trends in these gradients reveal environments influenced primarily by water mass interactions and not OAE2. The BZ is not a reliable biostratigraphic marker for the onset of OAE2 in the WIS, because it does not have a consistent, recognizable faunal composition and is not a time-equivalent ecosystem response. We also find that BZ samples reflect locality-dependent changes in the style of productivity, rather than basin-wide changes in oxygenation, and thus do not yield a consistent “marker” assemblage. Using a taxon-free metric preserves these ecological signals such that morphotypes are a good substitute when species-level identifications are difficult. Co-analyzing assemblages from multiple stratigraphic records revealed previously unrecognized environmental heterogeneity across OAE2, and similar approaches in other systems can help refine paleoenvironmental histories.
Acknowledgments
We thank L. Robinson and J. Whiteside for providing carbon isotope data for chemostratigraphy. We thank K. Elderbak and C. Lowery for providing detailed foraminiferal abundance data. R.B. also thanks the Texas A&M College of Geosciences Future Faculty Postdoctoral Fellowship Program.
Declaration of Competing Interests
The authors declare no competing interests.
Data Availability Statement
All data and supplementary material is available here: https://doi.org/10.5061/dryad.4b8gthtgn.
Open access
