Environment and sample selection

Samples are from the W.J. Zinsmeister collection, formerly housed at Purdue University, and now located at the Paleontological Research Institute in Ithaca, NY, USA. Specimens were collected from Seymour Island from the López de Bertodano and Sobral Formations, covering the Late Cretaceous to early Paleocene. The stratigraphic framework used in this study is described by Zinsmeister et al.50 and has a stratigraphic resolution of ±3.5 m. Samples come from Units 9 and 10 of the López de Bertodano Formation, part of the upper beds (Units 7–10, Molluscan Units), as defined by Macellari8. The Molluscan Units contain diverse and abundant bivalve, gastropod, and ammonite macrofauna. Since the Cretaceous, these sediments are thought to have experienced only mild diagenetic conditions51. The depositional environment is interpreted to be inner to outer shelf, a semi-enclosed or quiet environment, with water depths less than 200 m, on the basis of faunal assemblages, grain size and sorting analysis, diversity of macrofauna, and abundance of glauconite8. Interpretations of water depth variation through the López de Bertodano Fm. disagree, with some suggesting a shallowing-upward52 and others a deepening-upward8 trend through the section. Overall, the homogeneous nature of the Molluscan Units 7–9 suggests8,42,52, at most, minor variation in water depth over this interval. Although there is some evidence of methane seeps on Seymour and the surrounding islands, it is restricted to lower in the stratigraphy (Units 1–6)53. There is no indication of methane seeps in the upper López de Bertodano Formation on Seymour Island43.

Bivalve shells (116) were selected from the Zinsmeister collections at Purdue for use in the PhD thesis of author Andrea Dutton54 with the aim of including as many stratigraphic horizons as possible for maximum temporal resolution (the ‘full sample set’). Samples were selected from six species (L. larseni, C. antarctica, C. ellioti, E. regina, D. drygalskiana and N. nordenskjoldi) to provide internal consistency between horizons. All these species were thought to live in shallow waters and be infaunal suspension feeders, with the exception of Nordenskjoldia, which is thought to be an epibyssate suspension feeder8.

A subset of 29 shells were selected from the full sample set for clumped isotope analysis, representing five of the six species in the original study (eliminating D. drygalskiana; the ‘clumped isotope sample set’). These were again selected for maximum stratigraphic coverage and species overlap, but were filtered to only include the best-preserved samples, as indicated by low trace element concentrations (see below). A list of samples in the full sample set and those sub-selected for clumped isotope analysis can be found in Supplementary Data 1.

Sample preparation and extraction

Bulk samples of 50–100 μg were drilled from the umbo, mid-shell and ventral margin of all shells in the full sample set for stable isotope and trace element analysis. Samples were drilled to average over multiple growth bands to represent typical mean annual conditions. The umbo position was centred near the cusp of the umbo. The ventral margin (‘back’) was within 1 cm of the edge of the shell. The mid-shell position was half way between the umbo and ventral margin. In the case of broken, partial shells, it was not always possible to sample the ventral margin and/or mid-shell position. Broken and partial shells were not selected for the clumped isotope sample set.

The 29 shells selected for clumped isotope analysis were re-drilled in two positions, the umbo and ventral margin. Sampling locations were directly adjacent to earlier drill locations. Sample size requirements for clumped isotope analysis are much larger than for stable isotope analysis11, so 20–30 mg of sample material was collected at each position and homogenized. During drilling for clumped isotope analysis, the drill was operated at the lowest possible rotation speed (1,000 r.p.m.) and pressure was only applied for a few seconds at a time to minimize frictional heating. Based on tests with other samples in our lab, at this drill speed, we do not expect to see any resetting effects. Most shells were drilled from the exterior at a position that was clear of external encrustation or fractures. In a few cases, shells were cut in half and were drilled on the cut interior surface. No correlation was seen between measured values and interior versus exterior drilling. Six shells of the genus Eselaevitrigonia were only sampled in the umbo position because invariant stable isotopic compositions between positions suggested that identical environmental conditions were being recorded at all points on the shell (see below), and some shells were broken such that the ventral margin was no longer present.

Trace element analysis

Elemental analysis of Sr, Mg, Ca, Fe and Mn was carried out on powdered splits of all bulk samples using either an inductively coupled plasma optical emission spectrophotometer (ICP-OES) or an inductively coupled plasma mass spectrometer (HRICP-MS) depending upon the amount of powdered carbonate available for analysis. Analysis was performed in the Keck Elemental Geochemistry Laboratory at the University of Michigan (Director: K.C. Lohmann). Smaller samples were run on HRICP-MS owing to the lower detection limits attainable on this instrument. By comparison with lab standards and of replicate analyses, precision of both types of analyses was maintained at better than 3%. Because these samples were too small to allow for efficient weight determination, calculation of cation concentrations within the carbonate was achieved by assuming that Sr, Mg, Ca, Fe, Mn are carbonate-bound cations and that these five elements account for 100% of the carbonate-bound cations. Trace element data can be found in Supplementary Data 1.

Assessment of shell preservation

Many lines of evidence (lack of smectite-illite transition, mineralogical and vitrinite-reflectance data) suggest a maximum burial depth of less than 1 km and a maximum burial temperatures of less than 80 °C for the López de Bertodano Formation51,55. Visual inspection shows shell preservation is excellent, with many shells still displaying pristine mother-of-pearl sheen, indicating the presence of primary aragonite. This was confirmed by conducting X-ray diffraction of powdered material from seven shells representing each of the six species in the full sample set, taken from across the full range of the stratigraphic section. These were shell numbers C1556E (C. antarctica), L776E (L. larseni), D1110I (D. dryganskalia), E1109E (E. regina), N1614C (N. nordenskjoldia), L1161A2 (L. larseni) and C1430 (C. ellioti). In all cases, the mineralogy was completely aragonitic and there was no indication of partial conversion to calcite.

Thick or thin sections were made for nine shells (sample numbers C757C2, L757B, L1480B, L1480E, C1467A, L1529A, C915A2, L1430A1, and L1609). Primary textures, including excellent preservation of growth band structures, were observed in shell carbonate from these sections (Supplementary Fig. 9). These sections were also viewed under cathodoluminescence, which can detect elevated amounts of Mn2+ that may become incorporated in diagenetic carbonate due to post-depositional water-rock interaction. Shells were non-luminescent, except in rare instances where fractures cut across the shell. This observation is particularly notable because it indicates that pore- fluids in the host sediment were reducing and enriched in Mn2+, yet they did not affect bulk shell chemistry. Therefore, alteration is likely only a concern in areas where the shells are fractured. Fractured areas were avoided in sample extraction.

Shell preservation was also assessed through trace element composition. Acceptable levels of Fe and Mn for a well-preserved shell have been determined differently by different authors17,56. Conservative trace element thresholds of 500 p.p.m. for Fe and 200 p.p.m. for Mn were set, and samples with concentrations above these were deemed ‘potentially altered’ and were excluded from further analysis. Of 339 total bulk samples, 42 individual bulk samples (representing 32 unique bivalve shells) had elevated trace element concentrations, indicating overall excellent preservation of samples at Seymour Island.

Shells selected for clumped isotope analysis all had Fe and Mn concentrations below the defined thresholds in the umbo and ventral margin areas. In three cases (N1416B, C772A and C1109A), the trace element composition was only measured at one of the positions on the shell but both positions were sampled for clumped isotopes. In three other cases (C757C-umbo, C915A1-umbo and C915A1-back), the trace element measurements came from the opposite half of the bisected shell, in a mirror image position to where the shell was sampled for clumped isotope analysis. For two shells (L757 and L1430A1), the trace element samples were taken from the outside of the shell, whereas the clumped isotope samples were taken from the inside after the shell was cut in half.

Stable and clumped isotope analysis

Stable carbon and oxygen isotope data were generated for all bulk samples from the full sample set using an automated Kiel device coupled to a Finnigan MAT 251 isotope ratio mass spectrometer in the University of Michigan Stable Isotope Laboratory (Director: K.C. Lohmann). Before reaction with phosphoric acid in the Kiel device, individual powdered carbonate samples were first roasted in vacuo at 200 °C for 1 h to remove volatile contaminants. Precision of the data was maintained at better than 0.1‰ by calibration to carbonate standards (both internal lab standards and NBS-19), measured daily. All carbonate stable isotope compositions are reported relative to Vienna Peedee Belemnite (VPDB). Stable isotope data for bulk samples can be found in Supplementary Data 1 and are shown in Supplementary Figs 5–7.

Samples in the ‘clumped isotope sample set’ were measured for clumped isotopic composition at the University of Michigan Stable Isotope Laboratory. Carbonate powders were prepared on an offline sample preparation device described by Defliese et al.57. Samples measured in 2015 were cleaned using the ‘WarmPPQ’ configuration, whereas samples measured in 2014 used the ‘ColdPPQ’ configuration, and were corrected accordingly for fractionations introduced by the Porapak trap58.

Standard gases heated to 1,000 °C and equilibrated with water at 25 °C were measured alongside samples and were used to define the absolute reference frame59. An acid fractionation factor of 0.067‰ was used, chosen to reflect the 75 °C reaction temperature57. Two carbonate standard materials, Carrara Marble and Joulter’s Cay Ooids, were measured alongside samples. Accepted Δ 47 values for these two standards are 0.414‰ and 0.704‰, respectively57,58. Over the three measurement sessions presented here, the mean values for Carrara and Ooids were 0.423±0.005‰ (1 s.e., n=15) and 0.712±0.007‰ (1 s.e., n=13). Raw voltages are converted into delta values following Huntington et al.60, including the same defined values for λ and for the isotopic composition of VPDB and VSMOW, and Δ 47 values were converted to temperature using equation 6 of Defliese et al.57.

The stable isotopic composition of the powders in the ‘clumped isotope sample set’ was measured simultaneously with the measurement of the clumped isotopic composition. Oxygen isotope data were corrected for fractionation during acid digestion using an acid fractionation factor of 1.00836 for aragonitic samples, determined empirically using the aragonitic Ooids carbonate standard, separately calibrated to NBS 18 and NBS-19 (ref. 58). Oxygen isotopic composition of water was calculated using the aragonite-water fractionation of Kim et al.12, and is reported relative to VSMOW. Measured clumped isotopic and stable isotopic compositions of ‘clumped isotope sample set’ shells can be found in Supplementary Data 2, 3, 4, 5. Carbonate and gas standard data can be found in Supplementary Data 2, 3, 4.

Calculation of sample and horizon averages

Average δ13C and δ18O values are calculated as the mean of many (n=2–5) replicates, with error taken as 1 s.d. Average values for Δ 47 , temperature, and δ18O w are taken as the mean of n replicates with the error taken as 1 s.e. on the mean, as is typical in the clumped isotope community (Supplementary Data 5). δ18O w is first calculated individually for each replicate from each δ18O and Δ 47 -derived temperature pair, then mean δ18O w is taken as the mean of many replicates, with 1 s.e. error. In some cases, the calculated 1 s.e. is smaller than the long-term performance of our carbonate standards. We determine ‘external error’ on Δ 47 for a given sample to be the larger of (1) the calculated 1 s.e. of the mean of n replicates or (2) the standard deviation of all replicates of Carrara (0.019‰) divided by the square-root of n replicates (Supplementary Data 5). For two to five replicates, this gives a 1 s.e. of 0.013‰, 0.011‰, 0.009‰ and 0.008‰. The ‘external error’ on temperature is calculated as half of T(mean Δ 47 −extSE.)−T(mean Δ 47 +extSE), where T(x) is the Δ 47 −Temperature calibration function and ‘extSE’ is the external s.e. replacing the calculated s.e. Based on typical 1 s.e. values on δ18O w for a sample with a Δ 47 -error of 0.013‰, 0.011‰, 0.009‰ and 0.008‰, the ‘external error’ on δ18O w is assigned to be 0.92‰, 0.78‰, 0.64‰ and 0.57‰ for n=2–5 replicates, respectively. External error is used in all figures and calculations for Supplementary Data 6, 7, 8.

Most (11 of 18) horizons are represented by two samples taken from different positions on a single shell. Where more than one shell is combined (6 of 18 horizons), all samples are treated equally, regardless of position or shell. Only one horizon is represented by a single position on a single shell (E. regina at 66.9 Ma). Due to consistently observed patterns within species, we interpret differences between shell positions to reflect environmental aliasing, not scatter around a mean value in a noisy data set. Therefore, the error calculated on a horizon average does not represent certainty on the average value of a noisy data set, but instead represents the spread in data from different shell positions, possibly an indicator of seasonality, combined with some degree of measurement noise. We do not know which species or shell position best represents mean annual temperature, so we felt combining all points equally was the fairest way to treat the data.

Horizon means (Supplementary Data 6, 7, 8) are calculated using equation (1) as the inverse variance weighted mean of all samples from a given horizon, combining samples within the uncertainty of the stratigraphy (±3.5 m)50.

Error on the horizon mean is calculated using equation (2) as error on an inverse variance weighted mean, to maintain consistency.

Inverse variance weighted mean was selected to account for the variable reproducibility of samples and not overweight the mean towards a poorly reproducing sample (for example, Sample C1109A at 67.7 Ma). Supplementary Fig. 2 compares multiple methods for calculating horizon mean (for example, inverse variance weighted mean, traditional mean, average of shell averages). The use of traditional mean versus inverse variance weighted mean does not change the long-term pattern, interpretations, or conclusions. The largest difference is seen around 67.7 Ma where there are two poorly reproducing samples.

Taking the average of positions within each shell before combining shells is mathematically equivalent to averaging all points if using the inverse variance weighted mean. Averaging within-shell first only affects four horizons mean minimally at that, if using the traditional mean. We chose to combine adjacent samples based on stratigraphic position, but a composite record could also have been created by combining samples into equal-time bins. Equal-time bins combine samples from before and after the warming spike, resulting in a ‘smoothed’ temperature record. If the bin size is small enough (50 kyr), two warming pulses are still observed.

Age model construction

To construct an age model using the recently published magnetostratigraphic data from Seymour Island10, the stratigraphic heights of chron reversals were converted to the Zinsmeister stratigraphic framework50 based on the established position of the KPg boundary in each (1,059 m for this study, 865 for Tobin et al.10) and assuming equal unit thickness. The stratigraphic height of a point in the Zinsmeister framework is therefore calculated as the height in the Tobin framework plus the difference in the KPg boundary positions (1,059–865=194 m). This assumption is supported by the position of the Unit 8–Unit 9 boundary, which is at 652 m in the Tobin framework10. The age model conversion would predict that this Unit boundary should be found at 846 m in the Zinsmeister framework, and it is described as ∼850 m by Zinsmeister50. Based on this close agreement, when comparing extinction events observed in the Tobin framework with climatic events from this study, any events occurring between the base of Unit 8 and the KPB should be in agreement to within ∼4 m.

The ages of the relevant chron reversals were taken from the most recent Geological Time Scale (2012)61, and the age the KPg boundary was taken as 66.043±0.043 Ma (ref. 34). The position of chron reversals was taken as the average of the stratigraphic heights of the first sample listed above and below the reversal position10. This is different from the age model used by Tobin et al.10, which used the 2004 edition of the Geologic Time Scale with the age of the KPg boundary at 65.5 Ma. To plot the δ18O-derived temperature record from Tobin et al.10 (Fig. 2), ages for chron reversals were updated to the 2012 Geological Time Scale61.

Sample ages were linearly interpolated between age model tie points, listed in Supplementary Table 1. There is inherent uncertainty in magnetostratigraphy-based age models due to this interpolation and uncertainties in the construction of the Geological Time Scale61. Although not explicitly accounted for, this uncertainty should be considered when comparing events in our age model with absolute U-Pb or 40Ar/39Ar ages for the Deccan Traps7,33.

Strontium isotope analysis

Ten shells were analysed for strontium isotopes for use in corroborating the age model construction. Analyses were performed in the Biogeochemistry and Environmental Isotope Geochemistry Laboratory at the University of Michigan (Director: J.D. Blum). Powdered samples of bivalve aragonite were dissolved in nitric acid and strontium was subsequently separated using ion-exchange chromatography with Eichrom Sr-specific resin. The effluent was dried and loaded onto a tungsten filament for isotope ratio analysis using a Finnigan MAT 262 thermal ionization mass spectrometer. The 86Sr/88Sr ratio was normalized to 0.1194 and measured 87Sr/86Sr ratios were corrected for the difference between the long-term laboratory average of NIST 987 (0.710252±0.000026 (2σ)) and the accepted value of 0.710248. Analytical uncertainty was calculated using the standard deviation of NIST 987 over a 12-month period (2σ=0.000026).

One sample (L1480B) produced a strontium isotopic composition that is not possible to match to the LOWESS curve at any point between 65 and 69 Ma. L1480B also displays anomalously high Sr/Ca values (10–20 mmol mol−1), suggesting that it may have been diagenetically altered. Another sample (L1480E) from the same locality produced a reasonable strontium isotope value, so L1480B was ignored. L1480B was also avoided for clumped isotope analysis. Strontium isotope data measured in this study can be found in Supplementary Table 2.

Corroboration of age model with Sr-isotope data

Strontium isotopic data was not used in the creation of the age model described above. These data can, therefore, be used to assess the validity of the age model constructed based on magnetostratigraphy and stratigraphic correlation. In addition to new strontium isotope measurements made on bivalve samples from this study (Supplementary Table 2), strontium isotope data from other Maastrichtian bivalves from the Zinsmeister collection have been previously published (Supplementary Table 3)13. Because they come from the same collection, their stratigraphic positions can be directly compared. Ages were calculated for all published samples using the age model described above based on their stratigraphic position in the Zinsmeister stratigraphic framework50. Supplementary Fig. 10 shows a comparison of measured 87Sr/86Sr values and calculated ages with the most recent seawater strontium isotope curve (LOWESS 5)62. There is good agreement between the bivalve data and the seawater curve, lending strong credence to the magnetostratigraphic age model and the stratigraphic conversion between the Zinsmeister framework and the Tobin framework.

Stable isotope trends in the full sample set

In the paper we discuss systematic differences in δ13C, temperature, and δ18O w between the umbo and back position seen in samples in the clumped isotope sample set (Figs 1 and 2; Supplementary Figs 3,4 and 11). Supplementary Fig. 11 also shows the relationship between position and shell δ18O. Systematic position-specific behaviour in δ18O is complicated due to the competing influences of δ18O w and temperature, which vary together in a way to oppositely influence δ18O (Supplementary Fig. 8). Despite this, Cucullaea still consistently shows higher δ18O in the ventral margin position, and Lahillia generally shows higher δ18O in the umbo position.

Although we cannot look at trends in temperature and δ18O w in the full sample set, we can still look at variation in δ13C and δ18O. In the full sample set, Lahillia and Cucullaea maintain the position-specific relationships in δ13C and δ18O seen in the limited clumped isotope sample set (Supplementary Figs 5–7). Cucullaea has the highest δ13C and δ18O values in the ventral margin, whereas Lahillia has the highest δ13C and δ18O values in the umbo. The consistency in position-specific relationships in the larger data set suggests that the same position-specific behaviour in temperature and δ18O w is likely present in all samples, even if not measured for clumped isotopes.

Nordenskjoldia and Dozyia show similar position-specific relationships in δ13C and δ18O to Cucullaea and Lahillia, respectively, but of lower total magnitude. Eselaevitrigonia shows a much smaller variation in δ18O than Cucullaea or Lahillia and δ18O does not correlate to δ13C as in the other species. This lower variability in stable isotopes suggests consistent environmental and growth conditions between early and late life for this species. For this reason, Eselaevitrigonia was only sampled in the umbo position. This decision is supported by the small differences in temperature and δ18O w seen in Nordenskjolia, another species showing small variations in stable isotopes (Fig. 1).

Calculations of CO 2 emissions from Deccan Traps volcanism

Richards et al.6 tentatively place the KPB between the Bushe and Poladpur formations in the Deccan volcanic sequence due to evidence for a hiatus in volcanism and fracturing and faulting in the lower Bushe formation that did not extend into the higher Poladpur formation. This placement is consistent with existing dating of the Deccan Traps, which unfortunately has a sampling gap in the critical region preventing direct placement of the KPB7,33. Assuming this stratigraphic placement of the KPB, we determined an estimate of the minimum and maximum volume of lava erupted between the onset of Chron 29R Deccan volcanism and the KPB, which includes all formations between (and including) Jawhar and Bushe (Kalsubai and Lonavala subgroups).

Estimates of the volume of lava erupted in this interval range from 150 × 103 km3 (ref. 6) to 500 × 103 km3 (ref. 35). Using an emission rate of 14 Tg CO 2 km−3 lava34, a modern atmospheric mass of 5.15 × 1021 g and an atmospheric density of 28.966 g mol−1 air, these lava volumes equate to 270–900 p.p.m. CO 2 emitted to the atmosphere. Compare this to end-Maastrichtian atmospheric CO 2 levels, which are thought to be ∼360–380 p.p.m. CO 2 (ref. 36). The exact timing of lava extrusion during this interval is unknown, but a lack of red boles in the older portion of the section implies very few hiatuses7,63. The transient timing of marine carbonate dissolution suggests the main phase of degassing occurred beginning shortly after the onset of Chron 29R and ending within 200 kyr (ref. 37).

We also estimate CO 2 emissions during the interval between the KPB and the end of Chron 29R volcanism, which includes the Poladpur to Mahabaleshwar formations (Wai subgroup) based on our hypothetical placement of the KPB. The post-KPB portion of Deccan volcanism emitted another 461 × 103 km3 (ref. 6) to 500 × 103 km3 (ref. 35) of lava, equivalent to 830–900 p.p.m. CO 2 added to the atmosphere. The Wai subgroup eruptions contain more red boles7,63, suggesting hiatuses between eruptive events. Based on published dates7,33, the duration of the post-KPB eruptions is roughly twice the length of the pre-KPB portion (∼500 kyr versus 250 kyr), although this does not account for potential hiatuses within the eruptive sequence. Although the total volume of lava erupted (and therefore CO 2 emitted, based on our assumptions) was larger post-KPB, the climatological impacts might have been smaller due to higher background CO 2 levels and a longer duration of eruption. Additionally, the emission rate of 14 Tg CO 2 km−3 lava35 is an average value. Changes in the character of eruptions and the geochemical signature of the lava before and after the Bushe/Poladpur contact (the hypothesized KPB)6 could coincide with changes in the lava’s CO 2 content and therefore emission rate.

Data availability

The data shown and discussed in this paper is presented in full in the Supplementary Material.