The impact of global climate events on local ecosystems can vary spatially. Understanding this potential heterogeneity can illuminate which environments will be most impacted and the proximal drivers of ecosystem responses. Cenomanian–Turonian marine deposits of the Western Interior Seaway (WIS) record paleoceanographic changes associated with the Greenhorn transgression and the onset of Oceanic Anoxic Event 2 (OAE2). They provide an ideal setting to study basin-wide paleoecological responses during a global perturbation. Here, we integrate benthic foraminiferal assemblages from before, during, and after OAE2 via multivariate ordination analysis to examine spatial patterns in faunal responses across the western United States on a common scale and to interrogate a previously defined faunal marker commonly used for basin-wide correlation, the Benthonic Zone (BZ). We identify oxygenation and organic matter quality as primary and secondary controls of faunal variation across the 10 stratigraphic records and use this variation to infer paleoenvironmental changes. Stratigraphic trends reveal, in contrast to previous studies, deoxygenation at the onset of OAE2. They also reveal temporal patterns in oxygenation and productivity consistent with the gradual northward migration of a southern water mass into the WIS. This spatial heterogeneity hinders the use of the BZ as a temporal marker, because assemblages change in response to diachronous environmental change, and thus timing of the BZ with respect to OAE2 is not consistent across the basin. Our study demonstrates that regional processes can overshadow ecosystem responses to global events and underscores the importance of considering how changes in the position of water masses impact the expression of global biogeochemical perturbations.

Global environmental changes impact ecosystems at global to regional scales, but biotic responses often vary based on local conditions (e.g., Leckie et al. 1998; Seddon et al. 2014; Danise et al. 2019; Piazza et al. 2020). 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. 2014). 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. 2001; Miller et al. 2001; Webber 2004). 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 2022). 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.

OAE2 is associated with biotic turnover in the marine realm (e.g., Elder 1989; Leckie et al. 1998; Friedrich et al. 2006; Corbett and Watkins 2013; Coccioni et al. 2016; Elderbak and Leckie 2016; van Helmond et al. 2016; Boudinot et al. 2020) and is generally recognized by a positive δ13C excursion in the latest Cenomanian (e.g., Scholle and Arthur 1980; Pratt and Threlkeld 1984; Erbacher et al. 2005; Caron et al. 2006; Sageman et al. 2006; Elrick et al. 2009; Jenkyns 2010; Laurin et al. 2019; Jones et al. 2021). 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 1999; Fisher 2003; Elderbak and Leckie 2016; Eldrett et al. 2017; Lowery et al. 2018; Bryant et al. 2021). 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. 2010; Barclay et al. 2010; Monteiro et al. 2012; Raven et al. 2019; Boudinot et al. 2020). 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 1986; Gooday 1993; Jorissen et al. 1995). 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 1970; Leckie 1985; Leckie et al. 1998; Friedrich et al. 2006). Due to its stratigraphic association with OAE2, the BZ is often used as an informal biostratigraphic marker (e.g., Eicher and Worstell 1970; Elderbak et al. 2014; Lowery et al. 2014; Bryant et al. 2021). 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 1970; Frush and Eicher 1975; Leckie 1985; Leckie et al. 1991; West et al. 1998; Friedrich et al. 2006; Denne et al. 2014; Elderbak et al. 2014; Lowery et al. 2014; Parker 2016; Bryant et al. 2021). 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.

The WIS was a shallow (<200 m) epicontinental sea that occupied an asymmetrical foreland basin in western North America through the Late Cretaceous (Kauffman 1984; Caldwell and Kauffman 1993). 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. 1994; Sageman and Arthur 1994). 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 2003; Polyak 2003; Lockshin et al. 2017; Lowery et al. 2018). In the southern WIS (New Mexico and Texas), the water masses are, at times, influenced by tropical waters from farther south (Eicher and Diner 1985; Lowery et al. 2014, 2018; Fig. 1). The basin was shallower here because of the Comanche Platform (Donovan et al. 2012), 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 2003; Lowery et al. 2018; Bryant et al. 2021). 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. 1998; Fisher 2003; Elderbak et al. 2014; Lockshin et al. 2017).

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 1984; Caron et al. 2006; Sageman et al. 2006). Marker Bentonites A–D and limestone beds of Elder (1985) and Cobban and Scott (1972) are observed in association with the CIE at Pueblo and other WIS sites (e.g., Bowman and Bralower 2005; Caron et al. 2006; Cobban et al. 2006; Sageman et al. 2006; Jones et al. 2019; Bryant et al. 2021), 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. 1994; Leckie et al. 1998; Bowman and Bralower 2005; Lowery et al. 2014; Jones et al. 2021).

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. 2014; Lowery et al. 2014; Bryant et al. 2021). 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 1970; 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 1970); 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 1986; Corliss and Chen 1988; Murray et al. 2011; Table 1). We used morphotype definitions from Patarroyo and Martinez (2021), Koutsoukos and Hart (1990), and Corliss and Chen (1988) 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 (1970) and Tappan and Loeblich (1988). We chose not to differentiate between agglutinated taxa and calcareous taxa, although they are often assigned separate morphotypes (e.g., Reolid et al. 2008), 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. 1998; Lowery et al. 2018), however, in other studies, increased abundances of agglutinated taxa are interpreted as a taphonomic effect from increased carbonate dissolution (Murray and Alve 1999; Matsumoto et al. 2020).

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 2008; Patzkowsky and Holland 2012). 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. 2001; Miller et al. 2001; Webber 2004; Hendy 2013; Belanger and Garcia 2014; Sharon and Belanger 2022). 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. 2001; Webber 2004). Seven of the datasets used herein are derived from a single reference (Eicher and Worstell 1970) and thus are internally taxonomically standardized. The remaining three datasets (CNM, LC, CK) were created within a single laboratory group (Elderbak et al. 2014; Lowery et al. 2014; Bryant et al. 2021) using species concepts from Eicher and Worstell (1970). 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. 1982). 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. 2017). 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 1970; Elderbak et al. 2014; Lowery et al. 2014; Bryant et al. 2021). 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. 2017). 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 2021; 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 2005; Sageman et al. 2006), and southern (Lowery et al. 2014; Bryant et al. 2021) 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 2005; Sageman et al. 2006; Lowery et al. 2014; Bryant 2021; Bryant et al. 2021). 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 1971, 1985; Elder et al. 1994; Jones et al. 2019). 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 (1989) 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.

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.

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).

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 1990; Jorissen et al. 1995; 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 1986; Koutsoukos and Hart 1990; Jorissen et al. 1995), 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 1986; Koutsoukos and Hart 1990) 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.

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).

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 2004; Middelburg and Levin 2009). 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. 1998; West et al. 1998; Lowery et al. 2018; Boudinot et al. 2020; Bryant et al. 2021). 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 2016; Fortiz 2017; Boudinot et al. 2020). Similarly, when N. albertensis dominates at CNM, geochemical records reconstruct elevated C/N ratios that are likely driven by dysoxic conditions (Bryant et al. 2021). 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. 2006; Bernhard et al. 2012; Suokhrie et al. 2020). 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 2006; Friedrich et al. 2006) and with increased seafloor oxygenation in the WIS (Eicher and Diner 1985; Leckie et al. 1998). 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 1995; Holbourn and Kuhnt 2002; Gustafsson et al. 2003). Gavelinella dakotensis feeds at the sediment–water interface and has a preference for labile organic matter associated with episodic productivity (Bernhard 1986; Koutsoukos and Hart 1990; Gooday 1993; Thomas and Gooday 1996; Leckie et al. 1998; West et al. 1998), whereas T. laciniosa, which lives deeper in the sediment, has a preference for more refractory organic matter (e.g., Bernhard 1986; Koutsoukos and Hart 1990; Jorissen et al. 1995; Gustafsson et al. 2003). 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 2016; Boudinot et al. 2020; Bryant et al. 2021), 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. 2021), 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. 2020). 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. 1998; West et al. 1998; Lowery et al. 2014; Elderbak and Leckie 2016; Parker 2016; Bryant et al. 2021). 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. 2017). 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 1983; Bernhard 1986; Koutsoukos and Hart 1990), 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. 1990; Coccioni and Galeotti 1993; Gooday 1993; Ortiz et al. 2011; Orabi and Khalil 2014; Deprez et al. 2017; Arreguín-Rodriguez et al. 2021). 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 1986; Gooday 1993; Jorissen et al. 1995).

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. 2006; Elderbak et al. 2014; Lowery et al. 2014; Elderbak and Leckie 2016).

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 2005; Denne et al. 2014; Lowery et al. 2014; Bryant et al. 2021). 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 1971; Pratt 1984; Arthur et al. 1985; Leckie 1985; Denne et al. 2014). 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 1985; Leckie et al. 1998; Fisher and Arthur 2002; Fisher 2003; Lowery et al. 2018). 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. 2021). 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. 2014; Lowery et al. 2014). 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 1970; Leckie 1985; Eicher and Diner 1985; Leckie et al. 1998). 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 2005; Lowery et al. 2018; Bryant et al. 2021). 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 1971; Pratt and Threlkeld 1984; Denne et al. 2014) 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 1970; Leckie 1985; Lowery et al. 2018), 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. 1993; Fisher 2003; Elderbak and Leckie 2016; Lowery et al. 2018; Bryant et al. 2021), 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 1975; Leckie et al. 1998; Fisher 2003; Lowery et al. 2018). 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 2003; Arthur and Sageman 2005; Elderbak and Leckie 2016; Lowery et al. 2018; Bryant et al. 2021), 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 1999; Fisher and Arthur 2002) and export of organic matter (Stukel et al 2017). 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. 1998; West et al. 1998; Charbonnier et al. 2020; Bryant et al. 2021); 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 1970; Frush and Eicher 1975; Eicher and Diner 1985; Leckie 1985; Leckie et al. 1998). 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 2003; Gooday et al. 2009; Enge et al. 2016); 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. 2021). 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.

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.

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.

The authors declare no competing interests.

All data and supplementary material is available here:

This is an Open Access article, distributed under the terms of the Creative Commons Attribution-NonCommercial-NoDerivatives licence (, which permits non-commercial re-use, distribution, and reproduction in any medium, provided that no alterations are made and the original article is properly cited. The written permission of Cambridge University Press must be obtained prior to any commercial use and/or adaptation of the article.