Significance Debate lingers over what caused the last mass extinction 66 million years ago, with intense volcanism and extraterrestrial impact the most widely supported hypotheses. However, without empirical evidence for either’s exact environmental effects, it is difficult to discern which was most important in driving extinction. It is also unclear why recovery of biodiversity and carbon cycling in the oceans was so slow after an apparently sudden extinction event. In this paper, we show (using boron isotopes and Earth system modeling) that the impact caused rapid ocean acidification, and that the resulting ecological collapse in the oceans had long-lasting effects for global carbon cycling and climate. Our data suggest that impact, not volcanism, was key in driving end-Cretaceous mass extinction.

Abstract Mass extinction at the Cretaceous–Paleogene (K-Pg) boundary coincides with the Chicxulub bolide impact and also falls within the broader time frame of Deccan trap emplacement. Critically, though, empirical evidence as to how either of these factors could have driven observed extinction patterns and carbon cycle perturbations is still lacking. Here, using boron isotopes in foraminifera, we document a geologically rapid surface-ocean pH drop following the Chicxulub impact, supporting impact-induced ocean acidification as a mechanism for ecological collapse in the marine realm. Subsequently, surface water pH rebounded sharply with the extinction of marine calcifiers and the associated imbalance in the global carbon cycle. Our reconstructed water-column pH gradients, combined with Earth system modeling, indicate that a partial ∼50% reduction in global marine primary productivity is sufficient to explain observed marine carbon isotope patterns at the K-Pg, due to the underlying action of the solubility pump. While primary productivity recovered within a few tens of thousands of years, inefficiency in carbon export to the deep sea lasted much longer. This phased recovery scenario reconciles competing hypotheses previously put forward to explain the K-Pg carbon isotope records, and explains both spatially variable patterns of change in marine productivity across the event and a lack of extinction at the deep sea floor. In sum, we provide insights into the drivers of the last mass extinction, the recovery of marine carbon cycling in a postextinction world, and the way in which marine life imprints its isotopic signal onto the geological record.

There is abundant evidence for a massive bolide impact at the Cretaceous–Paleogene (K-Pg) boundary (66.04 Ma) that coincides in those sections with sufficient temporal resolution with a mass extinction horizon (e.g., ref. 1). Geological records document many effects from impact, including massive tsunamis, earthquake-driven gravity flows, and molten ejecta fallout (e.g., ref. 1). However, the mechanism (or mechanisms) by which impact drove global-scale ecosystem turnover and mass extinction is less certain. Among the most prominent hypotheses are global darkness and associated primary productivity loss leading to food chain collapse, acid rain, impact winter, and flash ocean acidification (2). Some of these mechanisms are supported by modeling work (e.g., ref. 3) but, critically, they generally lack empirical validation. Furthermore, the issue is complicated by the possibility of contributing effects from ongoing or intensified Deccan flood basalt volcanism (4, 5) and the short timescales of at least some of the effects of impact compared to the resolution of the geological record.

Here we apply the boron isotope-pH proxy to planktic and benthic foraminifera to examine whether ocean acidification occurred at the K-Pg boundary and to test competing hypotheses that have been proposed to explain changes in the marine carbon cycle in the aftermath of the mass extinction. We analyzed material from multiple sites where sufficient well-preserved foraminifera from the mixed layer of the surface ocean could be obtained, with >7,000 small, thin-walled postextinction foraminifera required per analysis in some cases. By combining measurements from continental shelf sites in communication with the open ocean [Geulhemmerberg Cave, The Netherlands (6), Owl Creek, Mississippi (7), and Brazos River, Texas (8)] and deep-sea drilling sites (Deep-Sea Drilling Program Site 465 and Ocean Drilling Program [ODP] Site 1209 in the Pacific and ODP Site 1049 in the Atlantic), we construct a composite global record spanning ∼800 ky (SI Appendix). Importantly, one shallow marine sample site studied here (Geulhemmerberg Cave) is thought to preserve sediments from the first 102 to 103 y after bolide impact (6), allowing us to pinpoint a short-lived signal of impact-induced ocean acidification that cannot be temporally resolved in deep marine records. Samples from ODP Site 1209 allow us to investigate longer-term changes in the ocean’s biological carbon pump as expressed in vertical pH gradients, as this site is one of few deep marine K-Pg boundary sites that preserves sufficient numbers of well-preserved epifaunal benthic and surface-dwelling planktic foraminifera for boron isotope analysis. The broad range of in situ temperatures, hydrostatic pressures, and environmental conditions represented in our sample set allows us to constrain bulk seawater δ11B (δ11B sw ; see SI Appendix), one of the main uncertainties in calculating pH from deep-time δ11B measurements. Assuming feasible limits on the calcium carbonate saturation state (Ω CaCO3 ) and oxygen utilization at each site (see SI Appendix, section 5C), we constrain K-Pg δ11B sw to within 39.05 to 39.85‰.

Our data suggest ocean pH was stable over the last 100 ky of the Cretaceous, despite suggestions of a major pulse of Deccan volcanism within this time interval (5). However, our data from the Geulhemmerberg Cave boundary clay indicate a marked ∼0.25 pH unit surface ocean acidification event within a thousand years (6) of the bolide impact in the Geulhemmerberg Cave boundary clay (Fig. 1C). These data provide geochemical evidence of rapid acidification in the immediate aftermath of the Chicxulub impact. Combined with model-derived estimates of alkalinity (SI Appendix), this change in pH corresponds to a rise in atmospheric partial pressure of CO 2 (pCO 2 ) from ∼900 ppm in the latest Maastrichtian to ∼1,600 ppm in the immediate aftermath of bolide impact (Fig. 1D). This measured pH drop may be a conservative estimate of the full magnitude of impact-induced acidification, as the earliest Ir-rich ejecta layer at Geulhemmerberg is thought to have been eroded within ∼100 y of the impact (6). A geologically rapid 0.2 to 0.3 pH unit change would have disadvantaged calcifying plankton vs. noncalcifiers (9) and could therefore explain the selective extinction of calcifying pelagic organisms during the K-Pg mass extinction (10⇓⇓–13) compared to noncalcareous groups such as dinoflagellates (e.g., ref. 14) and radiolarians (e.g., ref. 15). Similarly, since coccolithophores and planktic foraminifera from coastal areas may be better adapted to fluctuations in pH (16), and smaller foraminifera may better maintain calcification under low pH (17), acidification could also help to explain on-shelf/off-shelf and size-selective patterns of extinction among calcifiers (18, 19). The observed magnitude of acidification is compatible with the rainout of SO 2 and NO x released by the impactor (10, 20, 21), depending on the timing of sediment deposition at Geulhemmerberg. The effect of rainout products was transient (e.g., ref. 21), and if Geulhemmerberg clays were deposited on timescales greater than a few hundred years, a contribution to acidification from CO 2 release would have been required. This additional CO 2 could be sourced from carbonate target rocks (ref. 1 and references therein), the biosphere [through wildfires and decay (ref. 1 and references therein)], a transient weakening of the biological carbon pump (ref. 22 and references therein and SI Appendix), and/or a contribution from a pulse of volcanism (5).

Fig. 1. Records of surface ocean foraminiferal δ11B (B), calculated pH (C), and pCO 2 (D) across the K-Pg boundary, with high-resolution foraminiferal diversity counts from ref. 39 at the K-Pg Global Boundary Stratotype Section and Point (El Kef) plotted for context (A). pH is calculated assuming our best estimate of K-Pg δ11B sw , 39.45 ± 0.4‰. pCO 2 is calculated from pH along with total alkalinity estimates at each site from a GENIE late Maastrichtian simulation, adjusted for dynamic changes in alkalinity across the K-Pg using LOSCAR simulations from ref. 22 that match observed patterns of carbonate burial. Gray shaded areas are 1-sigma uncertainties, with thin lines representing 1,000 Monte Carlo simulations from the 10,000 that were run. For clarity, we only plot those samples that represent the surface mixed layer, which should be approximately in equilibrium with the atmosphere. Additional data from deep-sea benthic and thermocline-dwelling planktic foraminifera that do not reflect atmospheric pCO 2 can be seen in Fig. 2 and in Dataset S1. For details of the age models used, the estimation of δ11B sw , carbonate system calculations, and uncertainty propagation see SI Appendix.

Loss of plankton functional groups with the K-Pg impact has profound implications for marine carbon cycling (22), manifest in one of the earliest and most striking geochemical observations of the K-Pg boundary in open ocean sediments: a convergence in the δ13C values of surface- and deep-dwelling calcifiers (e.g., ref. 23). Paired measurements of δ11B (and hence pH) from surface and deep water after the K-Pg boundary allow us to address a long-running, unresolved debate about the cause of this convergence, and the long-term effect of the K-Pg mass extinction on marine carbon cycling. Surface ocean δ13C is typically maintained at heavy values relative to the deep sea through the biological pump, which exports isotopically light organic carbon from the surface ocean to depth. The sharp drop in surface ocean δ13C at the K-Pg was therefore initially interpreted as a near-total loss of marine new primary productivity after impact, in the “Strangelove Ocean” hypothesis (24). Later, the “Living Ocean” hypothesis (25) posited that such a reduction of vertical δ13C gradients arose instead from a decrease in the efficiency of the biological pump, that is, continued productivity but reduced sinking efficiency leading to more remineralizaton of organic carbon in the upper water column. Evidence for export to the deep ocean persisting in at least some regions led to the “Heterogeneous Ocean” hypothesis (11, 26), while it has also been argued that part of the reduction in δ13C gradients stems from the fact that carbon isotope vital effects differ in postextinction carbonate producers compared to the preextinction ones (11, 27). These latter scenarios appear more compatible with low levels of benthic foraminiferal extinction at the K-Pg (11, 28) than a near-total loss of primary or export production. Another alternative hypothesis not linked to any marine biogeochemical response is that extremely high atmospheric pCO 2 (thousands of parts per million) could have driven convergence in deep- and surface-ocean δ13C inorganically (29). Such a scenario might arise if outgassing from the Deccan Traps was accelerated following the impact (4).

Our records show that initial surface acidification was followed by a rapid rebound and overshoot in surface ocean pH (∆pH ≈ 0.5) within 40 ky of the K-Pg boundary, returning to preevent values after an additional ∼80 ky (Fig. 1A). This is consistent with model predictions of ocean pH rise due to the extinction of calcifying organisms and a consequent transient reduction in the marine alkalinity sink (22, 30). Such a period of imbalance in marine alkalinity fluxes in the earliest Danian is supported by observations of improved deep-sea carbonate preservation (22), the occurrence of heavily calcified deep-sea benthic foraminifera (31), and the deposition of abiogenic calcite crystals in the oceans (32). The absence of any prolonged (>50 ky) surface ocean acidification after the K-Pg suggests, therefore, that if there were a state shift toward significantly more active Deccan CO 2 degassing in response to impact (4), its effects were overwhelmed by the biogeochemical effect of extinction. Other independent pCO 2 proxy records (e.g., ref. 33) likewise do not support a substantial rise in CO 2 , and there is no evidence for deep-sea warming at this time (34). Thus, we discount high pCO 2 as the driver of reduced vertical δ13C gradients (29) and infer that the convergence in surface and deep-ocean δ13C arose from a change in oceanic biogeochemical cycling, combined to some extent with changes in carbonate vital effects.

To explore how changes in the carbon cycle would manifest in vertical pH and δ13C gradients, we used GENIE (Grid ENabled Integrated Earth System Model) (35), with a late Maastrichtian bathymetry and continental configuration (SI Appendix, Fig. S17 and SI Appendix). We simulate the response of marine δ13C and pH given either reduced or zero primary productivity (a “Strangelove Ocean”-type scenario) or less efficient biological pump (i.e., enhanced remineralization in the upper water column, “Living Ocean” scenario; full details are given in SI Appendix). For each simulation, we show depth transects for δ13C and pH at Shatsky Rise (ODP Site 1209) (Fig. 2 B and C). Pacific basin δ13C zonal means for 3 illustrative states of biological pumping are shown in Fig. 3. We use the observed patterns of boron and carbon isotope change to distinguish between the competing hypotheses for postextinction biogeochemical change. The fingerprint of the “Strangelove” scenario with a complete loss of productivity is a reversed vertical δ13C gradient with surface ocean δ13C values >1‰ lighter than benthic (Fig. 2B) and a muted but positive surface-deep pH gradient (Fig. 2C). The “Living Ocean” scenario with shallower organic matter remineralization leaves a muted but positive surface-deep δ13C gradient, but with a flattened pH gradient. Only with a 50% reduction in global new primary productivity can our model produce a flattened vertical δ13C gradient together with a positive pH gradient.

Fig. 2. Paired benthic–planktic foraminiferal δ11B measurements from 4 time slices following the K-Pg boundary at ODP Site 1209 (Shatsky Rise) are used to calculate the difference in surface and deep ocean pH (∆pH; red diamonds, A). These data show a trend from “normal” patterns of higher δ11B/pH in the surface ocean relative to the deep until 65.92 Ma, where this gradient disappears (see also SI Appendix, Fig. S12). Shown for context is the convergence in δ13C from surface waters (bulk CaCO 3 , in green) and deep ocean waters (benthic foraminifera, navy) (A). Water column δ13C and pH profiles in the cGENIE model the location of Shatsky Rise (see also Fig. 3) are shown in B and C, respectively. A weakened or reversed δ13C gradient can arise in a number of ways (B), but only a shallow remineralization (or “Living Ocean,” red dot-dashed line) model can produce a flat gradient in pH between the surface 80 m and 2,500-m water depth (C). Carbonate data from ODP 1209 include new data and benthic data from ref. 34. See SI Appendix for more details.

Fig. 3. Pacific zonal means of the simulated depth distribution of δ13C DIC in the Maastrichtian ocean for (A) the default (full-strength) biological pump, (B) global export production reduced to 50% of the default, and (C) no biological pump. The approximate paleo-latitude of ODP Site 1209 (Shatsky Rise) is marked by the vertical dashed white line. The action of the solubility pump, which absorbs atmospheric CO 2 at high latitudes due to increased solubility of surface waters with respect to CO 2 as they cool, imparts an inverse vertical δ13C DIC gradient at low latitudes in the absence of biological productivity (C; see the main text for explanation).

The “Strangelove” scenario produces a pronounced >1‰ reversed vertical δ13C gradient (Fig. 2B) rather than the observed approximately flat gradient. This phenomenon arises because of the solubility pump, that is, the increased solubility of CO 2 (and hence CO 2 uptake) in cooling surface waters as they reach high latitudes, where they ultimately sink to form deep water. Because the fractionation between dissolved CO 2(aq) (and thus atmospheric pCO 2 ) and aqueous bicarbonate (and thus dissolved inorganic carbon, DIC) increases by about 0.1‰ per °C cooling, the surface δ13C DIC at colder, higher latitudes is driven higher with respect to the atmosphere compared to the warmer, low-latitude ocean surface (Fig. 3C). Although this effect of the solubility pump is always there, typically the action of the biological pump overprints the pattern of the solubility pump to give the observed “normal” sign of δ13C, with most positive values at the surface (Fig. 3A and SI Appendix). Somewhere between a “Strangelove”-style shutdown and a full-strength biological pump, there must exist a weaker state of the biological pump that would impart a δ13C gradient opposite in sign and equal in magnitude to the solubility pump, canceling out the difference in δ13C between surface and deep ocean. We find this occurs with a 50% reduction in global export production for all but the highest-latitude Northern Pacific, where downwelling waters overlay an older (and lower δ13C) lens of southern-sourced water (Figs. 2B and 3B). This partially maintained productivity is consistent with the survival of benthic foraminifera at the K-Pg (11, 28).

The combination of our cGENIE simulations and the observed sequence of changes in vertical pH and δ13C gradients at Shatsky Rise (Fig. 2A) allow us to propose a hypothesis for the behavior of the marine carbon cycle after the K-Pg impact. During the first tens of thousands of years following impact, reduced vertical δ13C gradients without an associated collapse in pH gradients (Fig. 2A) are consistent with a substantial reduction in new primary productivity globally following mass extinction (Fig. 2 B and C and SI Appendix). The magnitude of productivity loss required depends on changes in carbon isotope vital effects in pelagic CaCO 3 producers, which likely account for some (but not all) of the observed collapse in δ13C gradients (11, 27, 36). Reduced global export productivity can, however, align with published observations of a “Heterogenous Ocean” with stable or increased productivity at certain sites [e.g., Shatsky Rise (11, 26)]. A global reduction in new primary productivity may have led to a buildup of nutrients in the surface ocean, and thus some locations (e.g., highly oligotrophic gyres) may have supported increased new/export primary productivity in the Danian relative to their preimpact state, even while export productivity globally was reduced (SI Appendix, Figs. S23 and S24).

After this period of globally reduced, but not collapsed, export productivity, we see a convergence in shallow and deep-ocean pH in our final time-slice ∼120 ky post-K-Pg (Fig. 2A). In our model scenarios, only by remineralizing organic carbon shallower in the water column can one produce similar pH values in both the deep sea and the subsurface (Fig. 2C). Shallow water remineralization alone (i.e., a “Living Ocean” scenario) does not fully collapse the δ13C gradient in our model, but combined with known carbon-isotope vital effects that drive surface-derived carbonates lighter than ambient δ13C DIC (27, 36, 37) the expression of such a scenario in the fossil record could be a full collapse. We therefore suggest that after an initial period (up to 40 ky) of globally partially weakened new primary productivity globally, it progressively began to recover, but without full restoration of pelagic ecosystem structure and function. In such a scenario, persistently ineffective ballasting of marine organic matter or other factors could have then resulted in more organic material being remineralized shallower in the water column [i.e., a “Living Ocean” scenario (25)], thereby reducing δ13C and pH gradients between the surface and deep sea.

Reduced vertical δ13C CaCO3 gradients persisted for over 1 My after the K-Pg impact (e.g., ref. 27), which suggests that the efficiency of organic carbon export to the deep sea, a key pelagic ecosystem function, was reestablished slowly. Our δ11B data do not support prolonged surface ocean acidification, making it difficult to attribute this delayed biogeochemical recovery to acidification from sustained (or enhanced) CO 2 degassing from the Deccan Traps. Rather, there may be intrinsic constraints on the time required to recover normal marine ecosystem function after such severe global perturbations, despite the short generation times that should make marine plankton ideally suited to rapid evolutionary radiation (38). In this way, the K-Pg event shows that even geologically rapid ocean acidification events can have prolonged and profound biotic ramifications.

Acknowledgments We thank the Integrated Ocean Drilling Program for supplying samples for this study. We thank Jan Smit and the organizers of the 12th International Conference on Paleoceanography Maastricht Excursion for facilitating sampling at Geulhemmerberg. M.J.H. thanks Donald Penman, David Evans, Simon D’Haenens, and Gavin Foster for helpful discussion; Dan Asael, Leanne Elder, Charlotte Bryan, Volkan Özen, and Myles Henehan for assistance in the laboratory; Jocelyn Sessa for provision of sample material; and Friedhelm von Blanckenburg for comments on an earlier draft. We thank the handling editor and 2 anonymous reviewers for their constructive feedback. M.J.H. acknowledges funding from the Yale Peabody Museum. A.R. acknowledges support from the European Research Council (ERC) as part of the “PALEOGENiE” project (ERC-2013-CoG-617313). L.A. acknowledges financial support by the Spanish Ministry of Economy and Competitiveness and European Regional Development Funds (project CGL2017-84693-R). D.N.S. acknowledges support from the Royal Society in the form of a Wolfson Merit award. J.W.B.R. was supported by ERC Starting Grant 805246 OldCO2NewArchives. J.D.W. received support from a Lerner-Gray Fellowship at the Richard Gilder Graduate School and the American Museum of Natural History (AMNH). J.D.W. and N.H.L. also gratefully acknowledge support from the Newell Fund (AMNH). S.E.G. was supported by Natural Environment Research Council (NERC) Independent Research Fellowship NE/L011050/1 and NERC large grant NE/P01903X/1 while working on this paper.

Footnotes Author contributions: M.J.H., A.R., E.T., D.N.S, and P.M.H. designed research; M.J.H., A.R., S.Z., L.A., J.W.B.R., and J.R.S. performed research; M.J.H. analyzed data; M.J.H. and S.Z. picked forams; J.D.W. and N.H.L. contributed samples and stratigraphic context; S.E.G. honed the model; B.T.H. trained taxonomy; N.J.P. and P.M.H. hosted the laboratory work; and M.J.H., N.J.P., and P.M.H. wrote the paper with contributions from all authors.

The authors declare no competing interest.

This article is a PNAS Direct Submission.

This article contains supporting information online at www.pnas.org/lookup/suppl/doi:10.1073/pnas.1905989116/-/DCSupplemental.