Introduction

As the greatest source of water vapor and heat on our planet, the western Pacific warm pool (WPWP) plays an important role in the global climate system1. Modern observations show that the thermocline depth of the WPWP is closely associated with the El Niño/Southern Oscillation (ENSO) evolution2 and large-scale changes in extratropical atmospheric circulation1. In turn, subsurface waters of the western equatorial Pacific (WEP) originate from the subduction of southern Pacific subtropical waters3,4,5,6,7,8, presenting extratropical control mainly through oceanic tunnels9. Therefore, historical records of subsurface waters could provide clues on the past ENSO-like process and the linkage of extratropical climate change to the tropical Pacific10.

However, sparse and limited paleoceanographic reconstructions have presented diverse histories of subsurface waters in the WPWP on orbital timescales. Two precession-dominated subsurface seawater temperature (subT) records off southern Sumatra11 and Mindanao Island12 in the far western part of the WPWP were controlled by local monsoon-induced oceanic upwelling processes responding to local summer precessional insolation changes11. Another newly published subT record from the WEP showed a dominant half-precession cycle under bihemispheric influences13. A subT record off Papua New Guinea demonstrated clear obliquity cycles since 110 ka, with a warming subT caused by weakened subtropical gyre circulation due to high obliquity14. Overall, most existing records of subsurface waters of the WPWP have focused on either the last glacial period10 or regions far from the center of the WPWP12,14 and have shown nonuniform information.

Here, we present new paired 360,000-year temperature and salinity records of subsurface water from the thermocline-dwelling planktonic foraminifera species Neogloboquadrina dutertrei in the central WPWP (core KX97322-4 in Fig. 1) combined with our previously published sea surface temperature (SST) and salinity (SSS) records from the mixed-layer species Globigerinoides ruber15, enabling further investigation of the orbital-scale upper ocean thermal/salinity structures of the WPWP and their linkage to the southern extratropical oceans and ENSO-like processes over the last four glacial–interglacial cycles.

Fig. 1: Core location of KX97322-4 and the distributions of modern temperature and salinity in the equatorial Pacific.
figure 1

The main figure shows temperature (color shading) and salinity (black contour) profiles from 40°S to 40°N along 160°E for 0–400 m water depth. Dashed arrows indicate Subtropical Lower Water (SLW), Subtropical Mode Water (STMW), and Subantarctic Mode Water (SAMW). The inset map shows the modern temperature distribution at 150 m water depth (color shading) based on the World Ocean Atlas 2013 dataset20 using Ocean Data View software25. The dashed 28 °C SST isotherm indicates the extent of the western Pacific warm pool (WPWP). The locations of core KX97322-4 (00°02′S, 159°15′E) from this study and previously published cores MD06-3018 (23°00′S, 166°09′E)28 and MD97-2125 (22°34′S, 161°44′E)27 from the southwestern Pacific and TR163-19 (2°15.5′N, 90°57.1′W)67 and ODP 1240 (0°01.3′N, 86°27.8′W)55 from the eastern equatorial Pacific (EEP) are shown.

Full size image

Results and discussion

Subsurface temperatures and salinities in the central WPWP over the past 360,000 years

Multiple plankton tows have demonstrated that the apparent calcification depth of N. dutertrei is 130 m on average off Papua New Guinea16 and 140–150 m in the eastern part of the WPWP within the upper thermocline17. Rongstad et al.18 suggested that the N. dutertrei distribution depth in the WEP has a mean of 150 m because the much deeper mixed layer and lower productivity allow warm water and light to penetrate deeper depths in the ocean. Thus, a water depth of 150 m is adopted as the calcification depth of N. dutertrei for relative simulation experiments in this study.

The shell Mg/Ca record of N. dutertrei fluctuates between 1.2 and 2.6 mmol mol−1 (Supplementary Fig. 1b), corresponding to subT variations between 18.2 and 26.0 °C since 360 ka (Fig. 2b), which is reasonable in terms of the temperature tolerance (15–32 °C) of N. dutertrei19. The average core-top Mg/Ca-derived subT of 23.0 ± 1.1 °C also conforms to modern observations20 at the calcification depth (~150 m) of N. dutertrei (Fig. 1). Furthermore, the difference in Mg/Ca of G. ruber between the Holocene and Last Glacial Maximum is ~1.1 mmol mol−1 (Supplementary Fig. 1a), comparable to the other relatively high sedimentation rate records (1.1–1.3 mmol mol−1) in the WEP21, implying that our temperature signals from core KX97322-4 are not considerably dampened by the effect of bioturbation.

Fig. 2: Variations in surface and subsurface temperatures and δ18Osw of core KX97322-4 from this study and previously published surface temperature and δ18Osw records from the southwestern Pacific over the last 360 kyr.
figure 2

a Mg/Ca-SST of G. ruber; b Mg/Ca-subT of N. dutertrei; c the vertical temperature difference between SST and subT (ΔTG-N); d G. ruber Mg/Ca-SST of cores MD97-2125 (ref. 27) (pink) and MD06-3018 (ref. 28) (red); e δ18Osw records of G. ruber of MD97-212527 (sky blue) and MD06-301828 (blue); f surface δ18Osw record from G. ruber; g subsurface δ18Osw record from N. dutertrei; h the salinity difference between subsurface and surface water (ΔSN-G); shaded blue bars indicate glacial intervals and dashed gray curves indicate polynomial fitting of glacial–interglacial cycles; i spectral analysis of SST and subT of core KX97322-4; j spectral analysis of SSS and subS; dashed curves represent the 95% confidence level. Shaded bars represent the 100-kyr and 23-kyr signals. All data are shown in Supplementary Data 1.

Full size image

As shown in Fig. 2, the amplitude of subT (~8 °C) is much larger than the range of SST (25.5–30.5 °C) from G. ruber over the past 360 kyr (Fig. 2a, b) and exhibits 100-kyr eccentricity and 23 (19)-kyr precession signals (Fig. 2i). In contrast, the SST record of core KX97322-4 is dominated by 100- and 30-kyr cycles, followed by notable 42- and 23-kyr signals (Fig. 2i). This discrepancy suggests different controls on the SST and subT in the central WPWP. Here, orbital-scale SST variations are forced by both local signals (i.e., low-latitude insolation) and global background climate signals (i.e., ice sheets, greenhouse gas concentrations)22. However, subsurface waters in the WEP primarily record low-latitude signals because of their austral origin via the shallow meridional overturning cell (SMOC)3,23. Although the boreal subtropical gyre could also influence the equatorial thermocline through the western boundary current, it is not discussed here due to its very small contribution attributed to the blocking effect of the Intertropical Convergence Zone (ITCZ)8,24.

The subsurface seawater δ18O (δ18Osw) deciphered from N. dutertrei shares similar features with the surface δ18Osw in terms of prominent glacial–interglacial variations, although the amplitudes of the subsurface δ18Osw (average of 2.31‰) are much larger than those of the surface δ18Osw (average of 1.42‰) (Fig. 2f, g). Furthermore, the average core-top SSS and subsurface seawater salinity (subS) calculated based on δ18Osw-iv values (see “Methods”) are 35.15 ± 0.22 and 35.18 ± 0.32 (Supplementary Fig. 1e, f), respectively, similar to the in situ salinities of this region at 25 m (34.78) and 150 m (35.10)25 depths. The spectral analyses show that subS exhibits 100- and 23 (19)-kyr signals (Fig. 2j), while SSS exhibits other cycles similar to SST.

Connection between WPWP subsurface waters and the subtropical southern Pacific Ocean

The fluctuations in subT and subS dominate the variations in ΔTG-N (SSTG. ruber − subTN. dutertrei) and ΔSN-G (subSN. dutertrei − SSSG. ruber), highlighting the prominent roles of subsurface waters in the WPWP on the upper water column structure. Modern observation3,4,8 and modeling studies5,6,7 support that the more saline subsurface water in the modern WEP mainly originates from the subduction of eastern6,7 and western4,7 south Pacific subtropical waters transported by the SMOC (Fig. 1). Based on an analogy to modern climate and some simulations on orbital-scale10,13,26, changes in the southern extratropical Pacific may account for the large subT and subS variations in the WPWP over the past 360 kyr on orbital timescales.

To identify the orbital-scale extratropical control on the WPWP, we compared our subsurface temperature and δ18Osw records with published surface temperature and δ18Osw records from sites MD97-2125 (ref. 27) and MD06-3018 (ref. 28) in the southwestern Pacific. δ18Osw represents the salinity signals as a passive tracer for the SMOC5. The two selected cores are located near the northern edge of the subtropical subduction region in the southwestern Pacific29 and on the route of the SMOC4 (Fig. 1), possessing the only available surface water records suitable for comparison. Additionally, we simulated the subT variation (150 m) at our core site as well as the SST variations at two main subduction zones in the southeastern Pacific (20°00′S, 120°00′E) and southwestern Pacific (23°00′S, 166°00′W) over the past 300 kyr (Supplementary Fig. 2) using a community Earth system model (CESM). Here, we present the simulated temperatures in October because subduction is stronger in austral late winter30. Our results show that the variations in subT of core KX97322-4 are quite similar to those in the G. ruber Mg/Ca-based SST (~22–27 °C) from sites MD97-2125 (ref. 27) and MD06-3018 (ref. 28) (Fig. 2b, d), although the SST from the southern subtropics exhibited only modest precession cycles28, probably due to their low resolution and/or other complex factors influencing the surface water in that region27,28. However, the simulated results reveal that the SSTs of the two main subduction regions are dominated by precession cycles and exhibit almost synchronous variations with our core over the past 300 kyr (Supplementary Fig. 2). Moreover, there are coherent patterns and similar amplitudes between the N. dutertrei δ18Osw record of KX97322-4 and the G. ruber δ18Osw records from MD97-2125 and MD06-3018 over the last 360 kyr (Fig. 2e, g). Overall, the results reinforce that the subT and subS signals in the WPWP probably propagated from the subtropical surface to the tropical subsurface through the lower branch of the SMOC.

We also characterize the temperatures and δ18Osw-iv-calculated salinities from cores KX97322-4, MD06-3018, and MD97-2125 in glacial (precession maximum) and interglacial (precession minimum) T–S plots (Fig. 3). The data generally fall between the 23.5 and 24.5 kg m−3 isopycnals, and we conclude that the subtropical surface water subducted along isopycnal lines to the equatorial subsurface and consistently influenced the variations in equatorial waters over the last 360 kyr. Additionally, the source of waters that were ventilated in the WPWP included the Subantarctic Mode Water, Subtropical Mode Water, and Subtropical Lower Water31 (Fig. 1). These waters passed mainly through the interior channel32 (with perhaps a small part through the western boundary) and mixed with the waters along the route before reaching the equatorial subsurface4,31. Based on the predominant range of the water densities (23.5–24.5 kg m−3) (Fig. 3), we argue that the Subtropical Lower Water (with a density between 24.4 and 24.8 kg m−3) and Subtropical Mode Water (with a density between 24.5 and 25.8 kg m−3)31 are the primary source waters contributing to the subsurface waters around our core site.

Fig. 3: T–S plot from the southwestern Pacific to the western equatorial Pacific.
figure 3

a Interglacial periods and b glacial periods. Glacial–interglacial (GIG)-1, GIG-2, GIG-3, and GIG-4 represent marine isotope stage (MIS) 1/2, MIS 5/6, MIS 7/8, and MIS 9/10, respectively. The four interglacials are defined as 0–5, 125–130, 240–243, and 333–337 kyr in accordance with the periods of precession maxima, and the four glacials are defined as 19–23, 132–139, 252–256, and 338–343 kyr in accordance with the periods of precession minima. Diagonals represent isopycnal surfaces (kg m−3). The colored shading represents the main distribution range of the water temperature–salinity data at each core site, and the closed contour lines are defined as the median of two boundary points.

Full size image

Connection between the WPWP subsurface waters and ENSO-like processes in the precession band

Modern observations demonstrate a strong linkage between ΔT25–150m (vertical temperature difference between 25 and 150 m water depth) at site KX97322-4 and the Niño 3.4 index, with a smaller (4.27 ± 1.27 °C) ΔT25–150m during a La Niña event and larger (8.15 ± 1.25 °C) ΔT25–150m during an El Niño event (Supplementary Fig. 3b) at decadal timescales, which is controlled by the deeper (shallower) thermocline and warmer (cooler) subsurface temperature at the core site, while seawater temperature at 25 m depth changes little (Supplementary Fig. 3b). Since the changes in downcore foraminiferal Mg/Ca-derived ΔTG-N (corresponds to ΔT25-150m) can be applied to reflect the orbital-scale variations in the depth of the thermocline assuming a constant calcification depth of foraminifera33, a lower (higher) ΔTG-N represents a deeper (shallower) thermocline and thus weaker (stronger) upper-ocean stratification. After the background glacial–interglacial (GIG) climate signal (dashed gray line in Fig. 2c) is removed, we compare the GIG-detrended ΔTG-N data (Fig. 4e) of core KX97322-4 with the CESM34 simulated Niño 3 SST of December–January–February (DJF) (Fig. 4i) since 300 ka, and the SST of DJF could sufficiently reflect the changes in paleo-ENSO in the precession band35. Our results show a high positive correlation between the GIG-detrended ΔTG-N of core KX97322-4 and the model-simulated Niño 3 DJF SST in the precession band (Fig. 4e, i). The minimum (maximum) ΔTG-N and lower (higher) Niño 3 SST both indicate a La Niña (El Niño)-like state, in accordance with the modern observed relationship36 (Supplementary Fig. 3b), corresponding to warmer subT (Fig. 4d and Supplementary Fig. 3b) during La Niña (and La Niña-like) conditions and vice versa, which is also supported by our simulated subT result (150 m depth) over the past 300 kyr (Fig. 4c). Moreover, cross-spectrum analyses show high correlations among GIG-detrended subT, ΔTG-N, and ΔSN-G of our core and simulated subT data and Niño 3 SST (Supplementary Fig. 4) in the precession band.

Fig. 4: Comparisons of the upper water column in the central WPWP with the isolation and ENSO-like processes in the precession band.
figure 4

a Climatic precession and insolation at 23°S on June 21; b insolation difference between 0° and 23°S, as well as between the austral summer and winter; c CESM-simulated ocean temperature at 150 m water depth in October at the core site of KX97322-4; d glacial-interglacial (GIG) detrended subT from N. dutertrei Mg/Ca from this study; e GIG-detrended vertical ΔTG-N and f ΔSN-G from this study; g SST difference between the WEP (core KX97322-415) and EEP (core TR163−19 (ref. 67)); h comparison of the thermocline depth between the western equatorial (KX97322-4) and eastern equatorial Pacific (ODP 1240 (ref. 55)); the original age model of ODP 1240 was tuned to LR04 for consistency using the Match 2.3.1 software57; i CESM-simulated Niño3 SST of DJF. The bold curves in c–i are the corresponding bandpass filters in the precession band. Upper blue bars indicate glacial intervals, and green bars indicate the periods when precession is less than zero. All data are shown in Supplementary Data 1.

Full size image

In addition to the temperatures, modern seawater salinities of upper water in the study area also respond to ENSO activity. The ΔS150–25m (salinity difference between 150 and 25 m water depth) at site KX97322-4 decreased (increased) during La Niña (El Niño) periods on decadal timescales, ascribed to the increased (decreased) SSS, while S150m remained nearly stable (Supplementary Fig. 3a). An in situ study showed that the SSS at our core site decreased during El Niño events owing to zonal salt advection and net freshwater flux37. However, in the precession band, the GIG-detrended δ18Osw difference between subsurface and surface δ18Osw values (Δδ18Osw) deduced from N. dutertrei and G. ruber behaved oppositely from the modern ΔS150–25m, showing increased ΔS during La Niña-like states (Fig. 4f, i), which probably resulted from the more saline subsurface water sourced from the southern Pacific, implying that the SMOC-influenced variations in subS contributed more than the precipitation-induced SSS changes at KX97322-4 site in the precession band. At the precession minimum (Pmin), enhanced evaporation during the austral winter led to increased salinity in the subtropical surface water38, thereby reinforcing subduction and propagating equatorward.

In addition, the hydrological changes in the Indo-Pacific warm pool also respond to ITCZ and monsoon variations, which also have precession cycles39,40. However, the ITCZ mainly works on the surface water cycle (i.e., SSS) north of the equator24 and does not have a potent effect on austral tropical-subtropical interactions. Long-term precipitation records in the central WPWP were more influenced by obliquity than by precession-paced ITCZ and associated Australian summer monsoon variations40. Coupled with the thick mixed layer in the central WPWP prohibiting vertical exchange with surface water3, we argue that the subsurface water at our core site is mainly influenced by the subduction of extratropical surface water by the SMOC.

Linkage of WEP subsurface water to global climate in the precession band

In the modern climate, the subtropical gyre has a substantial influence on the ENSO process6,7. On longer timescales, the subtropical gyre may also regulate ENSO-like variations in the tropics and further influence global climate change. A numerical simulation has shown that the modulation of the meridional radiation gradient related to changes in the Earth’s obliquity leads to modification of the subtropical gyre41. In addition, a subT record off Papua New Guinea demonstrated an obliquity cycle with subT warming caused by diminished subduction of relatively cool waters by the weakened subtropical gyre circulation during high obliquity14. The different signals from ours might be due to the different information documented by different core locations. Another subT record in the WEP exhibited a strong half-precession cycle ascribed to the bihemispheric influences at its location13. Our results agree with other numerical studies that have proposed that Pmin corresponds to enhanced Hadley circulation caused by the occurrence of aphelion near the boreal summer solstice42,43. Precession-paced solar insolation dominates low-latitude climate change via coupled ocean–atmosphere processes44,45. Precipitation records in the Indo-Pacific warm pool generally display dominant precession cycles, as well as productivity variations in the equatorial Pacific and Indian oceans, which have been deduced to be controlled by ENSO-like oscillations at the precession band46,47,48.

The thermal energy of the tropical and subtropical Pacific is primarily obtained from local insolation44. The local insolation at 23°S reaches its maximum in austral winter, corresponding to the climatic precession minimum (Fig. 4a). On the other hand, the minimum insolation differences between summer and winter at 23°S (Fig. 4b) indicate prolonged summers and warm winters in the Southern Hemisphere at Pmin. Both conditions should induce a warm anomaly in subtropical surface waters. This scenario is further supported by the simulated result showing the temperature and salinity differences between Pmin and Pmax (precession maximum) (Fig. 5a), the meridional cross-section in the WPWP show that a warm and salty anomaly could be traced from the subsurface in the WPWP to its austral source during Pmin. In line with the modeling results for the early Holocene (Pmin), warm subtropical water heats the region of Subantarctic Mode Water formation, and this warm anomaly is then transported to the tropical subsurface49, determined by local austral late winter insolation forcing26.

Fig. 5: Thermal state and hypothetical meridional and zonal ocean–atmospheric circulation in the Pacific under precession minimum (Pmin) and maximum (Pmax) conditions.
figure 5

a Cross-section of CESM-simulated annual mean temperature (color shading, °C) and salinity (contour lines) differences between precession minimum and precession maximum in the WPWP along 160°E; b composite difference (°C) of temperature along the 24–25 σ isopycnal surface during Pmin relative to the 300-kyr climatological mean from CESM simulation; c and d show the Pacific ocean–atmospheric circulation during Pmin and Pmax, respectively. The ENSO-like settings refer to the La Niña event in 2010 (c) and the El Niño event in 2015 (d), which were collected and made freely available by the International Argo Program and the national programs that contributed to it (http://www.argo.ucsd.edu, http://argo.jcommops.org). SEC and EUC represent the southern equatorial current and equatorial undercurrent, respectively. Color shading represents water temperature (°C), and the thickness of lines represents the strength of circulations.

Full size image

In addition, the largest insolation gradients between the equator and 23°S occur in the austral winter at Pmin50 (Fig. 4b), strengthening the Hadley circulation and intensifying the SMOC by enhancing the surface wind strength (Fig. 5c) and then reinforcing the transfer of warm and salty subtropical water to the equatorial subsurface through isopycnals (Fig. 5b). Furthermore, the enhanced ΔS at Pmin expedites surface warming by preventing downward heat transfer through the thermal barrier effect of a relatively shallow halocline3, causing heat storage in the WEP10. This process enhances the SST difference between the WEP and eastern equatorial Pacific (EEP) (Fig. 4g) and leads to a much deeper thermocline in the WEP than in the EEP (Fig. 4h, i), promoting a prevailing La Niña-like state (Supplementary Fig. 5a, b). Conversely, Pmax is dominated by an El Niño-like state (Supplementary Fig. 5c, d). Then, the heated WEP again delivers heat to the subtropical surface through the western boundary current, accomplishing positive feedback. Moreover, WEP warming strengthens the poleward Hadley circulation51 and encourages more heat to be transferred to high latitudes, thus accelerating the melting of sea ice52 and continental ice sheets through atmospheric bridges53. The associated shoaling of the thermocline in the EEP during a La Niña-like state additionally facilitates CO2 release into the atmosphere54. Phase analysis demonstrates that the peak subT in the WPWP is almost in phase with the minimum climatic precession (−8.1 ± 9.1°) and precedes the peak western equatorial SST (38.2 ± 1.0°), as well as the ENSO-related SST (101.3 ± 24.7°) and thermocline gradient between the WEP and EEP (94.1 ± 11.5°) in the precession band (Supplementary Fig. 6), emphasizing the pivotal role of subsurface waters in the WPWP.

It is worth mentioning that the last four deglaciations are invariably accompanied by a La Niña-like state (Fig. 4) in this study, which is consistent with the reconstructed EEP upwelling intensity55 and equatorial Pacific Ocean productivity records48. During deglaciations, as insolation increases, low-latitude warming is reinforced by interactions between the equatorial and Southern Hemisphere subtropical regions. The orbital configuration of an obliquity maximum and a precession minimum would accelerate the melting of high-latitude glaciers. Weaker stratification in the equatorial Pacific and sea ice retreat would also induce a greater release of CO2 to the atmosphere, providing positive feedback to global warming.

Concluding remarks

Our paleoceanographic and modeling results demonstrate that the subsurface temperature and salinity variations in the WPWP over the last 360,000 years are governed by precession-paced insolation that modulated subtropical temperature signals and SMOC variations. As an intermediary between the equatorial and Southern Hemisphere subtropical regions, the subsurface waters in the WPWP superimpose precession-driven subtropical climate variations on the equatorial region, generating a strong influence on ENSO-like processes. During precession minimum (deglaciation) phases, warming of the equatorial subsurface and a shoaling halocline amplify the WPWP warming, inducing a La Niña-like state, and the shoaled thermocline in the EEP would promote increased degassing of CO2 to the atmosphere, thus accelerating global warming.

Methods

Core KX97322-4 (00°01.73′S, 159°14.66′E, 2362 m) was retrieved from the Ontong Java Plateau in the central WPWP using a giant piston corer on the Science-1 vessel during the Warm Pool Subject Cruise executed in 2008. The sedimentation rate is 0.39–4.95 cm kyr−1 (Supplementary Fig. 1j), with an average time resolution of 0.57 kyr cm−1.

Age model

The age model was established based on downcore stable oxygen isotope measurements on the benthic foraminifer Cibicidoides wuellerstorfi (>500 μm)15 correlated with the reference benthic stack LR04 (ref. 56) using Match 2.3.1 software57, combined with five accelerator mass spectrometry radiocarbon (14C) dates on the planktonic foraminifer Trilobatus sacculifer (with a sac) (350–500 μm) measured at the National Ocean Sciences Accelerator Mass Spectrometry facility, Woods Hole Oceanographic Institution, USA (Supplementary Fig. 1g, h).

SubT and subS calculations

Approximately 30 tests of N. dutertrei (300–355 μm) were picked and crushed and then split into two parts for Mg/Ca and δ18O analyses for each sample. Samples for Mg/Ca analysis were cleaned following the treatment procedures without a reductive step58 because the core was located far from the continent15. Mg/Ca ratios were analyzed by inductively coupled plasma-optical emission spectrometry (ICP-OES, Thermo Fisher, iCAP6300 radial) at the Institute of Oceanology, Chinese Academy of Sciences with an analytical precision of 0.44% (1σ, relative standard deviation). We chose the calibration equation that was established for N. dutertrei from the WPWP with an average standard error of 0.5 °C16 to calculate subT as follows:

$${{{{{\rm{Mg}}}}}}/{{{{{\rm{Ca}}}}}}({{{{{{{{\rm{mmol}}}}}}\;{{{{{\rm{mol}}}}}}}}}^{-1}) = 0.21\,{{\exp }}[0.097 \times {{{{{\rm{subT}}}}}} ({\,}^\circ {{{{{\rm{C}}}}}})]$$
(1)

To better compare with the Mg/Ca ratio of G. ruber from our published work on the same core15, we also recalculated the SST using the equation for G. ruber developed by Hollstein et al.16 as follows:

$${{{{{\rm{Mg}}}}}}/{{{{{\rm{Ca}}}}}}({{{{{{{{\rm{mmol}}}}}}\;{{{{{\rm{mol}}}}}}}}}^{-1}) = 0.26\,\exp [0.097 \times {{{{{\rm{SST}}}}}}({\,}^\circ {{{{{\rm{C}}}}}})]$$
(2)

The Fe/Ca, Al/Ca, and Mn/Ca analyses indicated no substantial contamination (Supplementary Figs. 8 and 9). In addition, we measured the test weights of G. ruber and N. dutertrei to assess the dissolution effect on Mg/Ca values since Termination II (Supplementary Figs. 10 and 11). Moreover, the influence of salinity on the Mg/Ca ratio of the foraminiferal shell can be ignored in the WPWP16.

The tests for δ18O measurements were bathed in 3–6% H2O2 for 30 min followed by a 30 s sonication and rinsed with acetone and ultrapure water. Then, the δ18O of N. dutertrei was measured on a GV IsoPrime mass spectrometer at IOCAS calibrated to the Vienna Pee Dee Belemnite via the NBS18 carbonate standard, with an analytical precision of 0.06‰ (1σ). Then, the subsurface δ18Osw value was calculated with the equation as follows:

$${{{\delta }}}^{18}{{{{{{\rm{O}}}}}}}_{{{{{{\rm{sw}}}}}}}({{{{{\rm{SMOW}}}}}})=({{{T}}}_{{{{{{\rm{Mg}}}}}}/{{{{{\rm{Ca}}}}}}}-10.5+6.58\times {{{\delta }}}^{18}{{{{{\rm{O}}}}}})/6.58+0.27$$
(3)

which established for N. dutertrei obtained from the Indo-Pacific warm pool based on plankton tow and surface sediment samples and laboratory culture results59. Subsurface salinity values were calculated by the relationship as follows:

$${{{{{{\rm{\delta }}}}}}}^{18}{{{{{{\rm{O}}}}}}}_{{{{{{\rm{sw}}}}}}-{{{{{\rm{iv}}}}}}}({{{{{\rm{SMOW}}}}}})=0.74\times {{{{{\rm{Salinity}}}}}}-25.4({R}^{2}=0.9)$$
(4)

which generated for WEP intermediate water with an uncertainty of the absolute salinity better than ±0.69 (ref. 60). We selected this equation to avoid the potential influence of variable surface waters in other equations, and δ18Osw-iv is the ice-volume corrected δ18Osw value using the mean ocean water δ18O value61.

Spectral and phase analysis

We used redfit spectral analysis with Past 3.20 software62 to investigate the potential mechanism for subsurface temperature and salinity variations over the past 360 kyr. We also used Redfit-X_1.1 (ref. 63) to analyze the phase relations between the precession index and other proxies.

Detrended fluctuation and filtering analysis

We used a ninth-order polynomial fit to obtain the glacial-interglacial trend of the proxies and then subtracted the trend from the original data to obtain localized signals. We also used fast Fourier transform filters to perform bandpass filtering of the proxy data with a cutoff frequency between 0.04 and 0.046 to emphasize cyclic variations in the 23-kyr precession band.

Numerical simulation

To model tropical and subtropical temperature variations, we used the CESM with T31_gx3v7 resolution (3.75° × 3.75° resolution for the atmosphere and nominal 3° resolution for the ocean) to extract the Niño 3 SST in DJF and the upper water temperature of the equatorial and austral subtropical regions. As a spin-up, the CESM was run for 200 model years under fixed orbital parameters of 300 kyr bp, with all other boundary conditions set to their values in 1950 ad. Then, the model was integrated for 3000 model years with the transient orbital insolation forcing of the past 300,000 years, in which orbital parameters were advanced by 100 years at the end of each model year64. This transient experiment was called CESM_orb34, and its outputs for the last 3000 model years were used in our analysis. This coarse-resolution model can reproduce the modern variability in the tropical ENSO, including its two different spatial types (eastern Pacific-type and central Pacific-type), in a realistic manner due to its improved convection scheme35 (Supplementary Fig. 7). It is reasonable to select this model for our analysis in view of the successful application of coarse-resolution models to paleo-ENSO changes under orbital forcing13,65,66.