Ocean temperatures through the Phanerozoic reassessed

The oxygen isotope compositions of carbonate and phosphatic fossils hold the key to understanding Earth-system evolution during the last 500 million years. Unfortunately, the validity and interpretation of this record remain unsettled. Our comprehensive compilation of Phanerozoic δ18O data for carbonate and phosphate fossils and microfossils (totaling 22,332 and 4615 analyses, respectively) shows rapid shifts best explained by temperature change. In calculating paleotemperatures, we apply a constant hydrosphere δ18O, correct seawater δ18O for ice volume and paleolatitude, and correct belemnite δ18O values for 18O enrichment. Similar paleotemperature trends for carbonates and phosphates confirm retention of original isotopic signatures. Average low-latitude (30° S–30° N) paleotemperatures for shallow environments decline from 42.0 ± 3.1 °C in the Early-to-Middle Ordovician to 35.6 ± 2.4 °C for the Late Ordovician through the Devonian, then fluctuate around 25.1 ± 3.5 °C from the Mississippian to today. The Early Triassic and Middle Cretaceous stand out as hothouse intervals. Correlations between atmospheric CO2 forcing and paleotemperature support CO2’s role as a climate driver in the Paleozoic.

The 18 O/ 16 O ratios of biogenic carbonate and apatite are an essential proxy for paleotemperature and seawater δ 18 O. For decades, however, Earth scientists have debated the climate implications of the Phanerozoic 18 O/ 16 O record (reported as δ 18 O), which shows increasing values with decreasing age (e.g. Refs. [1][2][3][4]. The increase in biomineral δ 18 O toward the present has been interpreted three ways: (1) cooling of the Earth's oceans through time (e.g. Ref. 5 ), (2) evolving crustal cycling manifested in an increase in seawater δ 18 O (δ 18 O sw ; e.g. Refs. 2,3 ), and (3) progressive sample diagenesis with age (e.g. Ref. 6 ). Correct interpretation of this record is fundamental to understanding climate drivers, the limits of Earth's climate, and importantly, the evolution and thermal limits of metazoan life on geologic timescales.
To characterize and evaluate the Phanerozoic trend in biomineral δ 18 O and sea surface temperature (SST), we present a compilation of low-latitude (30° S-30° N) δ 18 O paleotemperatures for Phanerozoic carbonate and phosphate fossils and microfossils that, in contrast to earlier published records, corrects for the influence of ice volume and salinity (via paleolatitude relationships) on the δ 18 O sw of surface waters. Furthermore, for the first time we consider the impact of fractionation differences in belemnites 7,8 to develop a 500-million-year paleo temperature curve that resolves the problem of anomalously cold paleotemperatures previously obtained for intervals in the Jurassic and Cretaceous. The compiled data reveal similar trends for carbonate and phosphate δ 18 O, providing evidence for signal preservation, and argue for extreme warmth in the early Paleozoic (30 to > 40 °C), comparable to that experienced at the end-Permian and earliest Triassic 9,10 . are not shown. These are almost exclusively Paleozoic brachiopod shells that fail quality control (QC) standards (i.e., not "select"). Locfit regression lines (α = 0.05) for tropical and subtropical samples are shown with ± 95% confidence limits (CL). ForamP = planktonic foraminifera, Biv-gast = bivalves and gastropods, Belem = belemnites, Brach = brachiopods, and select = data that have passed QC standards. Belemnite values were adjusted by −1.5‰ to correct for 18 O enrichment. Letters in the bar between figures refers to geologic periods Cambrian, Ordovician, Silurian, Devonian, Mississippian, Pennsylvanian, Permian, Triassic, Jurassic, Cretaceous, Paleogene, and Neogene. Some symbols in key are enlarged for clarity.  14,15 . Paleotemperature differences also occur because of differences in the distribution of sample ages. Whereas brachiopods provide a robust record of the Hirnantian Cool Event, conodont data are essentially absent. In contrast, the hothouse at the end-Permian and Early Triassic is well represented in conodont data but not in brachiopod data 9,10 . Simple linear regression of carbonate and phosphate δ 18 O paleotemperatures averaged for Paleozoic and Triassic stages (Fig. 2B, Supplementary Tables 1 and S1) yields an equation: The strong correlation (R 2 = 0.6615) and slope near 1 are evidence that these materials retained their original isotopic signature. The ~ 4 °C average offset from the 1:1 line suggests higher paleotemperatures for the phosphate samples. This difference could reflect differences in aridity and the spatial and temporal distributions of specimens as mentioned above (e.g. Ref. 18 ). Another contributing factor may be depth habitat. At some locations, nektonic conodonts may have lived shallower in the water column than benthic brachiopods. Lastly, differences in paleotemperature may indicate uncertainties in the paleotemperature equations. Applying the equation of Friedman and O'Neil 19 for calcite decreases the temperature difference between Paleozoic stages from 3.6 to 2.0 °C. However, we use the Kim and O'Neil 12 equation because it is more widely applied by the paleoclimate community. Combining the results for conodonts and brachiopods provides a continuous and comprehensive paleotemperature record for low latitudes (30° S-30° N; Fig. 4). Late Cambrian (Furongian) and Early to Middle Ordovician low-latitudes experienced the highest SSTs, 47 °C and between 35 and 47 °C, respectively. Late Ordovician to Devonian low-latitude SSTs were considerably cooler (32-40 °C). Cooler still were Carboniferous to modern low-latitude SSTs, which varied between 19 and 35 °C with relatively cool "icehouse" temperatures observed in the Carboniferous to Middle Permian (19-24 °C), Middle Jurassic to Early Cretaceous (18-27 °C), and late Eocene to modern times (23-25 °C). The late Permian to Middle Jurassic and Late Cretaceous to Eocene are recognized as warm climatic intervals (25-30 °C) with the Permian-Triassic transition, Early Triassic, and Cenomanian-Turonian standing out as hothouse intervals. The Paleocene-Eocene Thermal Maximum (PETM), another very warm period, is not well represented in our oxygen isotope dataset due to insufficient data for shallow-dwelling organisms.

Discussion
The decrease in marine fossil δ 18 O with age, the trend expected with oxygen isotope exchange with meteoric water, has led some to question the preservation and efficacy of the original marine δ 18 O signal, especially with regard to the Paleozoic record. We address these concerns by restricting Paleozoic samples to the best-preserved material, brachiopod shell calcite and conodont apatite. Phosphate oxygen is more tightly bound in the crystal lattice than is calcite oxygen; in fact, apatite phosphate survives dissolution and reprecipitation as trisilverphosphate during sample preparation without altering its oxygen isotope composition. That both minerals yield similar δ 18 25 and for Devonian conodont apatite experiencing minimal versus extensive heating 26 .
As with previous studies, our δ 18 O results show an increase from very low late Cambrian through Middle Ordovician values toward higher δ 18 O values in the modern (Fig. 5). Not surprisingly, our paleotemperatures show many of the same relative changes as in previous studies, an expectation considering the overlap in data sources; however, substantial differences occur in absolute temperature with higher Early Paleozoic temperatures in our study and much lower Phanerozoic temperatures in other studies. These offsets reflect differences in (1) sample material (carbonates versus combined carbonates and phosphates), (2) screening approaches, (3) methods of estimating seawater δ 18 O, and (4) paleotemperature equations 1,4,[27][28][29][30] .
We start the discussion with Veizer and Prokoph's 4 seminal work because several studies adopt their data and interpretation (e.g., Refs. 22,33 ). With regard to sample material, the Veizer and Prokoph 4 δ 18 O curve is based on data only from brachiopods and planktonic foraminifera and excludes data from phosphates (e.g., conodonts), belemnites, bivalves, and gastropods. Importantly, planktonic foraminifera can be readily recrystallized on the cold sea floor; specimens not characterized as glassy or excellently preserved can yield δ 18 29 . This curve yields lower temperatures in the Paleozoic, early Mesozoic, and Cenozoic compared with ours (Fig. 5). These authors make use of phosphate δ 18 29 also includes samples for higher latitudes (up to ± 40°), which could lower paleotemperatures. However, in StabisoDB, phosphate samples from 30 to 40° paleolatitude only represent 7.5% of the samples between 0 and 40° paleolatitude, so this effect should be minor. Scotese et al. 's 30 Phanerozoic temperature curve shows less extreme temperatures in the Early Paleozoic and late Cenozoic (Fig. 5). The curve uses a combination of isotopic data 29 and paleotemperatures based on paleo-Köppen climate belts calibrated with modern Köppen belt temperature relations. While novel, this hybrid paleotemperature curve implicitly assumes modern temperature relations and thus may inject a uniformitarian bias into quantification of Earth's temperature history. As discussed earlier, our Phanerozoic curve is the first to correct δ 18 O sw for paleolatitude and for 18 O-enrichment in belemnites, which accounts for warmer proxy temperatures for the Jurassic compared with other isotopic studies.
Accepting that the Phanerozoic δ 18 O trend is not an artifact of diagenesis, the debate as to whether the trend reflects δ 18 O sw or temperature change distills to two endmember assumptions: (1) no long-term trend in seawater δ 18 O (e.g. Ref. 36 ) and (2) no long-term trend in low-latitude SST 3,4,27 . In theory, the oxygen isotope history of seawater can be calculated from the rates and temperatures of oxygen exchange between the mantle, crustal reservoirs, and the ocean. Mass balance models have been used to argue for near-constant (e.g. Refs. [36][37][38] 40 ). One limitation of crustal exchange models is the large uncertainty in hydrothermal fluid fluxes. These fluxes are estimated based on the difference between modeled and measured (conductive) ocean heat flux 41 . However, uncertainty in the temperature of off-ridge hydrothermal fluids can lead to greater than ± 50% uncertainty in hydrothermal flux and large uncertainty in the mineral-water 18 O fractionation. A recent hypothesis is that Snowball-Earth sequestration of glacial ice might have resulted in 18 O-enriched residual seawater that would re-equilibrate with ocean crust, lowering δ 18 O sw toward 0‰, followed by melting of 18 ; Fig. 4). Even higher TEX 86 SSTs of 40 to 45 °C are reconstructed using the TEX 86 BAYSPAR calibration 16,53 . Thus, low-latitude SSTs of 35 °C and warmer are not unique to the Early Paleozoic but occurred also during Mesozoic and Cenozoic warm intervals. Furthermore, these temperatures are close to tropical SST projected for the year 2100 if modern tropical seas (25-30 °C) warm by more than 4 °C 54,55 , as predicted by the RCP 8.5 scenario.
Comparing low-latitude temperature to seafloor accretion rate 22 allows us to examine the link between plate tectonics and climate. Seafloor accretion rate correlates significantly with low latitude temperature (R 2 = 0.47, p < 0.0001; Supplementary Fig. S1, Supplementary Table S2), with warmer temperatures associated with faster spreading rates. This linkage is attributed to high rates of volcanic CO 2 degassing with higher rates of seafloor spreading and subduction 30 .
Isotopic temperatures of > 40 °C for the late Cambrian and Early Ordovician, which extend beyond the temperature tolerances of most modern multi-cellular Eukarya (e.g. Ref. 56 ), present a conundrum. Is our understanding of the physiology and behavior of early animals incomplete? If Early Ordovician fauna were limited only to taxa able to tolerate unusually high temperatures, this would further strengthen Ordovician cooling as an explanation for the dramatic increase in Ordovician diversity, the Great Ordovician Diversification Event (GOBE) 57 . Furthermore, such cooling would raise the solubility of oxygen in seawater 29 and, along with a posited increase in atmospheric oxygen levels, permit greater metabolic activity and predation 58 , leading to paleoecological reorganization 59 .
Our results can be used to examine the sensitivity of low-latitude temperatures to changes in pCO 2 (lowlatitude Earth system sensitivity or ESS) Overall, paleotemperature exhibits a significant correlation (R 2 = 0.234, p = 0.0014) with pCO 2 doubling based on the proxy record of Foster et al. 21 ; Supplementary Table S3; based on 10-myr Locfit averages). This relationship indicates a role of pCO 2 in controlling low-latitude Phanerozoic SST. For Phanerozoic climate, changes in solar radiation must be considered in addition to the radiative forcing controlled by pCO 2 . To examine the relationship between changes in low-latitude temperature (ΔT LL ) and changes in radiative forcing of pCO 2 and solar radiation (ΔS [CO2,SOL] ), we convert pCO 2 doubling to radiative forcing by multiplying by 3.7 W m −2 , and correct for increasing solar radiation with time using the equation: ΔF sol = −23 8 W m −2 × 0.0000665 age (myr) 60 . ΔS [CO2,SOL] estimates yield significant correlations (p < 0.0001) with ΔT LL for the Paleozoic but not for the Mesozoic and Cenozoic (Fig. 6). Royer 61 also found a poor relationship between surface temperature and combined CO 2 and solar forcing for the Cenozoic and Mesozoic, noting that temperature data for the Mesozoic are sparse. Though not seen in our results based on limited planktonic foraminifera and macrofossil data, a strong relationship between Cenozoic temperatures and pCO 2 is seen in studies based on benthic foraminiferal δ 18  www.nature.com/scientificreports/ DMII regression was chosen because of comparable uncertainty in both X and Y. This equation yields an Earth system sensitivity (ESS) value for Paleozoic low latitudes of 2.9 K W −1 m 2 or 10.7 K per CO 2 doubling. This value is high compared with the range for 35-150 Ma based on an ensemble of climate model simulations (3.5 and 5.5 K 63 ), especially considering that (1) the 0°-30° latitude band accounts for only half of Earth's surface and (2) SST change will underestimate change in global surface temperature, especially during icehouse climate 61,64 . On the other hand, the value is within those calculated for discrete time intervals in the Pliocene 61 . In comparing radiative forcing and temperature change, other studies of Earth system sensitivity have used simple linear regression (SLR), which considers only uncertainty in Y 61,65 . SLR of our results yields the following equation for the Paleozoic: This equation provides an ESS value of 1.7 W −1 m 2 or 6.3 K per CO 2 doubling, similar to climate sensitivities determined for glaciated time 64 . Justification for using simple linear regression instead of Model II regression is discussed in Smith 66 and centers on the objectives of the study. Since our objective is to "define some mutual, codependent "law" underlying the interaction between X [ΔpCO 2 radiative forcing] and Y [ΔSST]", and since the slope "will be used to interpret the pattern of change" (Smith 66 , p. 482), we favor Model II regression. On the other hand, Smith 66 notes that "when an equation is used for prediction", SLR is the "method of choice". More detailed examination of choice of regression model is beyond the scope of this paper but clearly merits future consideration.
Note that this treatment does not account for changes in paleogeography, sea level, land ice area, vegetation, non-CO 2 greenhouse gases, and aerosols, some of which also serve as feedback mechanisms; however, it does provide an estimate of ESS for the highest sustained CO 2 levels in the past 420 myr (> 2000 ppm) and the worst-case scenario levels for a couple of centuries into the future 55 . Lastly, our finding that ESS was high in the Paleozoic compared with the Cenozoic supports studies suggesting higher ESS with higher CO 2 levels (e.g. Ref. 67 ).

Materials and methods
Details regarding samples and methods appear in Ref. 1 . The compilation builds upon previous efforts (e.g. Refs. 4,27 ) and focuses on carbonate and phosphate fossils and microfossils that (1) are widely distributed in the sedimentary record, (2) are precipitated with quantitative δ 18 O fractionation relative to temperature, and (3) exhibit excellent preservation. Samples include mollusks, brachiopods, planktonic foraminifera, fish teeth, and conodonts.  www.nature.com/scientificreports/ The late Cambrian through Triassic record is based on brachiopod calcite and conodont phosphate. Thick brachiopods from cratons tend to show the best preservation and least scatter in their δ 18 O values (e.g. Refs. 15,68,69 ). Targeting best preserved shells with petrographic and cathodoluminescence microscopy, combined with analyses of microsamples (< 100 µg), further reduces variability and allows for multiple analyses from a single shell. While our compilation includes all data, regressions for temporal trends and paleotemperature estimates only consider data from brachiopod shell that (1) is non-luminescent, (2) contains manganese contents < 250 ppm, or (3) thick-shelled and from areas known for excellent preservation (e.g., Moscow Basin; see Ref. 1 for additional details). Biogenic apatite is less prone to diagenetic overprinting; however, the unknown habitat of conodonts may represent an uncertainty for interpretation of δ 18 O values. While brachiopods are benthic organisms, conodonts were active swimmers that could have lived in warm surface waters or deeper and thus colder parts of the water column. The comparison of oxygen isotope values of conodont taxa from sediments of different water depths gave equivocal results (e.g. Refs. 70,71 ).
Belemnite rostra, calcite deposits in the cephalopod's posterior, are the most common material analyzed from Jurassic and Cretaceous sediments. These fossils typically have δ 18 O values higher than those of co-occurring bivalves, confounding paleotemperature studies. However, recent clumped isotope studies have revealed that belemnites precipitated in warmer waters than their δ 18 O values indicate, prompting researchers to conclude that belemnite guards are precipitated in true equilibrium with seawater, and that δ 18 O paleotemperature relations for other biominerals and laboratory precipitates do not represent true equilibrium 7,8 . In our dataset, belemnites average 1.7 ± 0.5‰ (N = 19) and 1.1 ± 1.2‰ (N = 13) higher in δ 18 O compared with brachiopods and bivalves (Appendix 4). To account for this effect, we have applied a −1.5‰ correction to belemnite data as suggested by Vickers et al. 8 .
Planktonic foraminiferal data provide paleotemperatures for Cretaceous through Cenozoic climate. Because planktonic foraminiferal tests commonly recrystallize with burial on the sea floor 34 , we only use data for planktonic foraminifera that exhibit exceptional preservation (e.g., "glassy", "excellent" 16,34 70 and presented in detail in the papers from which the data are derived. Briefly, carbonates of 0.05 to several milligrams are acidified with concentrated phosphoric acid and the CO 2 evolved is analyzed on an isotope ratio mass spectrometer. Isotopic data are reported in delta (δ) notation and reported versus PDB (Peedee belemnite) or VPDB (Vienna PDB). The latter refers to calibration to PDB using the NBS-19 calcite standard (δ 18 O = − 2.20‰ versus PDB 73,74 or the new carbonate standard, IAEA-603 (δ 18 O = − 2.37‰). The precision for oxygen isotope analyses of CaCO 3 is typically ± 0.05 to 0.10‰, which equates to roughly ± 0.25 to ± 0.5 °C at 25 °C 12 . Oxygen isotope analyses of biogenic apatite are either measured by (1) TC-EA IRMS (thermally coupled elemental analyzer-isotope ratio mass spectrometry) on trisilverphosphate precipitated after dissolving calcium fluorapatite or (2) in situ by secondary ion mass spectrometry (SIMS). Whereas phosphate-bound oxygen is analyzed by TC-EA IRMS, total oxygen including phosphate-, carbonate-and hydroxyl-bound oxygen is measured by SIMS with the δ 18 O offset between these methodologies not well constrained. We applied a correction of −0.6‰ to all SIMS δ 18 Tables S4-S6). Ice volumes through time are binned into simple categories of ice-free, low, moderate, and high based on studies of glacial sediments and sea level (e.g. Refs. [76][77][78][79]. For the δ 18 O of ice, we assume the δ 18 O values for the West Antarctica (−41‰) and Greenland ice sheets (−34‰) for "moderate" and "low" ice volumes respectively (Supplementary Table S4-S6). The calculated values for mean δ 18 O sw range from −1.08‰ for the ice-free state to 0.45‰ for high ice volume (Pleistocene average; Supplementary Table S7). Lastly, δ 18 O sw was averaged for 1-myr steps using a 2-myr window to smooth the impact of assigned ice volume changes.
Paleolatitudinal correction. Paleolatitudes were reconstructed using the GPlates software 80 Table S8) derived from gridded data modeled in LeGrande and Schmidt 83 . Using data for the Southern Hemisphere, which are less influenced by landmasses than Northern Hemisphere data, yields the relationship for 0-60° latitude: For hothouse climates, we generate the δ 18 O sw -latitude relation using the isotope-enabled ocean-atmosphere general circulation model (GCM) results of Ref. 82 (Supplementary Table S8) for the early Paleogene. The modeled δ 18 O sw values for 0-60° latitude in the southern hemisphere (from Fig. 1 of Ref. 82   www.nature.com/scientificreports/ where δ 18 O sw,iceV is the ice volume correction for global seawater and 0.45‰ is the icehouse endmember (average between glacial and interglacial states). We bin our data into the following climate zones by paleolatitude: tropical (± 10°), tropical-subtropical (10°-30°), temperate (30°-50°), and subpolar-polar (50°-90°) using the maps from Scotese and Wright 81 . These bins were selected based on the temperature gradients for the latest Cretaceous through Recent reported in Zhang et al. 84 . Over the time interval studied, spanning greenhouse and icehouse climates, paleotemperatures within 10° N or S are invariant with paleolatitude. Inflections in latitudinal temperature gradients at 30° and 50° define boundaries for the next two bins. Lastly, stage-averaged paleotemperatures for the ± 10° and 10°-30° bins were found to be statistically similar (ΔT(10°-30° minus 0°-10°) (°C) = −1.0 ± 2.1(2SE) °C) (N = 25) for carbonates and −2.0 ± 2.2 (N = 25) for phosphates) and thus were combined. At latitudes higher than 30°, data become sparse. Further, the increased latitudinal temperature gradient at higher latitudes along with greater influence of 18 O-depleted fresh water 83 increases the uncertainty in paleotemperature determinations. Lastly, sample ages are based on the GTS2020 timescale 85 .

Data availability
All data and Locfit regression tables used in this study are available in the Supplementary tables and auxiliary data files. The data are also available on the StabisoDB online database (http:// stabi soDB. org).