Precession cycles of the El Niño/Southern oscillation-like system controlled by Pacific upper-ocean stratification

Modern observations have presented linkages between subsurface waters of the western Pacific warm pool and both El Niño/Southern Oscillation-related and extratropic-controlled upper-ocean stratification on interannual timescales. Moreover, studies have showed that such controls may operate on orbital cycles, although the details remain unclear. Here we present paired temperature and salinity reconstructions for the surface and thermocline waters in the central western Pacific warm pool over the past 360,000 years, as well as transit modeling results from an Earth system model. Our results show that variations in subsurface temperature and salinity in the western Pacific warm pool have consistently correlated with the shallow meridional overturning cell over the past four glacial-interglacial cycles, and they vary on eccentricity and precession cycles. The shallow meridional overturning cell regulates subsurface waters of the western Pacific warm pool by changing subtropical surface water density and thus equatorial upper-ocean stratification, acting as an El Niño/Southern Oscillation-like process in the precession band. Therefore, the western Pacific warm pool is critical in connecting the austral shallow meridional overturning cell to the Earth’s climate system on orbital timescales. Subsurface temperature and salinity in the Western Pacific Warm Pool were linked to the shallow overturning circulation and varied on orbital timescales over the past 360,000 years in El Niño/Southern Oscillation-like processes, suggest foraminifera proxy records and numerical modelling.

A s 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 system 1 . Modern observations show that the thermocline depth of the WPWP is closely associated with the El Niño/Southern Oscillation (ENSO) evolution 2 and large-scale changes in extratropical atmospheric circulation 1 . In turn, subsurface waters of the western equatorial Pacific (WEP) originate from the subduction of southern Pacific subtropical waters [3][4][5][6][7][8] , presenting extratropical control mainly through oceanic tunnels 9 . 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 Pacific 10 .
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 Sumatra 11 and Mindanao Island 12 in the far western part of the WPWP were controlled by local monsoon-induced oceanic upwelling processes responding to local summer precessional insolation changes 11 . Another newly published subT record from the WEP showed a dominant half-precession cycle under bihemispheric influences 13 . 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 obliquity 14 . Overall, most existing records of subsurface waters of the WPWP have focused on either the last glacial period 10 or regions far from the center of the WPWP 12,14 and have shown nonuniform information.
Here, we present new paired 360,000-year temperature and salinity records of subsurface water from the thermoclinedwelling 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 ruber 15 , 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.

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 Guinea 16 and 140-150 m in the eastern part of the WPWP within the upper thermocline 17 . 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. dutertrei 19 . The average core-top Mg/Ca-derived subT of 23.0 ± 1.1°C also conforms to modern observations 20 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 WEP 21 , implying that our temperature signals from core KX97322-4 are not considerably dampened by the effect of bioturbation.
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 lowlatitude 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 δ 18 O (δ 18 O sw ) deciphered from N. dutertrei shares similar features with the surface δ 18 O sw in terms of prominent glacial-interglacial variations, although the amplitudes of the subsurface δ 18 O sw (average of 2.31‰) are much larger than those of the surface δ 18 O sw (average of 1.42‰) (Fig. 2f, g). Furthermore, the average core-top SSS and subsurface seawater salinity (subS) calculated based on δ 18 O sw-iv values (see  [5][6][7] support that the more saline subsurface water in the modern WEP mainly originates from the subduction of eastern 6,7 and western 4,7 south Pacific subtropical waters transported by the SMOC (Fig. 1). Based on an analogy to modern climate and some simulations on orbitalscale 10,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 δ 18 O sw records with published surface temperature and δ 18 O sw records from sites MD97-2125 (ref. 27 ) and MD06-3018 (ref. 28 ) in the southwestern Pacific. δ 18 O sw represents the salinity signals as a passive tracer for the SMOC 5 . The two selected cores are located near the northern edge of the subtropical subduction region in the southwestern Pacific 29 and on the route of the SMOC 4 ( 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 winter 30 . 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 cycles 28 , probably due to their low resolution and/or other complex factors influencing the surface water in that region 27,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 δ 18 O sw record of KX97322-4 and the G. ruber δ 18 O sw 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 δ 18 O sw-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 Water 31 (Fig. 1). These waters passed mainly through the interior channel 32 (with perhaps a small part through the western boundary) and mixed with the waters along the route before reaching the equatorial subsurface 4,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.
Connection between the WPWP subsurface waters and ENSOlike processes in the precession band. Modern observations demonstrate a strong linkage between ΔT 25-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) ΔT 25-150m during a La Niña event and larger (8.15 ± 1.25°C) ΔT 25-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 ΔT G-N (corresponds to ΔT 25-150m ) can be applied to reflect the orbital-scale variations in the depth of the thermocline assuming a constant calcification depth of foraminifera 33 , a lower (higher) ΔT G-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 ΔT G-N data ( Fig. 4e) of core KX97322-4 with the CESM 34 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 band 35 . Our results show a high positive correlation between the GIG-detrended ΔT G-N of core KX97322-4 and the model-simulated Niño 3 DJF SST in the precession band (Fig. 4e, i). The minimum (maximum) ΔT G-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 relationship 36 (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, ΔT G-N , and ΔS N-G of our core and simulated subT data and Niño 3 SST (Supplementary Fig. 4) in the precession band.
In addition to the temperatures, modern seawater salinities of upper water in the study area also respond to ENSO activity. The ΔS 150-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 S 150m 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 flux 37 . However, in the precession band, the GIG-detrended δ 18 O sw difference between subsurface and surface δ 18 O sw values (Δδ 18 O sw ) deduced from N. dutertrei and G. ruber behaved oppositely from the modern ΔS 150-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 precipitationinduced SSS changes at KX97322-4 site in the precession band. At the precession minimum (P min ), enhanced evaporation during the austral winter led to increased salinity in the subtropical surface water 38 , 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 cycles 39,40 . However, the ITCZ mainly works on the surface water cycle (i.e., SSS) north of the equator 24 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 variations 40 . Coupled with the thick mixed layer in the central WPWP prohibiting vertical exchange with surface water 3 , 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 process 6,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 gyre 41 . 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 obliquity 14 . 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 location 13 . Our results agree with other numerical studies that have proposed that P min corresponds to enhanced Hadley circulation caused by the occurrence of aphelion near the boreal summer solstice 42,43 . Precession-paced solar insolation dominates low-latitude climate change via coupled ocean-atmosphere processes 44,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 band [46][47][48] .
The thermal energy of the tropical and subtropical Pacific is primarily obtained from local insolation 44 . 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 P min . 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 P min and P max (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 P min . In line with the modeling results for the early Holocene (P min ), warm subtropical water heats the region of Subantarctic Mode Water formation, and this warm anomaly is then transported to the tropical subsurface 49 , determined by local austral late winter insolation forcing 26 .
In addition, the largest insolation gradients between the equator and 23°S occur in the austral winter at P min 50 (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 P min expedites surface warming by preventing downward heat transfer through the thermal barrier effect of a relatively shallow halocline 3 , causing heat storage in the WEP 10 . 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, P max 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 circulation 51 and encourages more heat to be transferred to high latitudes, thus accelerating the melting of sea ice 52 and continental ice sheets through atmospheric bridges 53 . The associated shoaling of the thermocline in the EEP during a La Niña-like state additionally facilitates CO 2 release into the atmosphere 54 . 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 intensity 55 and equatorial Pacific Ocean productivity records 48 . 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 CO 2 to the atmosphere, providing positive feedback to global warming. Thermal state and hypothetical meridional and zonal ocean-atmospheric circulation in the Pacific under precession minimum (P min ) and maximum (P max ) conditions. 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 P min relative to the 300-kyr climatological mean from CESM simulation; c and d show the Pacific ocean-atmospheric circulation during P min and P max , 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.
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 CO 2 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 software 57 , combined with five accelerator mass spectrometry radiocarbon ( 14 C) 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 δ 18 O analyses for each sample. Samples for Mg/Ca analysis were cleaned following the treatment procedures without a reductive step 58 because the core was located far from the continent 15 . 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°C 16 to calculate subT as follows: Mg=Caðmmol mol À1 Þ ¼ 0:21 exp½0:097 subTð CÞ ð1Þ To better compare with the Mg/Ca ratio of G. ruber from our published work on the same core 15 , we also recalculated the SST using the equation for G. ruber developed by Hollstein et al. 16 as follows: Mg=Caðmmol mol À1 Þ ¼ 0:26 exp½0:097 SSTð 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 WPWP 16 .
The tests for δ 18 60 ). We selected this equation to avoid the potential influence of variable surface waters in other equations, and δ 18 O sw-iv is the icevolume corrected δ 18 O sw value using the mean ocean water δ 18 O value 61 .
Spectral and phase analysis. We used redfit spectral analysis with Past 3.20 software 62 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 year 64 . This transient experiment was called CESM_orb 34 , 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 scheme 35 (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 forcing 13,65,66 .

Data availability
All data in the paper are present in the paper and the supplementary materials. All the data named as "Sable oxygen isotope and Mg/Ca ratios of planktonic foraminifera from KX97322-4 (KX22-4)" can also be obtained from PANGAEA (https://www.pangaea.de/) or NOAA (https://www.ncei.noaa.gov) with the reference ID: UKRLJ3 or from the authors.