Cretaceous sea-surface temperature evolution: Constraints from TEX 86 and planktonic foraminiferal oxygen isotopes

many similarities with trends present in individual records of Cretaceous climate change. For example, both SST proxies and benthic foraminiferal δ 18 O records indicate maximum warmth in the Cenomanian – Turonian interval. Our reconstruction of the evolution of latitudinal temperature gradients (low,< ±30°, minus higher,> ±48°, palaeolatitudes) reveals temporal changes. In the Valanginian – Aptian, the low-to-higher mid-latitudinal temperature gradient was weak (decreasing from ~10 – 17 °C in the Valanginian, to ~3 – 5 °C in the Aptian, based on TEX 86 -SSTs). In the Cenomanian – Santonian, reconstructed latitudinal temperature contrasts are also small relative to modern (< 14 °C, based on low-latitude TEX 86 and δ 18 O pl SSTs minus higher latitude δ 18 O pl SSTs, compared with ~20 °C for the modern). In the mid-Campanian to end-Maastrichtian, latitudinal temperature gradients strengthened (~19 – 21 °C, based on low-latitude TEX 86 and δ 18 O pl SSTs minus higher latitude δ 18 O pl SSTs), with cooling occurring at low-, middle-and higher palaeolatitude sites, implying global surface-ocean cooling and/or changes in ocean heat transport in the Late Cretaceous. These reconstructed long-term trends are resilient, regardless of the choice of proxy (TEX 86 or δ 18 O pl ) or calibration. This new Cretaceous SST synthesis provides an up-to-date target for modelling studies investigating the mechanics of extreme climates.


A B S T R A C T
It is well established that greenhouse conditions prevailed during the Cretaceous Period (~145-66 Ma).Determining the exact nature of the greenhouse-gas forcing, climatic warming and climate sensitivity remains, however, an active topic of research.Quantitative and qualitative geochemical and palaeontological proxies provide valuable observational constraints on Cretaceous climate.In particular, reconstructions of Cretaceous sea-surface temperatures (SSTs) have been revolutionised firstly by the recognition that clay-rich sequences can host exceptionally preserved planktonic foraminifera allowing for reliable oxygen-isotope analyses and, secondly by the development of the organic palaeothermometer TEX 86, based on the distribution of marine archaeal membrane lipids.Here we provide a new compilation and synthesis of available planktonic foraminiferal δ 18 O (δ 18 O pl ) and TEX 86 -SST proxy data for almost the entire Cretaceous Period.The compilation uses SSTs recalculated from published raw data, allowing examination of the sensitivity of each proxy to the calculation method (e.g., choice of calibration) and places all data on a common timescale.Overall, the compilation shows 1. Introduction
Understanding global temperatures in the Cretaceous, both in terms of absolute magnitude and spatial and temporal variation, is vital to elucidating the exact nature of greenhouse-gas forcing, climate sensitivity and ocean circulation at this time (e.g., Friedrich et al., 2012;PALEOSENS Project Members, 2012;Royer et al., 2012;Wang et al., 2014).Terrestrial MAT estimates from, for example, plant macrofossils, palynomorphs and δ 18 O of vertebrate tooth enamel (e.g., Askin, 1989;Wolf, 1990;Herman and Spicer, 1996;Amiot et al., 2004;Poole et al., 2005;Upchurch et al., 2015), provide important constraints on temperatures of continental interiors and high palaeolatitudes and in some cases information on seasonality (e.g., Herman and Spicer, 1996;Poole et al., 2005).In the marine realm, Cretaceous ocean sediments are host to valuable palaeotemperature proxy archives, offering considerable spatial and temporal coverage.Recent δ 18 O b compilations (e.g.Cramer et al., 2009;Friedrich et al., 2012) have provided important insights into the evolution of intermediate to deep-ocean temperature conditions, often assumed to reflect variations in global temperature, during the middle to Late Cretaceous.However, intermediate to deep-water temperature trends can differ between ocean basins depending on the influences of source waters with different temperatures, such that the relationship between surface and benthic temperatures depends on ocean circulation.Upper mixed-layer temperatures can provide more precise constraints on climate because of the direct contact between the atmosphere and the surface ocean.Thus, estimating the global evolution of Cretaceous sea surface temperature remains a prime objective of palaeoclimate research.
Both δ 18 O pl palaeothermometry and the TEX 86 palaeothermometer can provide estimates of SST, but each technique is subject to several proxy-specific caveats (see Sections 3.2 and 3.4 for detailed discussion).For example, δ 18 O pl values can be compromised by preservation and/or diagenetic alteration (e.g., Schrag et al., 1995;Pearson et al., 2001), the carbonate ion effect (e.g., Spero et al., 1997;Zeebe, 2001), and uncertainties related to the δ 18 O of seawater in which the foraminifer calcified (e.g., Poulsen et al., 1999;Zhou et al., 2008).Indeed, earlier research using δ 18 O pl values to reconstruct Cretaceous surface-ocean conditions suggested that tropical SSTs were similar or cooler than modern tropical ocean conditions (the cool tropics paradox; Sellwood et al., 1994;D'Hondt and Arthur, 1996) even though δ 18 O pl values from polar regions indicated sub-tropical-like SSTs.It has since been demonstrated that these inferred cool tropical temperatures reflect biased δ 18 O pl values derived from diagenetically altered, cool-biased planktonic foraminifera (Schrag et al., 1995;Pearson et al., 2001), indicating the importance of selecting only well-preserved foraminifera, i.e. "glassy" foraminifera found in clay-rich sediments.This realization essentially eliminated all published δ 18 O pl data from carbonate-rich deep ocean settings as reliable indicators for SST.
For the TEX 86 proxy, the exact mechanism(s) that relates sedimentary GDGT distributions to those produced by Thaumarchaeota in surface waters and therefore to SSTs are not fully understood (Pearson and Ingalls, 2013;Schouten et al., 2013b;Taylor et al., 2013;Tierney, 2014).For example, recent work has suggested that non-temperature factors such as oxygen concentration (Qin et al., 2015), growth phase (Elling et al., 2014), ammonium oxidation rate (Hurley et al., 2016) and other physiological, environmental, and ecological factors (Elling et al., 2015) may play an important role in governing GDGT distributions.Moreover, other factors will govern how the signal generated in surface waters is preferentially exported to sediments, with some work suggesting that TEX 86 ratios are in fact an integrated signal from the surface and shallow-subsurface (0 to ~250 m; Wakeham et al., 2003;Wuchter et al., 2005;Taylor et al., 2013;Hernández-Sánchez et al., 2014).Recent work has suggested that the TEX 86 signal is predominantly exported from deeper waters (Ho and Laepple, 2016).However, this assumption is likely flawed because it does not recognize modern mechanistic studies on TEX 86 -SST sensitivity and fails to predict SSTs in shallow settings (Tierney et al., 2017).
For both δ 18 O pl and TEX 86 palaeothermometry, a number of different calibrations exist for converting proxy data values to estimates of SST (e.g., Erez and Luz, 1983;Bemis et al., 1998).Consequently not all published absolute estimates of Cretaceous SSTs for a given proxy are directly comparable.Furthermore, since publication of earlier Cretaceous TEX 86 data, constraints on the application of the proxy have appeared, e.g., the Branched and Isoprenoid Tetraether (BIT) index (Hopmans et al., 2004).Hence, data rejection criteria have changed over time (see later discussion, Section 2.1.3).Aside from temperature proxy developments, the most recent Cretaceous time scale has undergone revisions since the publication of several datasets (Gradstein et al., 2004;Gradstein et al., 2012).Thus, given the current body of Cretaceous δ 18 O pl and TEX 86 data, it is timely to re-evaluate critically all currently available Cretaceous SST data derived from both of these techniques, assessing their quality as well as their ages.Here we present a global SST compilation for the Cretaceous Period in order to provide new constraints on long-term SST evolution and uncertainty, and to provide a target for modelling studies.

Data compilation
Original published marine palaeotemperature proxy data, raw GDGT data and planktonic δ 18 O data, were collated from all available locations where these determinations are accepted as primarily reflecting Cretaceous SSTs, as interpreted by the original authors and in subsequent publications.We then apply additional tests to ensure that all data have had the same scrutiny and to incorporate new understanding on proxy biases.While δ 18 O pl data are reported relative to an international standard, Vienna Pee Dee Belemnite (VPDB), indicating that it is valid to compare data produced in different laboratories, no such standard currently exists for the TEX 86 proxy.The most recent TEX 86 inter-laboratory comparison study (Schouten et al., 2013a) of variations in measured TEX 86 values across different laboratories indicated a range of 0.023-0.053(equivalent to 1.5-3.5 °C in the TEX 86 -SST calibration of Schouten et al., 2002), which should be considered when comparing different records generated in different laboratories.However, this variation is typically less than the errors in the temperature calibrations (2.5 to 4.0 °C; Kim et al., 2010) and additionally, ~55% of the GDGT data presented here were generated in one laboratory, namely at the Royal Netherlands Institute for Sea Research (NIOZ; Table 1).For comparison intra-laboratory precision (repeatability) is typically < 1 °C for TEX 86 (Schouten et al., 2013a) and < 0.08‰ for δ 18 O pl , equivalent to < 0.4 °C (Ravelo and Hillaire-Marcel, 2007).

Compilation of raw GDGT data
Fractional abundances of isoprenoid GDGTs (see structures given in Supplementary Fig. 1) were compiled for locations where published TEX 86 palaeotemperature proxy data have been generated from Cretaceous sediments.Where possible, we compute and report the fractional abundance of all individual GDGTs.For some datasets, original raw GDGT data were unavailable, and for these locations we rely solely on published TEX 86 values.In addition, where available, Branched and Isoprenoid Tetraether (BIT) indices, a proxy for input of soil-derived organic matter (Hopmans et al., 2004;Weijers et al., 2006), were either collected or determined (see below).Locations of Cretaceous sediments for which raw GDGT data/TEX 86 indices were obtained are shown in Fig. 1a; data references are given in Table 1.

GDGT-based SST indices: TEX 86 , TEX 86 H , and BAYSPAR
We calculate TEX 86 values using the original definition of Schouten et al. (2002): where GDGT-1, GDGT-2, GDGT-3 and Cren′ refer to the structures shown in Supplementary Fig. 1.Several global core-top based calibrations exist for converting TEX 86 estimates to SSTs.Originally the TEX 86 -SST relationship was described using a linear calibration (Schouten et al., 2002;Kim et al., 2008).Subsequently, a 1/TEX 86 expression (Liu et al., 2009) and then a logarithmic relationship (Kim et al., 2010) were applied to better fit the TEX 86 -SST relationship.Kim et al. (2010) provided two logarithmic calibrations, TEX 86 H and TEX 86 L , with the former recommended for reconstructing SSTs > 15 °C, such as those characteristic of the Cretaceous Period.More recently, a Bayesian model approach has been developed to predict SSTs from TEX 86 values, BAYSPAR (Tierney and Tingley, 2014;Tierney and Tingley, 2015).Motivation for this model arose from observations that the TEX 86temperature relationship appears to vary for different ocean regions and environments (e.g., Trommer et al., 2009;Ho et al., 2014) and that previous calibration models feature structured residuals indicating strong spatially varying trends in the response of TEX 86 to temperature (Tierney and Tingley, 2014).BAYSPAR model SST predictions are derived using an online graphical use interface (GUI; www.whoi.edu/bayspar) or directly via the corresponding MatLab code (Tierney, 2014;Tierney and Tingley, 2015).For pre-Quaternary applications, the BAYSPAR model searches modern core-top data (including data from the Red Sea) for TEX 86 values that are similar to the measured TEX 86 values in a given dataset and derives linear-regression parameters from these modern 'analogues'.This analogue approach is not geographically dependent (unlike the Quaternary BAYSPAR approach) and relies only on the assumption that it is reasonable to compare modern and ancient TEX 86 values, which is little different to assumptions incorporated in traditional TEX 86 -temperature calibrations (Tierney and Tingley, 2014).Another feature of the BAYSPAR approach versus traditional TEX 86 -temperature regression models is that the BAYSPAR model fully propagates uncertainties in the core-top data into resulting temperature predictions (Tierney and Tingley, 2014).In addition to TEX 86 -SST calibrations, there are now TEX 86 calibrations for deriving shallow subsurface temperatures, 0-200 m water depth, in certain settings (e.g., Kim et al., 2015;Tierney and Tingley, 2015).
Regardless of the choice of TEX 86 -SST calibration, application of the TEX 86 proxy in the Cretaceous Period, where TEX 86 values are frequently high (> 0.8), requires extrapolation of TEX 86 -SST calibrations above the upper limit of the modern range reflected in the core-top datasets, ~0.72 (excluding data from the Red Sea; Kim et al., 2008;Kim et al., 2010).Thaumarchaeota mesocosm studies at high temperature, 30-46 °C (e.g., Wuchter et al., 2004;Schouten et al., 2007;Pitcher et al., 2010) provide support for extrapolation of temperatures above the modern calibration datasets, but do not yield an alternative calibration due to an unusually low abundance of the crenarchaeol regioisomer.However, Pitcher et al. (2010) record a TEX 86 value of 0.99 at an incubation temperature of 46 °C for a Thaumarchaeote suggesting an upper limit for the maximum SST estimate the TEX 86 proxy can yield.
Depending on the TEX 86 -SST relationship, different core-top-derived calibrations produce a different maximum SST estimate; for the 1/ TEX 86 calibration (Liu et al., 2009), a linear calibration (Kim et al., 2008) and the TEX 86 H calibration (Kim et al., 2010) the maximum SSTs that can be computed, i.e. when TEX 86 = 1, are 34.1 °C, 45.4 °C and 38.6 °C, respectively.For the BAYSPAR model, the nature of the regression is dependent on the core-top analogues selected, as such the corresponding SST estimate when TEX 86 = 1 is not fixed.The maximum computable TEX 86 -SST estimate derived using the 1/TEX 86 calibration of Liu et al. (2009) is significantly lower than for the other calibrations discussed here and is least consistent with other evidence of greenhouse conditions where maximum SSTs may well have been higher than this limit.As such, we herein convert TEX 86 values to SSTs using three calibrations: (1) the TEX 86 H calibration and (2) a new TEX 86-Linear calibration, both based on the global core-top dataset (see below; Kim et al., 2010), and (3) a BAYSPAR model approach (Tierney and Tingley, 2014;Tierney and Tingley, 2015).TEX 86 H -derived SSTs were computed using the temperature equations of Kim et al. (2010).Similar to the original TEX 86 relationship (Schouten et al., 2002), TEX 86 H is calculated from the fractional abundances of GDGT-1, GDGT-2, GDGT-3 and the regio-isomer of crenarchaeol, Cren′, and is defined as: TEX 86 H is correlated to SST using a global core-top calibration (Kim et al., 2010) A similarly derived high temperature TEX 86 calibration for the Cretaceous was previously presented in Schouten et al. (2003).This newly modified TEX 86 linear calibration excludes all data from the Red Sea and also all data for which satellite-derived mean annual SSTs are < 15 °C.The maximum SST that can be derived using this calibration (when TEX 86 = 1) is 42.7 °C.
We also compute Cretaceous SSTs from GDGT data using BAYSPAR (Tierney and Tingley, 2014;Tierney and Tingley, 2015).Here, we apply the default settings of the BAYSPAR "Deep-Time" model (Tierney and Tingley, 2014;Tierney and Tingley, 2015), inputting all Cretaceous TEX 86 data into the BAYSPAR model as one whole dataset.This approach allows us to generate a single, BAYesian-derived Global Regression (BAY GlobR ) in order to predict Cretaceous SSTs, BAY GlobR -SSTs.This single 'global' approach maintains the original inter-site distribution of the TEX 86 data, avoiding potentially artificial inter-site relative shifts introduced by generating a specific regression for each dataset.
We note that we do not apply the TEX 86 L calibration proposed by Kim et al. (2010).Designed for reconstructing SSTs across all temperature ranges (− 3-30 °C; Kim et al., 2010), the TEX 86 L calibration differs from other TEX 86 calibrations since it only employs GDGT-1, GDGT-2 and GDGT-3, removes GDGT-3 from the numerator and excludes the crenarchaeol regio-isomer entirely; moreover, it is not mathematically related to ring number, thereby removing the inferred physiological foundation of GDGT-based temperature proxies (Schouten et al., 2002).Application of TEX 86 L alongside other TEX 86 calibrations, namely TEX 86 H , in the Cretaceous and Palaeogene has highlighted significant offsets between the two calibrations in a range of different settings (Taylor et al., 2013;Inglis et al., 2015), including the western proto-North Atlantic (Early Cretaceous; Littler et al., 2014), the western North Atlantic Shelf (Late Cretaceous; Linnert et al., 2014), the New Jersey Coastal Plain (Palaeocene; Taylor et al., 2013) and the southwest Pacific (middle Palaeocene to middle Eocene; Hollis et al., 2012;Bijl et al., 2013).Recent work has illustrated that it is the TEX 86 L calibration that is particularly sensitive to GDGT export depth issues, due to the vertical variation of GDGT-2/GDGT-3 ratios in the water column (Hernández-Sánchez et al., 2014;Villanueva et al., 2014;Kim et al., 2015;Kim et al., 2016), rendering the calibration inaccurate (see Inglis et al., 2015).Indeed, the TEX 86 L calibration can underestimate SST in some settings, e.g., in tropical Eocene settings, yielding SSTs lower than modern despite greenhouse conditions (Sluijs et al., 2014;Inglis et al., 2015) and, even in subpolar and polar regions the TEX 86 L calibration does not perform substantially better than TEX 86 H (Ho et al., 2014).As such, we focus on other indices and calibrations, e.g., TEX 86 H , TEX 86-Linear and BAY GlobR , here.

GDGT preservation and secondary effects
GDGT distributions can be potentially influenced by secondary, non-thermal effects including oxic degradation, thermal alteration, terrestrial inputs and methanogenesis.Diagenetic alteration of GDGT distributions related to oxic degradation of GDGTs, either in the water column or in sediments, could potentially modify primary TEX 86 values.However, studies suggest that TEX 86 appears to neither be significantly nor systematically influenced by degree of oxidation/diagenesis (Schouten et al., 2004;Huguet et al., 2008;Huguet et al., 2009;Kim et al., 2009;Lengger et al., 2013).Van Helmond et al. (2015) found preferential preservation of brGDGTs over isoGDGTs relative to %Total Organic Carbon in Cretaceous sediments, similar to observations from organic matter-rich turbidites affected by post-depositional oxidation (e.g., Huguet et al., 2008;Lengger et al., 2013).This preferential preservation in Cretaceous sediments was also reflected by higher BIT values when TEX 86 values were lower (van Helmond et al., 2015), thus, in addition to identifying samples biased by inputs of terrestrially derived isoGDGTs, the BIT index (see later discussion, Section 2.1.4)may help identify samples biased by preferential preservation.Another concern is the potential for thermal alteration, particularly in older, i.e.Mesozoic sediments (e.g., Bottini et al., 2015) which can bias TEX 86 values to cooler SST estimates (Schouten et al., 2004).However, the degree of sediment thermal maturity can be relatively straightforwardly assessed by evaluating homohopanoid isomer compositions (e.g., van Duin et al., 1997;Schouten et al., 2004).There is also some evidence for small, consistent differences in TEX 86 from interbedded lithologies of similar Cretaceous age (Littler et al., 2014).These sedimentary TEX 86 offsets were interpreted to represent ecological differences in Thaumarchaeotal populations, spatial distribution or seasonality preferences, emphasizing the importance of careful sample selection (Littler et al., 2014).In addition to diagenesis and sampling biases, the primary marine TEX 86 signal can be modified by the introduction of additional GDGTs either from terrestrial (soil-derived) sources (Weijers et al., 2006) or synthesizing sedimentary methanogenic and methanotrophic Archaea (Blaga et al., 2009;Zhang et al., 2011;Sinninghe Damsté et al., 2012).Further details and methods to identify these additional GDGT sources are described below.
To examine the influence of soil-derived GDGT input on TEX 86 values, we apply the BIT index.The ratio of branched GDGTs to crenarchaeol in marine sediments is a function of soil and riverine organicmatter input (and to a minor extent in situ production of brGDGTs also; Peterse et al., 2009), and is expressed as the BIT index (Hopmans et al., 2004): where I, II, III and Cren refer to the structures shown in Supplementary Fig. 1.TEX 86 estimates associated with BIT indices > 0.3, indicative of a potential soil-derived GDGT signal influence, should not be used for SST reconstruction (Weijers et al., 2006).However, this value may be conservative in some settings (e.g., certain Eocene sediments, Inglis et al., 2015) or overly inclusive in others (e.g., Wunstorf core, northern Germany, van Helmond et al., 2015).In addition, measured BIT values vary dramatically across different laboratories (based on findings from interlaboratory comparison studies; Schouten et al., 2009, Schouten et al., 2013a) such that a BIT threshold may be reached by one laboratory but not at another.A way to circumvent this issue is to crosscorrelate each data set internally, i.e.BIT values with TEX 86 values, as an additional control (e.g., van Helmond et al., 2015).Sedimentary methanogenic Archaea can synthesize GDGT-0 and, in lesser quantities, GDGT-1, -2 and -3 (Koga et al., 1993;Weijers et al., 2006).To assess sedimentary GDGT production qualitatively we compute %GDGT-0 (Sinninghe Damsté et al., 2012), an expression of the contribution of sedimentary archaeal methanogen-synthesized GDGTs to the sedimentary GDGT pool: where %GDGT-0 exceeds a threshold value > 67 an additional source of GDGT-0, i.e. sedimentary GDGT production via methanogenesis, and by extension potentially also GDGT-1, -2 and -3, is implied.This threshold is based on the range of %GDGT-0 values observed in enrichment cultures of Thaumarchaeota (Blaga et al., 2009;Sinninghe Damsté et al., 2012;Elling et al., 2015).
In the modern realm, the TEX 86 -SST proxy has been found to deviate from the general TEX 86 -SST relationship in certain ocean regions including the Red Sea and the Mediterranean Sea (e.g., Kim et al., 2008;Kim et al., 2010;Kim et al., 2015), suggesting additional environmental controls on the TEX 86 proxy in these settings.TEX 86 values recorded in Red Sea core-top sediments translate into much warmer TEX 86 H SSTs than measured values by up to ~8 °C (Trommer et al., 2009).As a result, Kim et al. (2008Kim et al. ( , 2010) ) excluded GDGT data from the Red Sea in their global core-top calibration dataset.However, core-top GDGT data from the Red Sea are included in the BAYSPAR calibration dataset (Tierney and Tingley, 2014;Tierney and Tingley, 2015).Similarly, in both the Mediterranean Sea (also a restricted basin) and along the Portuguese continental margin, TEX 86 H values do not correlate well with annual mean SST (Kim et al., 2015;Kim et al., 2016).Kim et al. (2015) discovered a deep-water derived GDGT influence on TEX 86 H values recorded in Mediterranean Sea sediments, which had been previously argued to occur in various settings on the basis of GDGT-2/ GDGT-3 ratios (Taylor et al., 2013).Using both surface sediments and suspended particulate matter collected at < 1000 m water depth in the Mediterranean Sea, Kim et al. (2015) found that proportions of GDGT-2 and also crenarchaeol regio-isomer increased with water depth.This phenomenon resulted in warm-biased TEX 86 H -temperature estimates from sub-surface derived isoGDGTs.At greater water depths, > 1000 m, surface sediment TEX 86 H values no longer correlate with water depth, instead exhibiting a strong relationship with annual mean SSTs (Kim et al., 2015).Similarly, Kim et al. (2016) (Taylor et al., 2013;Kim et al., 2015;Kim et al., 2016), potentially due to the effect of deepwater Thaumarchaeotal communities on sedimentary isoGDGT distributions (Kim et al., 2016), and (2) a statistically different TEX 86 Htemperature correlation from the general global correlation in certain > 1000 m water-depth settings (Kim et al., 2015).
To explore anomalous GDGT distributions in Cretaceous sediments we evaluate variations in the relative abundance of the crenarchaeol regio-isomer in such deposits using a new ratio: Produced by non-methanotrophic marine Thaumarchaeota, a significant enhancement in the proportion of the crenarchaeol regioisomer relative to crenarchaeol, f Cren′:Cren′ + Cren , i.e. a substantial deviation from values observed in the modern core-top dataset (Kim et al., 2010) could be indicative of a non-temperature control, potentially water depth (Kim et al., 2015), on isoGDGT fractional abundances.This ratio does not include GDGT-0, GDGT-1, or GDGT-2, which could have additional sources from methanotrophic Euryarchaeota (see earlier discussion).
We also apply the Ring Index to help distinguish samples that could have been influenced by non-thermal factors and/or deviate from modern analogues (Zhang et al., 2016).The Ring Index (RI) is a weighted average of the ring numbers in GDGT compounds: where In the modern core-top dataset, RI is significantly correlated (R 2 = 0.87; n = 531) with TEX 86 (Zhang et al., 2016); this relationship is expressed by the following quadratic regression: (10) The TEX 86 -RI relationship appears insensitive to shifts in GDGT production (related to depth and/or seasonality), transportation and changes in archaeal community structure, provided that temperature remains the dominant control on GDGT distributions (Zhang et al., 2016).Zhang et al. (2016) suggest that geological samples cannot confidently be used for palaeothermometry if they deviate from the modern TEX 86 -RI relationship: ΔRI values > | 0.3 | are thought to represent samples for which GDGT distributions reside outside the modern TEX 86 -RI relationship, based on the 95% confidence interval (2σ) of the modern regression, reflecting properties distinct from the modern production of GDGTs (Zhang et al., 2016).Using this approach, the Ring Index can help determine samples influenced by either soil-derived isoprenoidal GDGT inputs or methanotrophic archaeal communities and identify samples with atypical GDGT distributions probably impacted by non-thermal factors, e.g., Mediterranean Sea samples from < 1000 m water depth (Kim et al., 2015;Zhang et al., 2016).

Compilation of raw planktonic δ 18 O data
Original published raw planktonic δ 18 O data were collected from locations where palaeotemperature proxy data were of sufficient number and quality (see below) to offer accurate information on Cretaceous climate.Locations of Cretaceous sediments for which raw planktonic δ 18 O data were obtained are shown in Fig. 1b; data references are given in Table 3.

Conversion of planktonic δ 18 O values to sea-surface temperature
Planktonic δ 18 O values were converted to palaeotemperatures using the equation of Bemis et al. (1998).Similar to the TEX 86 proxy, temperature calibrations for δ 18 O pl in modern seawater or culture studies are somewhat limited beyond a maximum of approximately 30 °C (for review see Pearson, 2012).However, synthetic calcite studies (e.g., Kim and O'Neil, 1997) provide support for extrapolation of δ 18 O pl culture studies beyond 30 °C, akin to SSTs indicated for the Cretaceous oceans.Assuming that foraminiferal calcite was precipitated in isotopic equilibrium with Cretaceous seawater, we use a δ 18 O sw value of −1.0‰ standard mean ocean water (VSMOW) to represent the mean isotopic composition of seawater in a non-glacial world (Shackleton and Kennett, 1975), which corresponds to a δ 18 O sw value of − 1.27‰ (PDB) in the temperature equation (Hut, 1987) A δ 18 O sw value of −1.0‰ (VSMOW) lies within the range of more recent estimates for isotopic ice-free conditions, − 1.11 ± 0.03‰ (Lhomme et al., 2005) and − 0.89 ± 0.02‰ (Cramer et al., 2011), implying an overall range in uncertainty of ~0.28‰ equating to ~1.4 °C, although this figure does not account for any temporal variation during the Cretaceous, i.e. deviation from ice-free conditions.Reconstructions of Cretaceous δ 18 O sw from clumped isotope measurements on marine macrofossils indicate a δ 18 O sw value of − 1.0‰ (VSMOW) in the Late Cretaceous (late Campanian-Maastrichtian, Dennis et al., 2013).For the Early Cretaceous, however, estimates of δ 18 O sw from two studies measuring clumped isotopes in belemnites suggest values of −0.1 to + 1.2‰ (VSMOW, Bernasconi et al., 2011) and − 1.1 to +0.1‰ (average − 0.7‰), − 1.8 to −0.4‰ (average − 1.0‰) and − 0.2 to + 0.9‰ (average − 0.4‰) for the Berriasian, early Valanginian and late Valanginian, respectively (Price and Passey, 2013).These δ 18 O sw reconstructions from clumped isotope measurements span a range of latitudes and hence δ 18 O sw values.If episodes of glaciation occurred during the Cretaceous (see later discussion, Section 4.4.3),δ 18 O sw would have been higher, leading to an underestimation of SSTs if a constant δ 18 O sw value of − 1.0‰ (VSMOW) were assumed.Underestimation of SSTs during periods of ice growth would result in an overestimation of the temporal variability between glaciated and nonglaciated climate phases.
In addition, we also convert δ 18 O pl values to temperature with an added adjustment for differences in the oxygen-isotopic composition of local seawater (δ 18 O sw ) owing to changes in site palaeolatitudethe "salinity effect".In the modern surface ocean, δ 18 O sw varies significantly (~− 0.5 to +1.25‰) because local variations in precipitation versus evaporation coupled with fractionation between water and water vapour result in an increase in δ 18 O sw and salinity with enhanced evaporation (see Pearson, 2012, for review).Average latitudinal variations in δ 18 O sw can be predicted for Cretaceous palaeolatitudes by analogy with modern oceans (Zachos et al., 1994).Adjusting Cretaceous δ 18 O pl values for differences in latitude and hence changes in precipitation versus evaporation is thought to reduce the error in Cretaceous δ 18 O pl -SST reconstructions, assuming that latitudinal variations in δ 18 O sw in the Cretaceous were not significantly different from modern.This assumption disagrees with model predictions of meridional δ 18 O sw for the mid-Cretaceous (e.g., Poulsen et al., 1999;Zhou et al., 2008).However, currently available model views of δ 18 O sw only provide information for certain time slices of the Cretaceous, e.g., the Cenomanian (Zhou et al., 2008).As such, we still apply the latitudinal δ 18 O sw correction of Zachos et al. (1994), based on analogy with modern latitudinal gradients in δ 18 O in the surface ocean: Estimates of palaeolatitude, L, for each site were extracted from global palaeorotations provided by Getech Plc.We use the model of Getech Plc to be consistent with modelling efforts (Inglis et al., 2015;Lunt et al., 2016), but we recognize that there are other available tools for generating estimates of Cretaceous palaeolatitudes (van Hinsbergen et al., 2015).

Data quality
In compiling all available planktonic δ 18 O data for the Cretaceous, we undertook an initial screening of the data.We only include published planktonic δ 18 O data from reportedly well-preserved, i.e. goodto-exceptionally preserved specimens.We acknowledge that, by focusing on δ 18 O pl data from well-preserved foraminifera present in clayrich lithologies, the resulting datasets are skewed to sites proximal to the continents.Thus, even sparse SST data obtained from well-preserved carbonate at open-ocean sites, e.g., Maastrichtian δ 18 O data from remarkably well-preserved metastable carbonate from the carbonate platform of Wodejebato Guyot in the western equatorial Pacific (Wilson and Opdyke, 1996), are of particular value.In addition, we have not included planktonic δ 18 O data from older sites that have since been re-drilled using newer coring systems and subsequently sampled with improved recovery.For example, for Demerara Rise we do not include published planktonic δ 18 O data from Deep Sea Drilling Project Site 144 which was spot-cored but instead include data from the recored site, Ocean Drilling Program Site 1257, which recovered a far more complete sample of the stratigraphic record (Shipboard Scientific Party, 2004).Likewise, for Blake Plateau we include published planktonic δ 18 O data from ODP Site 1049 but not DSDP Site 390 (Shipboard Science Party, 1998).
Raw GDGT data from sediments reported to be strongly affected by thermally mature allochthonous organic-matter input, which can lower SST estimates substantially, were also excluded from our compilation: approximately 50% of the sample set from the Cismon core in northern Italy (Bottini et al., 2015).Other tests to screen GDGT data, e.g., the BIT index and the Ring Index, are discussed in the 'Results' section.Owing to a great many differences between the TEX 86 and δ 18 O pl temperatures proxies, e.g.their very definitions and their duration in use (TEX 86 being relatively younger), our discussion of potential secondary effects on TEX 86 (Section 3.2) is somewhat more extensive than that of δ 18 O pl (Section 3.4) with Pearson, 2012 providing a detailed review of the latter topic.

Cretaceous TEX 86 values
Raw TEX 86 values for Cretaceous sediments range between 0.51 and 0.96 (n = 1146;Figs. 2a, 3a), with a substantial number of these data being higher than any value observed in modern ocean sediments outside of the Red Sea (0.72; Fig. 3b; Kim et al., 2010).There is a significant 'time gap' during the late Aptian-Albian, ca.114-105 Ma, for which no GDGT data currently exist.TEX 86 values > 0.9 are recorded for the time intervals 100-87 Ma (encompassing the Cenomanian-Turonian), ~125 Ma (earliest Aptian) and also 139-130 Ma (Valanginian-Hauterivian).These high TEX 86 values generally correspond to low palaeolatitude sites, although not during all time intervals: for example, in the Campanian-Maastrichtian, ~83-67 Ma, TEX 86 values from the Shuqualak-Evans borehole located in Mississippi at 34-36°N palaeolatitude are similar to or warmer than TEX 86 values from PAMA Quarry and the Aderet borehole 1, located in Israel at ~17-19°N palaeolatitude.For some time intervals, data exist only for a small range of palaeolatitudes, precluding observations of latitudinal variations in TEX 86 -derived SSTs.The lowest TEX 86 values are generally recorded at the highest palaeolatitudes (> ± 39°).However, data from high-palaeolatitude locations (> ± 48°) are only available for the Valanginian to Late Aptian, 138-114 Ma, and one Arctic Ocean site dated as early Maastrichtian, ca.71 Ma (Jenkyns et al., 2004).

Critical evaluation of GDGT data
Application of the TEX 86 -palaeothermometer is complicated by factors influencing GDGT distributions other than sea-surface temperature (e.g., Weijers et al., 2006;Zhang et al., 2011;Sinninghe Damsté et al., 2012;Zhang et al., 2016) and the fact that TEX 86 values for the Cretaceous Period commonly lie outside of the upper range of TEX 86 values for the modern core-top datasets (e.g., Kim et al., 2008;Kim et al., 2010).Here we employ a variety of GDGT distribution parameters (BIT, %GDGT-0, MI, f Cren′:Cren′ + Cren , ΔRI) to investigate potential secondary controls on Cretaceous GDGT data and identify any samples that are problematic.We note that this exercise is only possible for samples where BIT values (n = 540) and/or fractional abundances of individual GDGTs (GDGT-1, GDGT-2, GDGT-3, crenarchaeol and the crenarchaeol regio-isomer; n = 810; or n = 678 for samples where GDGT-0 was also included) are available, not just TEX 86 (n = 1143).We then explore the suitability of the TEX 86-Linear (this study), TEX 86 H , (Kim et al., 2010) and BAYSPAR (Tierney and Tingley, 2014;Tierney and Tingley, 2015) calibrations for reconstructing Cretaceous SSTs.

Terrestrial input
BIT indices were reported for approximately 47% of the Cretaceous TEX 86 dataset.For BIT indices > 0.3, we exclude the corresponding TEX 86 data points from our temperature reconstructions.In the case of the Wunstorf core, in line with the original published data (van Helmond et al., 2015), TEX 86 data with BIT indices of 0.15-0.3were also excluded due to additional concerns that these samples were affected by post-depositional oxidation.In total we exclude 29 data points, 2.5% (n = 1143) of the TEX 86 dataset.Given the range of depositional settings and the total number of datasets, it is surprising that so few BIT indices have values > 0.3.One issue is that the BIT threshold of 0.3 was based on a mixing model of the Congo River, an intermediate temperature site (Weijers et al., 2006).Under extreme warmth (such as the Cretaceous) an increase in the fractional abundance of crenarchaeol could produce lower BIT values, potentially shifting the threshold indicated by the Congo mixing model (Weijers et al., 2006).
Cretaceous %GDGT-0 values range between 2 and 79 (Fig. 4), with a mean value of 19 (n = 678, σ = 15.5)suggesting that methanogenic contributions of GDGT-0 are relatively minor.Only two data points, from the Aderet borehole 1, Israel (Alsenz et al., 2013), have %GDGT-0 values > 67, indicating that associated TEX 86 values are potentially compromised by an additional, potentially methanogenic, source of GDGT-0.Cretaceous %GDGT-0 values exhibit a greater range than observed for the modern core-top dataset, 9-65 (n = 426; σ = 12.5; Fig. 4), but are narrower in span than %GDGT-0 values for the Eocene, 5-97 (n = 641; σ = 17.3;Inglis et al., 2015).The average %GDGT-0  for the Cretaceous, 19, is substantially lower than the %GDGT-0 mean for the modern core-top data and the Eocene compilation, which possess similar values of 45 and 42, respectively.To a first order, a low fractional abundance is likely due to the much higher temperatures of the study interval, as derived from this and previous studies (Schouten et al., 2002).As such, for %GDGT-0 to exceed 67% would require a substantial methanogen contribution, and it is likely harder to resolve this influence in warm settings and time intervals, such as the Cretaceous and Eocene.At lower TEX 86 values, ~0.6-0.7,%GDGT-0 values are higher and more variable relative to the modern core-top dataset (Fig. 4b).This result derives from lower abundances of crenarchaeol when TEX 86 values are lower (temperatures are cooler) such that contributions from GDGT-0 are more dominant, while increased variability reflects the range of spatial and temporal depositional settings.
In combination, MI and %GDGT-0 values indicate that Cretaceous sediments are unlikely to have been subjected to appreciable sedimentary GDGT input via archaeal methanogenesis and/or archaeal methanotrophy sufficient to compromise TEX 86 values.

Anomalous GDGT distributions
We apply the f Cren′:Cren′ + Cren ratio and ΔRI to investigate anomalous GDGT distributions in ancient sediments (Fig. 5a).Cretaceous f Cren′:Cren′ + Cren values (n = 810) mostly range between 0.03 and 0.24, which is broadly similar to the range observed in the modern core-top sediments 0.00-0.16with the inclusion of f Cren′:Cren′ + Cren values from the modern Red Sea (Fig. 5a).Similar to the modern core-top dataset, f Cren′:Cren′ + Cren in Cretaceous sediments increases with TEX 86 , suggesting that most Cretaceous sediments have a similar GDGT distribution-temperature relationship as the modern.There are some exceptions to these patterns; some Cretaceous samples deviate from the overall f Cren′:Cren′ + Cren -TEX 86 relationship, displaying f Cren′:Cren′ + Cren values > 0.25 (Fig. 5a).These outliers are from Aderet borehole 1, Israel and DSDP 463, Mid-Pacific mountains.The enhanced contribution of the crenarchaeol regio-isomer, Cren′, in these few outlier sediments could indicate a potential warm bias in corresponding TEX 86 values and/or a different TEX 86 -temperature response; therefore, sediments with f Cren′:Cren′ + Cren values > 0.25 are excluded from our Cretaceous SST compilation.Overall though, the f Cren′:Cren′ + Cren values indicate that conditions in the Cretaceous were very warm.This exercise for Cretaceous sediments (Fig. 5a) suggests that this approach could be a valuable tool for identifying anomalous GDGT distributions in other investigations, i.e. increased proportions of the crenarchaeol  regio-isomer could indicate a potential depth influence, similar to that observed in modern Mediterranean Sea and Portuguese continental margin sediments (Kim et al., 2015;Kim et al., 2016).
ΔRI values calculated for Cretaceous GDGT data range between − 0.71 and 1.33 (Fig. 5b).Of the Cretaceous GDGT distributions evaluated (n = 678), 128 samples have ΔRI values > | 0.30 |, indicating a significant deviation from modern analogues and/or an influence from non-thermal factors.The Ring Index is designed to identify spurious data, including those that would likely be detected by computing MIs and/or BIT indices.The application of the Ring Index here identifies substantially more anomalous/problematic samples than all other indices -BIT index, MI, %GDGT-0 and f Cren′:Cren′ + Crencombined (33 samples total).However, the majority of the Cretaceous ΔRI values > | 0.30 |, 80 samples, are from one location, PAMA quarry, Israel.Of the remaining 48 samples with ΔRI values > | 0.3 | (1049, n = 1; 463, n = 1; 534, n = 3; FL533, n = 1; 1207, n = 2; 1258, n = 1; Shuqualak, n = 10; Aderet borehole 1, n = 6; Bass River, n = 2; Wunstorf, n = 4; DSDP 398, n = 4; Meirs Farm 1, n = 12; Brazos River, n = 1), 9 were also identified as potentially influenced by non-temperature factors based on one or more other indices.The observation that > 76% of the Cretaceous GDGT data from PAMA quarry (n = 105) are identified as spurious by ΔRI, suggests a strong influence from nonthermal factors undetected by the other indices employed here.Given that BIT indices and hopane isomers (maturity indicators) are unavailable for these sediments, which were deposited in the centre of an upwelling system (Edelman-Furstenberg, 2009;Ashckenazi-Polivoda et al., 2011;Schneider-Mor et al., 2012), we exclude all data from all samples with ΔRI values > | 0.30| from our SST compilation and interpret the rest of the data from the PAMA quarry with caution.We also exclude data from all other locations with ΔRI values > |0.30 |.In general, these data represent a small proportion of the total samples at any given site, suggesting that if ΔRI is detecting influences from nonthermal factors they were not persistent.

Cretaceous planktonic foraminifer δ 18 O data
The new planktonic oxygen-isotope (δ 18 O pl ) compilation (n = 3843) indicates that δ 18 O pl values for the Cretaceous (120-66 Ma) range between 1.2 and − 5.4‰ (Fig. 6a).Lowest δ 18 O pl values occur during the Cenomanian-Turonian (MacLeod et al., 2013), highest δ 18 O pl values are recorded in the Late Aptian and Late Maastrichtian, although δ 18 O pl data are unavailable prior to ca. 120 Ma owing to a lack of published records due to a combination of low planktonic foraminiferal abundances, small test sizes, and poor skeletal calcite preservation in Lower Cretaceous sediments.

Planktonic foraminifer species, depth of habitat and seasonality
The range in δ 18 O pl -values and consequently δ 18 O pl -SST estimates for any given site partly reflects the variety of planktonic species from which the data were generated.For most Cretaceous δ 18 O pl -SST datasets, δ 18 O pl values are derived from more than one planktonic foraminiferal species.Importantly, as is observed in the modern ocean, for most sites, e.g., North Atlantic ODP Site 1050 (Ando et al., 2009), planktonic species can represent a combination of annual and/or seasonal mixed-layer dwellers and potentially also some (sub)thermocline species, with peaks in seasonal abundance varying for different taxa.Further complicating temperature-depth-related signals, foraminifera can vary their depth habitat over a life cycle (e.g., Hemleben et al., 1989).In the modern ocean, many planktonic foraminiferal species migrate through the mixed layer (Bijma and Hemleben, 1994;Schiebel and Hemleben, 2005), typically secreting gametogenic calcite in the upper thermocline.For extinct Cretaceous species, foraminiferal calcification regimes are less well understood (Houston et al., 1999;Bornemann and Norris, 2007).
It is problematic to exclude Cretaceous planktonic δ 18 O pl data on the basis of seasonality or depth habitat because of the limited understanding of planktonic foraminiferal ecology in the Cretaceous Period.Paired δ 13 C pl and δ 18 O pl measurements can potentially be used to evaluate foraminiferal habitats and seasonal preferences since δ 13 C pl reflects δ 13 C of dissolved inorganic carbon in the ambient seawater, which will vary both vertically and seasonally.However, Cretaceous foraminiferal δ 18 O pl and δ 13 C pl gradients are commonly insufficient to separate planktonic species on account of depth, i.e. foraminiferal δ 13 C pl reflects multiple factors (e.g., Pearson et al., 2001), resulting in a potential bias towards inclusion of only the lowest (warmest) δ 18 O pl values.Hence, we include all Cretaceous δ 18 O pl data with the understanding that they represent upper-ocean conditions, but not exclusively the surface ocean.

Foraminiferal calcite preservation
Foraminiferal tests can undergo alteration via dissolution, recrystallization and/or the addition of calcite overgrowths, all of which have the potential to compromise δ 18 O pl values (e.g., Pearson et al., 2001).In general, our compilation presents published δ 18 O pl data from foraminifera of good-to-exceptional preservation, ideally demonstrated by 'glassy' tests.However, diagenetic micro-recrystallization, whereby modification occurs but the structure of the shell, e.g., wall pores, is maintained so as to appear unaltered under a scanning electron microscope, could explain some of the high δ 18 O pl values derived from some Cretaceous sediments (Pearson et al., 2001).Diagenetic recrystallization produces higher δ 18 O pl values in planktonic foraminiferal calcite equating to lower palaeotemperatures (Pearson et al., 2001;Sexton et al., 2006), a result of the precipitation of diagenetic calcite from relatively cold bottom waters or pore waters below the sea floor and the fast rate of carbonate recrystallization (Rudnicki et al., 2001;Schrag et al., 1995).This process has the potential to exert the most significant influence on δ 18 O pl values of planktonic foraminifera from low-latitude sites, where temperature differences between surface waters and pore waters are greatest.
Comments on planktonic foraminiferal preservation for each Cretaceous δ 18 O pl dataset are provided in Table 3.The majority of Cretaceous δ 18 O pl datasets report either exceptional 'glassy' preservation or excellent/very good/good preservation (Table 3).We interpret exceptional/excellent preservation to reflect foraminiferal tests with no diagenetic alteration and very good/good preservation to imply minimal diagenetic effects on planktonic δ 18 O pl values, i.e. foraminiferal wall microstructures are visible via SEM, suggesting that primary temperature trends are preserved (Pearson et al., 2001).Evidence of more substantial diagenetic alteration (and correspondingly high δ 18 O pl values) is reported for two sites, ODP Site 1049 and ODP Site 1050.Highest δ 18 O pl values (and lowest δ 18 O pl palaeotemperatures) are reconstructed for the late Aptian North Atlantic (Fig. 6a; ODP 1049: Huber et al., 2011).Huber et al. ( 2011) identified diagenetic overprinting of upper Aptian samples from ODP Site 1049 associated with high δ 18 O pl values, > 0‰, which in some cases were greater than corresponding benthic δ 18 O (δ 18 O b ) values from the same sample set.Manifestly, diagenesis provides one explanation for the high δ 18 O pl values, although early diagenesis cannot fully explain the highest δ 18 O pl values if the benthic values are taken as the diagenetic endmember.In Albian-Cenomanian sediments from North Atlantic ODP Site 1050, Ando et al. (2009Ando et al. ( , 2010) also found foraminifera from certain intervals had undergone some diagenetic recrystallization impacting δ O pl values (Huber et al., 2002;Petrizzo et al., 2008;Ando et al., 2009), although δ 18 O pl values are much lower than for the upper Aptian North Atlantic (ODP 1049; Huber et al., 2011).Based on our assessment of preservation, we exclude these problematic δ 18 O pl data from ODP 1049 and ODP 1050 from our Cretaceous SST compilation.

Changes in surface-water δ 18 O
Differences in surface-water δ 18 O patterns, δ 18 O sw , in the Cretaceous Period relative to the present day likely influenced δ 18 O pl values and corresponding palaeotemperatures on a regional scale.Although we adjust Cretaceous δ 18 O pl values for differences in latitude and hence changes in precipitation versus evaporation (Zachos et al., 1994), this approach assumes that latitudinal variations in δ 18 O sw in the Cretaceous were not significantly different from modern.This supposition is likely to hold true for some locations, such as low-latitude open-ocean sites.However, this approach is rather simplistic because it does not account for regional variations in δ 18 O sw due to both mixing and surface hydrological processes (e.g., Poulsen et al., 1999;Zhou et al., 2008), temporal variations in precipitation/evaporation patterns (Haq, 2014 and references therein), short-term freshening events including those linked with OAEs, e.g., surface water freshening in the north Atlantic during the onset of OAE 1b (Erbacher et al., 2001;Wagner et al., 2008), or uncertainties in the effective fractionation between seawater and exported water vapour, i.e. the slope of δ 18 O sw versus salinity (Zhou al., 2008;Huber et al., 2011).Indeed, a high evaporative fractionation factor (δ 18 O sw /salinity) producing higher δ 18 O sw at elevated salinities in the early Albian North Atlantic, i.e. an accelerated hydrological cycle (Ufnar et al., 2004), is used by Huber et al. (2011) to explain why δ 18 O pl values are exceptionally high at ODP 1049 in the early Albian and thereby yield unreasonably cool temperatures for the mid-Cretaceous (~10-16 °C cooler than modern; Fig. 6).

Kinetic effects -changes in carbonate ion chemistry
Cretaceous δ 18 O pl values may have a cool bias owing to relatively low activity of the carbonate ion in a high pCO 2 ocean (Spero et al., 1997;Zeebe, 2001;Tyrrell and Zeebe, 2004).Royer et al. (2004) suggested a 3 to 5 °C cool bias in shallow-water carbonate δ 18 O palaeotemperature estimates for the Cretaceous.However, subsequent studies (e.g., Beck et al., 2005;Uchikawa and Zeebe, 2010) indicate that the nature of oxygen-isotope fractionation within the carbonic acid system is less systematic than previously thought; for example, thermodynamic theory predicts a much greater change in the δ 18 O of the dissolved inorganic carbon species per unit pH between pH 6 and 8 than between pH 8 and 9 (Uchikawa and Zeebe, 2010), although it is not known if the same trend of δ 18 O versus pH applies to the δ 18 O of foraminiferal calcite (Zeebe, 1999;Uchikawa and Zeebe, 2010).Indeed, a culturing study by De Nooijer et al. (2009) demonstrated that planktonic foraminifera could regulate intracellular pH during calcification.Even for the modern day, comparisons between planktonic foraminifer δ 18 O values and local SSTs frequently reveal under-or over-estimations in reconstructed SSTs, presumably due to δ 18 O-disequilbrium related to physiological and ecological effects (Niebler et al., 1999;Mohtadi et al., 2011), despite parameters such as salinity, δ 18 O sw and the carbonate system itself (e.g., Hönisch et al., 2013).Given these uncertainties, we make no formal adjustment to δ 18 O pl data for kinetic effects associated with past changes in ocean pH but note that the effect of lowering  Shackleton and Kennett, 1975;Hut, 1987).(c) Same as previous but with the addition of a δ 18 O sw correction for changes in palaeolatitude (Zachos et al., 1994).Estimates of palaeolatitude are derived using a palaeorotational model provided by Getech Plc.In the case of site TDP, all data pertain to one sample and so instead the mean is plotted with full ranges of the dataset represented by horizontal and vertical bars.In panels (b) and (c) the black bar represents the paleotemperature calibration uncertainty ( ± 0.7 °C, Bemis et al., 1998).
Our evaluation of potential secondary controls on δ 18 O pl values suggests that most available Cretaceous δ 18 O pl data are supposed to reflect surface temperatures or underestimate SSTs owing to kinetic and/or diagenetic effects biasing δ 18 O pl values to higher values (this potential underestimation of SSTs equates to an uncertainty on the order of ~3 to 5 °C).There are also some notable data outliers; based on our assessment: for example, we exclude ODP 1049 Aptian-Albian data from our temperature compilation owing to exceptionally high δ 18 O pl values: potentially due to a combination of salinity influences and variable preservation and similarly exclude data from ODP 1050 in cases where samples were reported to suffer from preservation issues.

Comparison of TEX 86
H -, TEX 86-Linear -and BAY GlobR -SSTs for the Cretaceous Compiled TEX 86 -derived SSTs for the Cretaceous Period range between 18.8 ± 2.5 °C and 37.5 ± 2.5 °C (TEX 86 H ; Fig. 7a), or 18.8 ± 2.0 °C and 45.1 ± 2.0 °C (TEX 86-Linear ; Fig. 7b), or 16.3 °C and 44.8 °C ( ± ~7-9 °C 90% confidence; BAY GlobR ; Fig. 7c), depending on the choice of calibration.The differences in maximum values arise from the fact that a significant proportion of the data exceed the highest TEX 86 value in the calibration dataset (TEX 86 = 0.72 or TEX 86 = 0.89 if Red Sea data are included, which is relevant for BAY GlobR ; Kim et al., 2010, Tierney andTingley, 2015), requiring us to project the calibration outside of its modern range.Furthermore, such high TEX 86 values are near the limit of the calibration (TEX 86 = 1, TEX 86 H -SST = 38.6 °C; Fig. 7).
Moreover, the logarithmic TEX 86 H calibration has structured residuals at high temperatures (elements of variation which are unexplained by the fitted model) that make it particularly problematic when applied beyond the calibration range (Tierney and Tingley, 2014).However, maximum SSTs predicted by TEX 86-Linear and BAY GlobR are exceptionally warm, > 40 °C, to the extent these temperatures raise questions regarding the maximum heat stress Cretaceous plants and mammals could have tolerated (e.g., Hay and Floegel, 2012).Given the different strengths and weaknesses of the TEX 86 calibrations, we present and discuss compiled Cretaceous TEX 86 -SSTs derived using both the TEX 86 H -SST calibration and also a linear approach, TEX 86-Linear , the latter yielding SSTs very similar to those using BAY GlobR .

Cretaceous planktonic foraminifer δ 18 O-SSTs
Cretaceous δ 18 O pl -SSTs vary between 4.6 ± 0.7 °C and 36.6 ± 0.7 °C (Fig. 6b) and between 5.6 ± 0.7 °C and 40.0 ± 0.7 °C corrected for palaeolatitude (Fig. 6c).The two δ 18 O pl -SST reconstructions (Fig. 6b, c) exhibit very similar long-term trends with warmest surface-ocean conditions occurring during the Cenomanian-Turonian interval at both low palaeolatitude and higher mid-palaeolatitude locations.Application of a δ 18 O sw palaeolatitude correction (Eq.( 13)) results in greater offsets in δ 18 O pl -SSTs between different sites, e.g., between low-palaeolatitude and higher mid-palaeolatitude locations through the Albian to Coniacian (~6-9 °C versus ~3-5 °C offset) and between lower mid-latitude and higher mid-latitude sites in the Campanian-Maastrichtian (~6 °C versus ~3 °C offset, Fig. 6b, c).Aside from the long-term trends, δ 18 O pl -SST estimates are often highly variable at any given site (Figs.6b, c, 7).This variability is unlikely to be strongly driven by preservation since (a) we have only included data derived from good to excellently preserved planktonic foraminiferal specimens, and (b) this variability exists within individual core datasets.Instead, this likely reflects some combination of habit depth, species variety, seasonality, local and regional fluctuations in δ 18 O sw resulting from changes in oceanic circulation, source water, continental runoff and the hydrologic cycle, and short-term climate variability driven by orbital cyclicity.-SSTs indicate maximum Cretaceous SSTs, 37-38 ± 2.5 °C that are similar to maximum δ 18 O pl -SSTs (corrected for palaeolatitude) but record slightly higher minimum temperature estimates, 33-34 ± 2.5 °C compared with 29-32 ± 0.7 °C, although these temperature offsets are typically within associated proxy errors (Fig. 8a).By comparison, both maximum and minimum TEX 86-Linear -SSTs at ODP Sites 1258, 1259 and 1260 are warmer than corresponding δ 18 O pl -SSTs by ~7-8 °C (Fig. 8b).These observations are similar to those for the Eocene where (i) temporal trends in TEX 86 -, TEX 86 H -and δ 18 O pl -SSTs are similar but (ii) absolute δ 18 O pl -SST estimates are similar or lower than those of TEX 86 -and TEX 86 H -SSTs (Zachos et al., 2006;Pearson et al., 2007;Hollis et al., 2012).

Comparison of TEX
For the Cretaceous mid-latitudes, offsets between TEX 86 H /TEX 86- Linear -and δ 18 O pl -SSTs are more apparent (Fig. 8).Indeed, although δ 18 O pl data from DSDP Site 525 (~38°S), ODP 762 (~44-45°S), ODP 1049 (~30°N) and ODP 1050 (~30°N) represent four lower mid-palaeolatitude settings, these data give considerably lower SST estimates -SST calibration (Kim et al., 2010).(b) SST estimates derived using a linear TEX 86 -SST calibration, TEX 86-Linear , based on "warmer" ocean data, > 15 °C, from the global core-top calibration of Kim et al. (2010) and excluding data from the Red Sea (see Supplementary Fig. 2 and text for further details).(c) SST estimates derived using a Bayesian global core-top regression (BAY GlobR ; Tierney and Tingley, 2015).The red horizontal lines indicate the upper limits of the modern TEX 86 H -, TEX 86-Linear -and BAY GlobR -SST calibration datasets.The black bar in each panel represents the uncertainty in the corresponding paleotemperature calibration.In the case of BAY GlobR , since the calibration uncertainty is expressed as confidence intervals for each individual TEX 86 datapoint, we approximate the BAY GlobR uncertainty to the mean (n = 993) width of the 90% confidence interval.As in Fig. 2, for individual datasets where only a few datapoints exist and/or the data span a very narrow temporal window or age uncertainties are large, only the mean is plotted with full ranges of the datasets represented by horizontal and vertical bars.(For interpretation of the references to colour in this figure legend, the reader is referred to the web version of this article.)than TEX 86 H /TEX 86-Linear -derived SST estimates from other lower midlatitude locations for the Campanian-Maastrichtian, namely the Shuqualak-Evans Borehole (35-36°N) and Brazos River (~36°N).Similarly, δ 18 O pl -SST estimates from ODP Site 690 (~62-64°S palaeolatitude) give cooler SST estimates than TEX 86 H /TEX 86-Linear -SST values from a higher palaeolatitude location (FL533; ~80°N), although in this case the two proxy records are from different hemispheres.Cretaceous δ 18 O pl estimates offer better agreement with TEX 86 H /TEX 86-Linear -SST estimates when a palaeolatitude δ 18 O sw adjustment is applied.However, in several cases (including both lower mid-palaeolatitude and higher mid-palaeolatitude records), δ 18 O pl -SST estimates are still significantly lower than TEX 86 H /TEX 86-Linear -SST estimates, suggesting that either assumptions regarding changes in δ 18 O sw are inaccurate and/or other controls on foraminiferal δ 18 O pl are important; alternatively, the TEX 86 proxy may be overestimating SST, particularly at higher palaeolatitudes, as suggested for the early Eocene (Hollis et al., 2012).Unfortunately, current Cretaceous TEX 86 -and δ 18 O pl -SST data are insufficient to determine how well these proxies agree at high latitudes.
Data from several high-resolution events are included in this compilation (Figs. 6,7,8): namely, datasets for OAE 1a (early Aptian, ~125 Ma) and OAE 2 (Cenomanian-Turonian boundary, ~94 Ma) and also a high-resolution δ 18 O pl dataset from the Albian of the Lower Saxony Basin (~107 Ma).The high-resolution datasets indicate changes in SSTs, both δ 18 O pl -and TEX 86 -SSTs, over relatively short periods of time, on the order of 0.1-1 Myr.Again, the variability is greater for the δ 18 O pl proxy than for the TEX 86 proxy, such that the latter yields SSTs that are higher but less variable across the entire data set.These differences imply that either δ 18 O pl variability is driven by more than simple orbital cyclicity, e.g., seasonality, depth of habitat and variations in local δ 18 O sw , or that the TEX 86 proxy is less sensitive to these short-term variations.
4.3.1.Possible reasons for disagreement between TEX 86 and δ 18 O pl -SSTs One major difference between the δ 18 O pl and TEX 86 proxies is the way in which the water column represented by the proxy is treated and interpreted.The δ 18 O pl palaeothermometer is based on the calibration of δ 18 O pl values with growth temperature in laboratory culture studies (e.g., Bemis et al., 1998) or, alternatively, via modern core-top calibrations (e.g., Mulitza et al., 2004).However, as discussed earlier, the data likely reflect a range of different growth (and hence temperature) depths (e.g., Mohtadi et al., 2011;Hönisch et al., 2013).By contrast, the TEX 86 proxy is based on calibration of TEX 86 values in core-top sediments with overlying SST and therefore assumes that the depth of export of the TEX 86 signal is predominantly from the upper mixed-layer and that this was the same in the past as it is today (e.g., Schouten et al., 2013b).This assumption has been questioned (Hernández-Sánchez et al., 2014;Kim et al., 2015) but it is partially incorporated into the calibration uncertainty along with seasonality, unlike the δ 18 O pl proxy.However, what remains unclear is how these export processes vary spatially and temporally, resulting in deviations from the modern coretop calibrations (Taylor et al., 2013;Kim et al., 2015), i.e. to what extent signal depth and seasonality in the Cretaceous was analogous to modern.Another likely reason for disagreement between the δ 18 O pl and TEX 86 proxies is the uncertainty in δ 18 O pl values related to the combined effects of salinity and variations in carbonate ion concentration.It is likely that these factors are important, particularly as regards kinetic effects that can potentially lower δ 18 O pl -SST estimates on the order of ~3-5 °C (Crowley and Zachos, 2000;Royer et al., 2004;Uchikawa and Zeebe, 2010).Unfortunately, it remains difficult to assess the extent of these influences on δ 18 O pl values as these variables differ both temporally and spatially and are relatively poorly constrained for the Cretaceous Period.

Cretaceous climate maxima
For over a decade, climate reconstructions of the Cretaceous greenhouse have yielded palaeotemperature estimates that encompass the highest (> 35 °C) proxy data-based SST estimates of the last 150 Myr 'super greenhouse' conditions, with the most pronounced warming in the late Cenomanian-Turonian (e.g., Jenkyns et al., 1994;Huber et al., 1995;Norris and Wilson, 1998;Clarke and Jenkyns, 1999;Wilson et al., 2002;Schouten et al., 2003;Voigt et al., 2004;Forster et al., 2007a;Forster et al., 2007b;Bornemann et al., 2008;Jarvis et al., 2011;MacLeod et al., 2013;van Helmond et al., 2014).The specific timing of this thermal maximum, whether at the Cenomanian-Turonian stage boundary or in the early or late Turonian, is still debated (e.g., Voigt et al., 2004).This pronounced warming is also reflected in palaeotemperature estimates from brachiopods of temperate conditions in mid-latitude shelf seas (Voigt et al., 2004).Our data compilation (Fig. 8) indicates that during most of the Cretaceous Period, temperatures of the warmest surface waters were significantly greater than the maximum surface temperatures recorded in the modern ocean (~28-30 °C).Maximum SSTs did not drop below 30 °C until the late Santonian (~84 Ma, δ 18 O pl ) or early Campanian (~80 Ma, TEX 86 ) during global cooling that started in the Coniacian, ~88 Ma.Our compilation (Figs. 6,7,8)   with values ranging between 17 °C and 36 °C (Fig. 8).This trend agrees with Earth system model zonal mean SSTs that, along with SST proxy data from a range of sources (not only GDGTs and planktonic foraminifera but also fish tooth enamel, mollusks, bivalves, brachiopods and belemnite rosta), together confirm that Late Cretaceous SST cooling was widespread (Tabor et al., 2016).

Cretaceous climate minima
Mean global minimum Cretaceous surface-ocean temperatures are less well constrained than maximum global mean SSTs because of critical gaps in proxy temperature data from higher latitudes (> ca.48°p alaeolatitude; Fig. 8) relative to data from low latitudes (< ca.30°p alaeolatitude).In particular, SST estimates are provided solely by the TEX 86 proxy prior to 114 Ma (this limitation applies for the low-latitude data also), whereas after this time SST estimates for higher latitudes are almost exclusively provided by δ 18 O pl data, with the exception of a small dataset from the lower Maastrichtian of the Arctic Ocean (Jenkyns et al., 2004).Interpreting each proxy separately, TEX 86 H and TEX 86-Linear from higher mid-latitude locations (> ± 48°palaeolatitude) suggest minimum SSTs > 22 °C (TEX 86 H /TEX 86-Linear ) in the Early Cretaceous and SSTs of ~19 °C (TEX 86 H /TEX 86-Linear ) for the early Maastrichtian.δ 18 O pl data from higher mid-latitude locations (> ± 48°p alaeolatitude) suggest minimum SSTs > 13 °C in the Albian (δ 18 O pl ) and > 20 °C (δ 18 O pl ) in the Turonian, while later in the Cretaceous minimum SSTs were ~10 °C (δ 18 O pl ) and ~5 °C (δ 18 O pl ) in the early and late Maastrichtian, respectively.Clumped-isotope palaeothermometry estimates from high-latitude belemnites indicate marine temperatures of 10-20 °C in the sub-Arctic (60-65°N) for the Berriasian to late Valanginian, 145-134 Ma (Price and Passey, 2013).These temperatures derived from sub-Arctic nekto-benthonic organisms are broadly consistent with the lower palaeolatitude estimates of SSTs from TEX 86 (Price and Passey, 2013).
The lowest SSTs during the Cretaceous are not observed at the higher latitude sites.In the Early Cretaceous, ~133-124 Ma, the lowest SSTs are in fact recorded at lower mid-latitude sites plus one low-latitude palaeolocation.In the late Valanginian to early Aptian (135-124 Ma), TEX 86 proxy data from sites in England (Speeton), Germany (Gott and Moorberg) and Italy (Cismon) yield similar or slightly cooler minimum SSTs, ~21-22 °C, compared with higher latitude sites (~25-27 °C, ODP Site 766 and DSDP Site 511).Similarly, in the Late Cretaceous, Campanian to late Maastrichtian (83-69 Ma), lowpalaeolatitude TEX 86 data from the Southern Tethys margin (Aderet borehole 1 and PAMA quarry, Israel; only the datapoints which were not excluded after screening) indicate cooler temperatures than contemporaneous TEX 86 data from a lower mid-palaeolatitude site, the Shuqualak-Evans borehole (Fig. 7).This difference may reflect the Southern Tethys oceanographic setting, namely, a high-productivity upwelling system (Alsenz et al., 2013).Modern SST reconstructions based on TEX 86 values in suspended organic matter from the Santa Barbara Basin, an upwelling area, and also the Benguela upwelling system, were suggested to reflect mainly cooler, deeper water temperatures (Huguet et al., 2007;Lee et al., 2008;Seki et al., 2012).For the Cismon core, given that some of the sediments (excluded here) are known to have been strongly affected by thermally mature allochtonous organic matter input (Bottini et al., 2015), it is more likely that TEX 86 data from the lowest maturity Cismon sediments are also affected by allochtonous input, producing cooler SST estimates, rather than upwelling or cold surface currents.

Evidence for glaciation in the Cretaceous
The occurrence of cooler episodes and/or ephemeral ice sheets at times in the Cretaceous has been proposed and debated in numerous studies (e.g., Kemper and Schmitz, 1981;Stoll and Schrag, 1996;Price, 1999;Stoll and Schrag, 2000;Miller et al., 2005;Bornemann et al., 2008;Price and Passey, 2013;Ladant and Donnadieu, 2016).Below we discuss evidence for glaciation in the Cretaceous in light of our Cretaceous SST compilation.4.4.3.1.Evidence for glaciation in the Early Cretaceous.Evidence for cooler episodes during the Early Cretaceous was originally provided by a variety of mineralogical findings including the presence of glendonites (e.g., Kemper and Schmitz, 1981;Kemper, 1987;Price and Nunn, 2010;Rogov and Zakharov, 2010;Grasby et al., 2017), putative tillites and dropstones (Frakes and Francis, 1988;Frakes et al., 1995) and supposed glacial diamictite (Alley and Frakes, 2003) at various high-latitude locations, including Arctic Canada and Australia.Possible cold phases in the Early Cretaceous have also been inferred from δ 18 O values measured in fish teeth enamels (Pucéat et al., 2003;Barbarin et al., 2012), belemnites (Pirrie et al., 1995;Pirrie et al., 2004;Price and Mutterlose, 2004;McArthur et al., 2007;Bodin et al., 2015), carbonate concretions from fluvial sediments (Ferguson et al., 1999), reptile remains (Amiot et al., 2011) and hydrothermal zircon (Yang et al., 2013), as well as bipolar distribution patterns of calcareous nannofossils (Mutterlose et al., 2003;Mutterlose et al., 2009) and the presence of steryl alkyl ethers (Brassell, 2009).Although suggestive of (local) cooler phases, not all of these findings necessarily provide definitive evidence for glaciation, or at least the spatial extent of it.Glendonites, for example, can form in marine and continental bottom waters of up to 4 °C and 7 °C, respectively (De Lurio and Frakes, 1999, and references therein), but potentially also higher temperatures, in association with methane seeps (Teichert and Luppold, 2013).Dropstones can derive from a number of rafting agents including vegetation (e.g.driftwood), kelp, corals and pumice, vertebrates (via ingestion), as well as icebergs, sea and lake ice (Bennett and Doyle, 1996).Indeed, Lower Cretaceous glendonites and dropstones have been mainly found in epicontinental/shelf seas indicating a potential terrestrial origin (e.g., Kemper, 1987;Frakes and Francis, 1988).
Since most available Cretaceous SST proxy data derive from locations < ± 54°palaeolatitude, inferences about conditions at polar latitudes rely on extrapolation of these data.The persistence of very warm, ~30-35 °C, (sub)tropical surface waters in the Cretaceous does not preclude the presence of ice at the highest latitudes (and/or at high altitudes), but would require an exceptionally steep Pole-Equator thermal gradient at higher latitudes in order to achieve sustained freezing temperatures.Indeed, clumped-isotope palaeotemperature estimates from belemnites for the Berriasian to late Valanginian imply generally warm (10-20 °C) polar (60-65°N) conditions consistent with a greenhouse climate punctuated by cooler intervals when transient polar ice was possible, particularly at high altitudes (Price and Passey, 2013).Such equable conditions (at least for ± 0-54°palaeolatitude) in the Early Cretaceous likely occurred in combination with a hydrosphere significantly different from today, with implications for moisture transport and potential cryosphere development at high southern latitudes (Flögel et al., 2011).The current paucity of high-latitude, > 70°S, climate data for the Early Cretaceous along with a lack of detailed understanding of the elevation history of Antarctica implies that, for now, evidence for an Early Cretaceous icehouse from available 'cool temperature' proxy data remains equivocal (Jenkyns et al., 2012).4.4.3.2.Evidence for a Late Aptian 'cold snap'.On shorter timescales in the Early Cretaceous, micropalaeontological and sedimentological studies (Kemper, 1987;Mutterlose et al., 2009) along with TEX 86derived SSTs from the subtropical Atlantic (McAnena et al., 2013) suggest that episodic (< 2 Myr) interludes of global cooling may have occurred in the Aptian-Albian, in particular a late Aptian 'cold snap' ~114 Ma during prolonged cooling in the late Aptian (Bodin et al., 2015;Erba et al., 2015).In further support of this contention, Herrle et al. (2015) found glendonite beds in upper Aptian to lower Albian sediments of the Polar Sea, ~118-112 Ma, interpreted to reflect the presence of cool shelf waters in the High Arctic at this time, but not necessarily ice.4.4.3.3.Evidence for glaciation in the Middle-Late Cretaceous.In the mid-Cretaceous, an interval of glaciation has been reported for the mid-Turonian based on discrete parallel shifts in both planktonic and benthic oxygen-isotope records and sea-level fluctuations (Miller et al., 2004;Bornemann et al., 2008;Galeotti et al., 2009).However, recent δ 18 O temperature reconstructions from Turonian shelf sediments of Tanzania (~25°S palaeolatitude; MacLeod et al., 2013) indicate stable, hot temperatures, suggesting that the mid-Turonian excursion may not be a global feature (Stoll and Schrag, 2000).An effectively icefree mid-Cretaceous climate is consistent with our compiled estimates of exceptional low-palaeolatitude and lower mid-palaeolatitude warmth during this time (~15-40 °C, both TEX 86 and δ 18 O pl proxy estimates; Fig. 8).Some modelling approaches have indicated that ice growth is possible in the mid-Cretaceous under specific circumstances, i.e. when pCO 2 is < 800 ppm (e.g., Barron et al., 1995;Tabor et al., 2016).However, a recent mixed resolution modelling study suggests that the palaeogeography of the Cenomanian-Turonian renders the Earth System resilient to the inception of Antarctic glaciation under CO 2 concentrations as low as 420 ppm (Ladant and Donnadieu, 2016).Furthermore, at such low pCO 2 , modelling studies also indicate that SSTs in extra-tropical regions will not be > ~30 °C (e.g., Bice and Norris, 2002).In the Upper Cretaceous, planktonic and benthic oxygenisotope records show a discrete positive shift in the early Maastrichtian (e.g., Barrera and Savin, 1999;Huber et al., 2002;Friedrich et al., 2009;Friedrich et al., 2012).This excursion may reflect a short-lived glaciation event (e.g., Barrera and Savin, 1999;Miller et al., 1999) and/or cooling of ocean bottom waters associated with a reorganisation of intermediate-to deep-water sources (e.g., Barrera et al., 1997;Friedrich et al., 2009).
In summary, our SST compilation (Fig. 8) documents very warm global Cretaceous SSTs, particularly during the middle to late Aptian and the Cretaceous Thermal Maximum (Cenomanian-Turonian).Under such conditions of global warmth, episodes of sustained continental glaciation seem improbable, at least until the latest Cretaceous.

Latitudinal SST gradients during the Cretaceous
Available evidence suggests that latitudinal temperature gradients were lower in the Cretaceous Period compared with the present day (e.g., Barron, 1983;Huber et al., 1995;Barrera and Savin, 1999;Bice and Norris, 2002;Huber et al., 2002;Littler et al., 2011).This conventionally accepted view is confirmed by our Cretaceous SST compilation, that indicates, regardless of the choice of proxy, TEX 86 or δ 18 O pl , generally low latitudinal SST gradients between low-and higher midpalaeolatitude sites, with small SST differences (ΔSST = 3-17 °C prior to the late Campanian; Figs. 8,9a,b).This observation stands even if the most extreme calibration comparison is made, higher latitude δ 18 O pl -SSTs with low-latitude TEX 86-Linear or BAY GlobR -SSTs (Fig. 9b).
Our compilation allows examination of the temporal evolution of higher/higher mid-versus low-latitude SST conditions.We compute the SST gradient (low-higher mid, ΔSST) based on compiled TEX 86 H and δ 18 O pl higher mid-palaeolatitude (48-54°S palaeolatitude, with the addition of a small amount of data from Arctic Ocean Core FL533, ~80°N palaeolatitude) and low-palaeolatitude (< ± 30°) SST estimates, separately fitted with a LOESS smooth function, span = 0.5, using the PAST software package (Hammer et al., 2001;Fig. 9a).Note, in general there is little difference if the Cretaceous latitudinal SST gradient is calculated with TEX 86-Linear SST data rather than TEX 86 H SST data (Fig. 9b) or solely δ 18 O pl SST data (Supplementary Fig. 3), the resulting latitudinal SST gradient reconstructions often overlap.Uncertainty envelopes ( ± 95% confidence band) for the curves were calculated using an inbuilt bootstrap method based on 999 random replicates (Fig. 9a, b).LOESS curves were interpolated at 1.0 Myr resolution in order to compute the gradient.Limitations of the SST proxy data (e.g., signal depth; see earlier discussion) mean that temporal trends are likely more robust than absolute SST estimates using this approach.Importantly, we regard our ΔSST reconstruction for the Early Cretaceous (Fig. 9c) as tentative since a crucial caveat of the available data is that Early Cretaceous ΔSSTs are based on TEX 86 -SST estimates (light blue) and δ 18 O pl SST estimates, after interpolating at 1 Myr resolution.Dashed lines represents 'tentative' latitudinal SST temperature gradient reconstructions for the Early Cretaceous, calculated in the same way as for the middle-Late Cretaceous but only using TEX 86 H /TEX 86-Linear SST data.Note, for the Albian-Santonian interval the different latitudinal SST gradient reconstructions overlap.The modern mean annual latitudinal SST gradient between ± 0-30°and ± 48-54°is also indicated (Kim et al., 2010) (d) Global benthic oxygen-isotope, δ 18 O b , data (Friedrich et al., 2012); palaeotemperatures also indicated.In all cases, narrower shaded band represents 95% confidence interval for the LOESS fit, derived using a bootstrap technique.Wider shaded bands, panels (a, b, c), represent the calibration error associated with the TEX 86 -SST estimates; ± 2.5 °C for the TEX 86 H calibration (Kim et al., 2010) and ± 2.0 °C for TEX 86-Linear calibration.The calibration error for the δ 18 O pl -SST data is not indicated but is substantially smaller, ± 0.7 °C (Bemis et al., 1998) It should be noted that our reference to 'latitudinal SST gradient' does not represent the full Equator-Pole temperature differential for the Cretaceous Period but the difference between our compilation of SSTs for low-and higher mid-latitude sites.Specifically, all data from 'lowlatitude' sites encompass SST estimates for ± 5-30°palaeolatitude, whereas 'higher mid-latitude' sites included in this compilation, with the exception of the Arctic Ocean (Core FL533, ~80°N palaeolatitude) originate from sites located between ± 48-54°palaeolatitude in the Southern hemisphere.Nevertheless, the observed long-term shifts in the magnitude of the low-minus higher mid-latitude SST gradient, at least from the middle Cretaceous (weak) to the Late Cretaceous (less weak) likely resemble long-term trends in Equator-Pole temperature gradients.
Interestingly, the ΔSST gradient appears weakest in the Aptian, rather than at the Cretaceous Thermal Maximum during the Cenomanian-Turonian interval.This observation reflects, at least in part, the fact that the latitudinal gradient is generated from both higher mid-latitude and low-latitude TEX 86 data for the Early Cretaceous, and from a combination of low-palaeolatitude TEX 86 and δ 18 O pl data minus higher mid-palaeolatitude δ 18 O pl data for the middle to Late Cretaceous (Fig. 9a, b, c).This switch in the comparison of data from different proxies may influence temporal trends, particularly for the high palaeolatitudes where minimum δ 18 O pl -SST estimates are cooler than minimum TEX 86 -SSTs (Figs. 8,9a,b) resulting in stronger latitudinal SST contrasts.In addition, this observation may be biased by a combination of a lack of mid-Cretaceous TEX 86 data from high palaeolatitudes (> ± 48°) and potential secondary influences either from cool biases in some δ 18 O pl records, inaccuracies in assumptions for conversion of δ 18 O pl values to palaeotemperature estimates, and/or warm biases in TEX 86 -SST estimates.
The LOESS fit approach provides a broad picture of the evolution of SSTs and ΔSSTs during the Cretaceous and, within error, captures the full range of raw SST estimates when the TEX 86 H calibration is applied (Fig. 9a).On shorter timescales, certain details may, however, be lost.
During the Cretaceous Thermal Maximum, SSTs at higher mid-latitudes (> ± 48°palaeolatitude) reached values similar to the tropics, ~18-31 °C, suggesting a dramatic collapse of the latitudinal SST gradient.However, the LOESS model predicts higher mid-latitude SSTs of ~25 °C and therefore does not record this potentially significant decrease in the latitudinal SST gradient at this time (Fig. 9a, b, c, Supplementary Fig. 3), suggesting that the LOESS model may overestimate the ΔSST gradient during the Turonian.Overall, taking into account uncertainties in the LOESS model and associated proxy errors, our latitudinal ΔSST reconstruction tentatively indicates that the latitudinal thermal gradient was higher in the Early Cretaceous and then weakened around the time of the early Aptian OAE 1a event.The ΔSST gradient was likely relatively low in the middle Cretaceous and then increased in the Late Cretaceous (Campanian-Maastrichtian).These reconstructions suggest that latitudinal temperature gradients were in general very low when global temperatures were high, in accord with similar observations made for the Eocene (Bijl et al., 2009;Hollis et al., 2012;Inglis et al., 2015).This result remains valid for the Cretaceous Period regardless of the combination of proxies used (TEX 86 and δ 18 O pl ).Moreover, this weak Cretaceous latitudinal SST gradient, at least for the middle to early Late Cretaceous, offers little support for periods of large-scale ice growth.

Comparison with benthic δ 18 O records
Comparisons between our SST compilation and global benthic foraminifer oxygen-isotope records for the middle-Late Cretaceous climate (Fig. 9d) suggest that, in general, surface and intermediate to deepwater temperatures responded similarly to changes in climate during this time interval.Friedrich et al. (2012) separate their δ 18 O b compilation into four intervals; (1) increasing temperatures from 112 to 97 Ma, Albian to mid-Cenomanian; (2) the subsequent Cenomanian-Turonian hot greenhouse interval, 97 to 91 Ma; (3) a longlasting cooling trend between 91 and 78 Ma, late Turonian to mid-Campanian; and (4) the interval after 78 Ma, mid-Campanian to Maastrichtian, with small inter-ocean δ 18 O b values.Similar to the δ 18 O b compilation (Fig. 9d; Friedrich et al., 2012), our SST compilation also suggests increasing temperatures from the latest Aptian through Cenomanian interval, although the evidence for this trend is affected by the scarcity of SST data (particularly TEX 86 data) during the Albian.Our SST compilation clearly demonstrates the significance of the hot Cenomanian-Turonian greenhouse, particularly at low palaeolatitudes and lower mid-palaeolatitudes (Figs. 8,9a,b).The cooling recorded in the δ 18 O b compilation during the late Turonian to mid-Campanian is evident in SSTs from both low palaeolatitudes and lower mid-palaeolatitude locations, and likely also higher mid-palaeolatitudes, although the data again require cautious interpretation owing to their paucity.For the final interval, mid-Campanian to Maastrichtian, the δ 18 O b values of all ocean basins are similar.Such regionally similar bottom-water palaeotemperatures contrast with surface-ocean conditions, which display the strongest latitudinal differences and also the most significant cooling of the entire Cretaceous during this time.These differences likely reflect the strong influence of changes in inter-basin exchange on deep-water temperature signatures or, alternatively, changes in watermass formation.A full connection between all ocean basins following the complete opening of the Equatorial Atlantic Gateway would have filled the North Atlantic with cool high-latitude waters, homogenizing global δ 18 O b values for the latest Cretaceous and reorganizing atmospheric circulation and the hydrological cycle both regionally and globally.Likewise, more deep-and/or intermediate-water formation at high latitudes and a reduction of the proposed formation of warm and salty water masses in the subtropics could have resulted in small δ 18 O b interbasin gradients (Barrera et al., 1997;Friedrich et al., 2004;Friedrich et al., 2009).

Missing climate proxy data in time and space
Our Cretaceous SST compilation highlights several important 'gaps' in the data currently available that, were they to be filled, would significantly improve understanding of Cretaceous surface-ocean conditions.Temporally, there is a paucity of SST data for the Berriasian and Albian stages.For proxy comparison purposes, there are only five sites, ODP Sites 1049, 1258, 1259, 1260 and DSDP Site 511 for which both δ 18 O pl and TEX 86 data exist, although for DSDP 511 the δ 18 O pl and TEX 86 data are stratigraphically non-contemporaneous.Part of the challenge is finding Cretaceous sediments with sufficient organic matter where foraminifera are still well-preserved.One suggestion would be to target coastal shelf sitesproductive settings with high organic-carbon burial and shallow waters reducing dissolution of foraminiferal calcite.
Spatially, prior to ~125 Ma there are no data from the palaeo-Pacific Ocean, which constituted more than half of the world ocean during the Cretaceous (e.g., Hay et al., 1999;Scotese, 2004).The ability to sample particularly Early Cretaceous oceans is in part limited by the loss of the ancient oceanic record via subduction.Geographical coverage improves for the middle-Late Cretaceous relative to the Early Cretaceous, although there is a paucity of SST data from higher palaeolatitudes, particularly during the Cenomanian to Coniacian and during the early to middle Campanian, compounded by lack of TEX 86 data from high-palaeolatitude locations in general.Moreover, the data from the higher palaeolatitudes typically reflect conditions between ± 48-54°palaeolatitude in the Southern Hemisphere, and are therefore not wholly representative of Equator-Pole temperature trends.
These palaeoclimate data 'gaps' are not restricted to SST proxy data; similar to planktonic foraminifera, there is a lack of benthic foraminiferal δ 18 O data for the Early Cretaceous (Fig. 9d; Friedrich et al., 2012), and a lack of high-resolution pCO 2 proxy estimates, especially for the long quiet magnetic period in the middle Cretaceous, prevents a more comprehensive understanding of Cretaceous-CO 2 climate linkage (Li et al., 2014;Wang et al., 2014).A central motivation to improving our understanding of Cretaceous climate is to understand how sensitive the Earth's climate may be to much higher pCO 2 (e.g., PALEOSENS Project Members, 2012).Royer et al. (2012) estimate a Cretaceous Earth System Sensitivity of ~3 °C but up to 6 °C for the late Cenomanian-Coniacian (~95-85 Ma) based on available palaeo-reconstructions of pCO 2 and temperature circa 2012.Our compilation adds further constraints on (global) SST evolution during the Cretaceous.Ultimately, however, further studies are needed to elucidate the magnitude and variability of pCO 2, e.g., understanding episodic cyclic changes in pCO 2 during the middle-late Early Cretaceous (Li et al., 2014), as well as other second-order controls on climate such as palaeogeography.

Proxy data quality and interpretation
A recommendation regarding the quality of TEX 86 data is that workers report BIT values and, perhaps more importantly for Cretaceous sediments, the maturity of the samples (biomarker sterane and hopane ratios; Schouten et al., 2004) alongside TEX 86 values and individual GDGT relative abundances.Regarding the quality of Cretaceous planktonic δ 18 O data, it would be beneficial to have more values from exclusively 'glassy' foraminifera.In addition, the development of standardized criteria for assessing foraminiferal preservation and better reporting of such preservation for each sample/sample interval would improve comparisons drawing on data generated from a variety of localities by different workers (i.e. this study).In terms of improving proxy confidence, the application of additional approaches like clumped-isotope palaeothermometry would provide independent constraints on SST and the oxygen-isotope content of seawater.Both the TEX 86 and δ 18 O proxies suffer from extrapolation beyond the modern calibration range, while the clumped isotopes approach does not, although such measurements present different challenges, e.g., the requirement for large amount of sample (carbonate) material (e.g., Spencer and Kim, 2015).In addition, this approach will not resolve other questions important for interpretation of δ 18 O data, e.g., planktonic foraminiferal ecology.

Future climate modelling efforts
Our new Cretaceous SST synthesis provides an up-to-date target for future modelling studies investigating the mechanics of Cretaceous warmth.In particular, modelling of Cretaceous climate could help disentangle the relative importance of CO 2 and palaeogeography (and resulting changes in atmospheric and ocean circulation) for Cretaceous surface-ocean conditions.For example, climate modelling can be used to test hypotheses about the nature of Cretaceous latitudinal sea-surface temperature gradients under higher pCO 2 conditions akin to data-model comparisons already undertaken for intervals of the Late Cretaceous (Tabor et al., 2016;Petersen et al., 2016).Similarly, our compiled SSTs can provide climate constraints when modelling other Cretaceous climate forcing factors including the role of clouds (Sloan and Pollard, 1998;Kirk-Davidoff et al., 2002;Abbot et al., 2012) and biological cloud feedbacks (Kump Pollard, 2008;Upchurch et al., 2015), variations in the mixing ratios of non-CO 2 greenhouse gases (Beerling et al., 2011), long-term changes in the percentage of atmospheric oxygen (Poulsen et al., 2015), the possibility of multiple ocean steady states (Poulsen and Zhou, 2013) and polar forest vegetation-induced warming (Otto-Bliesner and Upchurch, 1997;Zhou et al., 2012) throughout the entire Cretaceous.Modelling can also be employed to explore the influence of changes in ocean-basin configuration (e.g., Poulsen et al., 2001;Zhou et al., 2008;Lunt et al., 2016), quantification of seawater δ 18 O gradients (e.g., Roche et al., 2006;Zhou et al., 2008) with the inclusion of isotopic tracers in GCMs (cf.Speelman et al., 2010) and/or surface and intermediate-to deep-water circulation on Cretaceous SSTs.Some such studies have already been undertaken for certain Cretaceous time intervals, e.g., short-term, ~2.5 Myr, biogeochemical modelling of the carbon cycle in the Late Aptian (McAnena et al., 2013). Donnadieuet al. (2016) modelled ocean circulation modes during two Late Cretaceous time slices in order to assess the role of changes in major continental configuration, comparing simulated ocean dynamics with existing neodymium-isotope data.Changes in Cretaceous ocean chemistry, in particular neodymium-isotope patterns (e.g., Pucéat et al., 2005;Robinson et al., 2010;Le Houedec et al., 2012;Robinson and Vance, 2012;Jung et al., 2013;Murphy and Thomas, 2013;Voigt et al., 2013;Zheng et al., 2013) indicate some interesting shifts in ocean circulation modes in the Late Cretaceous that may be important for SSTs.Finally, it is likely that some of the trends in our SST proxy data compilation are due to sampling sites moving in latitude as consequence of sea-floor spreading and continental drift, e.g., northward movement of the Pacific plate (Berger and Winterer, 1974).Future analyses might include a systematic point-by-point comparison with climate model simulations to quantify this effect, as well as longitudinal SST variations, which are obscured in gradient plots.

Summary
Using published SST proxy data, planktonic foraminiferal oxygen isotopes and GDGT distributions, we have generated a SST compilation for the entire Cretaceous Period.Our compilation uses SSTs recalculated from raw data, allowing examination of the sensitivity of each proxy to the calculation method, including choice of calibration, and places all data on a common timescale.In addition, we have investigated secondary controls on Cretaceous GDGT distributions through application of a range of GDGT indices -BIT, MI, f Cren′/ Cren′ + Cren , and ΔRI, compared together with modern GDGT distributions.After screening the raw GDGT data for problematic samples and considering the impacts of preservation and other nontemperature influences on Cretaceous δ 18 O pl records, all robust data were then used to generate a Cretaceous SST compilation.Overall, our compilation shows many stratigraphic similarities with other records of Cretaceous climate change, including benthic foraminiferal δ 18 O records, with both SST proxies indicating maximum warmth (SSTs > 30 °C at low and lower mid-latitudes) in the Cenomanian-Turonian interval (97-90 Ma).Similarly, both δ 18 O pl -and TEX 86 -SST estimates indicate prolonged cooling of the surface ocean and possible changes in ocean heat transport in the Late Cretaceous through the Coniacian-Santonian to the end-Maastrichtian interval.

Fig. 1 .
Fig. 1.Map showing the modern-day locations of sites with published Cretaceous (a) GDGT data and (b) planktonic δ 18 O data compiled in this paper.Sites are colour-coded by approximate Cretaceous palaeolatitude, estimated using a palaeorotational model provided by Getech Plc.Abbreviations: Bass R. = Bass River; Brazos R. = Brazos River; Sp. = Speeton, Als.= Alstätte; Moor.= Moorberg; Cis.= Cismon; Ader.= Aderet 1 borehole; PAMA = PAMA Quarry; Wun.= Wunstorf core; M. Farm 1 = Meirs Farm 1. (For interpretation of the references to colour in this figure legend, the reader is referred to the web version of this article.) observed a strong positive relationship between water depth and surface-sediment TEX 86 H values (and also suspended particulate matter TEX 86 H values) along the Portuguese continental margin.Phylogenetic analyses revealed that Thaumarchaeota populations identified at 1 m and 50 m water depth were different from those residing at 200 m and 1000 m water depth, leading to the suggestion that high surface-sediment TEX 86 H values could be due to the increasing contribution of isoGDGTs from the deepwater population of Thaumarchaeota.Together, these studies suggest (1) a water-depth control on TEX 86 H valuesand by extrapolation, other TEX-based indicesin certain settings

Fig. 2 .
Fig. 2. (a) Compiled published raw TEX 86 values for the Cretaceous (n = 1143).(b) Compiled published raw TEX 86 values for the Cretaceous after quality screening filtering for potentially problematic samples (n = 993); exclusion criteria are summarized in Table2.For individual datasets where only a few datapoints exist and/or the data span a very narrow temporal window (e.g., Brazos River, Meirs Farm 1 and FL-533), or age uncertainties are large (e.g., ODP 692, ODP 693 and DSDP 249) only the mean is plotted with full ranges of the datasets represented by horizontal and vertical bars.For Meirs Farm 1 and Brazos River horizontal bars are obscured since they lie within the areas of the respective mean datapoints.

Fig. 3 .
Fig. 3. Distribution of all (a) compiled published raw Cretaceous TEX 86 data (n = 1143) and (b) modern coretop TEX 86 data (n = 426; Kim et al., 2010).Cretaceous TEX 86 values are commonly high, > 0.7, in comparison to modern TEX 86 values and also exceed the maximum values recorded in the modern core-top data, > 0.9.

3. 4 .
Critical evaluation of planktonic δ 18 O data 3.4.1.Secondary controls on planktonic foraminifer δ 18 O-SSTs in the Cretaceous Although the primary controls on δ 18 O pl values are temperature and ice volume, several secondary factors influence Cretaceous δ 18 O pl values including depth of foraminiferal habitat, seasonality, test preservation, changes in surface-water δ 18 O (δ 18 O sw ) due to variations in the evaporation-precipitation balance, and carbonate ion concentration.

Fig. 6 .
Fig. 6.Compiled δ 18 O pl temperature reconstructions for the Cretaceous Period (n = 3789).(a) Published raw δ 18 O pl values.(b) δ 18 O pl -derived sea-surface temperature estimates using the palaeotemperature equation (1) of Bemis et al. (1998) and assuming ice-free conditions, (δw = −1‰ (VSMOW);Shackleton and Kennett, 1975;Hut, 1987).(c) Same as previous but with the addition of a δ 18 O sw correction for changes in palaeolatitude(Zachos et al., 1994).Estimates of palaeolatitude are derived using a palaeorotational model provided by Getech Plc.In the case of site TDP, all data pertain to one sample and so instead the mean is plotted with full ranges of the dataset represented by horizontal and vertical bars.In panels (b) and (c) the black bar represents the paleotemperature calibration uncertainty ( ± 0.7 °C,Bemis et al., 1998).

Fig. 8 .
Fig. 8. Compiled TEX 86 (n = 993) and δ 18 O pl (n = 3789) SST estimates for the Cretaceous Period derived from (a) TEX 86 H and δ 18 O pl , (b) TEX 86-Linear and δ 18 O pl .δ 18 O pl -SST estimates include a δ 18 O sw correction for changes in palaeolatitude (Zachos et al., 1994).The red horizontal lines indicate the upper limits of the modern TEX 86 H -and BAY GlobR -SST calibration datasets.The black bars in each panel represent the corresponding paleotemperature calibration uncertainties.As in Figs. 2 and 6, for individual datasets where only a few datapoints exist and/or data span a very narrow temporal window or age uncertainties are large, or all data pertain to the same sample, only the mean is plotted with full ranges of the datasets represented by horizontal and vertical bars.(For interpretation of the references to colour in this figure legend, the reader is referred to the web version of this article.)

Table 1
Sites from which Cretaceous raw GDGT data/TEX 86 indices were compiled.

Table 2
Summary of TEX 86 sample exclusion criteria.

Table 3
Sites from which Cretaceous planktonic δ 18 O data were compiled.
. Upper Aptian shells exhibit moderate to good preservation*.
indicates that warmest Cretaceous surfaceocean palaeotemperatures (> 35 °C) occurred from the late Cenomanian to Turonian in the equatorial Atlantic (δ 18 O pl & TEX 86 H / TEX 86-Linear ; ODP Sites 1258, 1259 and 1260, DSDP Sites 367 and 603), (Jenkyns et al., 2004)o Late Cretaceous are generated from low-palaeolatitude TEX 86 and δ 18 O pl data minus higher palaeolatitude δ 18 O pl proxy data, with the exception of a few TEX 86 values for the Arctic Ocean(Jenkyns et al., 2004).This comparison implies that a low-palaeolatitude minus higher palaeolatitude TEX 86 gradient would tend to give lower ΔSSTs compared with a TEX 86 and δ 18 O pl minus δ 18 O pl gradient, due to δ 18 O pl yielding cooler and more variable SSTs.Linear and δ 18 O pl minus higher mid-latitude δ 18 O pl ).This trend is almost identical if instead only low-latitude δ 18 O pl SSTs minus higher mid-latitude δ 18 O pl SSTs are considered (Supplementary Fig.3).For the early-mid Campanian there is another gap in our ΔSST reconstruction arising from a lack of higher latitude SST data.By the late Campanian the latitudinal SST gradient was ~16-18 °C (low-latitude TEX 86 H / TEX 86-Linear and δ 18 O pl minus higher mid-latitude δ 18 O pl ) and gradually strengthened as a result of greater cooling at higher latitudes, reaching a maximum, ~19-21 °C (low-latitude TEX 86 H /TEX 86-Linear and δ 18 O pl minus higher mid-latitude δ 18 O pl ), in the late Maastrichtian.
. Green shaded region represents the Cretaceous Thermal Maximum.(For interpretation of the references to colour in this figure legend, the reader is referred to the web version of this article.)