Atmospheric CO2 estimates for the Miocene to Pleistocene based on foraminiferal δ11B at Ocean Drilling Program Sites 806 and 807 in the Western Equatorial Pacific

Constraints on the evolution of atmospheric CO2 levels throughout Earth’s history are foundational to our understanding of past variations in climate. Despite considerable effort, records vary in their temporal and spatial coverage and estimates of past CO2 levels do not always converge, and therefore new records and proxies are valuable. Here we reconstruct atmospheric CO2 values across major climate transitions over the past 16 million years using the boron isotopic composition (δ11B) of planktic foraminifera from 89 samples obtained from two sites in the West Pacific Warm Pool, Ocean Drilling Program (ODP) Sites 806 and 807, measured using high-precision multi-collector inductively coupled plasma mass spectrometry. We compare our results to published data from ODP Site 872, also in the Western Equatorial Pacific, that goes back to 22 million years ago. These sites are in a region that today is near equilibrium with the atmosphere and are thought to have been in equilibrium with the atmosphere for the interval studied. We show that δ11B data from this region are consistent with other boron-based studies. The data show evidence for elevated pCO2 during the Middle Miocene and Early to Middle Pliocene, and reductions in pCO2 of ∼ 200 ppm during the Middle Miocene Climate Transition, ∼ 250 ppm during Pliocene Glacial Intensification and ∼ 50 ppm during the Mid-Pleistocene Climate Transition. During the MidPleistocene Transition there is a minimum pCO2 at marine isotopic stage (MIS) 30. Our results are consistent with a coupling between pCO2, temperature and ice sheet expansion from the Miocene to the late Quaternary. Highlights. In this study, we reconstruct atmospheric pCO2 using δ11B data from ODP Sites 806 and 807 and compare them with ice core data. We therefore apply the same framework to older samples from these sites to create a long-term pH and pCO2 reconstruction for the past 16 million years, including recalculating pCO2 for ODP Site 872 from 17 to 22 million years ago. We find that major increases in surface water pH and decreases in atmospheric pCO2 were associated with decreased temperature in the Western Equatorial Pacific and associated with major episodes of ice sheet expansion in the high latitudes, providing more robust quantitative constraints on the past coupling between pCO2, temperature and cryosphere stability.

spheric CO 2 and examine past relationships between CO 2 and temperature. Such data are not only critical for constraining Earth-system sensitivity (Lea, 2004;Lunt et al., 2010;Pagani et al., 2010;Hansen et al., 2012Hansen et al., , 2013Foster and Rohling, 2013;Schmittner et al., 2011;Tierney et al., 2020), but are also of broad interest for contextualizing the evolution of climate and geological systems throughout Earth's history (Tripati et al., 2011;Foster et al., 2017;Tripati and Darby, 2018). However, discrepancies between proxy reconstructions still exist, including for major climate transitions of the Cenozoic. In particular, there remains a pressing need for robust and higher-resolution atmospheric CO 2 records.
Each of the above proxy methods has sources of systematic errors that we do not attempt to exhaustively document as they have been discussed in-depth elsewhere (e.g., Pagani et al., 2005;Tripati et al., 2011;Guillermic et al., 2020). However, we note that significant developments in the boronbased proxies include improvements to the accuracy and precision of measurements using multi-collector inductively coupled mass spectrometry (MC-ICP-MS) compared to early work with negative thermal ionization mass spectrometry (N-TIMS), where there were large instrumental mass fractionations and challenges with laboratory intercomparison Farmer et al., 2016;Aggarwal and You, 2017). There was also the realization that temperature-dependent K D and B/Ca sensitivities reported from sediment trap, coretop and downcore studies (Yu et al., 2007;Tripati et al., 2009Tripati et al., , 2011Babila et al., 2010;Osborne et al., 2020) differ from inferences from foraminiferal culture experiments (Allen et al., 2011;Allen and Hönisch, 2012) and inorganic calcite (Mavromatis et al., 2015), which complicates the use of the B/Ca proxy, although this type of dis-crepancy has also been observed with other elemental proxies (e.g., Mg/Ca). Such differences may be due to differences in growth rates (Gabitov et al., 2014), ontogenetic changes, a correlation in the field between temperature and other hydrographic variables that obscure robust statistical determination of parameter relationships, culture conditions resulting in organisms being stressed, and/or other factors.
The marine CO 2 proxy that appears to be subject to the fewest systematic uncertainties, based on our current understanding, is the boron isotopic composition (δ 11 B) of planktic foraminifera as measured using MC-ICP-MS and N-TIMS . This proxy provides constraints on seawater pH, if temperature, salinity, seawater δ 11 B, and the appropriate mono-specific calibration between δ 11 B carbonate and δ 11 B borate are constrained (Pearson and Palmer, 2000;Sosdian et al., 2018;Raitzsch et al., 2018;Guillermic et al., 2020). Seawater pH can be used to calculate seawater pCO 2 if there are constraints on a second parameter of the carbonate system (e.g., alkalinity, DIC). Atmospheric pCO 2 can then be constrained if the site being examined is in airsea CO 2 equilibrium or if the disequilibrium is known and stable through time.
However, there are relatively few studies generating highprecision boron-based records over major climate transitions in the Cenozoic using recent analytical methods and that incorporate our current understanding of the proxy (e.g., Greenop et al., 2014;Martínez-Botí et al., 2015a;Dyez et al., 2018;Sosdian et al., 2018;de la Vega et al., 2020;Rae et al., 2021;Raitzsch et al., 2021). Furthermore, of the existing studies using boron-based proxies, an additional uncertainty frequently exists, namely the short time interval of study (e.g., emphasizing on a climate transition) (Martínez-Botí et al., 2015b;Chalk et al., 2017) and whether the study sites remain in air-sea CO 2 equilibrium (Martinez-Botí et al., 2015b). Moreover, although estimation of atmospheric pCO 2 from seawater pH using this proxy is relatively straightforward, reconstructions are still impacted by uncertainties, including the lack of robust constraints on a second parameter of the carbonate system and our limited understanding of secular variations in the δ 11 B of seawater (Tripati et al., 2011;Greenop et al., 2017;Sosdian et al., 2018;Rae et al., 2021). Therefore, to provide additional constraints on the evolution of atmospheric pCO 2 from the Miocene through the Pleistocene, we developed new records from the western tropical Pacific. We use foraminiferal δ 11 B and trace elements in the planktic foraminiferal species Trilobus sacculifer and Globigerinoides ruber to reconstruct past seawater pH and atmospheric CO 2 at Ocean Drilling Program (ODP) Sites 806 and 807 in the Western Equatorial Pacific (WEP) over the last 16 million years (Myr). The sites are located on the western border of the tropical Pacific Ocean, the largest open-ocean region on the globe, and the warmest open-ocean region at present. These two sites have been examined in other boron-based studies (Wara et al., 2003;Tripati et al., 2009Tripati et al., , 2011Shankle et al., 2021), as has the region more broadly (Pearson and Palmer, 2000;Sosdian et al., 2018), because it is understood to be in equilibrium with the atmosphere and have relative stable hydrography. The region experiences equatorial divergence but is not strongly affected by upwelling and has a current estimated annual air-sea CO 2 difference of +28 ppmv (Takahashi et al., 2014). The pre-industrial air-sea CO 2 difference is calculated to be +16 ppm (GLODAP database corrected from anthropogenic inputs), with a value of 298 ppm, compared to the ice core value of 282 ppm at 1.08 ka. This pCO 2 difference is similar to our pCO 2 uncertainty (an average of ∼ 17 ppm (2 SD) for the youngest samples). If trade winds were much stronger, and equatorial divergence greater, then this could have driven some disequilibrium in the past. However, a few lines of evidence suggest the region was in quasi-equilibrium in the past: (1) zonal temperatures are at a maximum in pre-industrial times and during the Pleistocene, and we are able to reconstruct atmospheric pCO 2 values from the ice cores, (2) and temperature proxies indicate the region is relatively stable with respect to temperature compared to other parts of the ocean and also indicate a weak and stable zonal temperature gradient during the Miocene and Pliocene which would support air-sea stable conditions and air-sea (dis-)equilibrium conditions (e.g., Nathan and Leckie, 2009;Zhang et al., 2014;Liu et al., 2019).
Thus, this study builds on prior low-resolution reconstructions for these sites (Wara et al., 2003;Tripati et al., 2009Tripati et al., , 2011Shankle et al., 2021), Site 872 in the tropical Pacific (Sosdian et al., 2018), and other published boron isotope work, to provide additional data to constrain past seawater pH and pCO 2 for the WEP using MC-ICP-MS, thereby providing a new perspective on reconstructing past atmospheric CO 2 via marine sediment archives. We explore various constraints on the second carbonate system parameter using a number of different scenarios, following on the systematic work done by Tripati et al. (2009Tripati et al. ( , 2011 for B/Ca. We interpret these data using recent constraints on seawater δ 11 B (Lemarchand et al., 2002;Raitzsch and Hönisch, 2013;Greenop et al., 2017). For temperature estimation, we utilize a multi-variable model for Mg/Ca correcting for salinity, pH and seawater Mg/Ca (Gray and Evans, 2019), that builds on prior work with clumped isotopes in planktic foraminifera for Site 806 and other WEP sites, demonstrating that for the Last Glacial Maximum to recent times, salinity-corrected Mg/Ca values are needed to yield convergent estimates of mixedlayer temperatures .

Materials and methods
Below we describe site locations, analytical methods used and principal figures. The supplemental methods section de-

Preservation
Microfossils in sediments at these sites, as with any sedimentary sequences, have the potential to be influenced by diagenesis. Despite evidence of authigenic carbonate formation, recent modeling work concluded that the influence of dissolution and reprecipitation at Sites 806 and 807 was relatively minor (Mitnick et al., 2018). Prior work has also found minimal impacts on the B/Ca ratio of Pliocene foraminifera from Site 806 (White and Ravelo, 2020) and on the Mg/Ca ratio of Miocene Dentoglobigerina altispera shells at Site 806 (Sosdian et al., 2020). The weight-to-shell ratio is commonly used to monitor dissolution, and the only published record at Site 806 for the Pliocene does not show a trend consistent with dissolution of T. sacculifer (Wara et al., 2005). We do note that while the "coccolith size-free dissolution" index reported in Si and Rosenthal (2019) indicates higher dissolution rates in the Miocene, their records were thought to be biased from changes in foraminifera assemblages as discussed in White and Ravelo (2020).
To further assess the potential impact of dissolution in our geochemical data, the weight-to-shell ratio was examined in our samples. The weight-to-shell data used to monitor dissolution does not exhibit any trend within the interval studied consistent with dissolution. Absolute weight-toshell is increasing in the Miocene, which is not consistent with dissolution influencing the record (Fig. 2e). Additionally, reconstructed pH and pCO 2 values also exhibit reasonable correspondence with the ice core data. Downcore δ 11 B values from Sites 806 and 807 are similar, despite evidence for higher authigenic carbonate at Site 807 relative to Site 806 (Mitnick et al., 2018). Further, despite different sedimentation rates, our δ 11 B and Mg/Ca results are consistent between Sites 806 and 807, and with data from Site 872 (Sos- dian et al., 2018), which implies that diagenesis is not a primary driver of the reconstructed trends. A comparison of raw data as well as derived parameters is shown in Figs. 2 and 7.

Age models
The age model for Site 806 from 0-1. 35 Ma is based on Medina-Elizalde and Lea (2005); calculated ages correspond well with ages from the Lisiecki and Raymo LR04 stack (Fig. 2a). The fourth polynomial regression-based biostratigraphy from Lear et al. (2015) was used for the rest of the record, following other work (Sosdian et al., 2020). Ages for Site 807 are based on published biostratigraphy (Berger et al., 1993) with additional constraints placed by Zhang et al. (2007) for the interval from 0-0.55 Ma. Benthic δ 18 O values from Sites 806 and 807 show good correspondence for the last 0.55 Myr, and the low-resolution benthic δ 18 O record for Site 806 (Lear et al., 2003(Lear et al., , 2015 is consistent with the stack from Lisiecki and Raymo (2005) for the period studied (Fig. 3).

Species and trace element cleaning
Samples were picked and cleaned to remove clays at UCLA (Los Angeles, CA) and the University of Western Brittany (Plouzané, France). A total of 50-100 foraminifera shells were picked from the 300-400 µm fraction size for T. sacculifer (without sacc) and from the 250-300 µm for G. ruber (white sensu stricto). Picked foraminifera were gently crushed, clays were removed, and they were checked for coarse-grained silicates. Samples were then cleaned using a full reductive and oxidative cleaning protocol following Barker et al. (2003). A final leach step with 0.001N HCl was done prior to dissolution in 1N HCl. Boron purification used a published microdistillation protocol (see Misra et al., 2014b;Guillermic et al., 2020, for more detailed methods).

Standards
Variations in B isotope ratios are expressed in conventional delta (δ) notation with δ 11 B values reported against the ref-   Lear et al. (2003Lear et al. ( , 2015. (b) pCO 2 values calculated from boron isotopes (colored symbols -this study) with data from the literature (open gray triangles -compilation B are data recalculated in Rae et al., 2021) and ice core pCO 2 (black line - Petit et al., 1999;Lüthi et al., 2008;Bereiter et al., 2015). (c) Cross plot for the last 0.8 Myr of pCO 2 δ 11 B from this study and pCO 2 ice core (from ice core compilation, Bereiter et al., 2015), grey line is a simple linear regression (p = 0.25, R 2 = 0.09), blue line is a Deming regression taking both x and y uncertainties into account (p = 0.25). Details of the regression parameters are in Table S6. Ice core CO 2 error was calculated based on 2 SD of reported values, and ±1 kyr for the age of sediment samples. Boron-based pCO 2 error is calculated based on error propagation described by Eq. (S17). Data compiled are from Foster (2008) Multiple analyses of external standards were performed to ensure data quality. For boron isotopic measurements, JC p−1 (Geological Survey of Japan, Tsukuba, Japan, Gutjahr et al., 2020) was used as a carbonate standard, and NEP, a Porites sp. coral from University of Western Australia and Australian National University was also used (McCulloch et al., 2014). A boron isotope liquid standard, ERM© AE121 (certified δ 11 B = 19.9 ± 0.6 ‰, SD), was used to monitor reproducibility and drift during each session (Vogl and Rosner, 2012;Misra et al., 2014b). For trace elements, external reproducibility was determined using the consistency standard Cam-Wuellerstorfi (University of Cambridge) (Misra et al., 2014a).

δ 11 B analyses
Samples measured for boron isotopes typically ranged in concentration from 10 (∼ 5 ng B) to 20 ppb B samples (∼ 10 ng B). Sensitivity was 10 mV/ ppb B (e.g., 100 mV for 10 ppb B) in wet plasma at 50 µL/min sample aspiration rate. The intensity of 11 B for a sample at 10 ppb B was typically 104 ± 15 mV (2 SD, typical session) and closely matched the 98 ± 6 mV (2 SD, typical session) of the standard. Procedural boron blanks ranged from 15 to 65 pg B (contributed to less than 1 % of the sample signal). The acid blank during analy-ses was measured at ≤ 1 mV on 11 B (which also is <1 % of the sample intensity), and no memory effect was seen within and across sessions.

Calculations
Detailed calculations can be found in the Supplement. Briefly, Mg/Ca was used to reconstruct sea surface temperature (SST) using the framework from Gray and Evans (2019) correcting for influences of pH, salinity and secular variation in seawater Mg/Ca. δ 11 B carbonate was corrected using an empirical δ 11 B carbonate weight-to-shell ratio relationship. δ 11 B borate was determined using species-dependent sensitivities of δ 11 B carbonate to δ 11 B borate (Guillermic et al., 2020). pH was calculated using δ 11 B borate with different scenarios of secular seawater δ 11 B changes (Lemarchand et al., 2002;Raitzsch and Hönisch, 2013;Greenop et al., 2017). pCO 2 was reconstructed using pH-based δ 11 B carbonate and different scenarios of alkalinity (Tyrrell and Zeebe, 2004;Ridgwell and Zeebe, 2005;Caves et al., 2016 andRae et al., 2021). Further details including equations are provided in the Supplement.

Geochemical results
Geochemical data used in this study are presented in Fig. 2. Mg/Ca data (Fig. 2c) are consistent with previously published Mg/Ca values for Site 806 on T. sacculifer (Wara et al., 2005;Tripati et al., 2009;Nathan and Leckie, 2009). Although the record we generated does not overlap with Site 872, they are 1 Myr apart (15.7 and 16.7 Ma); there is a good correspondence between our Mg/Ca data and the published Mg/Ca record from T. trilobus at Site 872 (Sosdian et al., 2018). Mg/Ca from a different species, D. altispera (Sosdian et al., 2020), is also plotted with an offset, for comparison.
Comparison with Site 872 data that are part of the compilation from Sosdian et al. (2018) shows that their δ 11 B data are in line with our dataset (Fig. 2b), and all sites examined in the WEP (Sites 806, 807, and 872) are above the lysocline (Shipboard Scientific Party, 1991). The δ 11 B data for T. sacculifer exhibit a significant increase (4.2 ‰) from the Miocene to the present. Figure 2b also compares the δ 11 B data used in this study with published data from other sites and shows that raw δ 11 B data for the WEP can be lower than values for other regions.

Reproducing pCO 2 from ice cores
We sought to assess whether there is evidence for air-sea equilibrium or disequilibrium in the WEP during the large amplitude late Pleistocene glacial-interglacial cycles, in order to validate our approach. We reconstructed pCO 2 for the last 800 kyr (n = 16, Fig. 3). For the last 800 kyr, reconstructed pCO 2 values for Sites 806 and 807 are in the range of ice cores (Fig. 3, Petit et al., 1999;Siegenthaler et al., 2005;Lüthi et al., 2008;compilation from Bereiter et al., 2015). The two critical diagnostics we used for method validation are (1) that the δ 11 B-based reconstruction of pCO 2 is consistent with ice core atmospheric CO 2 and (2) that the boron-based reconstruction empirically reproduces interglacial-glacial amplitudes from ice cores. Figure 3b shows that both of these criteria are met despite large scatter. We also created a cross plot comparing these two independent constraints on pCO 2 (Fig. 3c). Two regressions between ice core pCO 2 and boron-based pCO 2 are shown, a simple linear regression (grey line) and a Deming regression that takes into account error in variables (blue line). Bootstrapping was used to calculate uncertainties in the regression models (n = 1000, Fig. 3c, Table S6). While slopes and intercepts are not statistically different from a 1 : 1 line, the regressions do not reach a high significance level (p = 0.25); boosting the resolution of the record could help provide better constraints for this type of comparison. No significant difference in variability was observed at either site. The age models for the sites are based on comparisons of the benthic δ 18 O records for both Sites 806 and 807 (Fig. 3a, Zhang et al., 2007;Lear et al., 2003Lear et al., , 2015 to the published isotopic stack (Lisiecki and Raymo, 2005).
We also note that reconstructed pCO 2 uncertainties (both accuracy and precision) could potentially arise from Mg/Caderived estimates of temperature; these uncertainties could be reduced using independent temperature proxies for the WEP such as clumped isotope thermometry (Tripati et al., 2010, a technique that is not sensitive to the same sources of error as Mg/Ca thermometry and therefore is an area planned for future work. Other sources of uncertainty that have a larger effect on pCO 2 calculations are the weightshell correction, while the TA and seawater boron isotope composition have a minor effect over this time interval. Between MIS 7 and 6, our reconstructions exhibit a decrease in temperature ( T ) of 1.2 • C, an increase in pH ( pH) of 0.08 and a decrease in pCO 2 ( pCO 2 ) of 58 ppm. Between stage 3 and 1, we observed an increase in temperature of 2.0 • C, a decrease in pH of 0.13 and an increase in pCO 2 of 76 ppm. We also compare results with recent reconstructions in Figs. S1 and S2 (Sosdian et al., 2018;Rae et al., 2021). These results highlight that we are able to reproduce the range of atmospheric pCO 2 in the ice core record and reproduce the amplitude of changes between transitions, with uncertainties typical for this type of work (Hönisch et al., 2019).

Sea surface temperature in the WEP
Mg/Ca data are consistent at Site 806 (Wara et al., 2005;Tripati et al., 2009Tripati et al., , 2011Nathan and Leckie, 2009) and Site 872 (Sosdian et al., 2018) in the WEP. The Mg/Ca in T. sacculifer has to date not shown a pH dependency (Gray and Evans, 2019), but Mg/Ca of G. ruber does and was therefore corrected from this effect (see Supplement). Data for both species were corrected from salinity and seawater Mg/Ca changes. Mg/Ca temperatures for Site 872 were reconstructed using published data and the same framework we use here and are presented in Fig. 4. Recalculated values for Site 872 are from D. altispera, with an offset applied relative to T. sacculifer, and show similar variations to our record for the Miocene Climate Optimum and Middle Miocene Climate Transition (MCO-MMCT) periods (Sosdian et al., 2020). Temperatures from Tex 86 and U K 37 are plotted for comparison but those records are limited to the last 12 and 5 Myr, respectively (Zhang et al., 2014).
The Mg/Ca data support high temperatures of 35.2 ± 1.3 • C (2 SD, n = 11) for the early Miocene until the MMCT, with a relatively small (ca. 1 • C) change into the MCO, and larger changes out of the MCO. Similarly warm SSTs for the MCO were reconstructed in the North Atlantic at Site 608 from Tex 86 (Super et al., 2018). Despite a gap in our compilation from 11.5 to 9.5 Ma, there is a SST decrease of ∼ 6 • C from the MCO to ∼ 7 Ma where temperatures similar to present-day values are observed. A decline in temperature during the MMCT is coincident with the timing of the constriction of the Indonesian Seaway, the pre-closure of the trans-equatorial circulation and subsequent formation of a proto-warm pool (Nathan and Leckie, 2009;Sosdian et al., 2020). From 12 to 7 Ma, the Mg/Ca SST record diverges from Tex 86 and U K 37 -based reconstructions, with higher temperatures. At the same time, a record for the North Atlantic showed a decrease of ∼ 10 • C from the MCO to ∼ 9 Ma (Super et al., 2018). From 7 Ma to the present, the record from multiple proxies -Mg/Ca, Tex 86 and U K 37 -in the WEP agree.
3.4 Scenarios of seawater δ 11 B and alkalinity used for pCO 2 reconstructions Figures 5 and 6 show the different histories of seawater δ 11 B and alkalinity used in our calculations, respectively. Details of calculations are provided in the Supplement methods. Following the approach of Tripati et al. (2009Tripati et al. ( , 2014 and recent literature (Sosdian et al., 2018;Rae et al., 2021), we explored multiple scenarios for the evolution of seawater boron geochemistry ( Fig. 5) and alkalinity for calculations of pCO 2 (Figs. 6, S1 and S2). During the interval overlapping with the ice core record, we observe that the choice of model used does not make a significant difference in reconstructed values. During earlier time intervals, we see there is a greater divergence, reflecting larger uncertainties in seawater δ 11 B and alkalinity further back in Earth history. Prior to 10 Ma and during the early Pliocene (∼ 4.5 to 3.5 Ma), calculations of pCO 2 diverge from published values largely because of the different assumptions each study has used for past seawater δ 11 B (Fig. 5). However, we find that when the uncertainty in reconstructed pH is fully propagated, the differences in reconstructed pH values calculated using each of the δ 11 B seawater histories is not significantly different (Figs. 5 and 6; see also Hönisch et al., 2019). In contrast to the results from Greenop et al. (2017), the record from Raitzsch and Hönisch (2013) exhibits substantial variations on shorter timescales. Such variability is a challenge to reconcile with the Li isotope record of Misra and Froelich (2012), given that Li has a shorter residence time than boron while having similar sources and sinks. For the remainder of this study, we use the δ 11 B seawater history from Greenop et al. (2017) because it is in good agreement with seawater δ 7 Li (Misra and Froelich, 2012). The recent calculations of seawater pH (Sosdian et al., 2018;Rae et al., 2021) agree with values from our study when uncertainties are taken into account (Fig. 5).
The four alkalinity models used in this study diverge prior to 9 Ma, with a maximum difference at ∼ 13 Ma that is also reflected in reconstructed pCO 2 values (Fig. 6). However, all four models yield pCO 2 estimates that are within error of each other when the full uncertainty is considered. Uncertainty in the evolution of seawater alkalinity and seawater δ 11 B leads to differences in the absolute values of reconstructed pCO 2 (Fig. S2), and a divergence in reconstructed . Blue filled symbols are from Sites 806 and 807 with blue circles for T. sacculifer and triangles for G. ruber; filled gray squares are data from Site 872 (Sosdian et al., 2018). Open symbols are SST derived from Mg/Ca at Site 806 (Wara et al., 2005;Tripati et al., 2009;Nathan and Leckie, 2009). Tex 86 and U K 37 are also plotted for comparison (Zhang et al., 2014). Orange open circles are SST data calculated with our framework from the species D. altispera at ODP Site 806 (Sosdian et al., 2020) with an offset of +8 • C. Blue line is a smooth line (LOWESS) going through the data. . Different models for the evolution of the boron geochemistry explored as part of this work. Due to the 1 ‰ uncertainty propagated for δ 11 B seawater , all scenarios yield reconstructed seawater pH values that are within error of each other. Propagated uncertainties were calculated using Eq. (S14) (see Supplement). (a) Different models for δ 11 B seawater used for the reconstruction of pCO 2 in this study (blue -Lemarchand et al., 2002;green -Greenop et al., 2017;red -Raitzsch and Hönisch, 2013). (b) Reconstructed pH based on our measured δ 11 B carbonate values using different models for δ 11 B seawater (blue - Lemarchand et al., 2002;green -Greenop et al., 2017;red -Raitzsch and Hönisch, 2013)  . Different models for the evolution of a second carbonate (e.g., alkalinity) system parameter explored as part of this work. The propagated uncertainties were calculated using Eq. (S16) (see Supplement). (a) Different models for alkalinity used for the reconstruction of pCO 2 in this study (brown -constant alkalinity of 2330 µmol/kg, blue - Ridgwell and Zeebe, 2005;green -Tyrrell and Zeebe, 2004;violet -Caves et al., 2016). Colored symbols are reconstructed pCO 2 based on our measured δ 11 B carbonate values, alkalinity scenario and δ 11 B seawater from Greenop et al. (2017); open squares (compilation A) are the pCO 2 compilation from Sosdian et al. (2018), open triangles (compilation B) are from the compilation by Rae et al. (2021), black symbols are from site 872. (b) Reconstructed pCO 2 using constant alkalinity of 2330 µmol/kg and δ 11 B seawater from Greenop et al. (2017). (c) Reconstructed pCO 2 using the constant alkalinity scenario from Ridgwell and Zeebe (2005) and δ 11 B seawater from Greenop et al. (2017). (d) Reconstructed pCO 2 using constant alkalinity scenario from Tyrrell and Zeebe (2004) and δ 11 B seawater from Greenop et al. (2017). pCO 2 values that is largest in the Miocene. The two scenarios that produce the highest divergence in values are those calculated using constant alkalinity relative to those calculated using values from Caves et al. (2016), with a maximum difference at 15.06 Ma of up to 250 ppm CO 2 , and with the latter model producing lower values (Fig. 6b and e). Thus, for the MCO, alkalinity is a critical parameter in calculations of absolute pCO 2 values. For the Miocene and earlier intervals, improved constraints on past secular variations of seawater δ 11 B and alkalinity will yield more accurate reconstructions of pCO 2 .
For the remainder of this paper, we use the model of Caves et al. (2016) to estimate alkalinity and δ 11 B seawater determined by Greenop et al. (2017) (e.g., Fig. 6e). We note that two recent syntheses of boron isotope data have been published and compare our results to these findings (Figs. 8 and S2). Sosdian et al. (2018) report values that are in line with our results in the Miocene, but their study does not replicate results from ice cores. Rae et al. (2021) presents reconstructed values that are higher in the Miocene, due to the utilization of different scenarios of seawater δ 11 B and alkalinity compared to this work.

Miocene
The study of Miocene climate is thought to provide insights into drivers and impacts of global warming and melting of polar ice (Flower and Kennett, 1994). The Miocene epoch (23-5.3 Ma) is characterized by a warm interval, the Miocene Climate Optimum (∼ 17-14.7 Ma -MCO), and an abrupt cooling during the Middle Miocene Climate Transition (∼ 14-13 Ma -MMCT) that led to the expansion of ice on Antarctica and Greenland. Climate modeling supports a role for decreasing CO 2 in this transition (DeConto and Pollard, 2003). However, reconstructions for the Miocene are still relatively limited (Sosdian et al., 2018;Rae et al., 2021;Raitzsch et al., 2021). Boron isotope and alkenone-based pCO 2 reconstructions support higher pCO 2 during the MCO and a decrease over the MMCT (Sosdian et al., 2018;Stoll et al., 2019), consistent with what was previously inferred from B/Ca (Tripati et al., 2009(Tripati et al., , 2011Sosdian et al., 2020).
We applied the same framework we used for calculations at Sites 806 and 807 to published boron isotope data from Site 872 (Sosdian et al., 2018) in order to extend the WEP record to the early Miocene (Figs. 7, 8). The Miocene data between Sites 806 and 872 do not overlap as both are low in resolution but do show excellent correspondence in their trends in δ 11 B and reconstructed pH. For example, the closest data points in time at the two sites are at 15.6 Ma at Site 806 with a δ 11 B = 14.47 ± 0.21 ‰ and at 16.7 Ma at Site 872 with a δ 11 B = 15.12 ± 0.25 ‰. The pH values we reconstruct are within error of published estimates from Site 872 (Sosdian et al., 2018, Figs. 7d and 8d). Collectively, these data suggest that the early Miocene WEP was characterized by a mixed-layer pH of 8.1 ± 0.1 (2 SD, n = 4) between 19.4 and 21.8 Ma, which decreased to reach a minimum during the MCO of 7.7 (± 0.11 0.14 ). Given the sensitivity in absolute pCO 2 to assumptions about the second carbonate system parameter, a few scenarios were explored for the combined 806, 807 and 872 reconstructed pH values. For all alkalinity scenarios we used, reconstructed pCO 2 shows an increase from the Early Miocene to the MCO, with the highest values in the MCO. Recalculated pCO 2 for Site 872 between 19.4 and 21.8 Ma is 232 ± 92 ppm (2 SD, n = 4), lower but within error of the ones presented in Sosdian et al. (2018) and also within error of a constant alkalinity scenario (Fig. 8d). The main difference between our calculations and published reconstructions occurs between 19.4 and 21.8 Ma, when the same δ 11 B data for Site 872 from Sosdian et al. (2018) recalculated in Rae et al. (2021) yield higher pCO 2 , with an average value of 591 ± 238 ppm (2 SD, n = 4) because of the different assumptions used in their calculations. This difference is important because the assumptions from Rae et al. (2021) would imply a relatively high and stable pCO 2 from the early Miocene to MCO (Fig. S2), which would imply a decoupling between pCO 2 and temperature with no pCO 2 change during an interval of decreasing benthic δ 18 O. However, our reconstructed pCO 2 data increase towards the MCO is in line with the observed benthic δ 18 O decrease and δ 13 C increase and suggest a coupling between temperature and pCO 2 over this period. This highlights the critical need for the use of a common set of assumptions for studies. Assumptions may vary between studies depending on the timescales studied, but a common framework is needed. In addition, further constraints on the second carbonate system parameter and on secular changes in seawater δ 11 B will reduce uncertainties in reconstructed pCO 2 , with improved precision.
The highest pCO 2 values we reconstruct are found during the MCO (Fig. 6e). For the MCO, our estimates are 511 ± 201 ppm (2 SD, n = 3, Table 2). The middle Miocene values we reconstruct are in line with previous studies (Greenop et al., 2014;Sosdian et al., 2018). Published δ 11 B-based reconstructions also support higher pCO 2 for the MCO of ∼ 350-400 ppm (Foster et al., 2012) or 300-500 ppm (Greenop et al., 2014) that was recalculated by Sosdian et al. (2018) to be ∼ 470-630 ppm depending on the model of δ 11 B seawater chosen. During the MCO relative maxima in pCO 2 , our data support very warm sea surface temperatures in the WEP (35.6 ± 0.6 • C 2 SD, n = 3; Fig. 8c), that merits further examination in future studies. In fact, the highest temperatures recorded in our samples occur when there is a minimum in the global composite record of δ 18 O of benthic foraminifera (Zachos et al., 2001(Zachos et al., , 2008Tripati and Darby, 2018).
At the end of the MMCT, we find evidence for changes in pCO 2 and temperature in the WEP (Fig. 8). From 13.5 to 12.7 Ma, we reconstruct an increase in pH of ∼ 0.21 and  Sosdian et al. (2018), and open triangles (compilation B) are compilation data from Rae et al. (2021). (e) Reconstructed pCO 2 (ppm) using boron-based pH and alkalinity from Caves et al. (2016), data presented are from this study. Propagated uncertainties are given by Eq. (S17) for the dark blue envelope, while the light blue envelope shows the uncertainties calculated based on Eq. (S16) (taking into account uncertainty in δ 11 B seawater ). Crosses are original pCO 2 values calculated in Sosdian et al. (2018) at Site 872; asterisks are recalculated pCO 2 values at Site 872 by Rae et al. (2021).  4). (d) Reconstructed pCO 2 (ppm) from this study (blue symbols) using boron-based pH and alkalinity from Caves et al. (2016). Propagated uncertainties are given by Eq. (S17) for the dark blue envelope, while the light blue envelope reflects the uncertainties calculated based on Eq. (S16) (taking into account uncertainty on δ 11 B seawater ). Orange data points and envelope are calculated pCO 2 values and associated uncertainty from our study using our framework and a constant alkalinity scenario. Open squares (compilation A) are compilation data from Sosdian et al. (2018), open triangles are data from Raitzsch et al. (2021) at Site 1092. Crosses are original pCO 2 calculated in Sosdian et al. (2018) at Site 872; asterisks are recalculated pCO 2 at Site 872 by Rae et al. (2021); dark red triangles are from Site 1092 (Raitzsch et al., 2021). Data for compilation A are from Hönisch and Hemming (2009)  a major decrease in pCO 2 of ∼ 215 ppm during an interval highlighted by Flower and Kennett (1995), who observed changes in δ 18 O indicative of rapid East Antarctic Ice Sheet growth and enhanced organic carbon burial with a maximum δ 13 C reached at ∼ 13.6 Ma (Shevenell et al., 2004;Holbourn et al., 2007). As discussed in Sect. 3.4 the alkalinity model used for the calculations have an important impact during the Miocene which is likely responsible for the different absolute pCO 2 values over the MCO. In comparison, a scenario of constant alkalinity would lead to a pCO 2 during the MCO of 714 ± 313 ppm (2 SD, n = 3) and a decrease of ∼ 540 ppm during the MMCT. Both those reconstructions could simulate the large-scale advance and retreat of Antarctic ice with such low pCO 2 values (Gasson et al., 2016). At the same time, we find evidence for a decline in SST of 3.4 • C to minimum values of 33.3 • C. The synchronous shifts in the δ 13 C and δ 18 O of benthic foraminifera are consistent with increased carbon burial during colder periods, thus feeding back into decreasing atmospheric CO 2 and supporting the hypothesis that the drawdown of atmospheric CO 2 can in part be explained by enhanced export of organic carbon Kennett, 1993, 1995). However, given the limited sampling of this study, we are only able to resolve a pCO 2 decrease toward the end of the MMCT (∼ 13.5 Ma). The higher-resolution δ 11 B pCO 2 from Site 1092 for the MMCT (Raitzsch et al., 2021) reports eccentricity-scale pCO 2 variability; the authors reported that low pCO 2 during eccentricity maxima was consistent with an increase in weathering due to strengthened monsoonal circulation, which would increase nutrient delivery and support higher productivity that in turn would impact carbon drawdown and burial, in line with modeling from Ma et al. (2011). The resolution of our data during the late Miocene is low, with a data gap from 12.5 to 9.2 Ma and another gap between 6.5 and 5 Ma. We note the pCO 2 peak at ∼ 9 Ma observed by Sosdian et al. (2018) is not seen in our record, although this is likely due to the low resolution of our dataset. Between 9.5 and 7.1 Ma we find evidence for a decrease in atmo-spheric CO 2 of 100 ppm associated with a decrease in temperature of 1.3 • C. pCO 2 estimates derived from alkenones for Site 1088 (Tanner et al., 2020) do not show the same trend as boron-based reconstructions from the WEP or other regions (Fig. 6), which might be due to other controls on the alkenone proxy (Badger et al., 2019). A recent publication from Raitzsch et al. (2021) reports a δ 11 B reconstruction of pCO 2 that is within error of other δ 11 B isotope data from the Southern Ocean (Sosdian et al., 2018), although not for the same period as Tanner et al. (2020). pCO 2 differences between our reconstruction and that of Sosdian et al. (2018) and Raitzsch et al. (2021) (Fig. 8) likely reflect assumptions made for calculations (of δ 11 B, TA) and the specific monospecific calibrations used for each study, as well as potential geographic differences in air-sea pCO 2 . These differences do not invalidate the boron isotope proxy but illustrate the impact that specific seawater parameters and calibrations can have on reconstructed pCO 2 values, as well as potential inferences of air-sea disequilibrium.

Pliocene
Oxygen isotope data from a global benthic foraminiferal stack show that the Pliocene epoch (5.3-2.6 Ma) was initially characterized by warm conditions followed by the intensification of glaciation that occurred in several steps, including during MIS M2 (3.312-3.264 Ma), followed by the Middle Pliocene Warm Period (Lisiecki and Raymo, 2005). The Middle Pliocene Warm Period (mPWP -3.29-2.97 Ma) is considered a relevant geological analogue for future climate change, given ∼ 3 • C warmer global temperatures and sea levels that were ∼ 20 m higher than today (Dutton et al., 2015;Haywood et al., 2016), and is a target for model intercomparison projects, for which accurate paleo-atmospheric pCO 2 estimates are critical (Haywood et al., 2016).
We calculate high pCO 2 values of 419 ± 119 ppm (2 SD, n = 3, Table 2) between 4.7 to 4.5 Ma during the Early Pliocene warm interval (Fig. 9). The pCO 2 data we report provide a higher data density for the Early Pliocene and exhibit a trend that is in line with the reconstruction from Rae et al. (2021). Our data support values of 530 ± 110 ppm over the mPWP (2 SD, n = 4), higher than previously published data (Figs. 9, S2 and Table 2), although we acknowledge our low data density may not fully sample variability over this period. The similarity between our reconstructed values and those published for Site 871 in the Indian Ocean (Sosdian et al., 2018) suggests that changes in Indonesian through-flow do not induce substantial changes in air-sea exchange in the WEP.
The warmth and local pCO 2 maxima of the mPWP (mid-Pliocene Warm Period) was followed by a strong decrease in temperature in upwelling and high-latitude regions from 3.3 to 2.7 Ma, coincident with glacial intensification in the Northern Hemisphere. This climate transition was hypothesized to be driven by the closure of the Panama seaway the opening of the high latitudes and subsequent modifications of oceanic circulation (Haug and Tiedemann, 1998). However, modeling from Lunt et al. (2008) supports an additional major role for CO 2 in the glaciation. pCO 2 thresholds have been proposed to explain the intensification of Northern Hemisphere Glaciation, with values proposed ranging from 280 ppm (DeConto et al., 2008) to 200-400 ppm (Koenig et al., 2011).
The pCO 2 concentrations that we calculate indicate a reduction to 350 ppm by 2.7 Ma, ∼ 280 ppm by 2.6 Ma, and ∼ 210 ppm by 2.4 Ma, in several steps. These results support roughly a halving of CO 2 values when compared to values of ∼ 530 ppm at 3.3 Ma. These values are consistent with the pCO 2 thresholds proposed by both DeConto et al. (2008) and Koenig et al. (2011) for the intensification of Northern Hemisphere glaciation and the low atmospheric CO 2 (280 ppmv) scenario from Lunt et al. (2008). Mg/Ca SSTs decline from 30 to 26 • C, supporting an Earth-system sensitivity of ∼ 4 • C and doubling of CO 2 over this range, although given uncertainties, higher values of ∼ 6 • C and doubling of CO 2 that have recently been proposed (Tierney et al., 2020) cannot be excluded.
We speculate that at 4.42, 3.45 and 2.67 Ma, it is possible that the declines in CO 2 and ice growth, associated with Pliocene glacial intensification, in turn drove substantial changes in pole-to-Equator temperature gradients and winds, that in turn may have impacted iron cycling (Watson et al., 2000;Robinson et al., 2005;Martínez-Garcia et al., 2011), stratification (Toggweiler, 1999;Sigman et al., 2010) and other feedbacks that impact the amplitude of glacial-interglacial cycles and have been implicated as factors that could have contributed to Pliocene glacial intensification. Specifically, as the mean climate state of the planet became cooler, and glacial-interglacial cycles became larger in amplitude, enhanced windiness and dust transport and upwelling during glacials (Martínez-Botí et al., 2015b ) may have enhanced iron fertilization and subsequent carbon export (Martínez-Garcia et al., 2011). While data resolution is limited, we speculate this could explain why glacialinterglacial amplitudes in WEP pCO 2 values decrease from the mPWP towards the Pleistocene, whereas variations in δ 18 O are increasing -a speculation that could be tested with increased data resolution.

Pleistocene
During the Pleistocene (2.58-0.01 Ma), the climate system experienced a transition in glacial-interglacial (G/IG) variability from low-amplitude, higher-frequency and obliquitydominated oscillations (i.e., ∼ 41 kyr) of the late Pliocene to the high-amplitude, lower-frequency (∼ 100 kyr) cycles of the last 800 kyr. This transition is termed the Middle Pleistocene Transition (1.2-0.8 Ma -MPT). Questions have been raised about the role of atmospheric CO 2 during this transition, including using boron-based proxies (Hönisch et al.,  4). (d) Reconstructed pCO 2 (ppm) from this study (blue symbols) using boron-based pH and alkalinity from Caves et al. (2016). Propagated uncertainties are given by Eq. (S17) for the dark blue envelope, while the light blue envelope reflects the uncertainties calculated based on Eq. (S16) (taking into account uncertainty on δ 11 B seawater ). Open squares (compilation A) are pCO 2 compilation from Sosdian et al. (2018), open triangles (compilation B) are from the compilation by Rae et al. (2021). Data for compilation A are from Hönisch and Hemming (2009) 2009; Tripati et al., 2011;Chalk et al., 2017). Previous boron isotope studies for ODP Sites 668 and 999 in the tropical Atlantic Ocean have suggested that a decline in atmospheric CO 2 did occur during glacial periods in the MPT, but not during interglacials (Hönisch et al., 2009;Chalk et al., 2017;Dyez et al., 2018).
Our pCO 2 concentrations for Sites 806/807 reported here are in good agreement with those determined from ice cores from the early Pleistocene (Yan et al., 2019, Figs. 9 and 10), and with the boron-derived pCO 2 from a recent compilation (Rae et al., 2021). Results for the MPT are broadly in the range of values reported by Hönisch et al. (2009) andChalk et al. (2017). Although our data are relatively limited, we note they have greater resolution for the middle and later part of the transition than prior publications that have drawn conclusions about the MPT (Hönisch et al., 2009;Chalk et al., 2017;Dyez et al., 2018) (Fig. 10d), and therefore we explore their implications.
Taken alone, or when combined with the published data from Chalk et al. (2017) (that are also based on MC-ICP-MS), our results support a possible reduction of both glacial and interglacial pCO 2 values. We also find evidence that during the MPT, glacial pCO 2 declined rapidly from 189 ± 30 ppm at MIS 36  to reach a minimum of 170 (± 52 24 ) ppm during MIS 30. We note that pCO 2 concentrations are within error when uncertainty is fully propagated and then remained relatively stable until the end of the MPT, whereas interglacial pCO 2 values decrease gradually to reach post-MPT values.
In our record for the last 16 Myr, the lowest pCO 2 is recorded at MIS 30 during the MPT, with values of 164 (± 44 35 ) ppm, which supports an atmospheric CO 2 threshold that leads to large sheet generation. During this transition, the pCO 2 threshold needed to build sufficiently large ice sheets that were able to survive the critical orbital phase of rising obliquity to ultimately switch to a 100 kyr world was likely reached at MIS 30, but a higher pCO 2 resolution of the MPT is needed for confirmation. The multiple feedbacks resulting from stable ice sheets (iron fertilization, productivity, changes in albedo, changes in deep water formation) might have sustained larger mean global ice volumes over the subsequent 800 kyr. An asymmetrical decrease between pCO 2 values during interglacials relative to glacials, with glacials exhibiting the largest change across the MPT, would have led to increased sequestration of carbon during glacials in the 100 kyr world, as discussed by Chalk et al. (2017), with increased glacial dust input and iron fertilization.
3.6 Changes in volcanic activity and silicate weathering, and long-term pCO 2 On million-year timescales, atmospheric CO 2 is controlled by its input through mantle degassing in the form of subaerial and sub-aqueous volcanic activity and its removal by chemical weathering of continental silicate rocks. Over the last 16 Myr, two relative maxima in atmospheric pCO 2 are observed in our record, one during the MCO (at 15.67 Ma) and a second around the late Miocene-early Pliocene (beginning at 4.7 and 4.5 Ma) (Fig. 11), though the timing for the latter is not precise. The strong pCO 2 increase from the early Miocene to MCO occurs when there is increasing volcanic activity associated with the eruption of the Columbia River Flood Basalts (Hooper et al., 2002;Foster al., 2012;Kasbohm and Schoene, 2018), with recent geochronological evidence published supporting higher eruption activity between 16.7 and 15.9 Ma (Kasbohm and Schoene, 2018), reinforcing the idea of an episodic pCO 2 increase during the MCO due to volcanic activity. Underestimation of net CO 2 outgassing from specific continental flood basalt eruption is possible, as both sub-aqueous and sub-aerial flood basalts, under right climatic conditions, are prone to enhanced chemical weathering. For example, the 4 ‰-5 ‰ drop in δ 7 Li record at the Cretaceous-Paleogene (K-Pg) boundary (Misra and Froelich, 2012) is attributed to rapid quasi-congruent weathering of the Deccan Traps (Renne et al., 2015) during their eruption. Courtillot and Renne (2003) estimate that about 50 % of emitted CO 2 , roughly equivalent to the amount emitted by the eruption of a million cubic kilometers of Deccan Traps, may be missing due to chemical and physical weathering. Additionally, the early Eocene (at ∼ 50 Ma) 3 ‰-4 ‰ rise in seawater δ 7 Li at a time where there is not significant uplift of the Himalayas (Misra and Froelich, 2012) is also attributed to incongruent weathering of previously erupted Deccan Trap basalts as the Indian subcontinent moved from arid mid-latitudes to the wet low latitudes (Kent and Muttoni, 2008). Thus, a significant part of the outgassed CO 2 can be consumed by chemical weathering of freshly erupted hot basalts (Courtillot and Renne, 2003). However, the congruency of chemical weathering of basalts, depending on regional climatic conditions (warmwet vs. cold-arid), will determine the shape and position of inflection points in the seawater δ 7 Li record. The possible quantification of increased rates of silicate weathering inferred from δ 7 Li (mentioned below) can be utilized to determine total eruptive volume (missing + existing) and volatile emissions from the Columbia River Flood Basalts. At the same time as continental flood basalt emissions, enhanced seafloor production could also be a second possible source of CO 2 ; however, we note there is evidence that the rate of seafloor production has remained virtually invariant over the last 60 million years (Rowley, 2002;Müller et al., 2016). The second CO 2 peak can possibly be caused either by the observed increase in global volcanism during the earlymiddle Pliocene (Kennett and Thunell, 1977) and/or by a change in silicate weathering regime. Strontium and lithium isotopes ( 87/86 Sr and δ 7 Li) have been used as proxies for silicate weathering flux and congruency. Although the strontium isotope record exhibits a monotonous increase, lithium isotope data (Misra and Froelich, 2012) are more variable, with a transition from a period of increasing seawater δ 7 Li   Zachos et al., 2008). (b) Records from lithium isotopes (δ 7 Li, orange, Misra and Froelich, 2012) and strontium isotopes ( 87/86 Sr, grey, Hodell and Warnke, 1991;Farrel et al., 1995;Martin et al., 1999;Martin and Scher, 2004), both proxies for silicate weathering. Orange arrows represent the different weathering regimes as indicated by the δ 7 Li, black crosses indicates when changes in weathering regime occur. (c) Reconstructed pCO 2 (ppm) using boron-based pH and alkalinity from Caves et al. (2016), color indicates the site (filled light blue = 806, filled dark blue = 807), symbols represent the species (circle = T. sacculifer and triangle = G. ruber), filled grey squares (compilation A) are recalculated data based on Sosdian et al. (2018) Sosdian et al. (2018). Propagated uncertainties are given by Eq. (S17) for the dark blue envelope, while the light blue envelope are the uncertainties calculated based on Eq. (S16) (taking into account uncertainty on δ 11 B seawater ). Also shown is the timing of major events. The rose band and dark rose band indicate the eruption of the Columbia River flood basalts (Hooper et al., 2002) and time of maximum eruption (Kasbohm and Schoene, 2018), respectively.
It is interesting to note that the rise in δ 7 Li (Fig. 11b) from the early Miocene to the MCO is synchronous with the rise in pCO 2 . Before 18.5 Ma, the pCO 2 is relatively stable, and δ 7 Li is increasing, suggesting the non-steady-state and incongruent nature of continental chemical weathering. From 18.6 to 16.7 Ma, the δ 7 Li record decreases by ∼ 2 ‰, consistent with decreasing weathering rates and an associated increase in pCO 2 . Between 16.7 and 15.9 Ma, when the eruption of the Columbia River Flood Basalts is at a maximum, δ 7 Li increases, in line with higher weathering rates that could arise from higher atmospheric CO 2 and the presence of fresh basalts. The δ 7 Li record then decreases again until the end of the MCO at ∼ 14.7 Ma, in line with a decrease in the eruption rate, sustaining high atmospheric CO 2 . A constant increase in δ 7 Li is then observed, until the early Pliocene, where there is evidence for a shift to a steady-state weathering regime. This increase in δ 7 Li is also consistent with the decrease in pCO 2 observed until the early Pliocene.

Conclusions
We developed a reconstruction of atmospheric pCO 2 based on δ 11 B of planktic foraminifera from ODP Sites 806 and 807 located in the Western Equatorial Pacific for the past 16 million years and extended the record to 22 Ma by reprocessing data from Site 872 (Sosdian et al., 2018). We build on past efforts to reconstruct atmospheric pCO 2 using different proxies from this region, including from carbon isotopes in marine organic matter (Rayno and Horowitz, 1996) and alkenones (Pagani et al., 2010), as well as foraminiferal B/Ca ratios (Tripati et al., 2009(Tripati et al., , 2011, all of which have been shown to have a number of complexities and potential sources of systematic error (e.g., Tripati et al., 2011). It also builds on efforts to use boron isotopes in other regions using MC-ICP-MS (Seki et al., 2010;Foster et al., 2012Foster et al., , 2014Greenop et al., 2014;Martínez-Botí et al., 2015b;Stap et al., 2016;Chalk et al., 2017;Dyez et al., 2018;de la Vega et al., 2020), and our recent work constraining fractionation factors and measuring small samples of foraminifera (Guillermic et al., 2020).
Our study contributes a new long-term reconstruction of atmospheric pCO 2 for the Neogene derived from boron isotopes from the tropical Pacific Ocean. Although the record is not continuous, with variable resolution, it captures both long-term and short-term variability associated with several key transitions and demonstrates the utility of examining sites in the Western Equatorial Pacific for future higherresolution studies. Results for Sites 806 and 807 in the Western Equatorial Pacific reproduce the amplitude of late Pleistocene glacial-interglacial cycles in pCO 2 . These observations are consistent with the sites being in equilibrium with the atmosphere, although further work would be useful to explore sources of uncertainty and differences relative to ice core pCO 2 .
pCO 2 values increase from the early Miocene to the MCO with estimated MCO pCO 2 values of 511 ± 201 ppm (2 SD, n = 3). These elevated values are potentially linked to the eruption of the Columbia River Flood Basalts, with values declining into the early Pliocene, including during Pliocene glacial intensification. The changes in pCO 2 we observed are in line with changes in δ 7 Li, a proxy of silicate weathering, and future modeling of multiple proxy records should be insightful. Early Pliocene data for ∼ 4.7-4.5 Ma support high pCO 2 of 419 ± 119 ppm, and elevated values during the mid-Pliocene Warm Period of 530 ± 110 ppm for the time interval ∼ 3.3-3.0 Ma. These data are low in resolution, thereby not fully sampling orbital and millennial-scale variability. The higher-resolution record for the Pliocene glacial intensification supports a reduction in pCO 2 during several steps, with values at 2.7 Ma of 350 ppm, 2.6 Ma of ∼ 280 ppm and 2.4 Ma of ∼ 210 ppm. We find support for a larger reduction in glacial pCO 2 during the Mid-Pleistocene Transition compared to interglacial pCO 2 , and a minimum in pCO 2 during glacial MIS 30. These findings confirm a role for CO 2 in the transition from a 41 to a 100 kyr world.
Higher-resolution boron isotope records from the WEP would allow for further resolution of these changes. Additional constraints on temperature, such as from clumped isotopes (Tripati et al., 2010) in the WEP , could allow for uncertainties in pCO 2 estimates from boron isotopes to be reduced and for new constraints on Earth climate sensitivity. Future constraints on the vertical structure of the tropical Pacific (Shankle et al., 2021) during these transitions may also potentially be illuminating. Data availability. All data are available in the Supplement. Reconstructed climate parameters and proxy data will be archived at the NOAA's NCEI World Data Service for Paleoclimatology on acceptance at https://www.ncei.noaa.gov/products/paleoclimatology (NOAA, 2021).
Author contributions. AT developed the project and wrote the proposals that funded the work. All authors contributed to the experimental design. MG performed the measurements with assistance from SM. MG conducted data analysis with input from AT. MG drafted the paper, which was edited by all authors. Interpretation was led by MG and AT, with input from SM and RE.
Competing interests. The contact author has declared that neither they nor their co-authors have any competing interests Disclaimer. Publisher's note: Copernicus Publications remains neutral with regard to jurisdictional claims in published maps and institutional affiliations.