Plant taxonomic turnover and diversity across the Cretaceous/Paleogene boundary in northeastern Montana

Abstract

The Cretaceous/Paleogene (K/Pg) mass extinction was a pivotal event in Earth history, the latest among five mass extinctions that devastated marine and terrestrial life. Whereas much research has focused on the global demise of dominant vertebrate groups, less is known about changes among plant communities during the K/Pg mass extinction. This study investigates a suite of 11 floral assemblages leading up to and across the K/Pg boundary in northeastern Montana constrained within a well-resolved chronostratigraphic framework. We evaluate the impact of the mass extinction on local plant communities as well as the timing of post-K/Pg recovery. Our results indicate that taxonomic composition changed significantly from the Late Cretaceous to Paleocene; we estimate that 63% of latest-Cretaceous plant taxa disappeared across the K/Pg boundary, on par with other records from North America. Overall, taxonomic richness dropped by ~23–33% from the Late Cretaceous to the Paleocene, a moderate decline compared with other plant records from this time. However, richness returned to Late Cretaceous levels within 900 kyr after the K/Pg boundary, significantly faster than observed elsewhere. We find no evidence that these results are due to preservational bias (i.e., differences in depositional environment) and instead interpret a dramatic effect of the K/Pg mass extinction on plant diversity and ecology. Overall, plant communities experienced major restructuring, that is, changes in relative abundance and unseating of dominant groups during the K/Pg mass extinction, even though no major (e.g., family-level) plant groups went extinct and communities in Montana quickly recovered in terms of taxonomic diversity. These results have direct bearing on our understanding of vegetation change during diversity crises, the differing responses of plant groups (e.g., angiosperms vs. gymnosperms), and spatial variation in extinction and recovery timing.

Non-technical Summary

The Cretaceous/Paleogene (K/Pg) mass extinction is well known for the global extinction of non-avian dinosaurs and other dominant vertebrate groups. However, plants also experienced a massive turnover at this time, propelling the expansion of modern lineages into the Cenozoic. Our study of plants across the K/Pg boundary in northeastern Montana highlights the geographic heterogeneity of plant turnover during the K/Pg mass extinction and the differing effects of mass extinctions on various plant taxonomic groups. Overall, we find a high percentage of plant species disappeared across the K/Pg boundary, even though vegetation recovered relatively quickly in the earliest Paleocene in Montana. We compare this record with other studies of plant communities around the globe during the K/Pg interval.

Introduction

Mass extinctions are pivotal events in Earth history (Gould 1985) that vitally influenced the assembly of modern communities and inform our current biodiversity crisis (Barnosky et al. 2011; Ceballos et al. 2020). The most recent of the “big five” mass extinctions (Raup and Sepkoski 1982), the Cretaceous/Paleogene (K/Pg; boundary dated to 66.052 ± 0.043 Ma; Sprain et al. 2018), constituted biotic turnover on a global scale that severely affected both marine and terrestrial biotas (MacLeod et al. 1997; D'Hondt 2005). Proposed ultimate causes include the Chicxulub bolide impact (e.g., Alvarez et al. 1980; Schulte et al. 2010; Morgan et al. 2022), Deccan volcanism (e.g., Keller et al. 2009), or some combination thereof (e.g., Arens and West 2008; Renne et al. 2015). The environmental impacts of these events have been well studied, including shock waves, global wildfires, impact winter conditions lasting years, and CO₂ emissions leading to warming lasting millennia (Tobin et al. 2017; Hull et al. 2020; Morgan et al. 2022). Given its relative recency, the K/Pg mass extinction has particular relevance for modern times. It is documented by a greater quantity and quality of data than other mass extinction events, and it involved a greater proportion of taxa with living descendants whose ecologies are well understood (Jablonski et al. 2003). Detailed studies of the K/Pg event in terrestrial settings have largely derived from study areas in the western interior of North America (WINA; Nichols and Johnson 2008), with the Hell Creek study area of northeastern (NE) Montana among the most prominent (Wilson et al. 2014a). The extensive and well-constrained geologic and paleontological records of this study area document both the mass extinction and the recovery of terrestrial vertebrate faunas (e.g., Archibald and Bryant 1990; Holroyd et al. 2014; Wilson 2014); however, far less attention has been directed toward understanding the patterns among plants in NE Montana.

Data from past mass extinctions indicate that plants did not follow the same patterns as animal groups (Wing 2004; McElwain and Punyasena 2007; Green et al. 2011). During these mass extinction events, plants did not undergo family-level extinctions (McElwain and Punyasena 2007; Thompson and Ramírez-Barahona 2023; Wilf et al. 2023); instead, selective extirpation of species at the local or regional level unseated dominant plant taxa and allowed so-called disaster taxa to flourish until stable plant communities reemerged (McElwain and Punyasena 2007; Vajda and Bercovici 2014; Wilf et al. 2023). By comparison, terrestrial faunas experienced up to 62% extinction of global families during past mass extinctions (Benton 1995). Potential drivers of this difference include plants’ ability to survive short-term disruption to their reproductive cycles (i.e., through vegetative reproduction, seed dormancy), differences in resource acquisition (i.e., autotrophy in plants vs. specialized heterotrophy in animals), and ecological variation (i.e., plant families encompassing a range of ecological or life strategies) (Wing 2004; McElwain and Punyasena 2007; Green et al. 2011). Alternatively, or in parallel, apparent differences in the magnitude of extinction among plant versus animal taxa may be the result of differences in our ways of conceptualizing, measuring, and comparing biodiversity among plant versus animal communities, particularly in the fossil record (Wing 2004; Green et al. 2011; Cleal et al. 2021; Pardoe et al. 2021). Evaluating whether such differences influenced plant extinction and survival during the K/Pg event specifically requires temporally well-resolved local records of plant communities.

Studies focusing on the K/Pg event among plants have documented varying levels of extinction and diversity loss at the K/Pg boundary (KPB) and wide variability in the timing of recovery around the globe (see “Background on Plants across the KPB”; Nichols and Johnson 2008; Wilf et al. 2023). Within WINA, estimates of plant species-level extinction/extirpation are consistently ~50–75%, whereas recovery times post-KPB vary from 300 kyr (Lyson et al. 2019) to millions of years (Johnson 2002; Wilf and Johnson 2004; Peppe 2010). Previous authors have suggested environmental change as a driver (and/or potential obfuscator) of the recovery signal and heterogeneity (Johnson 2002; Lyson et al. 2019).

The well-dated plant fossil record of NE Montana presents an excellent opportunity to document local floral dynamics across the KPB and to assess potentially important regional variation in turnover patterns. Previous paleobotanical studies in NE Montana have investigated palynofloras (e.g., Hotton 2002; Arens et al. 2014b) and a few individual megafloras (e.g., Shoemaker 1966; Arens and Allen 2014; Wilson et al. 2021). Based on this previous work into the K/Pg event in WINA and elsewhere, we hypothesize that in NE Montana (H1) plant community composition was significantly altered by the K/Pg mass extinction event (and that this pattern was not driven by differences in depositional environment); and that (H2) plant species richness declined across the KPB, driven largely by elevated disappearances. We also hypothesize that (H3) Paleocene floras in NE Montana remained taxonomically depauperate compared with Cretaceous floras for more than 1 Myr after the KPB (as seen in North Dakota; Johnson 2002; Wilf and Johnson 2004). We test these hypotheses by documenting changes in taxonomic composition and richness in a composite stratigraphic sequence of 11 floras spanning ca. 2.3 Myr across the KPB (Figs. 1, 2). We then compare our results from NE Montana with previous studies of plants across the KPB in other areas of the world to understand the spatial pattern of vegetation diversity during the mass extinction.

Figure 1. Chronostratigraphic framework and sampling of floras. A, Magnetostratigraphy, lithostratigraphy, and radioisotopic dating based on Moore et al. (2014), LeCain et al. (2014), and Sprain et al. (2015, 2018). North American Land Mammal Ages (NALMAs) placed based on Sprain et al. (2018): Lancian, Puercan 1 (Pu1), inferred Puercan 2 (?Pu2), and Puercan 3 (Pu3). Floras from this study are shown in stratigraphic position relative to observed marker beds (HFZ, Hauso Flats; IrZ, Iridium Z). Floral assemblage names are abbreviated as follows: SS, Seafood Salad; F, Fisk I; SG, Smurphy's Guess II; BL, Bruce Leaf; TS, The Swamp; L, Lerbekmo N; N, New York; Y, Yabba Dabba Do; TH, Tharp's Market; BS, Biscuit Springs; and J, Jane's. On the symbol of each floral assemblage, the vertical bar indicates stratigraphic position error (see text for method of estimation), the symbol shade indicates the temporal interval (Hell Creek [HC], Tullock Member [TM1, TM2]), and the symbol shape indicates the sedimentary facies. Alternating white and gray shaded areas delineate 15 m stratigraphic intervals. B, The sample size of each assemblage is shown in the bar graph.

Figure 2. Map of Hell Creek study area. A, Map of the United States showing inset of Hell Creek area in Montana. B, Inset of Hell Creek area in NE Montana showing the 11 sites in this study. Flora symbols indicate the temporal interval (shade) and sedimentary facies (shape). Floral assemblage abbreviations as in Fig. 1. Geologic formations mapped across the study area indicating approximate location of outcrops of the Hell Creek Formation (HCF) and Tullock Member (TM) of the Fort Union Formation (FUF) as well as other geologic units (based on Raines and Johnson 1995).

Background on Plants across the KPB

Well-studied floras from the KPB have been documented in WINA (Montana, North Dakota, Colorado, Wyoming, and New Mexico), Colombia, and Argentina (Fig. 3). The pattern of plant turnover at the KPB varies among these regions and with the type of plant fossil record analyzed (i.e., palynoflora vs. megaflora) (Nichols and Johnson 2008; Wilf et al. 2023). Palynofloral studies of the K/Pg mass extinction are more numerous and sample a larger geographic extent (both in terms of global coverage and study area size) than megafloral studies (Spicer and Collinson 2014; Vajda and Bercovici 2014). Palynofloral studies generally document extirpation of 10–30% of regional palynotaxa across the KPB (range from studies around the globe; Spicer 1989; Vajda and Bercovici 2014). More locally, 30% of palynotaxa in North Dakota (Nichols 2002; Nichols and Johnson 2008) and 15–30% of palynotaxa in NE Montana (Hotton 2002) disappear across the KPB. Comparatively few studies of megaflora capture the lead-up, mass extinction, and recovery intervals. Recent work on Colombian pollen and leaf assemblages indicates that taxonomic composition changed across the KPB, but extinction magnitude among megafloras was not calculated (Carvalho et al. 2021). Megafloras from Argentinian Patagonia were greatly affected by the mass extinction, with a maximum estimated loss of 90% of plant taxa (Stiles et al. 2020). Analyses of megaflora from WINA show moderate KPB extinction or extirpation, on the order of 50–75% of species-level taxa (Wolfe and Upchurch 1986; Johnson 2002; Barclay et al. 2003; Wilf and Johnson 2004; Lyson et al. 2019).

Figure 3. Summary of global studies of megaflora from the Cretaceous/Paleogene boundary (KPB) interval. A, Temporal range and density of megafloral sampling from western interior of North America (WINA), including NE Montana (MT; this study and PDM flora described by Arens and Allen [2014]), North Dakota (ND; Wilf and Johnson 2004; Peppe 2010), Denver Basin (CO; Lyson et al. 2019), San Juan Basin (NM; data from Flynn and Peppe 2019), and Hanna Basin (WY; Dunn 2003). B, Temporal range and density of megafloral sampling from outside North America including Colombia (Carvalho et al. 2021) and Patagonian Argentina (Stiles et al. 2020) with inset box showing the extent of WINA sampling in A. Vertical bars show the estimated age range of the floral assemblages with the number of floras from that time range shown at right. The age of the K/Pg boundary and polarity chrons measured in NE Montana are taken from Swisher et al. (1993), LeCain et al. (2014), and Sprain et al. (2015, 2018).

Work on the recovery after the K/Pg mass extinction shows similar heterogeneity across the globe. Megafloras from North Dakota remained depauperate throughout the Paleocene (ca. 10 Myr) (Johnson 2002; Wilf and Johnson 2004; Peppe 2010); records from the Bighorn and Hanna Basins similarly point to species-poor Paleocene floras (Hickey 1980; Wing et al. 1995; Dunn 2003). By contrast, plant communities in the Denver Basin apparently recovered much faster (ca. 300 kyr post-KPB; Lyson et al. 2019); it has also been postulated that recovery in the Denver Basin may have varied with topographic relief (Johnson and Ellis 2002; Johnson et al. 2003). Extinction selectivity may have differed in lowland versus upland environments, or environmental heterogeneity in the uplands may have led to higher origination (Badgley et al. 2017; Antonelli et al. 2018). Farther south, Paleocene floras in the San Juan Basin exhibited higher species richness compared with the northern Great Plains, possibly the result of a latitudinal diversity gradient in plant communities linked to a warmer and wetter climate in southern WINA (Flynn and Peppe 2019). Additionally, proximity to the bolide impact site or variation in ecological stability associated with such latitudinal diversity gradients may have contributed to latitudinal selectivity during the mass extinction (Vilhena et al. 2013; Wilf et al. 2023). In comparison to WINA floras, Patagonian and Colombian Paleocene floras were relatively species rich, but Patagonian floras remained homogeneous through space and time during the Paleocene (Iglesias et al. 2007, 2021; Stiles et al. 2020) and Colombian floras did not regain pre-KPB species richness until the Eocene (Carvalho et al. 2021).

The highly variable patterns of K/Pg extinction and recovery among global plant communities may correlate with proximity to the bolide impact site, underlying ecosystem vulnerability, paleoclimate and topography, or biases in the fossil record. Given this spatial heterogeneity, high-resolution, local fossil records present the best opportunity to understand the effects of the K/Pg mass extinction and recovery around the globe. Furthermore, the nuanced understanding of community ecology garnered by local studies of mass extinctions is vital for better predicting patterns of ongoing and future biodiversity loss.

Geologic Setting of the Hell Creek Study Area

The Late Cretaceous Hell Creek Formation (HCF) and Paleogene Tullock Member (TM) of the Fort Union Formation (FUF) are exposed in the northern Great Plains region and preserve rich fossil assemblages dated to the latest Cretaceous and earliest Paleogene. Researchers have developed a high-resolution chronostratigraphic framework for the HCF and TM in Montana based on biostratigraphy, lithostratigraphy, geochemistry, magnetostratigraphy, and geochronology. Sanidine-bearing tephras within discrete lignite beds (henceforth referred to as “coals”) have been dated using ⁴⁰Ar/³⁹Ar age analysis (Renne et al. 2013; Sprain et al. 2015, 2018) and used to interpolate ages of the paleomagnetic reversals identified in local sections (Swisher et al. 1993; LeCain et al. 2014; Sprain et al. 2015, 2018). We use linear extrapolation based on stratigraphic distance to these dated marker beds to integrate our floras into this age model (Archibald et al. 1982; Swisher et al. 1993; Renne et al. 2013; LeCain et al. 2014; Moore et al. 2014; Sprain et al. 2015, 2018; Fig. 1).

The HCF is typically ~85–90 m thick (Hartman et al. 2014; Moore et al. 2014) and spans the magnetochron 30n/29r boundary (66.311 ± 0.067 Ma; Sprain et al. 2018). The upper contact of the HCF is usually placed at the base of the lowest, laterally continuous coal of the overlying FUF, referred to as the lower Z coal. In exposures where a thin (~2 cm) claystone immediately below the lower Z coal preserves an iridium anomaly and other signatures of the Chicxulub impact event, the lower Z coal is referred to as the Iridium Z or IrZ coal (Archibald et al. 1982; Smit and Van Der Kaars 1984; Rigby and Rigby 1990; Murphy et al. 2002; Moore et al. 2014). Sanidine-bearing tephras recovered from the IrZ coal have yielded a pooled age of 66.052 ± 0.043 Ma (Sprain et al. 2018). In central Garfield County, where most of our floras were collected, the HCF–FUF formational contact is often coincident with the KPB (identified by the presence of an impact claystone bed; Archibald et al. 1982; Moore et al. 2014). However, this formational contact is diachronous across the Williston Basin (Clemens 2002; Nichols and Johnson 2002). We note this diachroneity in our stratigraphic framework (Fig. 1) and account for it in our age determination of floras. In eastern Garfield and western McCone Counties, the coal present at the base of the FUF is referred to as the McGuire Creek Z (MCZ). Prior studies indicate that the KPB (identified by chemostratigraphy and palynostratigraphy; Lofgren 1995; Arens et al. 2014a) is some few meters below the FUF–HCF contact (and thus, the MCZ) in this area (Sprain et al. 2015; Smith et al. 2018), although this placement is uncertain (Tobin et al. 2021). In this study, we separate the TM into the older TM1 interval (bounded by the IrZ and Hauso Flats [HFZ] coals; 66.052 ± 0.043 to 65.973 ± 0.047 Ma; Sprain et al. 2015, 2018) and the younger TM2 interval (bounded by the HFZ and W coals; 65.973 ± 0.047 to 65.118 ± 0.048 Ma; Sprain et al. 2015, 2018). TM1 floras temporally correlate with Puercan 1 (Pu1), whereas TM2 floras likely correlate with Puercan 2 and 3 North American Land Mammal Ages (NALMAs; only Pu1 and Pu3 local faunas are observed in Montana, placement of Pu2 is inferred based on our age model; Fig. 1). Therefore, grouping floras into these TM intervals also allows us to compare floral and faunal trends. We examine dynamics across floras ca. 1.3 Myr to 110 kyr pre-KPB (HC interval), floras from within the first 80 kyr after the KPB (TM1 interval), and floras ca. 80 to 900 kyr after the KPB (TM2 interval).

The HCF and FUF record terrestrial deposition in floodplain environments to the west of the receding Western Interior Seaway (Archibald 1982; Archibald et al. 1982; Fastovsky and McSweeney 1987; Swisher et al. 1993; Johnson et al. 2002; Murphy et al. 2002; Flight 2004; Hartman et al. 2014; Moore et al. 2014; Fastovsky and Bercovici 2016; Fowler 2020; Weaver et al. 2022). Depositional environments in the HCF and FUF can largely be divided into channel and overbank (i.e., flood basin) settings. Channel deposits are typified by coarser grain size sediments, inclined strata, heterolithic composition, and a higher degree of oxidation; overbank deposits are characterized by finer grain size sediments, horizontal massive or tabular strata, and lower oxidation (Fastovsky 1987; Johnson 2002). In this study, these characterizations are the basis for assignment of our sites to either channel or overbank facies and depositional environments (following the descriptions of Johnson [2002]: table 3; see “Methods”).

Methods

Fossil Collection

We describe and analyze plant megafossils collected from 11 sites (3810 identifiable specimens from 37 quarries) between 2015 and 2019 by a team from the University of Washington. All quarries were sampled with unbiased, complete census collections of all identifiable material found, taken to the University of Washington Burke Museum of Natural History and Culture, Seattle, WA (UWBM), and counted. We received permissions from the Bureau of Land Management, the Charles M. Russell Wildlife Refuge administered by U.S. Fish and Wildlife and the Army Corps of Engineers, and private landowners for all collecting (data on our localities and quarries are provided in Supplementary Table S1). We excavated fossils using hand tools (i.e., pickaxe, chisel, and rock hammer) and collected primarily compression and impression fossils of vegetative and reproductive plant structures.

We consider a plant fossil assemblage (or flora) as all fossils collected from a single site; individual quarries at a site span no more than 40 m along strike and 3 m of stratigraphic height, and present similar lithologies. Each quarry was given a unique UWBM locality number to distinguish them for future studies and collecting (Supplementary Table S1).

Sedimentological Description and Stratigraphy

At each site, we used a hand level and a Jacob's staff as well as a high-precision (<1 m vertical and horizontal) Trimble R2 GNSS receiver to measure stratigraphic thicknesses and distance to identified marker beds (i.e., dated coal bed, formational contact). We also recorded detailed lithological descriptions of the fossiliferous horizon(s) and surrounding 2–7 m of stratigraphy (Supplementary Fig. S1). For each stratigraphic distance, we estimated uncertainty as some combination of instrument error (GPS accuracy, height of the fossiliferous horizon) and geologic context (height of channel body where appropriate). Stratigraphic error bars may be vertically asymmetrical because channel bodies incise older strata, thereby increasing the uncertainty that a flora preserved within a channel body is younger than its stratigraphic height indicates (Weaver et al. 2022; Fig. 1).

Of the 11 floras in this study, four are Late Cretaceous (HC floras) and seven are early Paleocene (TM floras). The youngest Late Cretaceous flora (Bruce Leaf) is constrained to ca. 110 kyr before the KPB based on linear extrapolation of the stratigraphic positions and ages of the C30n/29r reversal and the IrZ coal nearby (data from Thomas Ranch compiled by Sprain et al. [2018]). The oldest Paleocene flora (The Swamp) is 1.3 m below the MCZ coal (66.024 ± 0.044 Ma; Smith et al. 2018); this stratigraphic position makes The Swamp likely Paleocene in age given the diachroneity between the HCF–FUF contact and the KPB in this region (Lofgren 1995; Arens et al. 2014a; see “Geologic Setting of the Hell Creek Study Area”). In addition, characteristic taxa (e.g., MT037 Paranymphaea crassifolia) indicate a Paleocene age for The Swamp flora (Johnson 1989, 2002), and previous workers have noted palynological indicators of Paleocene age at this site (N. C. Arens personal notes 1997). Three TM floras are from the TM1 temporal interval (within 80 kyr after the KPB), and four TM floras are from the TM2 temporal interval (80 to 900 kyr after the KPB; Fig. 1). Not all dated marker beds are present at all our sites, but their stratigraphic positions are inferred through interpolation; therefore, these age constraints should be considered approximate.

Depositional setting may influence the spatial resolution and composition of preserved floras, potentially hindering direct comparisons between fossil plant assemblages. Channel deposits often reflect both plants growing locally on riverbanks and plants transported from farther away, whereas overbank deposits (e.g., crevasse splays, ponds, floodplain soils) preferentially preserve plants growing locally on the floodplain (e.g., Behrensmeyer and Hook 1992; Demko et al. 1998). Plant fossils derived from channel deposits also tend to be more abraded compared with those found in overbank sediments, which commonly preserve more anatomical detail (Ferguson 2005). Using lithology to interpret depositional environment (channel vs. overbank lithofacies, following Johnson [2002]: table 3; see discussion in “Geologic Setting of the Hell Creek Study Area”), we classified eight of the sites as channel deposits (three Cretaceous and five Paleocene) and three as overbank deposits (one Cretaceous and two Paleocene) (Fig. 1, Supplementary Figure S1, Supplementary Table S5).

Morphotype Identification and Taxon Assignment

We identified 3810 specimens to morphotype based on organ type, gross morphology, and leaf architecture. The number of identified specimens in each of our 11 floras ranges from 173 to 618 (Fig. 1). Angiosperm leaves were grouped and described based on leaf architecture following Ellis et al. (2009). Each morphotype was assigned a unique alphanumeric code (e.g., MT001; Table 1; see also Wilson et al. [2021] for a detailed description of a subset of these taxa).

Table 1. Summary of morphotype samples by major taxonomic group, calculated from abundance data (Supplementary Table S3). K/Pg, Cretaceous/Paleogene.

We recognize 122 plant fossil morphotypes in this study: 95 are non-monocotyledonous angiosperms (“dicot”), 9 are conifers, 7 are pteridophytes, 5 are plants of indeterminate affinity, 4 are monocotyledonous angiosperms, 1 is a cycadophyte, and 1 is a ginkgophyte. Of these, most are fossils of leaves or shoots with leaves (96), some are reproductive structures (e.g., fruits or cones; 20), and some are other vegetative structures (e.g., roots; 6) (Table 1; see Supplementary Tables S2 and S4 for character-state descriptions of each morphotype). Each individual specimen identified in our census count constitutes a single plant organ (e.g., leaf, root, fruit, or cone) for most morphotypes or a shoot with multiple leaves in the case of conifer vegetative morphotypes.

Data Analyses

Fossil Abundance Data

We tabulated the abundance of each morphotype in our fossil census data (Supplementary Tables S3, S6) and transformed our data in several ways to account for potential biases. Where appropriate, we excluded singletons to mitigate the impact of preservational or taphonomic biases (Foote 2000). Additionally, reproductive morphotypes may come from the same species as vegetative morphotypes and thereby inflate species richness estimates. To account for this, for some analyses, we included only leaf morphotypes (or, even more restrictively, only dicot leaf morphotypes) as a minimum estimate of taxonomic diversity (following Wilf and Johnson 2004). We compare taxonomic abundance or presence/absence in an individual flora as well as pooled abundance from the HC, TM1, and TM2 temporal intervals. See below and Supplementary Table S7 for which dataset was used in each analysis.

Taxonomic Composition

We first analyzed changes in community composition using an ordination and analysis of similarity (ANOSIM) of relative abundance data (considering all morphotypes) to test our first hypothesis (H1) that plant community composition was significantly altered by the K/Pg mass extinction event (and that this pattern was not driven by differences in depositional environments). Ordination analyses are commonly used in community ecology to plot communities in low-dimensional space based on their taxonomic similarity (Reyment 1963; Kovach 1989). Here we used an unconstrained nonmetric multidimensional scaling (NMDS) ordination, because this method is appropriate for species abundance data, avoids assumption of linearity among variables, and is relatively robust to the arch effect (Minchin 1987; Clarke 1993). Although there are drawbacks to NMDS (e.g., see discussion in Donohue et al. 2013), we contend that careful selection of ordination parameters (number of dimensions and distance metric) and vetting of stress for significance can limit these concerns.

We conducted our NMDS ordination using the Bray-Curtis distance metric (Bray and Curtis 1957), reducing to k = 2 dimensions, and running the analysis 100 times to find the minimum stress solution. To evaluate the fit of the NMDS, we ran a Monte Carlo simulation (see Supplementary Material for code; n = 100). To better understand the factors driving the ordination patterns, we also fit our site variables (facies and age) to the NMDS axes and calculated the weight of these variables to each axis (envfit function in the vegan package using 1000 permutations; Okansen et al. 2019). Finally, we quantified the differences in taxonomic composition between HC and TM floras through an ANOSIM, which utilizes rank dissimilarity based on Bray-Curtis distance to pair with our NMDS ordination (Clarke and Warwick 1994).

We investigated which taxa might be driving shifts in taxonomic composition from HC to TM floras through several species-centered analyses. First, we examined the distribution of taxa in the NMDS ordination. Second, we compared how conifers versus angiosperms changed from HC to TM floras through a Mantel test (Mantel 1967). This test is based on a Pearson's correlation of the Bray-Curtis distances (calculated from conifer or angiosperm vegetative morphotype relative abundance) between our floras.

Taxonomic Turnover and Diversity

We further examined our hypotheses that plant communities changed across the KPB in both taxonomic composition (H1) and diversity (H2) by examining patterns of turnover (first appearances and last appearances) as well as diversity (taxonomic richness and evenness) across this interval. We first calculated the percent disappearance from HC to TM floras as a measure of extinction magnitude, using range-through occurrences. We excluded singleton taxa (following Wilf and Johnson 2004) but chose to include all taxonomic groups to examine interclade extinction dynamics. We grouped our floras into 15 m stratigraphic bins to minimize sampling bias in our composite section. We then calculated 50% stratigraphic confidence intervals for each of the biostratigraphic ranges, following Wilf and Johnson's (2004) expansion of the Strauss and Sadler (1989) method as outlined by Wilson (2005, 2014). This method uses the abundance, range, and sample size of a given taxon to estimate the uncertainty in its record and to calculate upper and lower confidence intervals.

Previous researchers have calculated per capita origination and extinction rates to test hypotheses about patterns of turnover (e.g., Wilf and Johnson 2004). However, given that our floras are unevenly spaced through this time interval, per capita rates would likely have low accuracy (i.e., large time gaps over which we cannot infer rates). Instead, to test our second hypothesis (H2), we calculated proportional first and last appearances (e.g., Wilson 2014; Wilson et al. 2014b) to evaluate turnover through time: calculated as the percentage of taxa with their first or last occurrence in that flora, out of the total range through richness of that flora. We omitted the oldest and youngest floras in our composite sequence where appropriate to avoid edge effects (Foote 2000) and used only dicot leaf morphotype abundance, following Wilf and Johnson (2004) (see Supplementary Figs. S2 and S3 for analyses with all morphotypes and all dicot morphotypes). We also subsampled our dicot leaf abundance data to n = 55 specimens (our smallest sample size) for n = 50 replicates and calculated 95% confidence intervals. Note that this method potentially overemphasizes first and last appearances in depauperate floras; alternative methods for accounting for sample size (e.g., capture–mark–recapture) could be applied to address this bias in future work. Here, we have accounted for sample size by bootstrapping to estimate confidence intervals on our proportional appearances and richness estimates.

Finally, to test our third hypothesis (H3) that Paleocene floras remained depauperate through our section, we calculated taxonomic diversity through time using (1) taxonomic richness and evenness and (2) richness curves that employ both extrapolation and interpolation (rarefaction). We calculated Pielou's evenness (Pielou 1966) based on abundance of dicot leaf morphotypes, all leaf morphotypes, and all morphotypes. Taxonomic richness is commonly used as a metric for diversity in paleoecology, but it does not consider the contribution of relative abundance and is heavily impacted by the often-low sample size in fossil assemblages (Gotelli and Colwell 2001). To account for this bias, we used two measures of richness: (1) range-through (standing) richness and (2) rarefied richness. Range-through richness assumes that if a taxon is present in the preceding and succeeding time intervals (floras), then its absence from the intervening flora represents a sampling error and not a true absence. For fossil data, this assumption is considered reasonable as a means of estimating maximum diversity and to accommodate uneven temporal sampling effort (Foote 2000). Note that sampling biases result in low range-through richness at the beginning and end of our record (edge effect; Foote 2000). We used leaf morphotypes only and calculated 95% confidence intervals on range-through richness in the same manner as for proportional appearances (described earlier). Rarefied richness, which has become standard in paleoecology (Siegel and German 1982; Foote 2000), uses subsampling of abundance data to estimate taxonomic richness, thereby allowing for the direct comparison of fossil assemblages with different sample sizes. We rarefied our abundance data to the smallest sample size (Table 2) and calculated standard error. Rarefied richness and evenness (Pielou's) were calculated for all morphotypes, only leaf morphotypes, and only dicot leaf morphotypes to allow comparison with previous work. To quantify differences in diversity through time, we used a Mann-Whitney U-test to compare taxonomic richness and evenness (comparing floras grouped by age and by depositional environment). This test is nonparametric and therefore does not assume that samples are normally distributed (McKnight and Najab 2010).

Table 2. Comparison of average rarefied richness and evenness of floras by temporal interval and sedimentary facies. TM, Tullock Member; HC, Hell Creek.

In our second approach to measuring diversity changes, we calculated richness curves using an interpolative–extrapolative richness metric (R package iNEXT package; Hsieh et al. 2016). Richness curves allow for comparison of taxonomic richness across assemblages of variable sample size as well as observation of sample evenness (Gotelli and Colwell 2001). Although extrapolated richness has been criticized for its low accuracy when sample sizes are low (i.e., Melo et al. 2003), recent work has pointed to its utility, particularly because interpolative richness may compress richness ranges (Close et al. 2018). Using leaf morphotype abundances, we calculated curves for each floral assemblage as well as for floral assemblages grouped by temporal interval (HC, TM1, or TM2) with 95% confidence intervals.

All analyses were conducted in R v. 4.0.2 (R Core Team 2021) using the community ecology package vegan v. 2.5–7 (Okansen et al. 2019), the fossil diversity dynamics package divDyn (Kocsis et al. 2019), and the species diversity package iNEXT (Hsieh et al. 2016). See Supplementary Material for code.

Results

Taxonomic Composition

Taxonomic composition is significantly different between HC and TM floras (see Fig. 4 for images of select HC and TM taxa). Our NMDS ordination is a relatively good fit to the underlying data; stress is low (0.0815) and significant based on Monte Carlo simulation (p = 0.010), indicating that the ordination is a good representation of the variance in our data (Clarke and Warwick 1994). In the NMDS biplot (Fig. 5), HC floras cluster low on the first and second NMDS axes, whereas TM1 and TM2 floras cluster higher on both axes, with the TM2 floras intermediate between HC and TM1 floras. Using age as a predictor of distribution in the ordination, we found that age significantly correlates with the ordination results (envfit R ² = 0.4599, p = 0.0070), loading mostly on axis 1 (Fig. 5). Furthermore, the differences in taxonomic composition between HC and TM floras and among HC, TM1, and TM2 floras are significant (ANOSIM p = 0.003 and p = 0.001, respectively). In contrast, sedimentary facies does not significantly correlate with floral composition (envfit R ² = 0.1230, p = 0.2877), and there are no significant differences in taxonomic composition between channel and overbank floras (ANOSIM p = 0.362).

Figure 4. A–E, Examples of taxa sampled in the Hell Creek (HC) and Tullock Member (TM) floras of this study: A, MT103 Nilssonia comtula, UWBM PB 105816.1; B, MT007 Ginkgo adiantoides, UWBM PB 96540; C, MT034 likely Azolla sp., UWBM PB 97843; D, MT106 Limnobiophyllum scutatum, UWBM PB 99364; E, MT023 “Dryophyllum” subfalcatum, UWBM PB 103630; F, MT002 Metasequoia occidentalis, UWBM PB 103535.2. G–K, Examples of the five most abundant taxa that appear in the TM: G, MT053 unnamed dicot, UWBM PB 116050; H, MT037 Paranymphaea crassifolia, UWBM PB 103656; I, MT039 unnamed dicot, UWBM PB 104937; J, MT054 unnamed pteridophyte, UWBM PB 99797.1; K, MT090 unnamed dicot, UWBM PB 104855.1. Scale bar, 1 cm.

Figure 5. Biplot of nonmetric multidimensional scaling (NMDS) ordination of morphotype abundance among the 11 floras (all morphotypes included). Distance was measured using Bray-Curtis dissimilarity. Stress (0.0815) is relatively low and significant based on Monte Carlo simulation (p = 0.010). Symbols represent taxa (colors correspond to major plant group affinity) and floral assemblages (colors and shapes represent age and sedimentary facies). Age and sedimentary facies were fit to NMDS axes, and the centroids of these groupings are plotted as vectors. Age explains significant variation in the analysis (envfit p = 0.0070), whereas facies does not (envfit p = 0.2877). Abbreviations: TM, Tullock Member; HC, Hell Creek; CON, conifer; CYC, cycadophyte; DIC, dicot (non-monocotyledonous angiosperm); GIN, ginkgophyte; MON, monocotyledonous angiosperm; PTE, non-seed vascular plant (pteridophyte); IND, indeterminate; K, Cretaceous; Pg, Paleogene; ob, overbank; ch, channel. Floral assemblage abbreviations as in Fig. 1.

The NMDS analysis shows that distinct suites of mostly dicot angiosperm taxa plot with HC or TM floras. Conifer taxa are more likely to be shared across HC and TM floras, plotting between them in ordination space (Fig. 5). Furthermore, conifers and angiosperms are not significantly correlated in their distribution across the floras (Mantel statistic r = 0.1528, p = 0.17). However, whereas conifers are more likely to be shared between HC and TM floras, the relative abundance of conifer specimens dramatically declined from HC (35%) to TM1 (6%) floras and then increased again in TM2 floras (25%). Other non-angiosperms (e.g., the cycad Nilssonia comtula MT103, and Ginkgo adiantoides MT007) are rare in the sampled floras, and their distributions are therefore harder to characterize. Several monocotyledonous angiosperm and pteridophyte taxa plot close to the TM floras in ordination space; these represent common wetland or aquatic taxa such as Equisetum sp. (MT127), likely Azolla sp. (MT034; Fig. 4C) and Limnobiophyllum scutatum (MT106; Fig. 4D).

Turnover and Diversity

Analysis of taxonomic turnover shows that a large proportion of plant morphotypes disappeared across the KPB. Of the 63 HC morphotypes, 48 (76%) disappeared at the KPB, whereas 41 of 52 (79%) vegetative HC morphotypes disappeared. Most of the taxa that disappeared at the KPB are singletons; only 8 of 23 (35%) non-singleton morphotypes disappeared at the KPB (Fig. 6). When we restrict the analysis to the uppermost 15 m of the HCF (i.e., Bruce Leaf flora), dated to the last ca. 110 kyr of the Late Cretaceous, 17 of 27 (63%) morphotypes were lost at the KPB. By restricting our analysis to the uppermost HCF, we exclude any earlier disappearances that were not associated with the mass extinction event and reduce the chance of artificially inflating our estimate of extirpation at the KPB; however, we acknowledge that the resultant sample size is very small (a single HC flora), and therefore we may not be accurately reflecting the diversity of vegetation just before the KPB.

Figure 6. Stratigraphic ranges for non-singleton taxa in this study. The ranges (thick, vertical black bars) are based on occurrences within the 11 floras, which are grouped into 15 m bins. The stratigraphic section at the left shows the stratigraphic height relative to the Cretaceous/Paleogene (K/Pg) boundary, geologic formations, paleomagnetic polarity chrons, temporal intervals (Hell Creek [HC], Tullock Member [TM1, TM2]), and sites included in each 15 m bin (marked by alternating white and gray horizontal bands). Floral assemblage abbreviations as in Fig. 1. Light gray vertical lines indicate 50% confidence intervals. Note that taxa found in the uppermost or lowermost bin in our section may be found before or after, but we do not have records of their occurrences outside this section. Morphotype and major plant group affinity of the taxon shown at bottom. Abbreviations: CON, conifer; DIC, dicot (non-monocotyledonous angiosperm); GIN, ginkgophyte; MON, monocotyledonous angiosperm; PTE, pteridophyte. Number of singletons per 15 m bin shown at right.

In general, last appearances were high throughout the HC interval; whereas the overall loss of dicot leaf morphotypes from the HC to TM was 72%, proportional last appearances from one HC flora to the next ranged from 46% to 54% (Fig. 7A). Last appearances peaked in the latest HC flora and dropped just after the KPB (Fig. 7A).

Figure 7. Turnover and diversity through our study interval within each of the 11 floras. A, Last appearances is the percentage of taxa that end their stratigraphic range in that flora. B, First appearances is the percentage of taxa that begin their stratigraphic range in that flora. Because the proportion of appearances (A, B) is influenced by the start and end of our sequence, we have excluded the oldest and youngest floras where appropriate. C, Range-through richness is the number of taxa with range-through occurrences in that flora. D, Rarefied richness is subsampling of abundance data to the smallest sample size (N = 55). On the left, the formations are indicated along with the magnetostratigraphy (based on Sprain et al. 2015, 2018), temporal intervals used in this study (Hell Creek [HC], Tullock Member [TM1, TM2]), and names of floral assemblages (as in Fig. 1). Values are plotted at the stratigraphic height of the corresponding flora, marker shape represents sedimentary facies, vertical lines on markers represent stratigraphic uncertainty, and horizontal lines represent 95% confidence intervals based on bootstrapping of data to N = 55 specimens. Analyses were conducted using dicot leaf morphotypes only.

Our analysis further indicates that a high proportion of first appearances across the KPB also played a substantial role in turnover; proportional first appearances peaked in the youngest HC flora (67%; Fig. 7B), reflecting high Cretaceous turnover. TM1 first appearances remained moderate (37–50%; Fig. 7C) as new TM taxa originated or immigrated.

Rarefied taxonomic richness was, on average, 23–33% greater in HC floras than in TM floras (considering all morphotypes, only leaf morphotypes, or only dicot leaf morphotypes; Table 2), although the differences are not significant (U-test p > 0.05; Table 2). More specifically, rarefied richness dropped 39–50% from HC to TM1 floras (U-test p > 0.05; Table 2) and range-through richness also dropped immediately across the KPB (Fig. 7C). Rarefied richness then increased by 35–61% from TM1 to TM2 floras (U-test p > 0.15; Table 2); range-through richness exceeded HC levels during the TM2 interval (Fig. 7C). Although none of these differences are statistically significant, the changes in both range-through and rarefied richness through time are large (Fig. 7C,D). Evenness was similarly highest in the HC, lowest in the TM1, and intermediate in the TM2 interval, although none of these changes are statistically significant (Table 2). In contrast, there is no discernible difference in rarefied richness or evenness by sedimentary facies (U-test p > 0.60; Table 2).

Rarefaction and extrapolation curves for these floras also point to a large drop in taxonomic richness at the KPB (Fig. 8). Average richness is higher in the HC than in the TM; in addition, TM2 floras are slightly richer than TM1 floras, both when using abundance of all leaf morphotypes (Fig. 8A) and when using abundance of only dicot leaf morphotypes (Fig. 8C). However, richness is highly variable from one flora to another (Fig. 8B); this variability among individual floras likely contributes to the nonsignificance of our results comparing rarefied richness grouped by temporal interval (Table 2). One HC flora in particular has extremely low taxonomic richness (Smurphy's Guess rarefied richness is 10.1 ± 1.2 SE; Fig. 8B). When we exclude this anomalous HC flora, the difference in average rarefied richness between HC (17.9) and TM (11.6) floras is significant (U-test p = 0.0333; leaf morphotypes only). In addition, richness rebounds relatively quickly in the earliest Paleocene, and our youngest TM1 flora (New York) has similar richness to many TM2 floras (New York = 9.7, TM2 average = 8.2; Table 2 dicot leaf morphotypes only, Figs. 7, 8). Given this rapid increase in taxonomic richness within the TM1 interval, comparing average HC to average TM1 taxonomic richness may not accurately capture the interval of greatest change in richness across the KPB. Considering just the immediate change across the KPB, rarefied richness dropped by 76% (leaf morphotypes) or 85% (dicot leaf morphotypes) from the youngest HC flora to the oldest TM1 flora (Figs. 7, 8).

Figure 8. Richness curves for each time interval and for each individual flora. A, Combined rarefaction and extrapolation curves based on leaf morphotype abundance data for floras in the Hell Creek (HC) and Tullock Member (TM1, TM2) intervals with shading for 95% confidence intervals to a target sample size of 3000. B, Rarefaction and extrapolation curves for each flora (based on leaf morphotype abundance data) to a target sample size of 600. Marker colors correspond to temporal interval (HC, TM1, or TM2); floral assemblage abbreviations as in Fig. 1. C, Combined rarefaction and extrapolation curves based on dicot leaf morphotype abundance data for floras in the HC, TM1, and TM2 intervals with shading for 95% confidence intervals to a target sample size of 3000. D, Rarefaction curves for HC, TM1, and TM2 intervals as in C, as well as rarefaction curves from North Dakota (data from Wilf and Johnson [2004], dicot species abundance data from the uppermost 15 m of Cretaceous strata and all Paleogene strata, truncated to 1500 specimens) and Patagonia (data from Stiles et al. [2020], dicot species abundance data from Cretaceous [Lefipán floras] and Paleogene [P Loros 1, P Loros 2, and LF] floras as in that study).

Discussion

Pattern of Floral Turnover across the KPB

Our results show that the K/Pg mass extinction led to a distinct change in the taxonomic composition and vegetation diversity in the Hell Creek study area of NE Montana, as we hypothesized (H1). This change in community composition is independent of depositional environment (facies) and instead reflects changes in vegetation across the KPB. Cretaceous vegetation was characterized by a distinct suite of taxa consisting of mostly angiosperms as well as several conifers and pteridophytes (e.g., “ferns”). Floras from the last ca. 1.3 Myr of the Cretaceous were, on average, taxonomically rich (high alpha diversity) and heterogeneous (high beta diversity, high last appearances) either spatially, temporally, or both. Work on North Dakota megafloras and Montana palynofloras previously suggested that the Cretaceous vegetation was dynamic (experiencing high turnover, going through several vegetation regimes) (Hotton 2002; Johnson 2002; Wilf and Johnson 2004). Our results generally support these findings, although further sampling in the Late Cretaceous is needed to test how similar the vegetation was to that in North Dakota and to reduce sampling bias, which is likely increasing our rate of last appearances.

At the KPB, we find evidence for a substantial drop in taxonomic richness and a significant loss of plant taxa, either due to extirpation or extinction. Between our HC and TM floras, there was a 23–33% drop in rarefied richness (Fig. 7D). We also document the disappearance of ~63% of latest Cretaceous taxa (from ca. 110 kyr before the KPB; a best estimate of extinction magnitude). These changes in taxonomic composition observed across the KPB were driven largely by a high proportion of first and last appearances just before the KPB (Fig. 7). We find no evidence that depositional taphonomic biases affected the taxonomic composition of our fossil assemblages (comparison of overbank vs. channel facies; Fig. 5, Table 2) and instead interpret the changes observed through this succession of floras as representative of local vegetation change through time. These results provide mixed support for our hypothesis that high disappearances contributed to turnover (H2); a combination of disappearances and appearances across the KPB resulted in the overall observed taxonomic turnover.

Post-KPB Floras and Taxonomic Recovery in the Early Paleocene

Earliest Paleocene (TM1) plant communities in NE Montana were species poor (Figs. 7, 8) and distinct in terms of the overall taxonomic composition (Figs. 4–6). Though many conifers persisted across the KPB, the relative abundance of those conifer taxa declined from 35% in the Late Cretaceous (HC) to 6% in the early Paleocene (TM1), with angiosperms becoming increasingly abundant in the relatively low-diversity TM1 floras. TM1 vegetation also includes several weedy or aquatic angiosperm taxa (e.g., Limnobiophyllum scutatum) (Figs. 4, 5), suggesting local environmental change in the Paleocene and/or that the mass extinction led to selection for wet-adapted taxa. Johnson (2002) and Hotton (2002) also both found that swamp or wet-adapted taxa were commonly Paleocene associated. The change in lithology from the HCF to the FUF has been interpreted as having resulted from changes in the local hydraulic flux (see “Geologic Setting of the Hell Creek Study Area”); the abundance of wet-adapted taxa may therefore be an indication of local environmental pressures.

Although TM1 floras were low in taxonomic richness and evenness, the relatively high species richness of TM2 floras does not support our hypothesis (H3) that the Paleocene was depauperate and homogeneous for an extended period. Rather, this pattern suggests that, locally, the post-KPB taxonomic recovery occurred relatively early in the Paleocene. Within 80 to 900 kyr after the KPB (TM2 floras), there was a change in taxonomic composition (novel angiosperms as well as the relative abundance of conifers rebounding to 25%) along with a ~55% increase in taxonomic richness, returning to pre–mass extinction diversity levels. That said, the youngest TM1 flora was also relatively species rich, indicating that taxonomic recovery may have begun already during the TM1 interval (Figs. 7, 8). This resembles the pattern of recovery among vertebrate faunas from NE Montana. Mammalian taxonomic richness and heterogeneity began increasing within 320 kyr of the KPB (Smith et al. 2018), and modern lineages (e.g., plesiadapiforms) arrived in NE Montana within ca. 150 kyr post-KPB, becoming a major component of the “fully recovered” mammalian fauna by ca. 847 kyr post-KPB (Wilson Mantilla et al. 2021). Future work to tighten our chronostratigraphic constraints on the TM1 and TM2 intervals would enable a more precise estimate of the timing and rate of vegetation recovery, as well as its correlation to faunal recovery.

Taxonomic richness is only one measure of biotic recovery from the K/Pg mass extinction. Previous authors have investigated whether plant ecological diversity (studied through, e.g., leaf mass per area) increased during the early Paleocene to test for additional signals of recovery (Blonder et al. 2014; Butrim et al. 2022). Our ongoing research into these additional measures of plant diversity across the KPB in NE Montana might point to alternative signs of recovery among plant communities and help evaluate drivers of post-KPB vegetation change.

Comparison with Global K/Pg Records

Overall, these results indicate that plant communities in NE Montana experienced similar rates of taxonomic turnover during the K/Pg mass extinction compared with other studied regions of WINA (Table 3, Fig. 3). In the Williston Basin of North Dakota and the Denver Basin of Colorado, 46–57% of Cretaceous megafloral species disappeared (Wilf and Johnson 2004; Lyson et al. 2019; Table 3), compared with our 63% disappearance rate immediately across the KPB. Conversely, palynotaxa experienced a lower rate of extinction (15–30% in the Williston Basin; Hotton 2002; Nichols 2002; Nichols and Johnson 2008), likely owing to the lower taxonomic and spatial resolution of palynomorph records (Cleal et al. 2021). Even though our study has fewer floras from this latest Cretaceous interval, our median sample size per flora is greater than in North Dakota or Colorado (Table 3), and our youngest Cretaceous flora is particularly well sampled (Bruce Leaf; Fig. 8B). Therefore, our sample size should be sufficiently large for a robust comparison among north-central WINA study areas; however, our Cretaceous temporal resolution is not sufficient to speculate on the specific causal mechanism(s) of the mass extinction. Relative to WINA, vegetation in Patagonia experienced an exceptionally high disappearance rate (up to 90%; Stiles et al. 2020). However, the studies in Patagonia, as well as those in Colombia (Carvalho et al. 2021), are based on floras distributed across large geographic distances (ca. 300–400 km) and are only coarsely dated to the Maastrichtian/Late Cretaceous and Danian/early Paleocene (Fig. 3). Therefore, we cannot exclude the possibility that this difference in estimated extinction magnitude between WINA and South America is due to differences in the scale of spatial and temporal sampling.

Table 3. Summary of sampling, extinction magnitude, and change in species richness (where available) associated with the mass extinction across well-studied regions around the globe (see Fig. 3 for additional context). (1) Data from Wilf and Johnson (2004) and Peppe (2010); (2) data from Lyson et al. (2019); (3) data from Flynn and Peppe (2019); (4) approximated from Dunn (2003: fig. 5.3); (5) data from census sites studied by Carvalho et al. (2021), note that these authors estimated palynofloral extinction magnitude (46%) but not macrofloral extinction magnitude; (6) data from Stiles et al. (2020).

The drop in taxonomic richness we observe from the HC to TM1 temporal interval, although substantial (33% considering dicot leaf morphotypes only), was less severe than in Patagonia (40%; Stiles et al. 2020), North Dakota (~49% comparing Paleocene floras with the latest 15 m of Cretaceous strata; Wilf and Johnson 2004), and Colorado (~50%; Lyson et al. 2019). However, these studies vary in their temporal range and resolution when calculating this drop in diversity (Table 3, Fig. 3), so directly comparing species richness from one study to another is difficult. In general, Cretaceous floras in Montana were less diverse than those in Patagonia and North Dakota (Fig. 8D); this may in part explain the lower diversity loss across the KPB in Montana.

Our NE Montana floras also indicate differences in the timing of taxonomic recovery within the WINA. NE Montana floras returned to pre–mass extinction taxonomic richness levels relatively quickly, by 900 kyr after the KPB, on par with Colorado (Lyson et al. 2019). In contrast, North Dakota floras remained depauperate for 10 Myr after the KPB (Wilf and Johnson 2004), even as they changed taxonomically during the early and middle Paleocene (Peppe 2010). In other basins of north-central WINA (e.g., Bighorn and Hanna Basins) Paleocene megafloral diversity was generally low, similar to North Dakota (Hickey 1980; Wing et al. 1995; Dunn 2003). Differences in the pattern of first appearances also point to a more rapid re-diversification in the earliest Paleocene in Montana. In North Dakota, only 11 nonreproductive morphotypes have first appearances in the early Paleocene (Wilf and Johnson 2004), whereas in Montana, we find 46. These results highlight the significance of origination and/or immigration in the rapid re-diversification of Montana vegetation post-KPB. In contrast to our megafloral results, Montana palynofloral richness remained depressed throughout studied Paleocene sections (Hotton 2002). This discrepancy could reflect differences in taxonomic resolution: genus- and family-level extinction both were lower but also slower to recover (as documented by palynofloras), whereas species-level extinction was higher but also faster to recover (as documented by megafloras).

The differences in floral recovery across WINA may have been due to differences in environment (e.g., local climate) or biotic factors (e.g., variation in vegetation resilience). For example, recovery rates in Colorado have been linked to elevation (Johnson et al. 2003). In this scenario, higher-elevation areas may have served as refugia for Cretaceous taxa (e.g., Johnson and Ellis 2002) and/or topographic (hence ecological) variation in high-elevation areas may have promoted rapid speciation (see, e.g., Antonelli et al. 2018). Based on the relative position of the Western Interior Seaway at this time (although debated; see, e.g., Slattery et al. 2015), we suggest that NE Montana would have been at a slightly higher elevation than North Dakota. If true, heterogeneous microhabitats in the upland may have allowed some NE Montana plant taxa to persist during the mass extinction and promoted accelerated speciation in the early Paleogene. However, without quantitative estimates of the paleoelevation in these study areas, this inference remains speculative. Alternatively, some aspect of plant community diversity or composition (e.g., in functional groups) may have allowed for higher resilience or faster recovery in NE Montana as compared with North Dakota. For example, the abundance of persistent conifer taxa in NE Montana may have lent resiliency to the local vegetation. Although conifers in both NE Montana and North Dakota show declines in abundance across the KPB (from 17% in the HCF to 4% in the FUF in North Dakota; Wilf and Johnson 2004), conifers were consistently more abundant in NE Montana during this time period (35% of HC, 6% of TM1, and 25% of TM2 interval specimens). More detailed ecological information about these taxa is necessary to evaluate this hypothesis. NE Montana palynofloral records, which sample a broader spatial area and therefore may capture the upland environments, document a large increase in the proportion of conifer pollen after the KPB (Hotton 2002). Conifers were perhaps restricted to upland or other refugia during the TM1 interval before rebounding in abundance in TM2 megafloras. Thus, the K/Pg mass extinction may have resulted in different patterns of survival in upland versus lowland environments and among conifers versus angiosperms, although further investigation of these patterns is necessary.

Similar to WINA, floral records in South America show variation in recovery rates. In Patagonia, Paleocene floras from within ca. 5 Myr after the KPB were more morphologically diverse than Cretaceous floras, potentially a sign of increasing ecological diversity, even as they remained taxonomically homogeneous (Iglesias et al. 2007, 2021; Stiles et al. 2020). Plant communities in Colombia did not fully recover for ca. 6 Myr after the KPB (Carvalho et al. 2021). Our view of the K/Pg mass extinction remains limited by the few regions with continuous plant fossil records; therefore, further work to elucidate global patterns in extinction magnitude and recovery is required to understand what led to these widely differing floral responses to the K/Pg mass extinction.

Implications of the Floral Record from NE Montana for the Study of Mass Extinctions

Our study points to regional and global variability in the impact of mass extinctions on plant communities, with susceptibility to the mass extinction and recovery likely shaped by local extrinsic (e.g., climate) and intrinsic (e.g., plant ecological strategies) factors. Our results therefore suggest that modern biodiversity loss must be assessed and addressed at a local level to predict species loss as well as potential for future recovery.

Our results are also consistent with the framework described by McElwain and Punyasena (2007) and Green et al. (2011) in which plants respond to mass extinctions by shifts in relative abundance and species-level extinction or extirpation, rather than global-, family-, or even genus-level extinction. Specific ecological strategies (e.g., fast growth, seed dormancy), common among many plant clades, may have conferred resiliency for some, if not many, species. Furthermore, ecological and physiological differences between plant groups might affect how they respond to mass extinctions within a region. For example, we find that conifer taxa were more likely to survive the K/Pg mass extinction, even though their abundance decreased in the earliest Paleocene (TM1). Conifers, having longer times between regeneration as well as specific adaptations to disturbance (e.g., thicker bark in many Pinaceae and Cupressaceae taxa; Brodribb et al. 2012) may have preferentially survived during the immediate environmental stressors at the KPB, even though their populations declined in abundance (perhaps restricted to refugia). In contrast, the lower abundance (smaller population size, smaller geographic range) of Late Cretaceous angiosperm taxa may have led to their higher disappearance rate during the mass extinction (Jablonski 2005, 2008; Payne and Finnegan 2007; Heim and Peters 2011). Later, during the TM1 interval, angiosperms may have increased in abundance because of ecological advantages such as faster seedling growth rates (Bond 1989), rapid reproduction (Bond 1989; Doyle and Hickey 1976; Hickey and Doyle 1977; Verdu 2002), higher hydraulic capacity and vein density (Lusk et al. 2003; Boyce et al. 2009; Brodribb and Feild 2010; Feild et al. 2011), and the ability to spread vegetatively (Boyce and Leslie 2012). These characteristics, which are thought to have promoted the success of angiosperms generally during the Cretaceous (Bond 1989; Augusto et al. 2014), may have been crucial during the period just after the mass extinction event (e.g., as seen in Colombia; Carvalho et al. 2021). Previous authors (e.g., Wing and Tiffney 1987) have similarly suggested that the environmental crisis at the KPB disturbed the dominance of stress-tolerant and competitive taxonomic groups (i.e., conifers), allowing for ruderal groups (i.e., angiosperms) to flourish (Wing and Boucher 1998). These results suggest fundamental differences in the response of conifers and angiosperms to disturbance and biodiversity crises. Future research on K/Pg plant ecology will help to elucidate the potential for distinct impacts on different plant taxa or functional groups during mass extinctions, with important implications for our modern biodiversity crisis.

Conclusions

Our study finds that plant community composition changed markedly leading up to and across the KPB in NE Montana. Turnover was high during the Late Cretaceous, culminating in the disappearance of 63% of latest Cretaceous plant taxa and the appearance of a novel suite of Paleocene plant taxa. Floras from the first 80 kyr of the Paleocene were 23–33% less rich than Late Cretaceous floras on average, with an increase in the representation of wet-adapted taxa. Between 80 and 900 kyr after the mass extinction, plant species richness rebounded to Late Cretaceous levels again. This pattern of rapid recovery echoes what has been described for mammalian communities from the same region (e.g., Smith et al. 2018; Wilson Mantilla et al. 2021; Claytor et al. 2022).

We find broad accord in terms of extinction magnitude (ca. 46–63% disappearance) across WINA, but a lower drop in taxonomic richness at the KPB in NE Montana compared with elsewhere in the region. This variation in extinction magnitude may result from differences in sampling intensity or reflect local variation in vegetation susceptibility to the environmental changes associated with the mass extinction event. Our results also support a growing body of evidence that post-KPB recovery may have varied based on latitude, local paleoenvironment, or underlying vegetation dynamics (e.g., Wilf and Johnson 2004; Lyson et al. 2019; Stiles et al. 2020; Wilf et al. 2023). Recovery in Montana was notably much faster (within 900 kyr) than in nearby North Dakota (millions of years; Johnson 2002; Wilf and Johnson 2004); these results are not driven by local depositional environment and taphonomic biases. Future work to integrate analyses of plant community diversity with plant ecology, local environment, and vertebrate paleontological records from this study area will enable a more nuanced view of how terrestrial communities are impacted by mass extinction events on shorter and longer timescales.