Replicability of paleotemperature records in the northern Okinawa Trough and its implications for paleoceanographic reconstructions

Geochemical proxies are frequently utilized in the reconstruction of past ocean temperatures. Due to resource constraints, these reconstructions typically rely on a single sediment core, raising questions about the local and regional representativeness of paleotemperature records. To address this, we analyzed four sediment cores located within a 10-km radius in the northern Okinawa Trough (OT), which share the same climatic forcing and thus should reflect similar climate variations. We compiled published data and generated new paleotemperature estimates based on three widely used geochemical proxies (foraminiferal Mg/Ca, \({\text{U}}_{37}^{{{\text{K}}^{\prime}}}\), \({\text{TEX}}_{86}\)). Analysis of the mean absolute deviations for nearby records based on the same proxy revealed that \({\text{U}}_{37}^{{{\text{K}}^{\prime}}}\) has the highest reproducibility, followed by Mg/Ca and \({\text{TEX}}_{86}\). However, inconsistencies in inter-proxy offsets among nearby sites suggest the presence of noise in the proxy records, likely stemming from instrumental errors and sediment heterogeneity. Furthermore, the Mg/Ca and \({\text{U}}_{37}^{{{\text{K}}^{\prime}}}\) paleotemperature records agree within uncertainty when accounting for inter-site variability and calibration uncertainties, challenging previous interpretations of temperature signals from different seasons. All proxies indicate similar glacial-interglacial trends, albeit with varying magnitudes of temperature change. Both Mg/Ca and \({\text{U}}_{37}^{{{\text{K}}^{\prime}}}\) records suggest a glacial cooling of ~ 3 °C, whereas \({\text{TEX}}_{86}\) sea surface temperature (SST) data indicate a stronger glacial cooling of approximately ~ 6–8 °C. Modern observations indicate a subsurface \({\text{TEX}}_{86}\) recording depth of 50–100 m, coinciding with the thermocline. However, the \({\text{TEX}}_{86}\) subsurface temperature (subT) record does not resemble the Mg/Ca records of thermocline-dwelling foraminifera species. Instead, there is a better agreement with benthic foraminiferal Mg/Ca records of Uvigerina spp. (~ 700 m) and the intermediate temperature record derived from radiolarian assemblages (~ 500 m), pointing to a \({\text{TEX}}_{86}\) recording depth that is deeper than the thermocline. In summary, our findings show that proxy noise can impact inter-proxy comparisons of paleotemperature records, but not the direction of glacial-interglacial shifts. Future research should prioritize constraining the recording depth of paleotemperature proxies and reducing calibration uncertainty for more precise and reliable quantitative paleotemperature reconstruction.

Sea surface temperature (SST) is a crucial parameter in the climate system, as it reflects the air-sea heat exchange and surface circulation in the ocean. However, the available instrumental SST data only cover the past ~ 200 years. Understanding SST changes beyond the instrumental era therefore necessitates the use of indirect inferences. One approach is to use paleo SST estimates derived from marine sedimentary archives. These SST reconstructions typically rely on geochemical proxies, such as the alkenone-based \({\text{U}}_{37}^{{{\text{K}}^{\prime}}}\) index, the foraminifera-based Mg/Ca ratio, and the archaeal glycerol dialkyl glycerol tetraether (GDGT)-based \({\text{TEX}}_{86}\) index.

The recorded signals in paleotemperature proxies can be influenced by sedimentary factors. One such factor is the spatial heterogeneity of the proxy signal in sediments caused by sedimentation and biogeochemical processes. This heterogeneity can introduce noise into the proxy records and compromise their ability to accurately capture local climate variability. One way to assess proxy noise is by analyzing replicate cores taken within proximity. Their proximity to each other ensures that these replicate cores are influenced by the same climatic conditions. While studies of nearby replicate cores have been conducted on snow and ice cores in Antarctica (Münch et al. 2016; Wolff et al. 2005), there have been few studies on marine sediment cores. Nonetheless, the few existing studies based on the radiocarbon content of foraminifera tests have shown that the proxy signal in marine sediments is not as spatially homogeneous as previously thought (Dolman et al. 2021; Zuhr et al. 2022). Another study, based on a compilation of published datasets, found that Holocene paleotemperature records are not highly correlated even for nearby sites, plausibly due to the low signal-to-noise ratios in the paleotemperature records (Reschke et al. 2019). Of note, global or regional compilation studies treat all proxy data as the same regardless of the sample processing or instrumental setting employed. This is despite of inter-laboratory differences reported for commonly applied paleotemperature proxies, namely \({\text{U}}_{37}^{{{\text{K}}^{\prime}}}\) (Rosell-Melé et al. 2001), Mg/Ca (Greaves et al. 2008; Rosenthal et al. 2004), and \({\text{TEX}}_{86}\) (Schouten et al. 2013, 2009). Therefore, there is a need to further characterize the signal content in paleotemperature records to separate the climate signal from proxy noise.

In this study, we analyzed paleotemperature datasets based on four sediment cores that were retrieved within a radius of < 10 km in the northern Okinawa Trough (OT). The proximity of these sites suggests that they likely reflect the same oceanographic and climatic changes; therefore, we refer to them as replicate cores. By combining published and newly generated data, our objective was to assess the reproducibility of (1) paleotemperature records based on the same proxy type and (2) inter-proxy deviations at nearby sites, as well as to determine whether the inter-proxy deviation varies as a function of inter-site distance. In addition, we attempted to constrain the recording depth of \({\text{TEX}}_{86}\) by comparing \({\text{TEX}}_{86}\) records with Mg/Ca records from mixed layer-, thermocline- and bottom-dwelling foraminifera, as well as an intermediate water record based on radiolarian assemblage. Lastly, we discussed the implications of our findings for the reconstruction of the Kuroshio Current (KC) and East Asian Monsoon (EAM) in the northern OT.

2.1 Oceanographic setting: modern and paleo

The OT is a basin with a length of approximately 1100 km, oriented in a northeast-southwest direction. It is located along the rim of the East China Sea (ECS) (Fig. 1). The core sites are located in the Danjo Basin, which is in the northernmost part of the OT (Tada et al. 2015). The surface ocean in the northern OT is influenced by the KC and the EAM. The KC originates from the equatorial Pacific and flows northward along eastern Taiwan, passing through the Yonaguni Depression. It then enters the OT and continues to flow northward along the slope of the ECS. This western-boundary current provides warm and saline water to the region. Long-term observations over 28 years in the central OT along a transect known as "section PN" have shown that the maximum volume transport of the KC in the summer is 27.0 Sv (1 Sv = 10⁶ m³/s), and the minimum volume transport in autumn is 23.9 Sv (Ichikawa and Beardsley 2002). The main axis of the KC does not directly flow through the core sites (Fig. 1a). Instead, a branch of the KC, known as the Tsushima Warm Current (TSWC), flows northward toward the Tsushima Strait and directly influences the core sites. During the summer, increased precipitation due to the East Asian Summer Monsoon (EASM) leads to higher discharge of the Changjiang River. This, in turn, results in the formation of the cool and less saline Changjiang Diluted Water (CDW) which is transported northward (Ichikawa and Beardsley 2002). The annual mean SST at the core sites is ~ 23 °C, with the highest SST of ~ 28 °C occurring in August and the lowest SST in March at ~ 17 °C (Fig. 1c; World Ocean Atlas 2018; Locarnini et al. 2018). Consequently, this yields a seasonal SST range of ~ 11 °C.

Studies based on various microfossil groups, including planktic foraminifera (Ijiri et al. 2005), diatoms (Shirota et al. 2021) and radiolarians (Matsuzaki et al. 2019), all indicate that the northern OT experienced substantial glacial-interglacial shifts in hydrography. As a result of the weakened KC and lower sea level during the Last Glacial Maximum (LGM), the northern OT received less warm, less saline water advected by the KC and more freshwater via riverine input, resulting in a nutrient-rich coastal environment with low temperature and salinity (Ijiri et al. 2005; Matsuzaki et al. 2019; Shirota et al. 2021).

2.2 Regional interpretation of paleotemperature proxies

The interpretation of paleotemperature proxy records in the study region can be informed by modern proxy observations (Table 1). In the northern OT, the Mg/Ca ratio of mixed layer-dwelling foraminifera such as Globigerinoides ruber and the \({\text{U}}_{37}^{{{\text{K}}^{\prime}}}\) index are commonly used as indicators of SST, albeit with different seasonal signals. The \({\text{U}}_{37}^{{{\text{K}}^{\prime}}}\) index has been regarded as representing the annual mean SST because the temporal distribution of alkenones in suspended particulate matter does not exhibit a discernible seasonal bias (Ko et al. 2018). Conversely, plankton net sampling conducted during different months in the northern OT (Yamasaki et al. 2010) indicated that the production of mixed layer-dwelling G. ruber is disproportionately higher during the summer months, suggesting that the Mg/Ca-derived temperature of G. ruber may reflect an SST biased toward summer. In contrast, foraminifera species such as Neogloboquadrina dutertrei, which reside in the thermocline throughout the year, record mean annual thermocline conditions. The interpretation of \({\text{TEX}}_{86}\) records is more complex compared to the other two proxies due to uncertainties surrounding the depth origin of the sedimentary \({\text{TEX}}_{86}\) signal. In other regions, the depth origin of \({\text{TEX}}_{86}\) signal can range from SST to intermediate depths (Ho and Laepple 2016; van der Weijst et al. 2022), or even bottom water depths (Kim et al. 2015). In the northern OT, \({\text{TEX}}_{86}\) records have been interpreted as reflecting a winter SST signal (Nakanishi et al. 2012; Yamamoto et al. 2013). However, it has been observed that GDGTs in suspended particulate matter are more abundant at shallow subsurface depths (50–100 m) in this region (Ko et al. 2022; Nakanishi et al. 2012), suggesting that \({\text{TEX}}_{86}\) may actually reflect temperatures near the thermocline, which is located at ~ 75 m depth in the northern OT (Fig. 1c).

2.3 Reconstructions based on paleotemperature proxies

Numerous reconstructions have been conducted in the northern OT to understand past changes in SST and sea surface salinity (SSS) in response to glacial-interglacial variability of the KC and EAM. In the following section, we present a brief overview of these reconstructions.

Previous studies have employed different geochemical proxies to reconstruct paleo SSTs in the northern OT. For instance, Ijiri et al. (2005) estimated temperatures using the \({\text{U}}_{37}^{{{\text{K}}^{\prime}}}\) index. In combination with records of planktic foraminiferal fauna and stable oxygen isotopes of G. ruber, their findings suggest that the KC entered the OT during the LGM, which contradicts several other studies (Chang et al. 2008; Ujiié et al. 2016, 2003). Yamamoto et al. (2013) demonstrated a pronounced cooling during the LGM based on the \({\text{TEX}}_{86}\) index. This cooling is attributed to a weakened KC that did not fully reach the northern OT and an intensified East Asian Winter Monsoon. More recently, Lee et al. (2021) estimated the range of seasonal temperatures (referred to as seasonality) during glacial periods based on temperature records obtained from the \({\text{U}}_{37}^{{{\text{K}}^{\prime}}}\) index and Mg/Ca ratio. This study also emphasized the influence of global atmospheric pCO₂, rather than regional climate variability, on controlling the annual mean temperature in the northern OT.

The combination of paleo SST and stable oxygen isotopes (δ¹⁸O) is a commonly used approach to reconstruct past variations in the intensity of the EASM, as δ¹⁸O is influenced by riverine input and precipitation. Datasets from Kubota et al. (2010) demonstrate that variations in EASM precipitation inferred from marine proxies (Mg/Ca and δ¹⁸O) during the Holocene align with those reflected by the δ¹⁸O of stalagmites. Millennial-scale fluctuations in SST and SSS in the northern OT may be explained by changes in the mixing ratio of the KC and CDW. However, Clemens et al. (2018) later proposed that greenhouse forcing and global ice volume, rather than EASM variability, were the primary drivers of glacial-interglacial seawater δ¹⁸O variations in the northern OT. This finding suggests that EASM precipitation may manifest differently in marine and terrestrial archives.

3.1 Study setup and approach

We compiled published paleotemperature data and generated new paleotemperature estimates for four sediment cores (Table 2) located within a 10-km radius of each other in the northern OT (Fig. 1; Table 3). These cores, specifically MD98-2195, KY07-04–01, KR07-12 PC01 and IODP U1429, were collected using different coring equipment (Table 2) during various research cruises. Hereafter, we will refer to these cores as MD98, KY07, KR07, and U1429. The age models of these cores have been previously published (Table 2), and are based on radiocarbon dating of foraminifera, tephra chronology, and benthic δ¹⁸O stratigraphy (Fig. 2). Given the inter-site distances of less than 10 km and similar water depths (Table 2), it is assumed that all sites experience the same climatic and oceanographic changes. Therefore, the paleotemperature records from these sites are expected to be similar. This unique setup allows us to assess, for the first time, how representative a proxy record is of the local paleotemperature signal by examining the similarity of paleotemperature records from the same location but based on different sediment cores and analyzed in different laboratories.

Among the four sites, MD98 and KY07 are the closest, with a distance of less than 1 km between them (Table 3) and both characterized by comparable sedimentary age-depth relationship and mean sedimentation rate (~ 80 cm/kyr; Fig. 2). U1429 is ~ 5 km east from these two sites and ~ 7 km from KR07. KR07 is located ~ 9 km northeast from both MD98 and KY07. The mean sedimentation rates at these sites, i.e., ~ 30 cm/kyr at KR07 and ~ 50 cm/kyr at U1429, are lower than that at MD98 and KY07. The sedimentation rates at all four sites are relatively stable without any large sudden shifts, except a depth horizon in KY07 that is adjacent to the Kikai-Akahoya tephra layer (Fig. 2b). If the inter-site deviation of the proxy values varies as a function of the distance between records, we would expect the smallest deviation between MD98 and KY07, and the largest deviation between these two sites and KR07.

In addition, we re-analyzed the \({\text{TEX}}_{86}\) record of MD98 using a more recent protocol. To assess whether different laboratory protocols yield comparable results in Mg/Ca data, we generated replicate data for the Mg/Ca proxy using selected samples from core KR07. For this purpose, we picked planktic foraminifera G. ruber from both sensu stricto (s.s.) and sensu lato (s.l.) morphotypes; the former is also recognized as G. ruber albus and the latter Globigerinoides elongatus (Aurahs et al. 2011). From this point forward, we will refer to a combination of these two morphotypes as G. ruber. Furthermore, to enhance the understanding of the uncertain recording depth of \({\text{TEX}}_{86}\), we generated Mg/Ca-based temperature records for thermocline-dwelling foraminifera (N. dutertrei) and benthic foraminifera (Uvigerina spp.). This was done to address the modern regional studies that have suggested a subsurface origin for sedimentary GDGTs (Ko et al. 2022; Nakanishi et al. 2012).

3.2 Foraminiferal Mg/Ca analysis

We conducted Mg/Ca analysis of mixed-layer-dwelling G. ruber from core KR07, as well as thermocline-dwelling N. dutertrei and benthic Uvigerina spp. from cores KY07 and KR07. The cleaning procedures were adapted from Shen et al. (2001). Approximately 10–20 individuals of foraminifera (sized between 250–355 μm) were collected per sample and the shells were crushed between glass plates to expose the chambers for cleaning. Four reagents were sequentially used to eliminate clay minerals, organic matter, exchangeable ionic compounds, and metal oxides. The reagents were added in the following order: (1) ethanol with a 40-min ultrasonic bath, (2) 1% H₂O₂ heated at 70 °C for 2 h, (3) 1.0 M NH₄Cl heated at 70 °C for 20 min, (4) 0.01 M reduction reagent (NH₂OH + CH₃COONa) heated at 70 °C for 3 h, and (5) 5 ×ばつ 10^–4 N HNO₃ for 5 min. Two Milli-Q water rinses were performed between steps to remove any remaining chemical agents and impurities dissolved in the previous stage. After cleaning, the samples were dissolved in 5% nitric acid before analysis. The Mg/Ca ratio was analyzed using a Sector-Field Inductively Coupled Plasma Mass Spectrometer (SF-ICP-MS; Thermo Scientific Element 2) at the Department of Geosciences, National Taiwan University (Lo et al. 2014). The precision of the Mg/Ca analysis during the measurement of our samples was ± 0.8% (2RSD), corresponding to ± 0.02 mmol/mol and ± 0.2 °C in the temperature scale when using the calibration by Hasting et al. (2001). Temperature estimates were derived by converting the Mg/Ca values using species-specific calibrations listed in Table 4.

Eleven Mg/Ca data of G. ruber albus from core KR07, spanning the early to late Holocene, were generated following the method of Kubota et al. (2019). The foraminiferal tests were cleaned using the reductive method, as described by Boyle and Keigwin (1985), with slight modification. The cleaning steps were conducted in the following order: (1) 25–40 individuals were crushed in a 1.5-mL microtube using a thin plastic rod, (2) 1-min ultrasonication with Milli-Q water, repeated four times, (3) 1-min ultrasonication with ultrapure water, repeated twice, (4) 30-min heating bath (~ 80–90 °C) with 100 μL reduction reagent (0.25 M citric acid in 16 M ammonia, made up to 1 M in hydrazine (NH₂NH₂)), (5) 10-min heating bath (~ 80–90 °C) with 250 μL (2 ×ばつ 10^–5 M H₂O₂ in 0.1 M NaOH), repeated twice, (6) leaching with 250 μL 0.001 M HNO₃ for a few seconds in an ultrasonic bath. Two Milli-Q water rinses were performed between steps. After cleaning, the samples were dissolved in 0.075 M HNO₃ before analysis. Mg/Ca analysis was performed using a SF-ICP-MS Thermo Element 2 installed in a clean room (Class 10,000) at the Mutsu Institute for Oceanography, Japan Agency for Marine-Earth Science and Technology (JAMSTEC). Isotopes of four elements (²⁴Mg, ⁴⁴Ca, ⁴⁸Ca, and ⁵⁵Mn) were analyzed using Sc as the internal standard with 0.3 M HNO₃, following the method described by Uchida et al. (2008). Four working standards were prepared by successive dilutions of the stock standard solutions to match the concentrations of Ca (100 ppb, 500 ppb, 2 ppm, and 5 ppm) and Mg (0.05 ppb, 0.2 ppb, 1.0 ppb, and 5 ppb), respectively, which covered the Ca and Mg concentration ranges of all samples (Kubota et al. 2015, 2010). The analytical method and equipment used were the same as those used for KY07 samples reported by Kubota et al. (2015, 2010). The precision of the replicate analysis of the working standard for Mg/Ca was better than ± 0.09 mmol/mol, corresponding to ± 0.3 °C in the temperature scale. All measurements were taken during the interval July 2010 to July 2013. Note that the Mg/Ca analysis for U1429 was conducted in the same laboratory but with a SF-ICP-MS Thermo Scientific Element XR (Clemens et al. 2018), an updated version of Element 2. For the analysis on U1429, the accuracy of Mg/Ca ratios was confirmed by analyzing a CaCO₃ reference standard, coral Porites standard material JCp-1, whose Mg/Ca values were internationally determined (4.199 ± 0.065 mmol/mol) (Hathorne et al. 2013), but JCp-1 was not used for the analysis of KR07 samples.

3.3 Alkenone and GDGT analyses

We reported new biomarker data for cores KR07 and MD98. Both, alkenone and GDGT analyses were conducted on core KR07, while only GDGT analyses were carried out on core MD98 as \({\text{U}}_{37}^{{{\text{K}}^{\prime}}}\) data have been published in Ijiri et al. (2005).

3.3.1 KR07

For each sample, 2–4 g of freeze-dried sediment was homogenized using an agate mortar and pestle. Before extraction, internal standards 10-nonadecanone (C₁₉ ketone) and C₄₆ GTGT were added to the sample to facilitate the quantification of alkenones and GDGTs, respectively. The total lipids were then extracted by utilizing a microwave-assisted extraction system (ETHOS X) with a mixture of dichloromethane (DCM) and methanol (MeOH) (2:1, v/v). The initial conditions of the microwave oven were set at 20 °C and 40 bar, which were subsequently increased to 40 °C and 90 bar, respectively, in 10 min and then maintained for 20 min. The extraction process was repeated twice. The supernatant from each extraction was combined, subjected to centrifugation at 3500 rpm for 5 min, and then concentrated under a stream of nitrogen gas. The resultant total lipid extract was then passed through a column of Na₂SO₄ and fractionated into three fractions via silica gel column chromatography. Alkenones and GDGTs were eluted using DCM and MeOH, respectively.

Alkenone analysis was performed at the Institute of Oceanography, National Taiwan University employing a gas chromatography system (Thermo Scientific TRACE 1300) equipped with a GC capillary column (Restek Rtx-200MS; 105 m ×ばつ 250 μm ×ばつ 0.25 μm) and a flame ionization detector (FID). The oven program was as follows: the temperature was maintained at 50 °C for 2 min, then increased to 255 °C at a rate of 20 °C/min, and finally raised to 320 °C at a rate of 3 °C/min and kept constant for 25 min. Nitrogen was used as the carrier gas with a flow rate of 1.5 mL/min. The analytical precision, determined through replicate analysis of an in-house standard (GH06), is 0.010 \({\text{U}}_{37}^{{{\text{K}}^{\prime}}}\) units, which corresponds to 0.3 °C when using the global calibration by Prahl et al. (1988).

GDGTs extracted from core KR07 were analyzed at the Alfred Wegener Institute (Bremerhaven) using an Agilent 1260 Infinity II ultrahigh-performance liquid chromatography-mass spectrometry (UHPLC-MS) system, consisting of a G1712B binary pump, a G7129A vial sampler with integrated sample thermostat, a G7116A multicolumn thermostat, and a G6125C single quadrupole mass spectrometer with an atmospheric pressure chemical ionization (APCI) ion source. Chromatographic separation (including 5-and 6-Methyl isomers of branched GDGTs) of the GDGTs was achieved by coupling two UPLC silica columns (Waters Acquity BEH HILIC, 2.1 ×ばつ 150 mm, 1.7 μm) with a 2.1 ×ばつ 5 mm pre-column as in Hopmans et al. (2016), but with the following chromatographic modifications: Mobile phases A and B consisted of n-hexane/chloroform (99:1, v/v) and n-hexane/2-propanol/chloroform (89:10:1, v/v/v), respectively. The flow rate was set to 0.4 mL/min and the columns were heated to 50 °C, resulting in a maximum backpressure of 425 bar. Sample aliquots of 20 μL were injected with isocratic elution for 20 min using 86% A and 14% B, followed by a gradient to 30% A and 70% B within the next 20 min. After this, the mobile phase was set to 100% B, and the column was rinsed for 13 min, followed by a 7 min re-equilibration time with 86% A and 14% B before the next sample analysis. The total run time was 60 min. GDGTs were detected using positive ion APCI-MS and selective ion monitoring (SIM) of (M + H)⁺ ions (Schouten et al. 2007) with the following settings: nebulizer pressure 50 psi, vaporizer and drying gas temperature 350 °C, drying gas flow 5 L/min. The capillary voltage was 4 kV and the corona current + 5 μA. The detector was set for the following SIM ions: m/z 744 (C₄₆ standard), m/z 1302.3 (GDGT-0), m/z 1300.3 (GDGT-1), m/z 1298.3 (GDGT-2), m/z 1296.3 (GDGT-3), m/z 1292.3 (crenarchaeol and crenarchaeol isomer), m/z 1022 (GDGT-Ia), m/z 1020 (GDGT-Ib), m/z 1018 (GDGT-Ic), m/z 1036 (GDGT-IIa and IIa’), m/z 1034 (GDGT-IIb and IIb’), m/z 1032 (GDGT-IIc and IIc’), m/z 1050 (GDGT-IIIa and IIIa’), m/z 1048 (GDGT-IIIb and IIIb’), m/z 1046 (GDGT-IIIc and IIIc’). The resulting scan/dwell time was 66 ms. Repeated analyses of the in-house standard (GH06) indicate an analytical precision of 0.007 \({\text{TEX}}_{86}^{\text{H}}\) units, which corresponds to 0.4 °C when using the global calibration by Kim et al. (2010).

3.3.2 MD98

Lipids were extracted with three repetitions from ~ 3 g of dried sediment using a DIONEX Accelerated Solvent Extractor ASE-350 with dichloromethane/methanol (6:4, v/v). The total lipid extract was separated into four fractions by SiO₂ column chromatography (5 mm diameter, 45 mm high): F1 with 3 mL hexane, F2 with 3 mL hexane-toluene (3:1), F3 with 4 mL toluene, and F4 with 3 mL toluene–methanol (3:1).

For GDGT analysis, F4 was dissolved in hexane:2-propanol (99:1, v/v) and filtered through a short bed of Na₂SO₄. GDGTs were analyzed using high-performance liquid chromatography-mass spectrometry (HPLC–MS) with an Agilent 1260 HPLC system coupled to an Agilent 6130 Series single quadrupole MS. Separation was achieved using a Prevail Cyano column (2.1 ×ばつ 150 mm, 3 μm; Grace Discovery Science, USA) maintained at 30 °C following the method of Hopmans et al. (2000) and Schouten et al. (2007). Conditions were as follow: flow rate of 0.2 mL/min, isocratic with 99% hexane and 1% 2-propanol for the first 5 min, followed by a linear gradient to 1.8% 2-propanol over 45 min. Ionization was achieved using atmospheric pressure chemical ionization. The spectrometer was run in SIM mode (m/z 743.8, 1018, 1020, 1022, 1032, 1034, 1036, 1046, 1048, 1050, 1292.3, 1296.3, 1298.3, 1300.3, and 1302.3). Compounds were identified by comparing mass spectra and retention times with those in the literature (Hopmans et al. 2000). Quantification was achieved by integrating the peak area in the (M + H)⁺ chromatogram and comparing these with the peak area of an internal standard (C₄₆ GTGT). The analytical accuracy is 0.45 °C, determined using the global \({\text{TEX}}_{86}^{\text{H}}\) calibration by Kim et al. (2010).

3.3.3 Biomarker indices

The \({\text{U}}_{37}^{{{\text{K}}^{\prime}}}\) value was calculated using the equation proposed by Prahl and Wakeham (1987), and then converted into a temperature estimate using the calibration of Prahl et al. (1988) (Table 4). The \({\text{TEX}}_{86}\) value was calculated using the GDGT Index I (also known as \({\text{TEX}}_{86}^{\text{H}}\)) of Kim et al. (2010), and subsequently converted to an SST estimate using the global calibration of Kim et al. (2010) and subsurface temperature (subT) estimate using the global calibration of Kim et al. (2012) and the regional calibration of Ko et al. (2022). The equations are detailed in Table 4. GDGT -2/3 ratio was calculated using the equation of Taylor et al. (2013) while the Branched Isoprenoid Tetraether (BIT) index was calculated using the equation of Hopmans et al. (2004).

3.4 Statistical analyses

3.4.1 Calculation of mean paleotemperature and uncertainty for each time slice

For each core site, we computed the arithmetic mean over three time slices, i.e., Late Holocene (LH; 0–4 ka), Early Holocene (EH; 6–10 ka), and LGM (19–23 ka) for all the paleotemperature records. We reported the standard error of the mean (1 sigma) as the uncertainty for the average values in each time slice. The magnitude of glacial cooling was determined as the temperature difference between the LH and the LGM time slice unless otherwise specified. The errors assigned to glacial cooling were calculated as the root sum squares of the standard errors for the two time slices.

3.4.2 Averaging paleotemperature records

To create proxy-specific paleotemperature stacks (average of nearby records), we first binned all records of the same proxy type into selected time windows, in this case 1000-year and 2500-year. We calculated the mean of each window for each nearby record. We then averaged these core-specific window mean values to obtain the mean and the standard error of the mean for each window, assigning equal weight to each record. The error envelope for the paleotemperature stack was calculated by summing up the standard error of the mean (2 sigma) and the reported calibration error for each proxy, i.e., ± 1.2 °C for G. ruber Mg/Ca (Hastings et al. 2001) and ± 1.5 °C for \({\text{U}}_{37}^{{{\text{K}}^{\prime}}}\) (Müller et al. 1998).

3.4.3 Calculation of Mean/Median Absolute Deviation (MAD)

The mean and median absolute deviations (MADs) were used to describe the dispersion of nearby records from the same proxy type; similar records would result in a small MAD. For each proxy, all records were grouped into bins of 1000-year. Then, the MAD for each bin was calculated using the following formulas:

Here, X represents the entire dataset and X_i refers to each data point in the dataset.

Finally, the mean or median of the MADs within each time window was reported as an estimate of proxy replicability.

3.4.4 Mann–Kendall test

The statistical significance of the trend in the Holocene portion (defined here as 0–10 ka) of the continuous paleotemperature records was determined using the Mann–Kendall test. We performed this test in R using the "mk.test" function from the "trend" package (Pohlert 2023).

3.5 Simulated vertical temperature profiles

Vertical profiles of temperature along the water column at the location of our study site for the pre-industrial, mid-Holocene, and LGM periods were derived from the Stony Brook Parallel Ocean Model (sbPOM) simulations conducted by Yang et al. (2022). The data were interpolated in the vertical (z) dimension to obtain temperature estimates at 5-m depth intervals.

4.1 Comparison of nearby paleotemperature records from the same proxy type

To evaluate the similarity of nearby paleotemperature records from the same proxy type, three aspects of the proxy time series are assessed. Firstly, the temporal patterns of the time series are examined, which encompass both the Holocene (0–10 ka) trend and the general direction of glacial-interglacial change. Secondly, the mean temperature of three selected time slices is analyzed, including the LH (0–4 ka), the EH (6–10 ka), and the LGM (19–23 ka). Thirdly, the dispersion of data within the aforementioned time slices is examined, along with the temporal patterns of the absolute deviation between core pairs (reported as inter-site deviation).

We note that all nearby records of planktic foraminiferal Mg/Ca (Fig. 3a), \({\text{U}}_{37}^{{{\text{K}}^{\prime}}}\) (Fig. 3b), \({\text{TEX}}_{86}\)-SST (Fig. 3c) and \({\text{TEX}}_{86}\)-subT (Fig. 3d) indicate a warming trend across the glacial-interglacial transition.

4.1.1 Planktic foraminiferal Mg/Ca SST

Our compilation includes G. ruber Mg/Ca SST data for three out of four sites: U1429 (continuous dataset), KY07 (continuous dataset), and KR07 (time slices covering 0–10 ka and 15–25 ka). There is no significant trend in the Holocene portion of the time series at any of the sites (Table 5).

The inter-site deviations range from 0 to 4.5 °C, without any evident glacial-interglacial fluctuations (Fig. 4a). The deviations between KR07 and two other sites (KY07 and U1429) (~ 1–1.5 °C; inter-site distance of ~ 7–9 km) are, on average, larger than those between KY07 and U1429 (~ 0.5 °C; inter-site distance of ~ 6 km) (Fig. 4a). As such, the inter-site offsets appear to be somewhat related to the inter-site distances between the records (Fig. 5).

All sites agree on the mean LH temperature within their respective errors caused by temporal variability within the time slice, while an inter-site offset of ~ 1 °C exists for the EH and LGM time slices (Table 6 and Fig. 6a). The standard error (1 sigma) of each time slice ranges from 0.1 to 0.7 °C (Table 6). The Mg/Ca-inferred LGM cooling ranges from 2.6 to 4.6 °C, with the most pronounced cooling observed at site KR07 while the weakest cooling is observed at site U1429. The median and mean absolute deviation of Mg/Ca SST between all core pairs is 0.6 °C and 0.8 °C, respectively (Table 7).

4.1.2 \({\text{U}}_{37}^{{{\text{K}}^{\prime}}}\)

Our compilation includes \({\text{U}}_{37}^{{{\text{K}}^{\prime}}}\)-inferred temperature data from three sites: MD98 (continuous dataset), U1429 (continuous dataset) and KR07 (only time slice data). The inter-site deviations range from 0 to 1.2 °C, without any clear glacial-interglacial fluctuations (Fig. 4b), nor a strong relationship with inter-site distance (Fig. 5). On average, the temperature deviations between KR07 and U1429 (~ 0.6 °C; inter-site distance ~ 7 km) are larger than the other two core pairs (~ 0.2–0.4 °C; inter-site distance ~ 6–9 km) (Fig. 4b). The Holocene interval of the \({\text{U}}_{37}^{{{\text{K}}^{\prime}}}\) time series at both U1429 and MD98 sites is characterized by a significant warming trend (Table 5).

The \({\text{U}}_{37}^{{{\text{K}}^{\prime}}}\) records show remarkable agreement among sites in all aspects, including the mean temperature of LH, EH and LGM, as well as the overall patterns and LGM cooling (~ 3.5 °C) (Table 6). The standard errors of the time slices are < 0.2 °C for all sites (Table 6), while the median and mean absolute deviation of \({\text{U}}_{37}^{{{\text{K}}^{\prime}}}\) SST of all core pairs is 0.2 °C and 0.3 °C, respectively (Table 7).

4.1.3 \({\text{TEX}}_{86}\)

\({\text{TEX}}_{86}\) data are available at only two sites in our compilation: MD98 (continuous dataset) and KR07 (only time slice data). Here, we examine temperature estimates obtained using SST and subT calibrations. These calibrations yield paleotemperature estimates with different absolute values (SST > subT) but the same temporal pattern. There is a significant warming trend in the Holocene interval of the MD98 record (Table 5).

Replicate data only exist for the EH and LGM time slices. \({\text{TEX}}_{86}\) temperature data from core MD98 are 1–2 °C higher than those from KR07, but both cores show a comparable EH-LGM cooling of ~ 6 °C (Fig. 6 and Table 6). The standard error in the LGM time slice of MD98 is twice that of the other time slices in both cores. The standard error of all time slices ranges from 0.3 to 0.9 °C (Table 6). The absolute inter-site deviations between MD98 and KR07 range from 0 to 2.7 °C, with an average of ~ 2 °C for SST estimates (Fig. 4c) and 1.5 °C for subT estimates (Fig. 4d). The median and mean absolute deviations of the core pair is 0.7 °C and 0.9 °C, respectively (Table 7).

4.2 Comparison of multiproxy paleotemperature datasets at individual site

Multiproxy temperature estimates are available for three out of the four sites discussed in this study. To compare different proxy estimates at a given site, we evaluate the temporal changes during glacial-interglacial periods and then compare the paleotemperature averages for each time slice. We assume a calibration error of 1.2 °C for Mg/Ca (Hastings et al 2001), 1.5 °C for \({\text{U}}_{37}^{{{\text{K}}^{\prime}}}\) (Müller et al. 1998), and 2.5 °C for \({\text{TEX}}_{86}\)-SST (Kim et al. 2010).

All available proxy records at site KR07 (Fig. 7a), MD98 (Fig. 7b), and U1429 (Fig. 7c) display a consistent warming trend across the glacial-interglacial transition, characterized by lower values during the glacial period and higher values during the Holocene.

4.2.1 KR07

In the EH and LGM time slices where data from all proxies (G. ruber Mg/Ca, \({\text{U}}_{37}^{{{\text{K}}^{\prime}}}\) and \({\text{TEX}}_{86}\)-SST) are available, Mg/Ca temperatures are the highest, followed by \({\text{U}}_{37}^{{{\text{K}}^{\prime}}}\) and \({\text{TEX}}_{86}\)-SST (Fig. 6b and 7a; Table 6). The agreement between the proxies varies with time. Mg/Ca-\({\text{U}}_{37}^{{{\text{K}}^{\prime}}}\) agreement was better during the LGM than the Holocene, whereas the agreement between these two proxies and \({\text{TEX}}_{86}\)-SSTs was poorest during the LGM with a ~ 5 °C offset (Table 6). The EH temperature estimates from all proxies fall within the range of 22.2–25.6 °C (Table 6), agreeing within errors.

4.2.2 MD98

\({\text{TEX}}_{86}\)-SST temperature estimates were higher than \({\text{U}}_{37}^{{{\text{K}}^{\prime}}}\) temperatures during the Holocene and lower during the LGM, resulting in a ~ 5 °C difference in the LGM cooling between these proxies (8.5 °C for \({\text{TEX}}_{86}\)-SST and 3.4 °C for \({\text{U}}_{37}^{{{\text{K}}^{\prime}}}\) index; Table 6 ; Figs. 6b and 7b).

4.2.3 U1429

Mg/Ca temperature estimates are consistently higher than those of \({\text{U}}_{37}^{{{\text{K}}^{\prime}}}\) by ~ 3 °C and indicate no significant Holocene trend that is present in the \({\text{U}}_{37}^{{{\text{K}}^{\prime}}}\) record (Table 5). The LGM cooling estimates inferred from these two proxies are comparable, with ~ 3 °C for Mg/Ca and ~ 4 °C for \({\text{U}}_{37}^{{{\text{K}}^{\prime}}}\) (Table 6).

4.3 Inter-site deviation between proxy pairs

\({\text{U}}_{37}^{{{\text{K}}^{\prime}}}\)-Mg/Ca deviations for sites KR07 and U1429 fall within a similar range of 1–4 °C, but the average downcore deviation is larger at site U1429 (~ 3 °C) than at site KR07 (~ 1.5 °C) (Fig. 8a). The \({\text{U}}_{37}^{{{\text{K}}^{\prime}}}\)-Mg/Ca deviations do not exhibit a prominent glacial-interglacial shift and vary between sites, e.g., the LH decrease at site U1429 is not evident at site KR07. Meanwhile, Mg/Ca-\({\text{TEX}}_{86}\) deviations at site KR07 fall within the range of 2–7 °C with an average of ~ 5 °C (Fig. 8b), which is larger than the \({\text{U}}_{37}^{{{\text{K}}^{\prime}}}\)-Mg/Ca deviations (Fig. 8a). Mg/Ca-\({\text{TEX}}_{86}\) deviations are greater during the LGM than the EH (Fig. 8b). Lastly, \({\text{U}}_{37}^{{{\text{K}}^{\prime}}}\)-\({\text{TEX}}_{86}\) deviations fall within the range of 0–4 °C for site MD98 and 1–5 °C for site KR07 (Fig. 8c). As a whole, the \({\text{U}}_{37}^{{{\text{K}}^{\prime}}}\)-\({\text{TEX}}_{86}\) deviations are the smallest during 5–15 ka, which spans the EH time slice, and higher during the LH and LGM.

4.4 Subsurface temperature records based on Mg/Ca ratio

The thermocline temperature record derived from thermocline-dwelling N. dutertrei in core KY07 and KR07 exhibits a cooling trend from 18 to 5 ka, followed by a warming trend from 5 to 0 ka (Fig. 9a). The paleotemperatures fall within the range of ~ 12–20 °C and are characterized by large positive excursions during the Holocene interval. The coldest paleotemperatures occurred during the mid-Holocene (~ 5 ka). In contrast, the bottom water temperature record derived from benthic Uvigerina spp. shows lower temperatures during 15–20 ka and higher temperatures during the Holocene (Fig. 9b).

4.5 \({\text{TEX}}_{86}\)-subT and GDGT indices

\({\text{TEX}}_{86}\)-subT records exhibit the same temporal patterns as their counterpart \({\text{TEX}}_{86}\)-SST, but with lower values and smaller amplitudes of change (Figs. 3d, 5, 7 ; Table 6). The range of \({\text{TEX}}_{86}\)-subT is 12–22 °C in core MD98 and 12–18 °C in core KR07 (Fig. 3d), with a similar EH-LGM anomaly (~ 5 °C) at both sites (Table 6). Both cores show higher values for the Holocene than for the LGM.

GDGT-2/3 values of core MD98 fall within the range of 4–5, except for a single outlier at ~ 6 ka which reaches ~ 7 (Fig. 9c). BIT values of core MD98 fall within the range of 0.005–0.025.

The uncertainties of paleotemperature time series stem from two main aspects: proxy analysis and age model. In the following discussion, we primarily focus on the uncertainties associated with proxy analysis, while directing readers to existing literature that addresses the uncertainties related to the age model of sediment cores, including their replicability (e.g., Dolman et al. 2021; Lougheed et al. 2020; Zuhr et al. 2022). While age model related uncertainty must be considered when comparing high-resolution paleotemperature records (e.g., variability on century to decadal scale), it is less critical in our examination of inter-proxy comparison within a sediment core, which can be performed without a high-resolution age model, as well as in assessing proxy replicability concerning glacial-interglacial patterns and time slice means.

5.1 Climate signal and proxy noise in paleotemperature records of nearby sites

Due to limited resources, it is typical for paleoceanographic reconstructions to rely on the analysis of a single sediment core. This assumes that the sedimentary archive and the proxy signal contained within accurately reflect local oceanographic and climatic changes. In the northern OT region, paleotemperature records of the same proxy type from nearby sites generally align, suggesting that this straightforward assumption holds true for glacial-interglacial patterns (Figs. 3 and 4), the extent of glacial cooling (with the exception of \({\text{TEX}}_{86}\) records) (Table 6), and the Holocene trend (Table 5). This finding indicates that the relative change and average signal of large climate shifts, such as glacial-interglacial oscillations, are not substantively compromised by proxy noise. However, caution is advised when interpreting high-frequency or abrupt excursions in paleotemperatures (e.g., those occurring at ~ 9 ka and 12–14 ka in the \({\text{U}}_{37}^{{{\text{K}}^{\prime}}}\) record at site MD98) that exceed the mean or median absolute deviations (Table 7). These anomalies are not observed at nearby sites, and are therefore likely not indicative of a climate signal. These excursions are not associated with any substantial shifts in the sedimentation rate, suggesting that turbidites are unlikely to be the cause. Another possibility is instrumental issues; however, this cannot be confirmed without access to the raw instrumental data. Therefore, the cause for these abrupt excursions in the \({\text{U}}_{37}^{{{\text{K}}^{\prime}}}\) record at site MD98 remains unclear.

Among the three proxy types evaluated, \({\text{U}}_{37}^{{{\text{K}}^{\prime}}}\) demonstrates the highest level of reproducibility. \({\text{U}}_{37}^{{{\text{K}}^{\prime}}}\) data generated from sediment cores (MD98, U1429, KR07) obtained using different coring gears, sample work-up techniques, and instruments in various laboratories yield nearly identical temperature estimates, with a median absolute deviation of just ~ 0.2 °C between core pairs (Table 7) and standard errors of just ~ 0.2 °C for time slices (Table 6). The largest deviations between \({\text{U}}_{37}^{{{\text{K}}^{\prime}}}\) records analyzed in three different laboratories fall within ~ 2 °C, consistent with the result of a community-wide round-robin study (Rosell-Melé et al. 2001).

On average, the inter-site deviations in \({\text{U}}_{37}^{{{\text{K}}^{\prime}}}\) records are half as large as those observed in Mg/Ca and \({\text{TEX}}_{86}\)-SST records (Figs. 3 and 4 , Table 7). This result suggests minimal inter-laboratory offsets and implies that the spatial distribution of \({\text{U}}_{37}^{{{\text{K}}^{\prime}}}\) signal in marine sediments is highly homogeneous within a small basin such as the Danjo Basin. The \({\text{U}}_{37}^{{{\text{K}}^{\prime}}}\) index is derived from alkenones, a type of lipid biomarker produced by haptophyte algae. Alkenones are abundantly produced in the photic zone, particularly during phytoplankton blooms, and subsequently descend to the seafloor, becoming incorporated into the sediment. For \({\text{U}}_{37}^{{{\text{K}}^{\prime}}}\) data generation, the alkenones are extracted from ground sediment powder. Each measurement likely represents tens of thousands to millions of alkenone molecules, providing a reliable approximation of the sediment sample’s true mean value.

On the contrary, the sample work-up and analytical protocol for Mg/Ca differ fundamentally from that of \({\text{U}}_{37}^{{{\text{K}}^{\prime}}}\) measurement. The increased variability in Mg/Ca records is most likely attributable to the fact that the measurement is based on a relatively small number of foraminiferal tests (20–40 individuals), as opposed to a much larger number of molecules for alkenone measurement. The lifespan of mixed layer-dwelling planktic foraminifera is generally about a month. Therefore, the temperature signal recorded by individual foraminifera specimen within the same sediment sample (typically 1 cm slice) may differ by up to several degrees Celsius (Groeneveld et al. 2019). Additionally, the distribution of magnesium within foraminiferal tests is not uniform, with a temperature difference of up to 7 °C between chambers (Sadekov et al. 2008). Foraminiferal tests are cracked open to yield several fragments during the cleaning processes to improve efficacy, and this can result in sample loss and the loss of fragments. As a result, the measurement of Mg/Ca ratio may be biased either lower or higher compared to the mean value of the foraminiferal population in the sediment sample, introducing noise into the Mg/Ca record. A previous study has reported larger variability in Mg/Ca records compared to \({\text{U}}_{37}^{{{\text{K}}^{\prime}}}\) records, based on a systematic statistical analysis of global data compilation (Laepple and Huybers 2013), indicating that this finding is not unique to our study area.

Our results also show that the variations in nearby records (Fig. 3) and inter-site deviations (Fig. 4) do not systematically vary with time or inter-site distance (Fig. 5). Therefore, we suggest that the noise observed in the Mg/Ca proxy records is not related to climate variability or the presence of multiple small-scale local hydrographic processes; the latter would lead to larger deviations for core pairs that are further apart. Interestingly, the Mg/Ca record at site KR07 appears to be noisier than those at other sites. This can be observed in its relatively high variability (Fig. 3), and larger inter-site deviations that include this core (by ~ 1–1.5 °C; Fig. 4a). These differences exist despite similar sedimentation rates and mesoscale processes, such as the susceptibility to the intrusion of Chinese coastal water, at all three sites which are located within a radius of 10 km.

To examine the consistency of KR07 with other cores, we re-analyzed several selected samples from KR07 using a different cleaning protocol that includes a mild reductive cleaning step at a different laboratory in NTU (Lo et al. 2013; Shen et al. 2001). Although the newly generated data are few, these data points indicate reduced variability and align more closely with the average of three Mg/Ca records (Fig. 3). Through a one-by-one horizontal comparison of data generated at NTU and JAMSTEC, we found that, except for two data points, the data generated at NTU are in good agreement (± 0.5 °C) with the data generated at JAMSTEC (Table 8). This result suggests that the large Mg/Ca excursions in core KR07 may originate from sediment inhomogeneity due to bioturbation. The degree of bioturbation is determined by the activity of benthic organisms, whose distribution in marine sediments is highly patchy and localized. Therefore, it is possible that the degree of bioturbation is higher at site KR07 compared to other sites. This could in turn result in larger differences between KR07 and the other sites, as well as variations between some of the replicate data.

5.2 Mismatch in the direction and magnitude of change in inter-proxy offsets between sites

The multiproxy approach, which involves using multiple proxies to reconstruct the same climate variable, has become increasingly prevalent in paleoceanographic reconstructions. It is commonly understood that agreement among proxies enhances confidence in the reconstruction, while discrepancies between proxies may indicate potential limitations or variations in the proxy signal due to secondary factors. At our study sites, SST proxy records, either based on the same sediment core or those from other sites, consistently show agreement regarding the direction of glacial-interglacial change (i.e., the cold LGM and warm Holocene periods), as well as in paleotemperature values when considering calibration uncertainties in proxies (Fig. 7). These results imply that the different geochemical proxies provide a consistent view of the glacial-interglacial changes in upper ocean temperatures in the northern OT.

Nevertheless, spatially varying inter-proxy deviations (i.e., difference between proxies in a sediment core), both in terms of absolute values and temporal patterns, are observed (Figs. 7 and 8). This is evident in both the Mg/Ca-\({\text{U}}_{37}^{{{\text{K}}^{\prime}}}\) comparison between sites KR07 and U1429, as well as the \({\text{U}}_{37}^{{{\text{K}}^{\prime}}}\)-\({\text{TEX}}_{86}\)-SST comparison between sites KR07 and MD98. Such discrepancies suggest that seasonal biases in the proxy records are unlikely to be the primary cause of inter-proxy offsets as often assumed in existing literature. As discussed in Sect. 2.2 and outlined in Table 1, in this region Mg/Ca, \({\text{U}}_{37}^{{{\text{K}}^{\prime}}}\) and \({\text{TEX}}_{86}\)-SST have been interpreted as indicators of summer, annual mean, and winter SSTs, respectively. These interpretations are based on modern observations and seem reasonable when examining a single sediment core. However, if seasonal biases were present in these proxies, they should be consistent across all sites given their proximity. In the northern OT, the observed discrepancies in inter-proxy offsets among sites in terms of the direction and magnitude of change suggests that non-climatic factors, such as sample processing and analytical techniques (e.g., cleaning protocol for foraminiferal tests in Mg/Ca analysis), as well as the heterogeneity of the proxy signal in marine sediments, likely contribute to these discrepancies in inter-proxy deviations between sites.

5.3 Depth origin of \({\text{TEX}}_{86}\) signal in the northern OT

Despite the aforementioned issues (Sects. 5.1 and 5.2), one aspect of paleotemperature reconstruction that appears to be consistent is the magnitude of glacial cooling, which is reproducible across proxies, sites, and laboratories for both \({\text{U}}_{37}^{{{\text{K}}^{\prime}}}\) and Mg/Ca records. However, \({\text{TEX}}_{86}\)-SST records suggest a much more pronounced glacial cooling of 8.5 °C, approximately 2–3 times greater than that inferred from the other two SST proxies (Table 6). Similar observations have been previously documented for tropical oceans (Hertzberg et al. 2016; Ho and Laepple 2015) and across different time scales in the global ocean (Ho and Laepple 2016). These studies have postulated that the pronounced LGM cooling in the \({\text{TEX}}_{86}\) records could be due to a glacial shift in the recording depth of \({\text{TEX}}_{86}\) (Hertzberg et al. 2016) or the application of an SST calibration when the sedimentary \({\text{TEX}}_{86}\) signal originated from greater depths (Ho and Laepple 2015). These hypotheses will be further examined below.

Considering the bathymetry of the OT, it is plausible that the study sites were situated much closer to land during the LGM when the global sea level was ~ 120 m lower than present (Fig. 1). Nevertheless, the BIT (Hopmans et al. 2004) values have consistently remained low (< 0.03) over the past 25 ka (Fig. 9c), indicating minimal terrigenous influence on the GDGTs in the basin throughout the studied interval. The GDGT-2/3 ratio, which might reflect changes in the contribution of deep-dwelling archaea (Taylor et al. 2013), remained stable and lacked significant variations during the last glacial cycle (Fig. 9c), suggesting that the marked changes in the northern OT \({\text{TEX}}_{86}\)-SST records are not linked to shifts in the recording depth of \({\text{TEX}}_{86}\).

Next, we evaluate the suitability of a subsurface calibration for our study area. A recent study by Ko et al. (2022) has proposed a subsurface origin for \({\text{TEX}}_{86}\), identifying a shallow subsurface maximum range of 20–50 m. However, this subsurface range may be underestimated as the study only analyzed samples from depths ranging from 0–50 m. An earlier study, based on suspended particulate matter collected in May 2008 from the northern OT, suggested a deeper subsurface GDGT abundance maxima at around 75–100 m (Nakanishi et al. 2012). This depth range is greater than the 20–50 m alkenone abundance maxima (Ko et al. 2018) and the 50 m G. ruber abundance peak (Lee et al. 2021 and references therein). Thus, modern observations, albeit limited, support the interpretation of \({\text{TEX}}_{86}\) as being a subsurface temperature proxy.

Utilizing the regional subT calibration of Ko et al. (2022) and the widely used global calibration from Kim et al. (2012) reveals a typical glacial-interglacial temperature pattern, with colder conditions during the LGM (13–20 °C) and warmer temperatures in the Holocene (18–24 °C) (Fig. 9a). Interestingly, this glacial-interglacial temperature range is broadly consistent with that of subT records inferred from the Mg/Ca ratio of thermocline-dwelling N. dutertrei at sites KY07 and KR07, and the lower thermocline-dwelling Globoconella inflata at a neighboring site MD98-2196 (Ujiié et al. 2016) (Fig. 9a). However, these subsurface Mg/Ca records exhibit an inverse temporal trend compared to the \({\text{TEX}}_{86}\)-subT record, showing higher values during the LGM than in the Holocene, with the coldest period occurring during the mid-Holocene. This unexpected finding is notable as both N. dutertrei and GDGTs are found at similar water depths of 75–100 m in the northern OT (Table 1). In fact, the glacial-interglacial pattern of \({\text{TEX}}_{86}\)-subT records aligns more closely with that of the intermediate water temperature record derived from radiolarian assemblages (Matsuzaki et al. 2019) and the bottom water temperature record derived from benthic foraminifera Mg/Ca (Fig. 9b). This finding suggests that the depth at which \({\text{TEX}}_{86}\) signals are captured and exported in the northern OT could be below the thermocline. We note that low GDGT -2/3 values in the OT imply a shallow subsurface origin of the sedimentary GDGTs, although a quantitative interpretation of GDGT- 2/3 in this region is yet precluded by data scarcity. Therefore, our findings underscore the need for further modern studies of GDGTs and \({\text{TEX}}_{86}\) in the region using suspended particulate matter and surface sediments to refine the \({\text{TEX}}_{86}\) calibration and subsequent \({\text{TEX}}_{86}\)-based interpretations.

Alternatively, it is possible that the temporal evolution characterizing the Mg/Ca records of thermocline-dwelling planktic foraminifera is not solely influenced by the temperature evolution of the thermocline. In the northern OT, N. dutertrei thrives in the upper 200 m, and its relative vertical distribution varies with seasonal changes in the water column structure (Yamasaki et al. 2010). Based on Yamasaki et al. (2010)’s multinet study, the abundance of N. dutertrei peaks at ~ 50 m (upper thermocline) in certain months but may be more dispersed vertically (0–100 m; upper thermocline and mixed layer) when the mixed layer deepens. If we assume that this ecological behavior remains consistent at glacial-interglacial time scale and that the signal is preserved in the sedimentary record, one potential explanation is that during the LGM, N. dutertrei occupied a broader depth range, including the mixed layer, due to intensified monsoonal wind-induced mixing (Ijiri et al. 2005; Yamamoto et al. 2013). In contrast, during the Holocene, N. dutertrei may have been confined to a narrower vertical depth range, as the water column was more stratified. When the vertical distribution of N. dutertrei is narrower, one could anticipate a greater vertical temperature contrast between that recorded by N. dutertrei and a mixed layer dweller, since the N. dutertrei is capturing only thermocline temperature signal instead of a combined mixed-layer and thermocline temperature signal. This shift in habitat depth range could be a plausible factor contributing to the distinctive glacial-interglacial pattern observed in N. dutertrei. Modern observations in the northern OT suggest that G. inflata may also experience a similar shift in habitat depth (Yamasaki et al. 2010). Notably, a similar glacial-interglacial pattern can be observed in a Mg/Ca record of N. dutertrei from the North Pacific Subtropical Gyre at site ASM-5PC, which is located at a similar latitude to our study site, but outside of the OT (Sagawa et al. 2011). This finding suggests that these thermocline proxy records may reflect a common environmental or ecological signal in the region.

Our hypothesis regarding the changing recording depth of thermocline-dwelling foraminifera is further substantiated by the LGM and mid-Holocene temperature-depth profiles simulated using a high-resolution ocean model (Yang et al. 2022). The model output indicates that the temperature of the entire water column at our study site was lower during the LGM compared to the mid-Holocene and pre-industrial periods (Fig. 10 ), consistent with the glacial-interglacial shifts observed in \({\text{U}}_{37}^{{{\text{K}}^{\prime}}}\), \({\text{TEX}}_{86}\), radiolarian transfer function, and Mg/Ca ratios of surface- and bottom-dwelling foraminifera (Figs. 7 and 9 ). The relatively stable hydrographic conditions at our study site are reasonable, given the insensitivity of the flow path of the Kuroshio Current to changes in sea level, glacial-interglacial climate, and a hypothetical tectonic shoaling of the Yonaguni Depression (Lee et al. 2013; Vogt-Vincent & Mitarai 2020; Yang et al. 2022). Collectively, both climate models and other paleotemperature proxies suggest that the warming inferred from thermocline-dwelling foraminifera for the LGM likely does not reflect the glacial-interglacial temperature change at the thermocline, but rather a shoaling of its recording depth during the LGM. The shoaling of its recording depth might be explained by changing nutrient condition and resulting food availability (Kubota et al. 2017) or favorable temperature of metabolism for each species. Nonetheless, further research is warranted to enhance our understanding of the habitat depth of thermocline-dwelling foraminifera, particularly through repeated sampling at the same location to capture interannual variability in their vertical distribution. Ultimately, this information can be employed to simulate the spatio-temporal distribution of foraminifera in the past via an eco-physiological model, thereby improving the interpretation of foraminifera-based proxies.

5.4 Implications and recommendations for paleoceanographic reconstructions

As discussed in Sects. 5.1 and 5.2, the presence of noise in paleotemperature proxy records in the northern OT primarily arises from uncertainties related to sample processing and instrumental analysis, alongside with sediment heterogeneity. This finding may pertain specifically to our study’s setting, namely that the sediment cores were retrieved from the depocenter of a small basin with a relatively high sedimentation rate and a lack of strong currents. Future studies conducted in different settings would provide additional insights into the governing factors of proxy signal heterogeneity in marine sediments. Regardless of the sedimentation setting, the noise in foraminifera Mg/Ca data may be reduced by increasing the number of specimens utilized in each measurement. One way to achieve this is by analyzing an aliquot of a sample containing several tens to hundreds of specimens. Replicating the aliquot measurements will further increase the reliability of the paleotemperature estimate, particularly for samples from climate transitions (e.g., glacial termination) that are relatively susceptible to the influence of bioturbation which mixes sediments from distinct climate states.

While in certain instances inter-laboratory offsets may be consistent, this does not appear to be the case with the dataset under examination. Hence, one way to enhance the signal-to-noise ratio of paleotemperature data is through averaging or stacking nearby records of the same proxy type, which effectively minimizes random noise while preserving the shared climatic signal. A similar approach has been applied to the age modeling of marine sediments in the South China Sea based on radiocarbon data (Zuhr et al. 2022). We synthesized paleotemperature stacks of Mg/Ca and \({\text{U}}_{37}^{{{\text{K}}^{\prime}}}\) based on datasets featuring a minimum of two continuous records; \({\text{TEX}}_{86}\), for which only a singular continuous record was available, was therefore excluded from this analysis. Our approach of incorporating proxy records analyzed across different laboratories makes the derived uncertainty estimate particularly suitable for comparing regional proxy datasets and global compilations, which typically integrate data from numerous research institutions. A more robust delineation of proxy-specific uncertainty linked to spatial heterogeneity necessitates that proxy analyses (of the same proxy type) to be conducted uniformly in terms of sample processing methodologies and instrumental settings.

When considering both inter-site variability and calibration uncertainty resulting from regression as forms of proxy uncertainty, the Mg/Ca and \({\text{U}}_{37}^{{{\text{K}}^{\prime}}}\) records generally agree within error for most parts of the reconstructions (Fig. 11). This finding, coupled with the discrepancies in \({\text{U}}_{37}^{{{\text{K}}^{\prime}}}\)-Mg/Ca deviations between proximal sites (Sect. 5.2), suggests that it remains uncertain whether the difference between these two proxies truly reflects a different seasonal temperature, as posited in earlier studies. The consistent ~ 2 °C divergence between the \({\text{U}}_{37}^{{{\text{K}}^{\prime}}}\) and Mg/Ca stacks over time may as well be attributed to calibration choices (Supplementary Material). Furthermore, foraminiferal Mg/Ca ratio may also be influenced by non-thermal factors such as changes in salinity and pH (Gray and Evans 2019). Correcting Mg/Ca SST estimates to account for both glacial-interglacial variations in global surface salinity and pH increases the overall variability in the records, exacerbates inter-proxy and inter-site deviations, and yields core-top Mg/Ca temperatures that are close to the annual mean SST (Supplementary Material andFig. 1c). Further research is needed to quantify the impact of non-thermal factors on Mg/Ca ratios in this region and, crucially, how they have changed over time, in order to improve the accuracy of Mg/Ca-based paleotemperature estimates. In this regard, we note that the similar magnitude of Mg/Ca-based glacial cooling compared to that based on the \({\text{U}}_{37}^{{{\text{K}}^{\prime}}}\) index implies that the transition to a low SSS coastal environment during the LGM, as indicated by microfossil assemblages (Ijiri et al. 2005; Matsuzaki et al. 2019; Shirota et al. 2021), was not sufficiently substantial to bias the Mg/Ca SST estimates. Therefore, having multiple nearby cores may provide additional constraints for interpreting paleotemperature proxies that are independent from modern observations. While modern observations undeniably form the foundation of proxy development, modern field studies rarely extend beyond a year. Given that the proxy signals in these modern samples represent a temporal "snapshot", they may not be directly comparable to the substantially longer timescale captured by sedimentary archives (Tapia et al. 2022). We propose that conducting more surface sediment studies in the region could elucidate whether there is a seasonal bias in \({\text{U}}_{37}^{{{\text{K}}^{\prime}}}\) and Mg/Ca SST estimates.

Paleotemperature gradients are commonly computed in paleoceanographic investigations. Both vertical and horizontal temperature gradients have been used to reconstruct past changes in the KC or EAM (Kubota et al. 2017; Ujiié et al. 2016; Wang et al. 2016; Yamamoto 2009). Given the relatively large uncertainty (~ 2 °C) in \({\text{U}}_{37}^{{{\text{K}}^{\prime}}}\) and Mg/Ca estimates when inter-site variability and calibration uncertainty are considered, we recommend that the calculation of horizontal temperature gradients be based on the same type of proxy. This approach eliminates the need to account for calibration errors when interpreting relative changes in the temperature gradient. On the other hand, for vertical temperature gradient calculations, two different calibrations are inevitably utilized, such as for \({\text{U}}_{37}^{{{\text{K}}^{\prime}}}\) and \({\text{TEX}}_{86}\), or for multiple species of foraminifera, resulting in gradient estimates that are associated with larger total uncertainty. In this respect, the calibration error is generally equivalent to or greater than the uncertainty associated with spatial heterogeneity (Figs. 6 and 11 , Table 6). Consequently, besides characterizing the spatial representativeness of a paleotemperature proxy, addressing proxy calibration errors would constitute a significant advancement toward achieving a quantitative paleotemperature reconstruction.

Finally, as demonstrated in Sect. 5.3, climate models have the potential to serve as an additional constraint for the interpretation of proxy data, particularly in instances of proxy discrepancy. In recent years, climate model outputs have been utilized to test assumptions regarding the proxy recording process, including recording depth and season (e.g., Ho & Laepple 2015, 2016; Lohmann et al. 2013), as well as the effect of bioturbation (Bienzobas Montávez et al. 2024; Dolman & Laepple 2018). Neither climate models nor proxies perfectly represent the climate system. Therefore, similar to the advantages of a multiproxy approach, employing both proxy data and climate models in reconstructions has the potential to enhance the robustness compared to methods that rely on either one alone.

In this study, we analyzed multiproxy (\({\text{U}}_{37}^{{{\text{K}}^{\prime}}}\), foraminiferal Mg/Ca and \({\text{TEX}}_{86}\)) paleotemperature data obtained from four sediment cores collected within a 10-km radius in the northern OT. Both the \({\text{U}}_{37}^{{{\text{K}}^{\prime}}}\) and Mg/Ca data suggest a glacial cooling of ~ 3 °C in the surface ocean. In contrast, \({\text{TEX}}_{86}\)-SST data suggest a glacial cooling of ~ 8 °C or ~ 6 °C, either at surface or subsurface, depending on the interpretation of the proxy. The disparity in cooling observed in the \({\text{TEX}}_{86}\) data may stem from calibration issues. Among all proxies, \({\text{U}}_{37}^{{{\text{K}}^{\prime}}}\) exhibits the highest level of replicability at nearby sites, whereas \({\text{TEX}}_{86}\) exhibits the lowest. Inconsistencies in the Mg/Ca-\({\text{U}}_{37}^{{{\text{K}}^{\prime}}}\) discrepancies among nearby sites further complicate the interpretation of these proxies as indicators of different seasonal signals. Notably, the mean absolute deviations between neighboring cores within the same proxy type are less than 1 °C, indicating precision levels that surpass the corresponding proxy calibration error. Consequently, future research endeavors should prioritize efforts to minimize calibration errors.

The datasets supporting the conclusions of this article are included as an Excel file in the supplementary material of this manuscript.

BIT:

Branched Isoprenoid Tetraether

CDW:

Changjiang Diluted Water

EAM:

East Asian Monsoon

EASM:

East Asian Summer Monsoon

ECS:

East China Sea

EH:

Early Holocene

GDGT:

Glycerol dialkyl glycerol tetraether

HPLC–MS:

High-performance liquid chromatography-mass spectrometry

KC:

Kuroshio Current

LGM:

Last Glacial Maximum

LH:

Late Holocene

MAD:

Mean/median absolute deviation

OT:

Okinawa Trough

SF-ICP-MS:

Sector-field inductively coupled plasma mass spectrometer

SSS:

Sea surface salinity

SST:

Sea surface temperature

subT:

Subsurface temperature

TSWC:

Tsushima Warm Current

UHPLC-MS:

Ultrahigh-performance liquid chromatography-mass spectrometry

We are grateful to Xinyu Guo and Haiyan Yang from Ehime University, who generously provided the vertical temperature profiles from their previously published sbPOM simulations. We would also like to thank Wei-Zhu Chen for processing biomarker samples from KR07, Raul Tapia (Institute of Oceanography, NTU) and Elena Dominguez (Institute of Oceanography, NTU) for assistance with the taxonomic identification of planktic and benthic foraminifera, Pei-Ting Lee (Institute of Oceanography, NTU), Hsun-Ming Hu (Department of Geosciences, NTU), and Chun-Yuan Huang (Department of Geosciences, NTU) for technical assistance with trace element analysis of foraminifera. We also thank K. Kimoto, T. Omura, H. Yamamoto, M. Takada, M. Sato, and Y. Nakamura for their assistance in our experiments at JAMSTEC. The manuscript benefited from constructive criticisms from Alicia Hou (Institute of Oceanography, NTU), Raul Tapia (Institute of Oceanography, NTU) and Pierrick Fenies (Institute of Oceanography, NTU), and three anonymous reviewers solicited by PEPS. We also extend our gratitude to the handling editor and chief editors Yasufumi Iryu and Ryuji Tada for handling our manuscript.

This study was supported by grants awarded to SLH from the National Science and Technology Council of Taiwan (NSTC), namely 112–2611-M-002–024-, 111–2611-M-002–027- and 110–2611-M-002–028-. Trace element analysis of foraminifera conducted at National Taiwan University was partially supported by grants from the National Taiwan University (NTU) Core Consortiums Project (112L894202), Higher Education Sprout Project of the Ministry of Education (112L901001) and the projects under the NSTC (111–2116-M-002–022-MY3). Trace element analysis of foraminifera conducted at JAMSTEC was partially supported by JAMSTEC and by a project "Rebuilding of Micropaleontological Reference Center collections" at the National Museum of Nature and Science, Japan.

Authors and Affiliations

1. Institute of Oceanography, National Taiwan University, No. 1, Sec. 4, Roosevelt Road, Da’an District, Taipei, 10617, Taiwan R.O.C.

   Ru-Yun Tung & Sze Ling Ho

2. Department of Geology and Paleontology, National Museum of Nature and Science, 4-1-1, Amakubo, Ibaraki, Tsukuba, 305-0005, Japan

   Yoshimi Kubota

3. Faculty of Environmental Earth Science, Hokkaido University, Kita-10, Nishi-5, Kita-Ku, Sapporo, 060-0810, Japan

   Masanobu Yamamoto

4. Alfred-Wegener-Institute, Helmholtz-Center for Polar and Marine Research, Am Handelshafen 12, 27570, Bremerhaven, Germany

   Jens Hefter

5. High-Precision Mass Spectrometry and Environmental Change Laboratory (HISPEC), Department of Geoscience, National Taiwan University, No. 1, Sec. 4, Roosevelt Road, Da’an District, Taipei, 10617, Taiwan R.O.C.

   Chuan-Chou Shen

6. Research Center for Future Earth, National Taiwan University, No. 1, Sec. 4, Roosevelt Road, Da’an District, Taipei, 10617, Taiwan R.O.C.

   Chuan-Chou Shen

Contributions

SLH proposed, conceived and designed the study. R-YT, MY, JH, and YK analyzed samples and wrote the methodology descriptions of their respective analyses. C-CS provided access to ICP-MS at NTU and ensured quality control for trace element analysis. R-YT compiled the data and performed statistical analyses under the guidance of SLH. SLH and R-YT interpreted the findings and produced the initial draft. All authors reviewed, revised and approved the final manuscript.

Corresponding author

Correspondence to Sze Ling Ho.

Competing interests

The authors declare that they have no competing interest.

Publisher's Note

Springer Nature remains neutral with regard to jurisdictional claims in published maps and institutional affiliations.

Open Access This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if changes were made. The images or other third party material in this article are included in the article's Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article's Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this licence, visit http://creativecommons.org/licenses/by/4.0/.

Reprints and permissions