Core location and palaeogeography

The Shuqualak–Evans core is from Shuqualak, Mississippi, USA (32°58′49′′N, 88°34′8′′W) and was sampled for TEX 86 from a depth of 9.45 m down to 251.46 m, spanning the Santonian/Campanian boundary interval through to the uppermost Maastrichtian.

Figure 1 shows a model of the palaeogeographic evolution of North America during the latest Cretaceous. In these reconstructions and others (for example, ref. 20) the Shuqualak–Evans borehole was situated on a broad shelf bordering the North Atlantic Ocean and Gulf of Mexico during the latest Cretaceous. Surface ocean circulation reconstructions (summarized in ref. 31) suggest that this location was likely not influenced by waters of the Western Interior Seaway and this is supported by calcareous nannofossil assemblage components and abundances, which suggest an open ocean water mass.

Biostratigraphy and age-model

The age-model for the Shuqualak–Evans core is based on integrated calcareous nannofossil41 and planktonic foraminifera42 biostratigraphic datums (Supplementary Fig. 1 and Supplementary Table 1), with calibrated absolute ages taken from Gradstein et al.42 For the uppermost Campanian through Maastrichtian (25.91–9.45 m), the age model is calculated as a linear function, with tie points at base Lithraphidites quadratus (16.76 m, 69.18 Ma42) and base Micula prinsii (12.80 m, 67.30 Ma42). This suggests a low sedimentation rate of ~2.1 m per Myr ago. We assign a minimum age of 66 Ma to the shallowest sample (9.45 m), as the nannofossils indicate a Cretaceous age yet the age-model predicts a Palaeocene age and we cannot further refine the age-model above 12.80 m, which appears to be topmost Maastrichtian. For most of the Campanian (245.36–30.48 m), the age-model is a linear function, with tie points at base Broinsonia parca subsp. constricta (239.27 m, 81.38 Ma42) and base Uniplanarius sissinghii (134.11 m, 77.61 Ma42). This indicates a high sedimentation rate of ~27.9 m per Myr. For the lowermost sample analysed for geochemistry, the age model is constrained by base Broinsonia parca subsp. parca (245.36 m, 81.43 Ma42) and the presence of Arkhangelskiella cymbiformis (252.83 m, assigned maximum age of 83.20 Ma42), indicating a low sedimentation rate of ~4.2 m per Myr ago. Note that the co-occurrence of Dicarinella asymetrica and A. cymbiformis (at 251.46 m) suggests that the age of the lowermost sample analysed for TEX 86 is around the Santonian/Campanian boundary interval. The TEX 86 and SST data are given in relation to the sample depths in Supplementary Fig. 2 and Supplementary Table 2.

GDGT extraction and analysis

Samples were solvent extracted using the technique previously published by Schouten et al.13,43 Approximately 6 g powdered sample was ultrasonically extracted using one time methanol, three times dichloromethane (DCM)/methanol (1:1, v/v) and three times DCM. All extracts were combined and dried under a continuous N 2 flow at 40 °C. Any water remaining in the samples was removed by passing the extracts (dissolved in DCM/methanol (3:1, v/v)) over a column containing anhydrous Na 2 SO 4 . Extracts were split into polar and apolar fractions by column chromatography, using hexane/DCM (9:1, v/v) and DCM/methanol (1:1, v/v) sequentially as the eluents and Al 2 O 3 as the stationary phase. The polar extract containing the targeted GDGTs was dissolved in hexane/propanol (99:1, v/v) and then filtered through a PTFE (polytetrafluoroethylene) 0.45 μm filter. After drying down, the samples were redissolved in a certain volume of hexane/propanol (99:1, v/v), which depends on the weight of each polar fraction. All 48 samples were analysed in the Department of Earth Sciences at UCL on an Agilent 1200 series HPLC attached to a G6130A single-quadrupole mass spectrometer. The analytical protocol followed is as described in Schouten et al.43. The abundance of both isoprenoid and branched GDGTs was measured in selective ion monitoring mode. Ion peaks of the respective GDGTs were integrated to determine the relative abundance of each molecule in the sample. These abundances were then used to determine the TEX 86, TEX 86 H (GDGT-index 2), and TEX 86 L (GDGT-index 1) indices13,14. These indices are defined as follows:

where Cren' represents the crenarchaeol regioisomer. Our calculated TEX 86, TEX 86 H (GDGT-index 2) and TEX 86 L (GDGT-index 1) values are presented in Supplementary Table 2 and Supplementary Fig. 2.

SST calculations from GDGT abundances

To calculate SSTs, we used the following equations from Kim et al.14 which are based upon a comprehensive modern core-top dataset:

At temperatures >15 °C expected during greenhouse periods in Earth history, it has been recommended14 that the TEX 86 H index should be used, as both indices should yield the same estimate of SST according to the modern core-top calibration dataset, but TEX 86 H has an associated calibration error that is significantly lower (±2.5 °C) compared with that of TEX 86 L (±4 °C). However, application of TEX 86 L and TEX 86 H to datasets from Early Cenozoic sediments21 has revealed that they do not always yield the same estimates of SST, with TEX 86 L generally yielding lower temperatures. At high-latitude Palaeocene–Eocene sites in New Zealand, Hollis et al.21 found that TEX 86 L yielded SST estimates more comparable to inorganic proxies (δ18O, Mg/Ca) than TEX 86 H, but suggested that, at low latitudes and high TEX 86 values, TEX 86 H might be as appropriate. In their analysis of the ability of the different GDGT-based proxies to replicate inorganic SST estimates, Hollis et al.21 noted that, at high TEX 86 values above 0.70, the over-estimation of SST is less than 5 °C using TEX 86 H. Furthermore, at TEX 86 values above 0.75, they suggested that TEX 86 L underestimates SST. The oxygen-isotopic compositions of well-preserved planktonic foraminifera of Cenomanian–Santonian age from Demerara Rise (~5 °N palaeolatitude) yield SSTs typically in the range of 35 to >37 °C24, comparable to SSTs derived by TEX 86 H of 35 to 37 °C from the same site (recalculated by us from the published TEX 86 data). Therefore, through the consideration of previous suggestions and multiproxy records of low-latitude mid-Cretaceous SSTs, it would initially appear that the TEX 86 H calibration represents the best estimate of SSTs for the Shuqualak–Evans borehole. Furthermore, in order to be able to compare our new data with previously published Late Cretaceous TEX 86 data2,24, from which TEX 86 L values are not yet available, we have had to use the TEX 86 H index in Fig. 3. Nonetheless, in Supplementary Figs 2 and 3 we show both the TEX 86 H and TEX 86 L data (where possible) and the calculated SSTs to illustrate the potential range of values. We also present in Supplementary Fig. 2 SSTs calculated using the recently developed BAYSPAR approach19. The BAYSPAR model considers how the relationship between TEX 86 and temperature varies spatially and considers uncertainties in the modern SST-TEX 86 relationship. Critically, the stratigraphic trends for Shuqualak–Evans are near identical, irrespective of which proxy is used and, thus, our conclusions regarding the temporal evolution of the direction of Late Cretaceous climatic and latitudinal gradient change remains valid, even if absolute values are harder to constrain.

Recent work, based on an analysis of modern water-column GDGT abundance profiles, the core-top calibration dataset and a compiled Paleogene dataset, suggests that TEX 86 H may be less appropriate than TEX 86 L for use at sites where the water depth was approximately shallower than 1000, m (such as Shuqualak). This is likely due to variations in the export dynamics of individual GDGT compounds with depth, and an apparent temperature/water depth bias in the core-top calibration dataset18. In both the modern core-top calibration dataset and the Paleogene dataset, sites deposited in <1,000 m of water exhibit low GDGT-2/GDGT-3 ([2]/[3]) ratios and high offsets between SSTs calculated by TEX 86 H and TEX 86 L (ΔH-L). We have been able to obtain the raw GDGT data for the sites on Demerara Rise, which were thought to have been deposited at water depths of <1500, m44. The GDGTs from these sites exhibit low [2]/[3] ratios and high ΔH-L, which Taylor et al.18 suggest is characteristic of water depths <1,000 m. We therefore contend that it is appropriate in Supplementary Fig. 3 to compare Demerara Rise and Shuqualak–Evans using either TEX 86 H or TEX 86 L (as the water depth of all sites was likely about, or shallower than, 1000, m). The application of the TEX 86 L calibration to our data from Shuqualak–Evans suggests SSTs of ~28 °C for the earliest Campanian and ~20 °C for the Campanian–Maastrichtian transition, which are some ~7–8 °C lower than the estimates of the TEX 86 H model, and much closer to modern SSTs at comparable latitudes, which perhaps is surprising given that Late Cretaceous pCO 2 levels are thought to have been higher than present (~600 to 800 p.p.m. versus 280 p.p.m. for the preindustrial modern)40. Furthermore, the use of TEX 86 L suggests almost no temperature gradient between 5° and 60° absolute palaeolatitude during the Turonian–Santonian, which also seems unlikely. The issue of how best to calculate temperatures from GDGT data is ongoing and it may be that a calibration based on suspended particulate organic matter may overcome some of the issues described above16,18.

Repeated analysis of an in-house standard and selected samples suggest that analytical reproducibility of TEX 86 is better than ±0.009, in line with previous studies45 that suggest an analytical error for TEX 86 index of ±0.01. In our data, we estimate that the error on SST estimates associated with analytical precision is <±0.4 °C for the TEX 86 H calibration, and <±0.6 °C using TEX 86 L. Analytical error is far less than the standard error associated with the core-top calibrations, which for TEX 86 H is ±2.5 °C, and for the TEX 86 L calibration is ±4.0 °C21.

BIT and MI indices

The GDGTs analysed for the TEX 86 palaeotemperature proxy are mainly produced by marine thaumarchaeota. However, the same GDGTs are also produced by terrestrial soil organisms and methanotrophs.

GDGTs of terrestrial origin can be washed into the marine realm by rivers, potentially biasing SST reconstructions. Apart from isoprenoid GDGTs, which are used for the TEX 86 techniques, terrestrial organisms also produce branched GDGTs46. Branched GDGTs are typical of terrestrial organisms, but they do not occur among marine thaumarchaeota46. Thus, the branched GDGTs are used to quantify the terrestrial GDGT contamination in marine sediment samples by calculating the Branched and Isoprenoid Tetraether (BIT) index46. The BIT index is based on the ratio of branched GDGTs to the isoprenoid GDGT crenarchaeol46. In our study, the measurement of these branched GDGTs was included in the analytical protocol. There is no significant relationship between BIT and TEX 86 in our data (Supplementary Fig. 4) and the BIT index is between 0.05 and 0.15 (Supplementary Fig. 2). This is well below the recommended threshold of 0.2, above which soil microbial contamination may be problematic47. Thus the SST estimates made in this study are probably not altered by an influence of terrestrial GDGTs.

GDGTs produced by methanotrophs within the sediments can distort the TEX 86 signal and lead to erroneous estimates of palaeotemperature. The degree of influence of methanotropic archaea can be estimated using the Methane Index (MI)17. Normal marine sediments have values <0.3, whereas sediments influenced by high rates of methane production have values >0.5. The interval from 0.3 to 0.5 marks the transition between the two environments17. The MI values from Shuqualak vary from ~0.1 to 0.25 (Supplementary Fig. 2), suggesting normal marine conditions and a lack of GDGT production by methanotrophs. We therefore conclude that methanotrophic production of GDGTs has not impacted adversely on our TEX 86 records.

Palaeotemperatures from foraminiferal oxygen isotopes

In Fig. 3 and Supplementary Fig. 3, we show estimates of palaeotemperature based upon the oxygen-isotopic (δ18O) composition of a global compilation of benthic foraminiferal data4, selected low-latitude planktonic foraminiferal datasets12,25 and mixed-layer-dwelling planktonic foraminifera from southern high-latitude sites (DSDP Site 511 and ODP Site 690)6,7. Calcareous dinoflagellate cysts from ODP Site 690 have similar d18O values to planktic foraminifera48 but are not included as a suitable temperature calibration for Cretaceous dinoflagellates is not available. Note for the Turonian-age samples from Tanzania, we have recalculated the SST range using the typical range of planktonic δ18O values measured (−4 to −5 ‰)25. In order to apply a consistent approach to calculating temperatures from δ18O of foraminiferal calcite, we have recalculated temperatures from the original published oxygen-isotopic datasets. We used equation (6) to calculate temperatures49:

where T=temperature (°C), δ cc =δ18O of foraminiferal calcite (‰, VPDB) and δ w =δ18O of ambient seawater.

We used a δ w value of −1.27, which includes a correction of −0.27 for the conversion of Vienna Standard Mean Ocean Water to Vienna PeeDee Belemnite and an ice volume value of −1‰. For calculation of SSTs from planktonic foraminifera, we applied an additional latitudinal salinity correction50 to δ w , using the palaeolatitudes of DSDP Site 511 and ODP Site 690 in the Late Cretaceous (58°S and 67°S, respectively) and equation (7):

where x=absolute palaeolatitude between 0° and 70°. Palaeolatitudinal positions for each site were using the ODSN Plate reconstruction software ( http://www.odsn.de/odsn/services/paleomap/paleomap.html).