Lower oceanic δ 13 C during the last interglacial period compared to the Holocene

. The last time in Earth’s history when high latitudes were warmer than during pre-industrial times was the last interglacial period (LIG, 129–116 ka BP). Since the LIG is the most recent and best documented interglacial, it can provide insights into climate processes in a warmer world. However, some key features of the LIG are not well constrained, notably the oceanic circulation and the global carbon cycle. Here, we use a new database of LIG benthic δ 13 C to investigate these two aspects. We ﬁnd that the oceanic mean δ 13 C was ∼ 0.2 ‰ lower during the LIG (here deﬁned as 125–120 ka BP) when compared to the Holocene


Introduction
The most recent and well documented warm time period is the last interglacial period (LIG), which is roughly equivalent to Marine Isotope Stage (MIS) 5e (Past Interglacials Working Group of PAGES, 2016;Shackleton, 1969). The LIG began at the end of the penultimate deglaciation and ended with the last glacial inception (∼ 129-116 ka BP; Dutton and Lambeck, 2012;Govin et al., 2015;Masson-Delmotte et al., 2013;Menviel et al., 2019). The LIG was globally warmer than the pre-industrial period (PI, ∼ 1850(PI, ∼ -1900IPCC, 2013;Shackleton et al., 2020), with PI estimated to be ∼ 0.4 • C cooler than the peak of the Holocene (10-5 ka BP) (Marcott et al., 2013). Though not an exact analogue for future warming, the LIG may still help shed light on future climates. In particular, we seek to constrain the mean LIG ocean circulation and estimate the global oceanic mean δ 13 C.
As greenhouse gas concentrations were comparable to the Holocene, the LIG was most likely relatively warm because of the high boreal summer insolation (Laskar et al., 2004). During the LIG, the atmospheric CO 2 concentration was relatively stable around ∼ 280 ppm (Bereiter et al., 2015;Lüthi et al., 2008), while during the Holocene CO 2 first decreased by about 8 ppm starting at 11.7 ka BP before increasing by ∼ 17 to 277 ppm at ∼ 2 ka BP ( Fig. 1a)  . CH 4 and N 2 O peaked at ∼ 700 and ∼ 267 ppb, respectively, during both the LIG and the Holocene (Flückiger et al., 2002;508 S. A. Bengtson et al.: Lower oceanic δ 13 C during the last interglacial period compared to the Holocene Figure 1. LIG and Holocene time series of (a) CO 2 stack smoothed with a spline based on the age model AICC2012 , (b) sea surface temperatures (SSTs) determined from alkenones and aligned with oxygen isotopes from the Iberian Margin (MD01-2444, blue, Martrat et al., 2007b) and the North Atlantic (GIK23414-6, green, Candy and Alonso-Garcia, 2018), (c) EPICA Dome C ice core (EDC96) deuterium measurements (orange) and estimated surface air temperature anomaly relative to the mean of the last 1 kyr (red, Bazin et al., 2013;Jouzel et al., 2007) on the AICC2012 timescale, and (d) spline of atmospheric δ 13 C from EPICA Dome C and the Talos Dome ice cores (Holocene, Eggleston et al., 2016) and Monte Carlo average of three Antarctic ice cores atmospheric δ 13 C (LIG, Schneider et al., 2013) both based on the age model AICC2012. Shading around the lines indicates 1σ . Vertical grey shading indicates the periods of analysis in this paper. Vertical dotted grey lines indicate the commencement of the LIG and Holocene. Petit et al., 1999;Spahni et al., 2005). Global sea level was 6-9 m higher at the LIG compared to PI (Dutton et al., 2015;Kopp et al., 2009), thus indicating significant ice mass loss from both Antarctica and Greenland.
Strong polar warming is supported by terrestrial and marine temperature reconstructions. A global analysis of sea surface temperature (SST) records suggests that the mean surface ocean was 0.5 ± 0.3 • C warmer during the LIG compared to 1870-1889 (Hoffman et al., 2017), similar to another global estimate, which suggests SSTs were 0.7 ± 0.6 • C higher during the LIG compared to the late Holocene (McKay et al., 2011). However, there were differences in the timing of these SST peaks in different regions compared to the 1870-1889 mean: North Atlantic SSTs peaked at +0.6 ± 0.5 • C at 125 ka BP (e.g. Fig. 1b) and Southern Hemisphere extratropical SSTs peaked at +1.1 ± 0.5 • C at 129 ka BP (Hoffman et al., 2017). On land, proxy records from mid-latitudes to high latitudes indicate higher temperatures during the LIG compared to PI, particularly in North America (Anderson et al., 2014;Axford et al., 2011;Montero-Serrano et al., 2011). Similarly, the EPICA DOME C record suggests that the highest Antarctic temperatures from the last 800 kyr occurred during the LIG (Masson-Delmotte et al., 2010) (Fig. 1c).
Polar warming was also associated with significant changes in vegetation. Pollen records suggest a contraction of tundra and an expansion of boreal forests across the Arctic (CAPE, 2006), in Russia (Tarasov et al., 2005), and in North America (Govin et al., 2015;Muhs et al., 2001;de Vernal and Hillaire-Marcel, 2008). The few Saharan records suggest a green Sahara period during the LIG (Drake et al., 2011;Larrasoaña et al., 2013), consistent with a stronger West African monsoon (Otto-Bliesner et al., 2021). Although these reconstructions indicate changes in vegetation distribution during the LIG, the total amount of carbon stored on land remains poorly constrained.
Recent numerical experiments of the LIG as part of the Paleomodel Intercomparison Project Phase 4 (PMIP4) simulate significant warming over Alaska and Siberia in boreal summer, with mean annual temperature anomalies of close to zero, which is in good agreement with the proxy record (Otto-Bliesner et al., 2021). Despite this and other recent data compilations and modelling efforts (including Bakker et al., 2013), there are many open questions remaining about the LIG. In particular, stronger constraints are needed on the extent of Greenland and Antarctic ice sheets; on ocean circulation and the global carbon cycle, including CaCO 3 accumulation in shallow waters; and peat and permafrost carbon storage (Brovkin et al., 2016).
It is important to constrain the state of the Atlantic Meridional Overturning Circulation (AMOC) at the LIG given its significant role in modulating climate. Seven coupled climate models integrated with transient 130-115 ka BP boundary conditions simulate different AMOC trends, with some models producing a strengthening of the AMOC, while others simulate a weakening during the LIG (Bakker et al., 2013). Paleoproxy records suggest equally strong and deep North Atlantic Deep Water (NADW) during the LIG and the Holocene (e.g. Böhm et al., 2015;Lototskaya and Ganssen, 1999), with a possible southward expansion of the Arctic front related to changes in the strength of the subpolar gyre (Mokeddem et al., 2014), and AMOC weakening during a few multi centennial-scale events between 127 and 115 ka BP (e.g. Galaasen et al., 2014b;Helmens et al., 2015;Lehman et al., 2002;Mokeddem et al., 2014;Oppo et al., 2006;Rowe et al., 2019;Tzedakis et al., 2018).
Stable carbon isotopes are a powerful tool for investigating ocean circulation (e.g. Curry and Oppo, 2005;Eide et al., 2017) and the global carbon cycle (e.g. Menviel et al., 2017;Peterson et al., 2014). Since the largest carbon isotope fractionation occurs during photosynthesis, organic matter is enriched in 12 C (low δ 13 C), while atmospheric CO 2 and surface water dissolved inorganic carbon (DIC) become enriched in 13 C (high δ 13 C). Organic matter on land includes the terrestrial biosphere, as well as carbon stored in soils, such as in peats and permafrost. Different photosynthetic pathways (which differentiate C 3 and C 4 plants) fractionate carbon differently, producing typical signatures of about −37 ‰ to −20 ‰ for C 3 plants (Kohn, 2010) and around −13 ‰ for C 4 plants (Basu et al., 2015), though these values vary with a number of factors, including precipitation, atmospheric CO 2 concentration and δ 13 C, light, nutrient availability, and plant species (Cernusak et al., 2013;Diefendorf et al., 2010;Diefendorf and Freimuth, 2017;Farquhar, 1983;Farquhar et al., 1989;Keller et al., 2017;Leavitt, 1992;Schubert and Jahren, 2012). In the ocean, phytoplankton using the C 3 photosynthetic pathway are found to have fractionation during photosynthesis that depends on the concentration of dissolved CO 2 . Thus, atmospheric δ 13 CO 2 during the LIG (Fig. 1d) is influenced by the cycling of organic carbon within the ocean, changes in the amount of carbon stored in vegetation and soils, temperature-dependent airsea flux fractionation (Lynch-Stieglitz et al., 1995;Zhang et al., 1995), and on longer timescales by interactions with the lithosphere (Tschumi et al., 2011). The mean surface DIC is enriched by ∼ 8.5 ‰ compared to the atmosphere due to fractionation during air-sea gas exchange (Menviel et al., 2015;Schmittner et al., 2013).
NADW is characterised by low nutrients and high δ 13 C as a result of a high nutrient and carbon utilisation by marine biota and fractionation during air-sea gas exchange in the northern North Atlantic. Along its path through the Atlantic basin interior, organic matter remineralisation and mixing with southern source waters lowers δ 13 C, with δ 13 C values of ∼ 0.5 ‰ in the deep Southern Ocean.
The tight relationship between the water masses' apparent oxygen utilisation, nutrient content and δ 13 C allows δ 13 C to be used as a water mass ventilation tracer (e.g. Boyle and Keigwin, 1987;Curry and Oppo, 2005;Duplessy et al., 1988;Eide et al., 2017). The δ 13 C of benthic foraminifera shells, particularly of the species Cibicides wuellerstorfi, has been found to reliably represent the δ 13 C signature of DIC (Belanger et al., 1981;Duplessy et al., 1984;Zahn et al., 1986) and has therefore been used to better constrain the extent of different water masses. Mass balances of δ 13 C between the atmosphere, ocean, and land have been previously used to constrain changes in terrestrial carbon between the last glacial maximum (∼ 20 ka BP) and Holocene (e.g. Peterson et al., 2014). However, on longer timescales, exchanges with the lithosphere including volcanic outgassing (Hasenclever et al., 2017;Huybers and Langmuir, 2009), CaCO 3 burial in sediments and weathering, release of carbon from methane clathrates, and the net burial of organic carbon also influence the global mean δ 13 C. It has been estimated that the amount of carbon both entering and exiting the lithosphere due to weathering and burial of organic carbon fluxes could be from 0.274 to 0.344 Gt C yr −1 , though these vary through time (Hoogakker et al., 2006). Over timescales greater than 10 kyr, the influence of weathering and burial of carbon might dominate the δ 13 C signal (Jeltsch-Thömmes et al., 2019;Jeltsch-Thömmes and Joos, 2020), and thus a mass balance cannot be accurately applied to evaluate terrestrial carbon changes between the LIG and Holocene.
Here, we present a new compilation of benthic δ 13 C from Cibicides wuellerstorfi spanning the 130-118 ka BP time period. We use this data to compare the δ 13 C signal of the LIG with that of the Holocene and to determine the difference in average ocean δ 13 C between the two time periods. We then investigate the AMOC during the LIG with our new benthic δ 13 C database. Finally, we qualitatively explore the role of the various processes affecting the δ 13 C difference between the LIG and the Holocene.

Database
We present a new compilation of benthic δ 13 C covering the periods 130-118 and 8-2 ka BP. From these two sets of data, we select data pertaining to the LIG and compare it to data from the Holocene. Our database only includes measurements on Cibicides wuellerstorfi as no significant fractionation between the calcite shells and the surrounding DIC has been measured in this species (Belanger et al., 1981;Duplessy et al., 1984;Zahn et al., 1986).

Age models
Due to the lack of absolute age markers, such as tephra layers, the LIG age models mostly rely on alignment strategies that tie each record to a well-dated reference record. The age model tie points used in this study are taken from the original age model publications. The reference records (LS16; Lisiecki and Stern, 2016) consist of eight regional stacks (one for the intermediate and one for the deep ocean each for the North Atlantic, South Atlantic, Pacific, and Indian oceans) of benthic δ 18 O that were dated through alignment with other climatic archives such as ice-rafted debris records, synthetic ice core records, and speleothems. The use of regional stacks, rather than a single global stack, improved stratigraphic alignment targets and provided more robust age models. The estimated age model uncertainty (2σ ) for this group of cores is 2 kyr. Please refer to Lisiecki and Stern (2016) for further details. Oliver et al. (2010) defined their age tie points assuming that sea level minima and benthic δ 18 O maxima are synchronous. The benthic δ 18 O records Table 1. List of cores for the last interglacial period (LIG). Provided is the core name ("Core"), latitude ("Lat", • ), longitude ("Long", • ), depth ("Dep", m), the region, and the reference. Regions are abbreviated as follows: NEA: northeastern Atlantic; NWA: northwestern Atlantic; SWA: southwestern Atlantic; SEA: southeastern Atlantic; SA: southern Atlantic; NP: northern Pacific; SP: southern Pacific; I: Indian. Reference abbreviations are as follows: BW96: Bickert and Wefer (1996); CL82: Curry and Lohmann (1982); dA03: ; KJ8994: Jones (1989, 1994); KS02: Keigwin and Schlegel (2002); L99: Labeyrie et al. (1999); MB99: Mackensen and Bickert (1999); OH00: Oppo and Horowitz (2000); SH84: Shackleton and Hall (1984); SS0405: Shackleton (2004, 2005); VH02: Venz and Hodell (2002); V99: Venz et al. (1999); ZM1011: Mackensen (2010, 2011).  Table 2.
The published age models for the additional cores were determined using similar alignment techniques: SSTs were correlated to the NGRIP Greenland ice core for CH69-K09 and MD95-2042 (Govin et al., 2012). The age model for MD03-2664 was determined by correlating MD03-2664 δ 18 O with previously dated MD95-2042 δ 18 O (Galaasen et al., 2014b). ODP 1063 and U1304 δ 18 O were originally aligned to the LR04 stack (Lisiecki and Raymo, 2005). In order to align all of the records, adjustments to the age models of cores from Oliver et al. (2010) and the five additional cores (CH69-K09, MD95-2042, MD03-2664, ODP 1063, and U1304) were made by aligning the δ 18 O minima during the LIG to the corresponding δ 18 O minima of the nearest LS16 stack. The δ 18 O data before and after the alignment is given in Fig. S1 in the Supplement.
The Holocene age models are based on planktonic foraminifera radiocarbon dates Waelbroeck et al., 2001) that have been converted into calendar ages using IntCal13 and using reservoir ages based on modern observations (Key et al., 2004), which are assumed to have remained fairly stable across the Holocene. The age uncertainty associated with these Holocene radiocarbon-based age models is generally less than 0.5 kyr. However, it is important to note that Holocene age models from Oliver et al. (2010) were derived using the same method as their LIG age models, leading to larger age uncertainties of about 1-2.5 kyr for this set of Holocene records (four cores). The tie points were used to derive a full age-depth model assuming a constant sedimentation rate between tie points (i.e. linear interpolation).

Spatial coverage
The spatial distribution of the database for the Holocene and the LIG is shown in Fig. 2, and the depth distribution in each ocean basin is shown in Fig. 3. There are more data in the Atlantic Ocean (65 LIG, 118 Holocene) than in the Pacific (15 LIG, 19 Holocene) and Indian (3 LIG, 7 Holocene) oceans. We used this database to determine (1) if there is a significant difference in the average ocean δ 13 C signal at the LIG compared to the Holocene and (2) if ocean circulation patterns were comparable. Due to the sparsity of data in the Indian and Pacific oceans, our investigation is primarily focused on the Atlantic. Additionally, the temporal uncertainties (∼ 2 kyr) do not permit an investigation of centennialscale events, and therefore we restrict our analysis to mean LIG and Holocene conditions. Global distribution of benthic foraminifera δ 13 C covering the periods studied here: the Holocene (7-2 ka BP) (a) and LIG (125-120 ka BP) (b). Symbol size indicates the number of values per core, colour indicates average δ 13 C, and the triangle direction indicates the proxy depth (upward-pointing triangle: between 1000 and 2500 m depth; downward-pointing triangle: deeper than 2500 m). Four specific regions used in Sect. 3.1 are outlined: eastern equatorial Pacific (black, grey), equatorial Atlantic (yellow, green), southeastern Atlantic (cyan, blue), and northeastern Atlantic (magenta, red). Regional boundaries used to calculate the global volume-weighted mean δ 13 C (Sect. 3.2) are indicated by dotted black lines as defined in Peterson et al. (2014).

Results
The δ 13 C signal varies significantly regionally and with depth. The highest average δ 13 C values are associated with NADW and are generally found at depths between ∼ 1500 and 3000 m in the North Atlantic, with organic matter remineralisation and mixing with southern source waters leading to a δ 13 C decrease along the NADW path. The lowest δ 13 C values are in the deep South Atlantic (> 4000 m) because the Antarctic Bottom Water (AABW) end member is much lower than its NADW counterpart. Since the Indian and Pacific oceans are mostly ventilated from southern-sourced water masses, δ 13 C generally decreases northward in these two basins.
Since the number of cores is not consistent across the two time periods and given the high regional variability observed in δ 13 C, it is not possible to simply average all available data to determine the global mean δ 13 C. Furthermore, the spatial heterogeneity of the data density adds to the complexity of the problem. To address these points, we first analyse differences between the LIG and Holocene records for pre-defined small regions with high data density. We then calculate regional volume-weighted δ 13 C means for larger regions from which we estimate the global LIG-Holocene anomaly.
3.1 Regional reconstruction of δ 13 C We define regions with high densities of cores to reconstruct regional mean δ 13 C (Fig. 2). These regions need to be small enough to assume reasonably small spatial variability in the δ 13 C signal and yet still have enough data to establish a reliable statistical difference between the two time periods.
Based on these requirements, four regions are selected: the northeastern Atlantic, the equatorial Atlantic, a region off the Namibian Coast (southeastern Atlantic), and a region around the Galapagos Islands (eastern equatorial Pacific). The boundaries of each region are defined in Table 3.
We then define the time periods within the LIG and the Holocene to perform our analyses. For the Holocene, as most of the available data are dated prior to 2 ka BP, we define the end of our Holocene time period as 2 ka BP. To capture as much of the Holocene data as possible, we include data back to 7 ka BP, ensuring that we do not include instability associated with the 8.2 ka BP event (Alley and Ágústsdóttir, 2005;Thomas et al., 2007). This provides a time span of 5 kyr of data that we will consider for our analysis of the Holocene.
For the LIG, we seek to avoid data associated with the end of the penultimate deglaciation, which is characterised by a benthic δ 13 C increase in the Atlantic until ∼ 128 ka BP (Govin et al., 2015;Menviel et al., 2019;Oliver et al., 2010, Fig. 4). In addition, a millennial-scale event has been identified in the North Atlantic between ∼ 127 and 126 ka BP (Galaasen et al., 2014b;Tzedakis et al., 2018). Considering the typical dating uncertainties associated with the LIG data (2 kyr), we thus decide to start our LIG time period at Circular, coloured points connected by lines show each average δ 13 C value per core per time slice. Black symbols represent δ 13 C averages per slice. Each slice has a corresponding averaged depth (right y axis, m), with 1 standard deviation on either side shown on the bars. Slices with an average depth within ±300 m of the mean core depth of all slices are represented with a square point. Slices with an average depth 300 m shallower than the mean are shown with an upward triangle, and slices that are than 300 m deeper than the mean are shown with a downward triangle. Shading shows 1 standard deviation on either side of the mean for slices where more than one point exists. 125 ka BP. To ensure that the two time periods are of the same length (5 kyr), we define the LIG period for our analysis to be 125-120 ka BP. We note that our definition should also avoid data associated with the glacial inception (Govin et al., 2015;Past Interglacials Working Group of PAGES, 2016). We verify that the LIG time period has sufficient data across the selected four regions, noting that the highest density of data falls within the 125-120 ka BP time period -particularly in the equatorial Atlantic and southeastern Atlantic (Fig. 4b, c). Table 3. Regional summary of δ 13 C below 2500 m depth for the LIG (125-120 ka BP) and Holocene (7-2 ka BP) using a single value per core for each time slice. Shown are the non-volume-weighted means (δ 13 C, ‰), standard deviations (σ , ‰), and counts (N) for both time periods, along with the time period regional anomalies ( δ 13 C, ‰), propagated standard deviations for the anomaly (σ , ‰), and p values from two-sample t tests between the two time periods. We round the data to the nearest 1 kyr, find an average per 1 kyr, and refer to this as a time slice. We consider the LIG-Holocene anomaly across these two time periods for the four regions selected, and qualitatively consider the influence of changes in the average depth in which the proxies were recorded, as indicated by the direction of the black triangles in Fig. 4.
The average δ 13 C anomaly between the LIG and Holocene periods for cores deeper than 2500 m is consistent across the different regions despite their geographic separation, suggesting a significantly lower δ 13 C during the LIG than the Holocene, with differences ranging from −0.13 ‰ in the northeastern Atlantic to −0.20 ‰ in the equatorial Pacific (Table 3). The statistical significance between the two time periods is established using a two-tailed t test on data of one mean value per core and spans all time slices (125-120 and 7-2 ka BP). The t test shows that there is a statistically significant difference everywhere except in the equatorial Atlantic, with confidence intervals varying from 0.13 in the equatorial Atlantic to 0.04 in the northeastern Atlantic. When using a single tail t test instead, the difference becomes significant in the equatorial Atlantic, giving a new p value of 0.02. Figure 4 suggests that the difficulty in determining significance in this region for cores deeper than 2500 m might be due to a singular outlier measurement in the equatorial Atlantic: a value of −0.23 ‰ from GeoB1118 at ∼ 3.5 ka BP. If this value is excluded, then an anomaly of −0.22 with a p value less than 0.005 is observed.
We also compare the distribution of δ 13 C for cores deeper than 2500 m for three overlapping periods within the LIG (early LIG: 128-123 ka BP; LIG: 125-120 ka BP; late LIG: 123-118 ka BP). The results for the four regions are shown in Fig. 5. The statistical characteristics do not show much variation between the LIG and late LIG δ 13 C distributions. In the equatorial Pacific, the difference between the early LIG and the Holocene is smaller than between the LIG and Holocene, but this is countered with a larger difference in the equatorial Atlantic between the early LIG and Holocene. The spread in the data is generally larger during the Holocene than during the other time periods, which might be due to the greater number of measurements during the Holocene. The spread of data during the early LIG is slightly larger than during the LIG and late LIG in the equatorial and southeastern Atlantic. The equatorial Atlantic is the only region that displays significantly more points with lower δ 13 C during the early LIG. Overall, these distributions do not suggest that the LIG-Holocene anomalies that we have determined would be significantly impacted by slight variations in the selected time window. We perform an analysis of variance (ANOVA) on each region and post hoc tests on the data. We find that the Holocene data are significantly different from the three LIG periods in the northeastern Atlantic, the southeastern Atlantic, and the equatorial Pacific, while the three periods within the LIG are not significantly different from each other for any of the regions.
3.2 Volume-weighted regional δ 13 C The second approach we use to further constrain the LIG-Holocene δ 13 C anomaly is to estimate the volume-weighted regional δ 13 C. We define our regional boundaries based on the regions described in Peterson et al. (2014); however, we only include the regions where there are enough data to justify an analysis. For all the data in each of these regions, we calculate a mean value by taking the direct averages of all data. We divide the ocean basins into eight regions (Table 4, shown in Fig. 2) and calculate the volume-weighted averages δ 13 C for each of these regions. Since the Atlantic and Pacific oceans have more data than the Indian Ocean, there is greater confidence in the δ 13 C estimates for these regions. These regional averages are then used to calculate a global volumeweighted δ 13 C.
Results for the Atlantic and Pacific oceans are given in Fig. 6 and show a mean LIG-Holocene anomaly of −0.21 ‰ and −0.27 ‰, respectively, slightly higher than the range of estimates for the four regions selected in Sect. 3.1. There is a higher offset estimated in this definition of the southwestern Atlantic (−0.45 ‰) than in Sect. 3.1; however, there are only four cores available in this region during the Holocene and two during the LIG. Table 4. Regional breakdown of δ 13 C data for all depths during the Holocene (7-2 ka BP) and LIG (125-120 ka BP) averaged across the 1 kyr time slices. For each region: the average number of data points (labelled as "Points") and cores per time slice (labelled as "Cores"), the average standard deviation of δ 13 C per time slice (‰), the mean depth (m) across time slices, and the standard deviation of depth (m) between time slices (σ depth ). Regions are abbreviated as follows: NEA: northeastern Atlantic; NWA: northwestern Atlantic; SA: southern Atlantic; SEA: southeastern Atlantic; SWA: southwestern Atlantic; I: Indian, NP: northern Pacific; SP: southern Pacific.

Holocene LIG
Area δ 13 C Points Cores σ δ 13 C Mean depth σ depth δ 13 C Points Cores σ δ 13 C Mean depth σ depth (‰)  The estimated LIG-Holocene anomaly in the southern Pacific is relatively high at −0.39 ‰, giving a relatively large Pacific anomaly estimate of −0.27 ‰. This could be due in part to the deeper location of the LIG cores compared to the mean of the Holocene cores (439 m difference, Table 4). There is less confidence in the estimate of the Pa- Figure 6. Comparison of volume-weighted δ 13 C for the Atlantic (red) and Pacific (blue) for the LIG and Holocene, calculated using the regions from Peterson et al. (2014) from data covering all depths. Solid coloured lines indicate the mean volume-weighted δ 13 C, and the shading indicates the volume-weighted sum of square deviations from the mean. The horizontal bars indicate the mean of the stable period determined from the regional analysis as defined in Sect. 3.1 (LIG: 125-120 ka BP; Holocene: 7-2 ka BP), with the δ 13 C indicating the mean anomaly between these two averages and the standard deviation (σ (δ 13 C), ‰). cific volume-weighted mean since the proxy data are sparse, and the majority of cores are from the eastern equatorial Pacific as shown in Fig. 2. We also note that the average depths of cores from the Pacific Ocean (LIG: 2711 m; Holocene: 2131 m) and Indian Ocean (LIG: 2383 m; Holocene: 2303 m) are shallower than that of the Atlantic Ocean (LIG: 3531 m; Holocene: 3157 m; Fig. 3). However, as the vertical gradient below 2000 m depth in the Pacific Ocean is small (e.g. Eide et al., 2017), this might not significantly impact our results.
There is a small positive trend in the average Atlantic δ 13 C from 125 ka BP, reaching a maximum value at 118 ka BP (Fig. 6). The average core depth over the 125-120 ka BP time period does not suggest that a change in the mean depth could explain this variation. Fitting a linear regression over this period indicates an increase in δ 13 C of 0.03 ‰ kyr −1 in the Atlantic, with a p value of 0.01 and an R 2 of 0.85 (Fig. 4a). For the Pacific, there is a ∼ 0.13 ‰ increase in δ 13 C between 7 and 5 ka BP, which could be associated with the early Holocene terrestrial regrowth (Menviel and Joos, 2012).
For the Indian Ocean, we only include four cores, as these are the only ones spanning both the LIG and Holocene. An LIG anomaly of −0.13 ‰ in the Indian Ocean compared to the Holocene is therefore associated with higher uncertainties. The whole ocean mean LIG δ 13 C anomaly is −0.25 ‰, but it is associated with higher uncertainties than each region anomaly.
Both the regional analysis of our new database and our volume-weighted estimate indicate that the global mean δ 13 C was about 0.2 ‰ lower during the LIG than during the Holocene. We further test the robustness of this result in the next section.

Reconstruction of the LIG Atlantic Meridional
Overturning Circulation In this section, we analyse the spatial δ 13 C distribution in the Atlantic Ocean to assess potential changes in the penetration depth and southward expansion of NADW during the LIG, defined here as 125-120 ka BP, with respect to the Holocene. A change in NADW might influence our estimate of the mean δ 13 C, given that most of the available data are localised in the Atlantic Ocean. We use simple statistical regression models to reconstruct NADW and AABW separately with a quadratic-with-depth and linear-with-latitude equation following the method of Bengtson et al. (2019). For consistency, the regression algorithm only includes records from cores that span both the LIG and Holocene and uses a weighted least-squares approach, where the weighting equals the number of samples per core. The modelled region is defined between 40 • S and 60 • N as this is the region where we can expect to find both the NADW and AABW δ 13 C signals.
The results are shown in Fig. 7. We test the robustness of our statistical model using the jackknifing technique. We systematically exclude each individual core from the database one at a time, fit the parameters using this modified database, and compare the model prediction against the core, which was excluded. This produces small variations in the average mean response of the statistical models (the standard deviations were 0.01 ‰ for both the LIG and Holocene, respectively).
We calculate end member values based on proxies located near the water mass sources. These are taken as 0.79 ‰ and 1.02 ‰ for NADW for the LIG and Holocene, respec-tively, and −0.09 ‰ and 0.23 ‰ for AABW for the LIG and Holocene, respectively. The end member values are calculated as the average of cores shallower than 3000 m but deeper than 1000 m and located between 50 and 70 • N for NADW. The NADW end member cores have an average depth of 2043 m and a standard deviation of 478 m during the LIG. For the AABW end member, the only eligible core is ODP1089, which is at ∼ 41 • S and 4621 m.
The mean volume-weighted δ 13 C for the Atlantic Ocean between 40 • S and 60 • N based on this interpolation is 0.53 ‰ for the LIG and 0.70 ‰ for the Holocene (Fig. 7). This suggests a 0.17 ‰ lower Atlantic δ 13 C at the LIG than the Holocene. Our statistical reconstruction points to a very similar NADW depth (∼ 2600 m) for both time periods (Fig. 7). The NADW depth is defined here as the depth of maximum δ 13 C in the North Atlantic.
We also investigate the meridional gradient in δ 13 C in the Atlantic Ocean to determine whether the NADW southward penetration, transport, and remineralisation rates were significantly different during the LIG compared to the Holocene. We only consider cores that are located between depths of 1000 and 3000 m in order to stay within the main pathway of NADW (Fig. 8a). Though there is significant scatter, in accordance with our previous findings, a moving average through the Holocene and the LIG data shows that LIG δ 13 C is typically lower than the Holocene counterparts. However, the meridional δ 13 C statistical model gradients are not very different for the LIG (0.0036 ‰ per degree latitude) and the Holocene (0.0030 ‰ per degree latitude) (Fig. 8a), suggesting a similar southward penetration of NADW.
Using δ 13 C of the end members for NADW and AABW, we use a simple binary mixing model for all cores deeper than 1000 m to estimate changes in NADW penetration (Fig. 8b). The LIG and Holocene δ 13 C slopes in the Atlantic are similar, indicating similar southward penetration of NADW during both time periods. This suggests that the differences in δ 13 C between the two time periods is most likely due to change in end member values, while the mean Atlantic oceanic circulation was likely similar.
Based on our analysis, there appears to be no significant difference in the mean time-averaged AMOC between the LIG and the Holocene. Negative LIG-Holocene anomalies are found for each of the smaller regions selected (northeastern Atlantic, equatorial Atlantic, southeastern Atlantic, and eastern equatorial Pacific), with statistical significance seen in all regions except the equatorial Atlantic, where an unusual low δ 13 C value in one core is responsible for narrowing the difference between the two period means. Additionally, our volume-weighted mean δ 13 C estimates have similar anomalies in the Atlantic and Pacific oceans (−0.21 ‰ and −0.27 ‰, respectively, Fig. 6).

Figure 7.
Reconstructed Atlantic δ 13 C (‰) meridional section during the LIG (125-120 ka BP) and Holocene (7-2 ka BP). The circular points represent the proxy data, showing the average δ 13 C with colour and the number of points per core with size. The stars represent the proxy data that make up the end members. Background shading shows the reconstructed δ 13 C using a quadratic statistical regression of the proxy data following the method described in Bengtson et al. (2019).

Discussion
One of the goals of our study is to assess the mean change in oceanic δ 13 C between the LIG and the Holocene. Given the uncertainties in the chronologies, avoiding data that pertains to deglaciation, and capturing the same length of time during the LIG and the Holocene, we chose the periods 125 to 120 ka BP for the LIG and 7 to 2 ka BP for the Holocene. Using a similar geographical distribution of data points for both periods, we find that the oceanic δ 13 C was ∼ 0.2 ‰ lower during the LIG than the Holocene.
Our analysis of the δ 13 C signal suggests consistent LIG-Holocene δ 13 C anomalies in different regions of the Atlantic basins, as well as in the Pacific and Indian oceans, even if there are significant uncertainties with the later due to fewer available records. The δ 13 C distribution in the Atlantic Ocean suggests that there was no significant mean change in the southward penetration or depth of NADW during the LIG (125-120 ka BP) compared to the Holocene (7-2 ka BP). A statistical reconstruction of the early LIG (128-123 ka BP) δ 13 C compared to our 125-120 ka BP reconstruction does not reveal a significant difference in either the NADW core depth or NADW extent, as indicated by the meridional δ 13 C gradients (Fig. S2). The volume-weighted average δ 13 C during the early LIG is 0.06 ‰ lighter than during the LIG period considered here (125-120 ka BP). Since both time slices (128-123 and 125-120 ka BP) are 5 kyr averages and include dating uncertainties of ∼ 2 kyr, it is not possible to resolve potential centennial-scale oceanic circulation changes (e.g. Galaasen et al., 2014b;Tzedakis et al., 2018).
Explanations for the 0.2 ‰ lower δ 13 C anomaly in the ocean may include a redistribution between the oceanatmosphere system. Such a redistribution can result from a change in end member values (Fig. 8). As fractionation during air-sea gas exchange is temperature dependent, globally higher SSTs at the LIG could lead to a lower oceanic δ 13 C. However, the effect of this is likely small (Brovkin et al., 2002) and would also lead to a higher atmospheric δ 13 CO 2 at the LIG, which is inconsistent with Antarctic ice core measurements that suggest an anomaly of −0.3 ‰ . Lower nutrient utilisation in the North Atlantic would decrease surface ocean δ 13 C and thus the δ 13 C end members. However, this would also imply that less organic carbon would be remineralised at depth. Therefore, it is unlikely that the lower average oceanic mean δ 13 C results from a change in end members through lower surface ocean nutrient utilisation. Currently, there is still a lack of constraints on nutrient utilisation in these end member regions during the LIG compared to the Holocene. Therefore, the lower δ 13 C in the ocean-atmosphere system cannot be explained by a simple redistribution of δ 13 C between the atmosphere and the ocean.
An alternative explanation for the anomaly is a change in the terrestrial carbon storage, which has a typical signature of approximately −37 ‰ to −20 ‰ for C 3 -derived plant material (Kohn, 2010) and −13 ‰ for C 4 -derived plant material (Basu et al., 2015). The total land carbon content at the LIG is poorly constrained. Proxies generally suggest extensive vegetation during the LIG compared to the Holocene (CAPE, 2006;Govin et al., 2015;Larrasoaña et al., 2013;Muhs et al., 2001;Tarasov et al., 2005;de Vernal and Hillaire-Marcel, 2008), which would imply a greater land carbon store. However, other terrestrial carbon stores including peatlands and permafrost may also have differed during the LIG compared to the Holocene. With an estimated ∼ 550 Gt C stored in peats today (mean δ 13 C ∼ −28 ‰; Dioumaeva et al., 2002;Novák et al., 1999) and ∼ 1000 Gt C in the active layer in permafrost, which may have been partially thawed during the LIG (Reyes et al., 2010;Schuur et al., 2015;Stapel et al., 2018), less carbon stored in peat and permafrost at the LIG could have led to a lower total land carbon store compared to the Holocene. However, it is not possible to infer this total land carbon change from the oceanic and atmospheric δ 13 C anomalies because it cannot be assumed that the mass of carbon and 13 C is preserved within the ocean-atmosphere-land biosphere system on glacial-interglacial timescales. There is indeed continuous exchange of carbon and 13 C between the lithosphere and the coupled ocean, atmosphere, and land biosphere carbon reservoirs. Isotopic perturbations associated with changes in the terrestrial biosphere are communicated to the burial fluxes of organic carbon and CaCO 3 and are therefore attenuated on multi-millennial time scales (Jeltsch-Thömmes et al., 2019;Jeltsch-Thömmes and Joos, 2020). Nevertheless, when hypothetically neglecting any exchange with the lithosphere, we find that the change in terrestrial carbon needed to explain the difference in δ 13 C would be on the order of 310 ± 44 Gt C less during the LIG than the Holocene (Text S1).
In addition, due to the warmer conditions at the LIG than during the Holocene, there could have been a release of methane clathrates which would have added isotopically light carbon (δ 13 C: ∼ −47 ‰) to the ocean-atmosphere system. However, available evidence suggests that geological CH 4 sources are rather small (Bock et al., 2017;Dyonisius et al., 2020;Hmiel et al., 2020;Saunois et al., 2020) making this explanation unlikely, although we cannot completely exclude the possibility that the geological CH 4 source was larger at the LIG than the Holocene. Similarly, since δ 13 CO 2 from volcanic outgassing has a similar value to atmospheric δ 13 CO 2 (Brovkin et al., 2016) and modelling suggests volcanic outgassing likely only had a minor impact on δ 13 CO 2 (Roth and Joos, 2012), it is unlikely that volcanic outgassing of CO 2 played a significant role in influencing the mean oceanic δ 13 C.
While changes in terrestrial carbon could have impacted the oceanic δ 13 C at the LIG, the LIG-Holocene differences in the isotopic signal of both the atmosphere and ocean were most likely due to a long-term imbalance between the isotopic fluxes to and from the lithosphere, including the net burial (or redissolution) of organic carbon and CaCO 3 in deep-sea sediments, and changes in shallow water sedimentation and coral reef formation (Jeltsch-Thömmes and Joos, 2020).

Conclusions
We present a new compilation of benthic δ 13 C from 130 to 115 ka BP covering the LIG. Over this time period, benthic δ 13 C generally display a maximum value at ∼ 121 ka BP (±2 kyr), in phase with the maximum atmospheric δ 13 CO 2 (LIG value of −6.5 ‰ at ∼ 120 ka BP). As there are significant chronological uncertainties associated with LIG records, we analyse data between 125 and 120 ka BP to avoid data associated with millennial-scale events and the deglaciation. We compare this LIG benthic δ 13 C data to a similar database covering the Holocene (7-2 ka BP). We find that during these specific time periods, LIG oceanic δ 13 C was about 0.2 ‰ lower than during the Holocene. This anomaly is consistent across different regions in the Atlantic Ocean. Even though there are fewer records available, benthic δ 13 C data from the Pacific Ocean also support an anomaly of about 0.2 ‰.
An analysis of δ 13 C gradients across the Atlantic Ocean suggests that there were no significant changes in mean longterm ocean circulation between the two intervals. While reduced high northern latitude peat and permafrost caused by higher temperatures at the LIG than during the Holocene (Otto-Bliesner et al., 2021) could have lead to a lower atmospheric and oceanic δ 13 C, the most likely explanation for the lower LIG oceanic δ 13 C is a long-term imbalance in the weathering and burial of carbon. Additional studies are required to further constrain the LIG carbon balance. 520 S. A. Bengtson et al.: Lower oceanic δ 13 C during the last interglacial period compared to the Holocene Data availability. The data are published on Research Data Australia at DOI https://doi.org/10.26190/5efe841541f3b (Bengtson et al., 2020).
Author contributions. SAB, LCM, and KJM designed the research. CDP and LEL provided significant portions of the δ 13 C data. SAB, LCM, KJM, and LM analysed the data and developed the methodology. FJ assisted in the interpretation of the results. SAB prepared the manuscript with contributions from all coauthors.
Competing interests. The authors declare that they have no conflict of interest.
Acknowledgements. Shannon A. Bengtson acknowledges funding from the Australian Government Research Training Program Scholarship. Laurie C. Menviel and Katrin J. Meissner acknowledge funding from the Australian Research Council. Computational resources were provided by the NCI National Facility at the Australian National University through awards under the Merit Allocation Scheme, the Intersect allocation scheme, and the UNSW HPC at NCI Scheme. Fortunat Joos acknowledges funding from the Swiss National Science Foundation. This study was facilitated by the PAGES QUIGS working group. Review statement. This paper was edited by Marit-Solveig Seidenkrantz and reviewed by two anonymous referees.