Ice sheet and precession controlled subarctic Pacific productivity and upwelling over the last 550,000 years

Domain: earth_science
Published: 2024-04-25
Editorial outcome: major_revision

Abstract
The polar oceans play a vital role in regulating atmospheric CO2 concentrations (pCO2) during the Pleistocene glacial cycles. However, despite being the largest modern reservoir of respired carbon, the impact of the subarctic Pacific remains poorly understood due to limited records. Here, we present high-resolution, 230Th-normalized export productivity records from the subarctic northwestern Pacific covering the last five glacial cycles. Our records display pronounced, glacial-interglacial cyclicity superimposed with precessional-driven variability, with warm interglacial climate and high boreal summer insolation providing favorable conditions to sustain upwelling of nutrient-rich subsurface waters and hence increased export productivity. Our transient model simulations consistently show that ice sheets and to a lesser degree, precession are the main drivers that control the strength and latitudinal position of the westerlies. Enhanced upwelling of nutrient/carbon-rich water caused by the intensification and poleward migration of the northern westerlies during warmer climate intervals would have led to the release of previously sequestered CO2 from the subarctic Pacific to the atmosphere. Our results also highlight the significant role of the subarctic Pacific in modulating pCO2 changes during the Pleistocene climate cycles, especially on precession timescale (~20 kyr).

Introduction
Atmospheric carbon dioxide concentrations (pCO2) recorded in ice cores have varied cyclically with oscillations typically ranging from ~180 ppmv during glacials to ~280 ppmv during interglacials1. Yet, the mechanisms responsible for the variations in pCO2 remain elusive. Previous studies have proposed the Southern Ocean as a key player in modulating pCO2 fluctuations over glacial-interglacial cycles2–4. However, recent literature highlights a notable gap in our understanding, particularly concerning the subarctic Pacific which is regarded as the polar twin of the Southern Ocean5. The North Pacific is the largest reservoir of respired carbon, comprising nearly half of the global ocean inventory today6. This substantial size of reservoir raises the possibility that it may have substantially influenced the carbon cycle in the past as well7. Recent research highlight the substantial impact of the subarctic Pacific in regulating the rapid deglacial rise in pCO2 levels8–11, and imply that it may have been a crucial role over longer timescales in the Pleistocene12–14. Furthermore, recent research challenges conventional wisdom by proposing that, during the latter stages of the last deglacial period, the North Pacific, rather than the Southern Ocean, may have contributed the lion’s share of CO2 degassing during the second half of the last deglacial15.

The exchange of CO2 between the atmosphere and the deep ocean carbon reservoirs, driven by a combination of physical and biological processes2,8,10,16, is the primary modulator of glacial-interglacial pCO2 changes. In this context, upwelling and vertical mixing in the North Pacific have been viewed as a crucial mechanism that influence Pleistocene pCO2 fluctuations by regulating the release of deeply-sequestered CO2 to the atmosphere8–10,14. To date, however, the mechanisms driving upwelling and, more broadly, the role of the subarctic Pacific in regulating Pleistocene pCO2 variations remain poorly understood, due to the scarcity of high-resolution records spanning several glacial cycles. Indeed, only a few studies have addressed the mechanisms controlling the rate of upwelling in the subarctic Pacific, with a predominant focus on the last deglacial transition10,17,18.

Export production proxies have been widely used to reconstruct changes in the rate of upwelling and vertical mixing8,14. Based on the productivity proxy, an upwelling record covering the past 850 kyr was reconstructed in the Bering Sea14, yet a 230Th normalization method is necessary to evaluate the efficacy of the productivity proxy19. Moreover, with multiple factors influencing upwelling in the subarctic Pacific10,14,18, it is crucial to differentiate these elements and assess their respective contributions. This distinction is pivotal to reveal the mechanisms that govern upwelling. In this study, we present high-resolution, 230Th-normalized export production and upwelling records from the subarctic Pacific, spanning the last five glacial-interglacial cycles. Based on the reconstructed productivity and upwelling records, in conjunction with transient climate model simulations, we explore the mechanisms related in particular to orbital forcing, greenhouse gasses and ice-sheet extent, underpinning the variations in productivity and upwelling on orbital timescales.

Results
Productivity changes over the past 550 kyr
Sediment core LV76-16-1 (48.85°N, 168.46°E, 2,374 m water depth) was recovered from the Emperor Seamount chain in the subarctic Pacific (Fig. 1). The 7.6 m long core is texturally homogeneous, dominated by fine silt and clay. The age model of core LV76-16-1 was established by assigning the abrupt increases in Ca/Ti to glacial terminations20 (Supplementary Fig. 1a, b). This correlation is based on the observation that each deglacial transition is characterized by transient increases in the sedimentary accumulation of biogenic carbonate all across the North Pacific12,13,21. Indeed, the downcore variability in Ca/Ti ratio of core LV76-16-1 is comparable to variations in the Ca/Al ratio and Ca counts at nearby ODP Site 88213 and core MD241622 (Fig. 1; Supplementary Fig. 2), both of which have a well-constrained chronology. Specifically, the chronology for ODP Site 882 was established by correlating high-resolution X-ray fluorescence (XRF) scanning Ba/Al ratios with the millennial-suborbital variability of δD from Antarctica ice core12,22. The age model for MD2416 were developed by aligning XRF scanning Ca counts with those from ODP Site 88222. This demonstrates that the chronology of core LV76-16-1 has a high enough resolution to allow for discussions of suborbital-scale variations. The calculated sedimentary accumulation rates of core LV76-16-1 show a coherent downcore variability, typically ranging between ~1.2 and 1.7 cm/kyr (Supplementary Fig. 1c).

Paleoproductivity proxies such as biogenic barium (BioBa), biogenic silica (opal), and carbonate have been widely used to reconstruct past changes in biological productivity23. Calcium (Ca) normalized to Ti reflects the sedimentary concentration of biogenic carbonate (CaCO3), and is consistent with the CaCO3 content measured in discrete samples (Supplementary Fig. 1b). Similarly, Ba/Ti in marine sediments is often used to reflect the concentration of BioBa. Unlike opal and carbonate, BioBa is considered to reflect integrated export production from the photic zone23, and can thus be considered as the sum of biogenic opal and carbonate, as silicious (diatoms) and calcareous (coccoliths) phytoplankton constitutes most of the primary producers in the subarctic Pacific24. This is corroborated by the good correspondence between the BioBa record and the normalized sum of the percentages of opal and CaCO3 contents downcore (Supplementary Fig. 3a, b).

The sedimentary opal content and Ca/Ti and Ba/Ti ratios all display, coherent, cyclic variations throughout the past 550 kyr (Fig. 2). Comparison of their downcore variations with the benthic δ18O stack20 shows that variations in these three proxies track glacial-interglacial cycles, with high values during interglacials and generally lower values during glacials (Fig. 2). This indicates that the temporal variations of these productivity proxies are strongly linked to glacial-interglacial cyclicity. In addition, comparison of the evolution of Earth’s precession and obliquity parameters in the past with our export production records suggest that their variations are also controlled by orbital forcing. The Ca/Ti ratio correspond well and positively with obliquity (Fig. 2a), while the opal content and the Ba/Ti ratio correlate negatively with precession (Fig. 2b, c). It is noteworthy that the effect of orbital forcing was particularly pronounced during interglacials when the global ice volume was small, while it was generally suppressed during full glacial intervals (Fig. 2a–c). For example, the precession peaks during MIS 1, 5, 7, 9 and 13 are clearly visible in the opal and Ba/Ti records, yet they are more subdued during MIS 2, 6, 8, 10 and 12 (Fig. 2b, c). Further analysis shows that the precessional signal filtered out from the Ba/Ti record compares very well with precession throughout the last 550 kyr (Supplementary Fig. 4a), confirming the role of precession. Wavelet analysis shows that the variations of the Ba/Ti and upwelling records contain three major periodicities of ~100 kyr, ~40 kyr and ~20 kyr (Supplementary Fig. 4b, c). The power of the ~20 kyr cycle (so does the power of the ~40 kyr cycle) is varying in time mainly due to the amplitude modulation of precession by eccentricity. The weak precession and strong obliquity signals during MIS 11 in all productivity proxies may be related to the relatively low amplitude of precession variability, caused by low eccentricity and the large variation in obliquity during MIS 11, a particularity which has also been observed in other climate variables25,26.

Productivity proxies may be subject to differential preservation during diagenesis21, due to the effects of degradation and dissolution8,27, which may complicate their interpretation. However, the consistent variations between the Ba/Ti and the sum of normalized CaCO3 and opal contents, along with clearly expressed orbital signals, indicate that dissolution of these biogenic components in the water column and the sediments did not substantially obscure the sedimentary records. Thus, the downcore BioBa record can be interpreted as a robust indicator of changes in integrated export production. This is further supported by the overall good agreement between the variations in sedimentary opal, CaCO3 and BioBa concentrations with the 230Th-normalized fluxes of each proxy over the last glacial cycle (Supplementary Fig. 5), which can also be applied for the old portion of the record. The reconstructed productivity history of core LV76-16-1 is largely consistent with previous studies covering a large region of the subarctic Pacific, including the Bering Sea28, Okhotsk Sea29, NW12,30 and NE Pacific31 (Supplementary Fig. 6a–d), all of which show higher export production during warm periods.

Export productivity in the subarctic Pacific can be enhanced by the alleviation of light limitation32, enhanced supply of iron33 and/or upwelling10. Previous results at nearby ODP Site 882 preclude the role of sea-ice-driven light limitation to account for generally reduced export production during glacial times, because the site, analogous to the location of core LV76-16-1, is located southeast of the summer sea-ice limit even during glacial maxima12. Enhanced productivity can also be caused by an increased iron supply associated with ice-rafted detritus (IRD) input via the melting of sea ice/icebergs34, and/or aeolian transport of dust35. However, these factors can also be excluded, as the IRD variations in core LV76-16-1 are inversely correlated with changes in opal and BioBa concentrations (Supplementary Fig. 3c). On the other hand, previous study suggested that Fe input associated with increased dust supply during the glacials did not exert a substantial influence on export production21. Although the subarctic Pacific is a high-nutrient low-chlorophyll (HNLC) region (Fig. 1), biological productivity in pelagic ecosystems is primarily sustained by the supply of nutrients, including Fe, to the euphotic zone via upwelling. Therefore, biological productivity in this region may have been primarily driven by changes in the upwelling of nutrient-rich subsurface water. This inference is consistent with previous studies in the subarctic Pacific covering various timescale9,10,12,14. Our results suggest that warm interglacial conditions and high boreal summer insolation, resulting from high obliquity and low precession, conspire to sustain high biological productivity in the subarctic Pacific. Under such warm conditions, the biological pump can be enhanced by the intensification of nutrient-rich subsurface water upwelling, facilitated by reduced sea-ice coverage and elevated ocean surface temperatures14,36.

Mechanism driving the nutrient upwelling in the subarctic Pacific Ocean
As biological productivity in the study region is primarily driven by the supply of nutrients via upwelling, we define an upwelling index based on bulk sedimentary δ15N and productivity proxies (see the Methods for detailed information), following the approach outlined by ref. 14. Similar to productivity proxy of Ba/Ti, primary variations in the upwelling index closely correlate with glacial-interglacial cycles, while secondary variations exhibit a strong correspondence with precession forcing (Fig. 2d; Supplementary Fig. 4). Our reconstructed upwelling index is consistent with a recent study from the Bering Sea14 (Fig. 3c, d), revealing that warm interglacial conditions and low precession (equivalent to high boreal summer insolation) conduce to stronger upwelling (Fig. 2d; Supplementary Fig. 4). Similar to the productivity proxies, upwelling during MIS 11 displayed a more pronounced response to obliquity when compared to other interglacial periods.

Both proxies and model simulations10,37 demonstrate that wind stress plays a key role in controlling Ekman suction and nutrient upwelling in the subarctic Pacific, which is ultimately related to the changes in the position and intensity of the northern westerlies. The northern westerlies are located at ~40°N today, and are the prevailing westerly winds that flow from the subtropical high-pressure belts to the subpolar low-pressure belts in the Northern Hemisphere37. These winds resulted from baroclinicity caused by the meridional thermal gradient between the warm subtropical atmosphere and the colder polar atmosphere38,39, and they are ultimately fueled by the non-uniform distribution of annual solar radiation inputs to the global Earth-atmosphere system37. In contrast to the subarctic Pacific area (north of ~40°N), productivity in the subtropical Pacific (south of ~40°N) was higher during glacials, at least over the last 200 kyr40–43 (Supplementary Fig. 6e–g). We thus propose that the observed seesaw pattern in the export production and nutrient upwelling records spanning the subarctic and subtropical Pacific regions were caused by a meridional shift and/or changes in the strength of the northern westerlies which seem sensitive to changes in glacial conditions. A compilation of planktic foraminiferal δ18O records from the North Pacific and climate model results indicate a southward migration (~3°) of the westerlies and increased wind stress during the Last Glacial Maximum18.

We propose that upwelling in the subarctic Pacific was mainly controlled by changes in wind stress and consequently by the variations in the northern westerlies which are controlled by glacial-interglacial ice sheets and orbital forcing. To test our hypothesis and to investigate how the westerlies responded to changes in ice sheets (ICE), greenhouse gases concentrations (GHG), and orbital parameters (Orb), we analyzed the results of three transient simulations performed using the LOVECLIM1.3 model. These simulations cover the period of 133–75 ka, which includes three precession cycles, one-and-a-half obliquity cycle, and large variations in GHG and ICE, and thus are suitable for our research purpose. The first two simulations, which were used and described in ref. 44, consider only the effect of orbital forcing (Orb experiment) and the combined effects of orbital forcing and GHG (OrbGHG experiment), respectively. In the third simulation, the changes in Northern Hemisphere (NH) ice sheets were included (OrbGHGICE experiment), and the experimental setup related to ice sheets was described in ref. 45.

The Orb experiment allows us to investigate the effects of orbital forcing alone. The model outputs show that the intensity of the westerlies and wind stress across the subarctic Pacific show high consistency, and are mainly controlled by precession, with obliquity playing a secondary role (Fig. 4c, d, red line). The model results indicate that low precession (or high boreal summer insolation) leads to stronger westerlies across the subarctic Pacific and by inference stronger upwelling, which is consistent with the abovementioned observation that upwelling was stronger when precession is low. Comparison of two extreme precession cases (i.e., the precession minimum (Pmin) at 127 ka and precession maximum (Pmax) at 97 ka), with similar obliquity, suggests that the westerlies over the northern mid-latitude Pacific were intensified at Pmin (Fig. 4e, g). Further analysis of the model results shows that, compared to Pmax, the low pressure centered over the Kamchatka Peninsula is deepened and the high pressure over the subtropical Pacific is strengthened at Pmin, leading to a larger latitudinal pressure gradient and thus, stronger westerlies.

Comparison between Orb and OrbGHG simulations shows that GHG has little effect on the westerlies over the subarctic Pacific, compared to the effect of precession. However, the OrbGHGICE simulation shows that the NH ice sheets have a large effect on the intensity and position of the westerlies, an effect which is much larger than that of precession (Fig. 4c, g). The effect of ice sheets (expressed by the difference between OrbGHGICE and OrbGHG) explains a much larger amount of variance of the westerlies in OrbGHGICE (R=0.90) than the effect of orbital forcing (R=0.32). The much greater impact of ice sheets can also be observed in Fig. 4g, and by comparing Fig. 4e, f (note that the color scale is different). Figure 4g shows that in response to large NH ice sheets, such as those at 91 ka, not only is the intensity of the westerlies substantially reduced but also its position is shifted southwards. This finding supports the observation that upwelling was weaker during glacials and stronger during interglacials. It is evident that the westerlies respond much less to precession changes, being about three times weaker, and their position is also less affected compared to the response to changes in ice sheets (Fig. 4g). Further analysis shows that large NH ice sheets greatly weaken the low-pressure system centered over the Bering Strait and the surrounding regions (pressure increases) and they slightly weaken the subtropical high-pressure system (pressure decreases), leading to a reduced pressure gradient between the subtropical and subpolar zones, and finally to the weakening and southward shift of the westerlies over the northern Pacific.

Therefore, our model results and proxies consistently show that the westerlies and upwelling over the subarctic Pacific are strongly influenced by ice sheets and precession. Model results also show that the changes in ice sheets exert the primary influence on the westerlies over the subarctic Pacific, with precession playing a secondary role. This aligns with the observation that variations in export productivity and the reconstructed upwelling index primarily follow glacial-interglacial cycles, with precession having a comparatively weaker impact throughout the sequence (Fig. 2a–d; Supplementary Fig. 4). Overall, warm conditions, resulting from a low ice volume and low precession, promote stronger upwelling of nutrient-rich subsurface water.

Implications for changes in glacial-interglacial atmospheric pCO2
The cyclic changes in pCO2 on glacial-interglacial timescales are mostly attributed to changes in the upwelling of CO2-rich water in the deep ocean2,8,10, while the Southern Ocean is considered as the primary control on this pCO2 variability46. Furthermore, Jaccard et al. (2013) proposed a two-mode mechanism in the Southern Ocean that may control the amplitude and timing of pCO2 changes during glacial-interglacial cycles4. Peak interglacial pCO2 changes in the Southern Ocean’s Antarctic zone were driven by CO2-rich water upwelling, while the transition to glacial conditions involved increased remineralized carbon sequestration in the ocean interior through a stronger Sub-Antarctic zone biological carbon pump aided by iron fertilization4.

While the Southern Ocean plays a critical role in regulating the Pleistocene pCO2 variability, similar processes could occur in the subarctic Pacific14,46, considering that approximately 30 ppm of pCO2 was released into the atmosphere as a result of enhanced overturning in the subarctic region during the last deglaciation47. More importantly, the variations in the upwelling index are consistent with those of pCO2 and of the modelled northern westerlies since ~550 ka (Fig. 3a, d, f). This suggests that wind-driven upwelling exerts a primary control on changes in pCO2. Moreover, the CaCO3 peaks coincided with the strongest upwelling during the warmest interglacials (Fig. 3e) throughout the entire subarctic Pacific12,13,21. Besides the effects of increased carbonate production, this reflects abrupt releases of deep-sequestered CO2 from the ocean back into the atmosphere, leading to a higher calcite saturation state in the bottom water4,12.

The most prominent feature in our productivity and upwelling records is the clear precessional signal in the subarctic Pacific (Fig. 2b–d; Supplementary Fig. 4), which is not evident in the Southern Ocean (Fig. 3b; Supplementary Fig. 7d, e), either in the Antarctic zones4,48 or the Sub-Antarctic49, or in the Equatorial Pacific50. Most upwelling records inferred from productivity proxies in these regions are dominated by 100 kyr cyclicity with a subsidiary 40 kyr cycle4,50–52. Our analysis clearly shows that precession (~20 kyr cycle) plays a more important role in the upwelling in the subarctic Pacific than that in the Southern Ocean. Therefore, the relatively high-frequency (20-kyr) cyclic changes in atmospheric pCO2 (Fig. 3a; Supplementary Fig. 4d) could be related to wind-driven upwelling in the subarctic Pacific. A recent study provides partial support for our inference that CO2 degassing during the second half of the last deglacial occurred in the North Pacific and elsewhere, instead of in the Southern Ocean15.

In summary, our results shed light on the significant role of the subarctic Pacific in regulating atmospheric pCO2 concentrations, primarily through the dynamic interplay of ice sheets and insolation on the northern westerlies and wind-driven upwelling. In particular, we have introduced a mechanism to explain the variations of upwelling in subarctic Pacific, through which differentiate the forcing of ice sheets and insolation. We propose that the subarctic Pacific plays a crucial role in regulating atmospheric CO2 concentration variations on precession timescale, which would fill a critical gap in our knowledge of the subarctic Pacific’s significance in regulating global atmospheric CO2 changes across different timescales. Under warm climatic conditions, the strengthening of the northern westerlies and the enhancement of upwelling in the subarctic Pacific could lead to increased atmospheric pCO2 levels. With the intensification of global warming, the anticipated poleward shift of the northern westerly winds53 is poised to further amplify the upwelling rate in the subarctic Pacific. This, in turn, would release more CO2 into the atmosphere and further accelerate global warming through positive feedbacks.

Methods
Opal and CaCO3 content analysis
Concentration of biogenic opal of 766 samples (1 cm intervals, ~1 kyr resolution) from core LV76-16-1 was determined by alkaline extraction of silica54, and was measured by molybdate-blue spectrophotometry. The long-term accuracy of this method is ±0.5 wt%, as deduced from replicate and in-house standard measurements. A total of 383 samples were analyzed at 2 cm intervals (~2 kyr resolution) for measuring CaCO3 content by dissolution in HCl.

Elemental abundance determination
A total of 383 samples were analyzed for major and trace elements. All samples were first digested with HNO3-HF (1:1) in closed Teflon beakers, and analyzed with ICP-MS (Thermo Scientific XSERIES 2) for trace elements. Quality control were implemented throughout the entire experimental process, including the inclusion of a blank experiment, the use of GSD-9 standard material, and replicate measurements. The relative standard deviations of the major and trace element analyses are all less than 5%.

Concentrations of biogenic Ba (BioBa) can be expected based on the assumption that the composition of Ti in the terrigenous material remained constant in space and time. The elemental ratio from the upper continental crust (UCC)55 are used for normalization.

BioBa = Ba_sample - (Ba/Ti)UCC × Ti_sample (1)

XRF scanning elemental abundance
Elemental abundances (in counts per second; cps) in core LV76-16-1 were measured by X-ray fluorescence (XRF) core scanning at ~5 mm intervals (~0.5 kyr resolution) using the Itrax XRF Core Scanner, with 20 s count times, 30 kV X-ray voltage, and an X-ray current of 40–55 mA.

Sedimentary δ15N analysis
A total of 379 bulk samples were measured for nitrogen isotopic analysis on an Isoprime 100 isotopic ratio mass spectrometer (IRMS). Values for δ15N are reported in permil notation relative to atmospheric nitrogen gas. A standard sample (protein, δ15N = 6.0‰) was inserted for every 10 samples, and the measured nitrogen isotopes of standards yielded a precision of better than ±0.2‰. The reproducibility of replicate nitrogen isotopic analyses for duplicate samples was generally better than ±0.3‰.

Ice-rafted debris (IRD) analysis
Ice-rafted debris counts were carried out for 382 samples at ~2 cm interval using a standard reflected light binocular microscope. Each sample weighing approximately 5 g was wet-sieved through a 150 μm sieve. The IRD is defined as the grains (>150 μm) per gram of the dry bulk sediment, as this fraction was regarded as the lower size threshold for IRD categorization56.

Calculation of Upwelling index
The upwelling index is calculated following the approach in ref. 14 by using productivity proxy and sedimentary δ15N. Briefly, the normalized sum of CaCO3 and opal (proxy of productivity) were linearly interpolated at an interval of 2 kyr (mean resolution is 1.5 kyr), which was subtracted by the normalized Δδ15N to create a semi-quantitative proxy of the ‘upwelling index’, with higher values indicating increased nutrient upwelling. The sedimentary δ15N from ODP Site 1012 in the eastern tropical North Pacific57 is thought to be a site of nearly complete nitrate utilization throughout the Pleistocene14,22.

Upwelling index = (CaCO3 + opal)_normalized - (Δδ15N_LV76-1012)_normalized (2)

Analysis of U-Th isotopes
Uranium and Thorium isotopes of 27 samples were determined by laser ablation Multi-Collector Inductively Coupled Plasma Mass Spectrometry (MC-ICP-MS) on a Neptune plus. Sample preparation and analysis followed the methods described in ref. 58. Briefly, about 0.2–0.3 g of dried sample was homogenized, heated to remove volatile components, packed in a molybdenum capsule and isolated in a graphite tube to avoid oxidation, then heated at 1300 °C for 10–15 min before being quenched in water. By doing this, the samples were transformed into homogeneous glass58. Analysis of 230Th, 232Th, and 238U in glass samples were done by the standard sample bracketing method, on a Neptune plus MC-ICP-MS, coupled with a New Wave 193 Research Laser Ablation system in the State Key Laboratory for Mineral Deposits Research in Nanjing University. The LA-MC-ICPMS was optimized using NIST 612 glass to achieve the best 238U and 232Th intensities, resulting in a U/Th ratio of around 1.1–1.3. The UO+/U+ in our experiments was less than 1%. The 232Th tail at 230Th mass was ~120 ppb with RPQ on.

Estimates of fluxes using 230Th normalization
The application of 230Th normalization allows high-resolution sediment mass flux reconstructions over time which are insensitive to lateral sediment redistribution19. The 230Th-normalized flux of biogenic component (opal, CaCO3, and BioBa; Supplementary Fig. 5) were calculated following the method described in ref. 19:

Flux_i = (C_i × β230 × d) / 230Thex0 (3)

Where Flux_i represents the flux of a given constituent i (e.g., opal, CaCO3, BioBa), with a concentration C_i in the sediment deposited at a specific water depth (d). The production rate of 230Th (β230 = 2.56 × 10^5 dpm/cm^3/kyr) was used59. 230Th_ex0 (dpm/g) is corrected for radioactive decay since its time of deposition, the fraction supported by uranium within lithogenic material (238U/232Th = 0.7 ± 0.1 in the Pacific Ocean) and the fraction of the in situ 230Th produced by decay of authigenic 238U19.

Model and simulations
The model used in this study is LOVECLIM1.3, a three-dimension Earth system Model of Intermediate Complexity, with its atmosphere (ECBilt), ocean and sea ice (CLIO) and terrestrial biosphere (VECODE) components being interactively coupled. The model configuration is identical to that employed in ref. 44 and detailed description can be found therein.

In our study, two types of simulations are used. The first type is transient simulations covering the period of 133–75 ka without acceleration. It includes three simulations. The first two simulations, Orb and OrbGHG, were performed in ref. 44 and detailed description of experimental setup can be found therein. Here we only give some brief introduction. In order to isolate and understand the impact of orbital forcing, in the Orb simulation, only the change of orbital forcing60 was taken into account, with the GHG and ice sheets being fixed to their pre-Industrial condition. In the OrbGHG simulation, the change of GHG61–63 is additionally taken into account. In the third simulation, OrbGHGICE, the change of NH ice sheets64 was taken into account, but the Southern Hemisphere ice sheets remain fixed to the pre-Industrial condition. The initial conditions were provided by a 2000-year equilibrium experiment with the NH ice sheets, GHG concentrations and astronomical parameters at the starting date of the simulated period. In the presence of land ice, albedo, topography, vegetation and surface soil types corresponding to ice-covered condition were prescribed at corresponding model grids in LOVECLIM1.3. Detailed description of the ice sheet setup can be found in ref. 45.

A transient simulation covering the last 800 ka with 10x acceleration was conducted, with the results of the last 550 ka used. In the transient simulation, variations in orbital forcing and GHG were taken into account, and the ice sheets were fixed to their pre-Industrial condition. To consider the impact of the NH ice sheets on the westerlies throughout the last 550 ka, we apply the relationship between the NH ice volume and the effect of NH ice sheets on the westerlies obtained from the 133–75 ka simulations, which are shown highly linearly correlated (correlation coefficient R = 0.9), on the ice volume of the last 550 ka to estimate the westerly strength over the sub-arctic Pacific, as shown in Fig. 3f.

Peer review rounds
[{"round": 1, "reviewer_comments": {"reviewer_1": "Yao et al. present new high resolution reconstructions of 230Thorium-normalised export productivity in the subarctic northwest Pacific Ocean over the last 550 kyr derived from opal, calcium carbonate, major and trace elements, sedimentary d15N, and ice-rafted debris, using sediment core drilled from the Emperor Seamount chain. While the presentation of the data in the manuscript is overall good, these data were not sufficiently integrated with previously published data in the region, and I found that the major conclusions were not adequately supported by data analysis. I also would like to see better delineation of what is new by the authors, along with better integration of existing relevant data.\nAt this point, I recommend major revision. In the present form, I do not think that this publication is suitable for Nature Communications. Detailed comments to the authors are provided below, I hope these are helpful for revision.\nTitle:\nI think that ice sheet and precession should not be placed together in the same breath here. Major findings of this work show that ice sheet is the main driver and precession only contribute to a \"lesser\" extent. Consider revising the title to highlight major and novel findings of this work.\nAbstract:\nI was confused. The way the Abstract is written seems like this manuscript is about the Southern Ocean. I understand that the authors are trying to compare the study area to the Southern Ocean in terms of carbon reservoir. Consider rewriting the first two sentences to make a better transition to the subarctic Pacific.\nL33: Spell-out “NW”.\nL37: What kind of “biological production”? In L33, you used export productivity.\nL38: Delete “somewhat”.\nIntroduction:\nSame comment in the Abstract about the way the Southern Ocean is highlighted. Why start with the Southern Ocean? If you do it this way, make sure that you transition the focus to your studied region and discuss the connection of the Southern Ocean to the North Pacific.\nThe order of the Introduction does not logically flow from large scale picture to focused. We know that most of the studies focus on the role of the Southern Ocean in atmospheric carbon dioxide concentrations during glacial-interglacial cycles. The authors can improve the first paragraph by highlighting the gap in recent studies so that it focuses on the North Pacific. The third paragraph is very important as this is the novel aspect of this work, yet the way the authors introduced what they have done is muted.\nL63: Insert “an” before “important”.\nL71: Replace “their” with “this”.\nProductivity changes over the past 550 kyr:\nL81: Call out Extended Data Fig. 1a first before 1b, c, etc. Please rearrange the figure order.\nL85: Replace “with” with “to”.\nL86: Please see comment in L81.\nL87: Replace “between ~1.2-1.7” with “between ~1.2 and 1.7”.\nL111-113: The precession is not evident in the plot, especially prior to ~250 kyr. Only MIS 5 shows precession. I can see a distinct 100 kyr and 41 kyr cycles in all the records, similar to the LR04 stack. Since this is highlighted as a significant finding of this research, I suggest running a spectral analysis (simple periodogram, multitaper, and the likes). I would like to see the precession peak in the plot/plots. It would also help doing the wavelet analysis or filter out the precession signal in the Ba/Ti record, to see the temporal evolution of the cyclicities. There could be a shift in orbital frequencies through the 550 kyr record.\nL117-118: Add an explanation how astronomical forcing influence the surface waters (change in water mass characteristics) that led to high biological productivity. This section lacks a discussion on oceanography.\nL136: Replace “was” with “is”.\nL138: Add reference/s for “aeolian transport of dust”.\nL140: You can delete “on one hand”.\nL141: Which previous studies?\nMechanism driving the nutrient upwelling in the subarctic Pacific Ocean\nL152: What do you mean by “first-order variations”?\nL153: What do you mean by “second-order variations”?\nL153: I am not convinced about the precession signal in the upwelling index. Can you please do a spectral and/or wavelet analysis?\nL156: “…responded mainly to obliquity…” Not only during MIS 11 but all through the record. Precession can be observed from ~250 kyr but the amplitude is not comparable to the interglacial periods. Please delete dash in “MIS-11”.\nL161: Delete the second “s” in “Hemispheres”.\nL162: Replace “These winds are resulting…” to “These winds resulted from…”\nL166: Add “…glacials, at least over the last 200 kyr.” because the record is only until around this time.\nL198: Replace “Comparisons of…” to “Comparisons between…”\nL208-210: Sentence not clear. Please rewrite.\nL221: Please mention time interval when precession “exerts a weaker influence”. Is it through the 550 kyr record?\nImplications for changes in glacial-interglacial atmospheric pCO2\nL224: Please discuss the link between the Southern Ocean and the North Pacific. See comments previous comments on this.\nL241: Not convinced about the “clear” precessional signal unless a spectral plot is shown. Look the same in the Southern Ocean, which is dominated by obliquity and/or eccentricity.\nL249-256: This manuscript has some significant and novel findings and yet it looks to me that they are not adequately described and claimed. This last paragraph can be better written to highlight these findings.\nMethods:\nPlease mention respective resolution in kyr of sampling interval in each analysis type.\nLine 267: Replace “by” with “with”.\nLine 268-269: Please rephrase “The entire experimental process was under quality\n269 control using a blank experiment…”\nL276: Delete “the”.\nL282: Spell-out “N2”.\nL292: Add “…the approach in Ref13…”\nL294: Delete “/sample”.\nLine 328-329: Rephrase.\nLine 331: Delete “of”.\nLine 332-333: Rephrase.\nLine 334: Add short explanation why “In the Orb simulation, only the change of orbital forcing was taken into account…”.\nFigures:\nSome figure/s need/s to be rearrange based on the order it was called out in the manuscript. Please see comments above.\nMake sure all figures are called out in the manuscript.", "reviewer_2": "In this paper Yao et al. try to show that the upwelling in the subarctic North Pacific is especially pronounced during warm interglacial periods with high boreal summer insolation. I have to say that I am not convinced by this finding, especially the association between more intense upwelling and high boreal summer insolation.\nThis sort of analysis has been down before by Galbraith et al. (2008). According to my reading of the Galbraith et al. paper, periods of more intense upwelling in the western North Pacific take place during times with low boreal summer insolation. This is exactly the opposite conclusion that Yao et al. come to in this paper.\nFigure 6 in Galbraith et al., in particular, identifies three intervals with low levels of nitrate assimilation during marine isotope stage 5. These should be times of more intense upwelling. The intervals of low nitrate assimilation occur around 72,000, 96,000, and 120,000 years ago, times when the boreal summer insolation was reduced. Yao et al. even develop an upwelling index that is based on the same approach as Galbraith et al., yet they come to the opposite conclusion.\nThe problem comes down to the age model developed for core LV76-16-1 in relation to the age model developed for core MD2416 in Galbraith et al. (2008). Galbraith et al. at least had a benthic d18O record that they were able to align with LR04. Yao et al. are stuck using %CaCO3 and Ca/Ti ratios to align with LR04 instead. I don't see that Yao et al. refer to a benthic d18O record for LV76-16-1 anywhere in their paper. The approach used by Yao et al. is simply not as valid, especially when trying to align a core chronology with the ~23 kyr precessional cycle.\nTo be sure, precessional cycles were not the main focus of the Galbraith et al. paper. But the observations that they published speaks for themselves. Yao et al. have simply drawn the wrong conclusions about what drives the productivity and upwelling in the subarctic Pacific in the precessional band."}, "author_response": "Response to Reviewers' comments:\nReviewer #1\nMajor comments:\nComment 1:\nYao et al. present new high-resolution reconstructions of 230Thorium-normalised export\nproductivity in the subarctic northwest Pacific Ocean over the last 550 kyr derived from opal,\ncalcium carbonate, major and trace elements, sedimentary d15N, and ice-rafted debris, using\nsediment core drilled from the Emperor Seamount chain. While the presentation of the data in\nthe manuscript is overall good, these data were not sufficiently integrated with previously\npublished data in the region, and I found that the major conclusions were not adequately\nsupported by data analysis. I also would like to see better delineation of what is new by the\nauthors, along with better integration of existing relevant data.At this point, I recommend\nmajor revision. In the present form, I do not think that this publication is suitable for Nature\nCommunications. Detailed comments to the authors are provided below, I hope these are\nhelpful for revision.\nReply:\nThanks for the suggestions. We have added more information on the novelty of our study\nwith a better integration of existing relevant data. Please see Lines 45-59,70-76 and Lines 268-\n280 in the revised manuscript. In addition, integration of previous studies has also been made\nby comparison of our records with other records (see lines 136-139, 166-169). A summary of\nthe novelty of our study is given here below.\nUpwelling and vertical mixing in the North Pacific play a crucial role on the Pleistocene\npCO changes by regulating the release of deeply-sequestered CO into the atmosphere.\n2 2\nHowever, the mechanism driving the upwelling and more generally the role of the subarctic\nPacific in regulating Pleistocene pCO variations remain poorly understood. A few studies have\n2\nspecifically addressed the factors controlling the upwelling in the subarctic Pacific, but their\nrecords are too short so allow focusing only on the last deglacial transition (e.g., Gray et al.,\n2018; Du et al., 2018). High-resolution upwelling records spanning multiple glacial cycles in\nthe subarctic Pacific are very limited. There exists so far only one long upwelling record that\nwas obtained from productivity proxies and sedimentary δ15N over the past ~800 ka at the\nIODP Site U1343 in the Bering Sea (Worne et al., 2019). The variations of the reconstructed\nupwelling index at this site and at our site are highly consistent on orbital timescale (see Fig.\n3c, d in revised manuscript), showing the robustness of the upwelling index reconstructed in\nboth studies. In addition to the long, high-resolution upwelling index, our manuscript also\nincludes the following new aspects:\n(1) We have introduced a new mechanism to explain the variations in upwelling and\nproductivity in subarctic Pacific. Our reconstructed upwelling records and climate model\nresults consistently show the crucial role of precession and the Northern Hemisphere ice sheets\nin regulating the upwelling and export production in subarctic Pacific via their influences on\nthe westerlies. In particular, our model simulations allow to separate and quantify the effect of\ndifferent forcing factors, which can hardly be done by using proxy records alone.\n(2) Our productivity proxy is validated by 230Th normalization, which is critical for\ninterpreting export production in marine sediments(Francois et al., 2004) and thus upwelling.\nSuch validated records spanning multiple glacial cycles are rare in the subarctic Pacific region.\n(3) We propose that the subarctic Pacific plays a crucial role in regulating atmospheric\nCO concentration variations on precession timescale, which is completely original and has not\n2\nbeen addressed in other study. This finding would fill a critical gap in our knowledge of the\nsubarctic Pacific's significance in regulating global atmospheric CO changes across different\n2\ntimescales.\nComment 2:\nTitle:\nI think that ice sheet and precession should not be placed together in the same breath here.\nMajor findings of this work show that ice sheet is the main driver and precession only\ncontribute to a \"lesser\" extent. Consider revising the title to highlight major and novel findings\nof this work.\nReply:\nThank you for the suggestion. Indeed, both our proxy records and model simulations show\nthat ice sheet is the primary driver and precession is secondary. However, both drivers are\nneeded in explaining the variability of the productivity and upwelling in the subarctic Pacific,\nwith ice sheet contributing mainly to the variability of the 100-kyr glacial cycles and precession\nto the variability of ~20-kyr cycles. Moreover, the clear precession signal in our records, which\nis distinctly different from the Southern Ocean records, is one of the original findings in our\nstudy. Therefore, we would keep precession in the title, but slightly change the title to “Ice\nsheet and precession control on biological export production and upwelling in the subarctic\nPacific over the last 550,000 years”.\nComment 3:\nAbstract:\nI was confused. The way the Abstract is written seems like this manuscript is about the Southern\nOcean. I understand that the authors are trying to compare the study area to the Southern Ocean\nin terms of carbon reservoir. Consider rewriting the first two sentences to make a better\ntransition to the subarctic Pacific.\nReply: Thanks for the suggestion. We agree with the reviewer and have rewritten the first two\nsentences in the Abstract (Lines 29-31) as the following:\n“The polar oceans play a vital role in regulating atmospheric CO concentrations during\n2\nthe Pleistocene glacial cycles. However, despite being the largest modern reservoir of respired\ncarbon, the impact of the subarctic Pacific remains poorly understood due to limited records.”\nComment 4:\nIntroduction:\nSame comment in the Abstract about the way the Southern Ocean is highlighted. Why start\nwith the Southern Ocean? If you do it this way, make sure that you transition the focus to your\nstudied region and discuss the connection of the Southern Ocean to the North Pacific.\nThe order of the Introduction does not logically flow from large scale picture to focused. We\nknow that most of the studies focus on the role of the Southern Ocean in atmospheric carbon\ndioxide concentrations during glacial-interglacial cycles. The authors can improve the first\nparagraph by highlighting the gap in recent studies so that it focuses on the North Pacific. The\nthird paragraph is very important as this is the novel aspect of this work, yet the way the authors\nintroduced what they have done is muted.\nReply: As suggested by reviewer 1, we have rewritten the first and third paragraphs in the\nIntroduction section. Please see Lines 45-59 and 70-76.\nComment 5:\nL111-113: The precession is not evident in the plot, especially prior to ~250 kyr. Only MIS 5\nshows precession. I can see a distinct 100 kyr and 41 kyr cycles in all the records, similar to\nthe LR04 stack. Since this is highlighted as a significant finding of this research, I suggest\nrunning a spectral analysis (simple periodogram, multitaper, and the likes). I would like to see\nthe precession peak in the plot/plots. It would also help doing the wavelet analysis or filter out\nthe precession signal in the Ba/Ti record, to see the temporal evolution of the cyclicities. There\ncould be a shift in orbital frequencies through the 550 kyr record.\nL153: I am not convinced about the precession signal in the upwelling index. Can you please\ndo a spectral and/or wavelet analysis?\nL241: Not convinced about the “clear” precessional signal unless a spectral plot is shown. Look\nthe same in the Southern Ocean, which is dominated by obliquity and/or eccentricity.\nReply:\nWe thank the reviewer for the suggestion, and agree that the suggested analysis would\nmake our point more convincing. In our revised manuscript, we have added two figures as\nsupplementary materials (Extended Data Fig. 4 and Extended Data Fig. 7, also given below)\nto show the relevant analysis we have added and to strengthen our conclusion.\nExtended Data Fig. 4a clearly shows that the precessional signal filtered out from our\nBa/Ti record compares very well with precession throughout the last 550 kyr. The wavelet\nanalysis of our Ba/Ti and upwelling records (Extended Data Fig. 4b, c) shows three major\nperiodicities of ~100 ka, ~40 ka and ~20 ka. It is worth noting that the power of the ~20 ka\ncycle (so does the power of the ~40 ka cycle) is varying in time mainly due to the amplitude\nmodulation of precession by eccentricity. However, what is new and exciting in our study is\nthat the power of the ~20 ka periodicity in our proxies is equivalent to or even stronger than\nthe power of the 40-ka one over most time of the last 550 ka. This is distinctly different from\nthe productivity records (used to indicate upwelling) in the Southern Ocean (Extended Data\nFig. 7d, e) in which the precession signal is hardly visible, but it is clearly visible even in the\ncurves of our raw data (Extended Data Fig. 7a-c). Our analysis clearly show that precession\n(~20 ka cycle) plays a more important role in the upwelling in the subarctic Pacific than that\nin the Southern Ocean, which further suggests that the upwelling in the subarctic Pacific might\nhave contributed to the ~20 ka cycle in the variations of the atmospheric CO concentration.\n257-261).\nExtended Data Fig. 4. Variations in Ba/Ti, wavelet analysis of Ba/Ti and upwelling index\nfrom core LV76-16-1. (a) Comparison between Ba/Ti (dark line), 20-kyr filtering (blue line),\nPrecession (red line) and Eccentricity (green line) (Berger and Loutre, 1991). (b) Wavelet\nanalysis of Ba/Ti from core LV76-16-1. (c) Wavelet analysis of upwelling index from core\nLV76-16-1.\nExtended Data Fig. 7. Comparison between productivity, upwelling proxies of Core\nLV76-16-1 and records from the Southern Ocean. (a) Variations in upwelling index of\nLV76-16-1. (b) Variations in Ba/Ti of LV76-16-1. (c) Variations in opal of LV76-16-1. (d)\nVariations in opal of Core RC13-259 from the Southern Ocean (Charles et al., 1991). (e)\nVariations in Ba/Fe of ODP Site 1094 from the Southern Ocean (Jaccard et al., 2013).\nComment 6:\nL117-118: Add an explanation how astronomical forcing influence the surface waters (change\nin water mass characteristics) that led to high biological productivity. This section lacks a\ndiscussion on oceanography.\nReply:\nAn explanation regarding this issue has been incorporated into the revised manuscript\n(Lines 155-159).\n“Our results suggest that warm interglacial conditions and high boreal summer insolation,\nresulting from high obliquity and low precession, conspire to sustain high biological\nproductivity in the subarctic Pacific. Under such warm conditions, the biological pump can be\nenhanced by the intensification of nutrient-rich subsurface water upwelling, facilitated by\nreduced sea-ice coverage and elevated ocean surface temperatures (Worne et al., 2019; Kender\net al., 2018).”\nComment 7:\nL156: “…responded mainly to obliquity…” Not only during MIS 11 but all through the record.\nPrecession can be observed from ~250 kyr but the amplitude is not comparable to the\ninterglacial periods.\nReply:\nPlease see our reply to Comment 5. The precession signal is clear throughout the 550 kyr\nalthough its power is varying in time due to the amplitude modulation by eccentricity. The\npower of the ~20 kyr cycle is equivalent to or even stronger than the ~40 kyr cycle over most\ntime of the 550 kyr (Extended Data Fig. 4).\nComment 8:\nL224: Please discuss the link between the Southern Ocean and the North Pacific. See comments\nprevious comments on this.\nReply:\nWe have included relevant information regarding the link between the Southern Ocean\nand the North Pacific (Lines 247-250).\n“While the Southern Ocean plays a pivotal role in regulating the Pleistocene pCO\n2\nvariability, similar processes could occur in the subarctic Pacific (Worne et al., 2019; Sigman\net al., 2010), considering that approximately 30 ppm of pCO was released into the atmosphere\n2\nas a result of enhanced overturning in the subarctic region during the last deglaciation (Rae et\nal., 2014).”\nComment 9:\nL249-256: This manuscript has some significant and novel findings and yet it looks to me that\nthey are not adequately described and claimed. This last paragraph can be better written to\nhighlight these findings.\nReply:\nAs recommended by the reviewer, the final paragraph has been revised to elucidate the\nsignificance of these findings (Line 268-280).\n“In summary, our results shed light on the significant role of the subarctic Pacific in\nregulating atmospheric pCO concentrations, primarily through the dynamic interplay of ice\n2\nsheets and insolation on the northern westerlies and wind-driven upwelling. In particular, we\nhave introduced a new mechanism to explain the variations of upwelling in subarctic Pacific,\nthrough which differentiate the forcing of ice sheets and insolation. We propose that the\nsubarctic Pacific plays a crucial role in regulating atmospheric CO concentration variations\n2\non precession timescale, which would fill a critical gap in our knowledge of the subarctic\nPacific's significance in regulating global atmospheric CO changes across different timescales.\nenhancement of upwelling in the subarctic Pacific could lead to increased atmospheric pCO\n2\nlevels. With the intensification of global warming, the anticipated poleward shift of the\nnorthern westerly winds (Yang et al., 2020) is poised to further amplify the upwelling rate in\nthe subarctic Pacific. This, in turn, would release more CO into the atmosphere and further\n2\naccelerate global warming through positive feedbacks.”\nOther Comments:\nL33: Spell-out “NW”.\nReply: Done.\nL37: What kind of “biological production”? In L33, you used export productivity.\nReply: We have changed “biological production” to “export productivity” to maintain\nconsistency.\nL38: Delete “somewhat”.\nReply: Done.\nL63: Insert “an” before “important”.\nReply: Done.\nL71: Replace “their” with “this”.\nReply: Done.\nL81: Call out Extended Data Fig. 1a first before 1b, c, etc. Please rearrange the figure order.\nReply: Done.\nL85: Replace “with” with “to”.\nReply: Done.\nL86: Please see comment in L81.\nReply: Done.\nL87: Replace “between ~1.2-1.7” with “between ~1.2 and 1.7”.\nReply: Done.\nL136: Replace “was” with “is”.\nReply: Done.\nL138: Add reference/s for “aeolian transport of dust”.\nReply: Done.\nL140: You can delete “on one hand”.\nReply: Done.\nL141: Which previous studies?\nReply: We have revised this sentence to enhance clarity of meaning. “On the other hand,\nprevious study suggested that Fe input associated with increased dust supply during the glacials\ndid not exert a substantial influence on export production (Kohfeld and Chase, 2011).” Please\nsee Lines 148-150.\nL152: What do you mean by “first-order variations”?\nReply: We have replaced this phrase by “primary variations” to enhance clarity of\nmeaning.\nL153: What do you mean by “second-order variations”?\nReply: We have replaced this phrase by “secondary variations” to enhance clarity of\nmeaning.\nL161: Delete the second “s” in “Hemispheres”.\nReply: Done.\nL162: Replace “These winds are resulting…” to “These winds resulted from…”\nReply: Done.\nL166: Add “…glacials, at least over the last 200 kyr.” because the record is only until around\nthis time.\nReply: Done.\nL198: Replace “Comparisons of…” to “Comparisons between…”\nReply: Done.\nL208-210: Sentence not clear. Please rewrite.\nReply: This sentence has been changed to “It is evident that the westerlies respond much\nless to precession changes, being about three times weaker, and their position is also less\naffected compared to the response to changes in ice sheets.” Please see Lines 222-225.\nL221: Please mention time interval when precession “exerts a weaker influence”. Is it through\nthe 550 kyr record?\nReply: Please see our reply to Comment 5. The precession signal is clear throughout the\n550 kyr although its power is varying in time due to the amplitude modulation by eccentricity.\nThus, we rephrased this sentence. “This aligns with the observation that variations in export\nproductivity and the reconstructed upwelling index primarily follow glacial-interglacial cycles,\nwith precession having a comparatively weaker impact throughout the sequence.” Please see\nLines 233-236.\nPlease mention respective resolution in kyr of sampling interval in each analysis type.\nReply: Done.\nLine 267: Replace “by” with “with”.\nReply: Done.\nLine 268-269: Please rephrase “The entire experimental process was under quality control\nusing a blank experiment…”\nReply: This sentence has been changed to “Quality control were implemented throughout\nthe entire experimental process, including the inclusion of a blank experiment, the use of GSD-\n9 standard material, and replicate measurements.”\nL276: Delete “the”.\nReply: Done.\nL282: Spell-out “N2”.\nReply: Done.\nL292: Add “…the approach in Ref13…”\nReply: Done.\nL294: Delete “/sample”.\nReply: Done.\nLine 328-329: Rephrase.\nReply: This sentence has been changed to “The model configuration is identical to that\nemployed in ref. (Yin et al., 2021) and detailed description can be found therein.”\nLine 331: Delete “of”.\nReply: Done.\nLine 332-333: Rephrase.\nReply: This sentence has been changed to “The first two simulations, Orb and OrbGHG,\nwere performed in ref. (Yin et al., 2021) and detailed description of experimental setup can be\nfound therein.”\nLine 334: Add short explanation why “In the Orb simulation, only the change of orbital forcing\nwas taken into account…”.\nReply: This sentence has been changed to “In order to isolate and understand the impact\nof orbital forcing, in the Orb simulation, only the change of orbital forcing (Berger and Loutre,\n1991) was taken into account, with the GHG and ice sheets being fixed to their pre-Industrial\ncondition.”\nSome figure/s need/s to be rearrange based on the order it was called out in the manuscript.\nPlease see comments above.\nMake sure all figures are called out in the manuscript.\nReply: Done.\nReviewer #2\nComment 1:\nIn this paper Yao et al. try to show that the upwelling in the subarctic North Pacific is especially\npronounced during warm interglacial periods with high boreal summer insolation. I have to say\nthat I am not convinced by this finding, especially the association between more intense\nupwelling and high boreal summer insolation.\nThis sort of analysis has been down before by Galbraith et al. (2008). According to my reading\nof the Galbraith et al. paper, periods of more intense upwelling in the western North Pacific\ntake place during times with low boreal summer insolation. This is exactly the opposite\nconclusion that Yao et al. come to in this paper.\nFigure 6 in Galbraith et al., in particular, identifies three intervals with low levels of nitrate\nassimilation during marine isotope stage 5. These should be times of more intense upwelling.\nThe intervals of low nitrate assimilation occur around 72,000, 96,000, and 120,000 years ago,\ntimes when the boreal summer insolation was reduced. Yao et al. even develop an upwelling\nindex that is based on the same approach as Galbraith et al., yet they come to the opposite\nconclusion.\nReply:\nWe thank the reviewer for the comments, and appreciate to have the opportunity to clarify\non these points.\nAfter a thorough examination of the Galbraith et al. (2008) paper and a comparison of\ntheir data with insolation (which was not conducted in their paper), it appears that there is an\ninaccurate interpretation in the reviewer's assessment of Figure 6 in Galbraith et al. (2008). The\nfigure below is a plot that we have made to compare directly the record of Galbraith et al. (2008)\n(a) with summer insolation (b), our upwelling index (c) and an upwelling index in the Bering\nSea (d). It shows that the peaks of strong upwelling (indicated by low δ15N values) in Galbraith\net al. (2008) during MIS-5e and MIS-5c align well with summer insolation peaks. The peak of\nMIS-5a lags insolation peak by ~5 kyr, possibly related to chronological uncertainties. In our\nrecord (c) and the one from Bering Sea (d), the three strong upwelling peaks during MIS-5a, -\n5b and -5c all correspond well with high summer insolation. Therefore, all these three records\nsupport our conclusion that more intense upwelling is associated with high boreal summer\ninsolation. This relationship between upwelling and insolation can be extended to older periods\nof the past 600 kyr, as can be observed in Fig. 2d and Fig. 3c, d in the revised manuscript.\nMoreover, this relationship is confirmed by our model results as shown in our manuscript.\nComparison of upwelling index of core LV76-16-1 with other upwelling and climate\nrecords. (a) Sedimentary δ15N in core MD2416 in the NW Pacific (Galbraith et al., 2008). (b)\nSummer insolation at 65o in the northern Hemisphere (Berger and Loutre, 1991). (c)\nReconstructed upwelling index of core LV76-16-1 (this study). (d) Upwelling index at IODP\nSite U1343 in the Bering Sea (Worne et al., 2019). The vertical dashed lines represent peaks of\nboreal summer insolation.\nComment 2:\nThe problem comes down to the age model developed for core LV76-16-1 in relation to the\nage model developed for core MD2416 in Galbraith et al. (2008). Galbraith et al. at least had a\nbenthic d18O record that they were able to align with LR04. Yao et al. are stuck using %CaCO3\nand Ca/Ti ratios to align with LR04 instead. I don't see that Yao et al. refer to a benthic d18O\nrecord for LV76-16-1 anywhere in their paper. The approach used by Yao et al. is simply not\nas valid, especially when trying to align a core chronology with the ~23 kyr precessional cycle.\nReply:\nAs explained in our reply to Comment 1 of Reviewer 2, the problem claimed by the\nreviewer on the relationship between upwelling and insolation actually does not exist.\nNevertheless, here we would like to provide more clarification and evidence to demonstrate\nthe robustness of our age model.\nOur age model was established by associating abrupt increases in the Ca/Ti ratios with\nglacial terminations. The correlation is grounded in the observation that each deglacial\ntransition is marked by abrupt increases in the sedimentary accumulation of biogenic carbonate\nacross the entire North Pacific region (Jaccard et al., 2005; Kohfeld and Chase, 2011). This\nmethod of age model construction has been widely employed in the subarctic Pacific, as\ndemonstrated at sites such as ODP Site 882 (Jaccard et al., 2005) and core MD2416(Galbraith\net al., 2008). Furthermore, the downcore variations in the Ca/Ti ratio of core LV76-16-1 in our\nstudy closely resemble the variations in the Ca/Al ratio and Ca counts at the neighboring ODP\nSite 882 and core MD2416 (Extended Data Fig. 2), both of which have well-established\nchronology ((Jaccard et al., 2005; Galbraith et al., 2008). The consistency observed across\ndifferent cores confirms the robustness of our age model and underscores its suitability for\naddressing the issues raised in our manuscript.\nExtended Data Fig. 2. Comparison of Ca/Ti in core LV76-16-1 with other records. (a)\nCa/Ti ratio of core LV76-16-1 (this study). (b) XRF scanning Ca counts in core MD2416 in\nthe NW Pacific (Galbraith et al., 2008). (c) Ca/Al ratio at ODP Site 882 in the NW Pacific\n(Jaccard et al., 2005). (d) The LR04 benthic δ18O stack (Lisiecki and Raymo, 2005). The\nvertical yellow bars represent interglacial periods.\nReferences\nBerger, A. & Loutre, M-F. Insolation values for the climate of the last 10 million of years. Quat.\nSci. Rev. 10 (4), 297–317 (1991).\nCharles, C. D., Froelich, P. N., Zibello, M. A., Mortlock, R. A. & Morley, J. J. Biogenic opal\nin Southern Ocean sediments over the last 450,000 years: implications for surface water\nchemistry and circulation. Paleoceanography6, 697–728 (1991).\nDu, J., Haley, B. A., Mix, A. C., Walczak, M. H. & Praetorius, S. K. Flushing of the deep Pacific\nOcean and the deglacial rise of atmospheric CO concentrations. Nat. Geosci. 11, 749–\n2\n755(2018).\nFrancois, R. M., Frank, M. M. R., van der Loeff & Bacon, M. P. Th normalization: An essential\ntool for interpreting sedimentary fluxes during the late Quaternary. Paleoceanography19,\nPA1018 (2004).\nGalbraith, E. D. et al. Consistent relationship between global climate and surface nitrate\nutilization in the western subarctic Pacific throughout the last 500 ka. Paleoceanography\n23, 1–11 (2008).\nGray, W. R. et al. Deglacial upwelling, productivity and CO outgassing in the North Pacific\nJaccard, S. L. et al. Glacial/interglacial changes in subantarctic North Pacific stratification.\nScience308, 1003–1005 (2005).\nJaccard, S. L. et al. Two modes of change in Southern Ocean productivity over the past million\nyears. Science339, 1419–1423 (2013).\nKender, S. et al. Closure of the Bering strait caused mid-Pleistocene transition cooling. Nat.\nCommun. 9 (2018).\nKohfeld, K. E. & Chase, Z. Controls on deglacial changes in biogenic fluxes in the North\nPacific Ocean. Quat. Sci. Rev. 30, 3350–3363 (2011).\nLisiecki, L. E. & Raymo, M. E. A Pliocene-Pleistocene stack of 57 globally distributed benthic\nδ18O records. Paleoceanography20, 1–17 (2005).\nRae, J.W.B. et al. Deep water formation in the North Pacific and deglacial CO rise.\nSigman, D.M., Hain, M.P. & Haug, G.H. The polar ocean and glacial cycles in atmospheric\nCO concentration. Nature 466, 47–55 (2010).\nsubarctic Pacific nutrient upwelling over the last 850,000 years. Earth Planet. Sci. Lett.\n522, 87–97 (2019).\nYang, H. et al., Poleward shift of the major ocean gyres detected in a warming climate. Geophys.\nRes. Lett.47, e2019GL085868 (2020).\nYin Q. Z., Wu Z. P., Berger A., Goosse H. & Hodell D. Insolation triggered abrupt weakening\nof Atlantic circulation at the end of interglacials. Science373, 1035–1040 (2021)."}, {"round": 2, "reviewer_comments": {"reviewer_1": "Yao et al. have substantially revised their manuscript and addressed all my comments, especially my questions on cyclicities.\nThe manuscript is much more organized, readable and clear than the first version, and excellent in sections. With that said, I am accepting this version as ready for publication.\nI want to emphasize that this manuscript covers a lot of interesting and new grounds at the intersection of biological, physical, and climate changes of the recent glacial cycles in the Pacific, and I believe this should be published.", "reviewer_2": "The main complaint in my original review was about Yao et al.'s insistence that the strongest upwelling in the subarctic North Pacific occurs during maxima in the northern summer insolation. I identified work by Galbraith et al. (2008) that seemed to show the opposite tendency, i.e., stronger upwelling during minima in the northern summer insolation. I would agree that the upwelling in the subarctic North Pacific varies with the precessional forcing, but it is far from certain that the strongest upwelling occurs during the insolation maxima.\nMaxima and minima in the precessional forcing are separated by only 10-11,000 years. The choronologies of the cores in question, LV76-16-1 in this paper and MD2416 in the Galbraith et al. paper, are constrained by spikes in the Ca or Ca/Ti contents that occur every 100,000 years or so during the big deglaciations. The chronologies at all other times are based on an assumption that the sedimentation between the spikes is constant over time. It would not take much variability in the sedimentation between the deglaciations to obscure 10-11,000-yr differences in the timing of events.\nA third record has since come to my attention in Yao et al.'s response to Reviewer 1. This record is from IODP site U1343 in the Bering Sea, as described in Worne et al. (2019). Worne et al. see evidence of upwelling variations that are similar to the variations seen at the other sites. The sedimentation rate at site U1343 is higher than at the other locations, however, and the sediments were recovered from a shallower depth. The sediments also contain benthic forams from which a d18O record was constructed. This is important because Worne et al.'s d18O record provides time markers for the intervals between the deglaciations.\nAccording to the upwelling index in the red curve in Worne et al.'s Figure S3, the strongest upwelling during stages 5 and 7 in the Bering Sea occurred around 70, 95, 190, and 210 ka. As I was inclined to think in my original view, the strongest upwelling t site U1343 occurred when the when the northern summer insolation was relatively weak. Yao et al. use the same upwelling index, yet their results seem completely different. The differences, it seems to me, arise from the chronologies of the two cores. One is better resolved than the other.\nYao et al.'s entire paper is built around their certainty about the upwelling variations at site LV76-16-1. I don't see any justification for this level of certainty. My recommendation at this point is that the paper should be rejected.", "reviewer_3": "This is a comment to the manuscript by Yao et al. and specifically the reviewer’s comments, which show some opposite opinions after two rounds of review.\nReviewer 1 provides an extensive and constructive review, which is was answered in detail by Yao et al. Reviewer 2 has a problem with the age model presented by Yao et al. and argues that the paper by Galbraith et al. shows opposite results because of the age model problem. However, Yao et al. argue that the interpretation of the reviewer 2 was not correct and provide an explanation that appears valid.\nWhat additionally argues for me that the age model, and therewith the story, is valid, is that their age model and trends in productivity on themselves are not something new for the north Pacific. These productivity peaks have been identified at many locations following a similar mechanism, including reference to the Galbraith et al. paper. This is also where their age model is based on, i.e. comparing these peaks to the neighboring ODP Site 882 and align them to construct the age model.\nObviously, and also determined by setting, the age model would not be valid for high-resolution conclusions on leads and lags during the termination. But Yao et al. are also staying away from this.\nOne, as a minor side comment, thing that I noticed in the age model plot was that the variations in sedimentation rate are not very large due to the limited number of tie-points. Although difficult to constrain, these peaks in productivity probably increased sedimentation rates significantly during the time intervals that these conditions prevailed. But this would not change the correlation of the peaks with those at other locations in the north Pacific as these would have experienced the same changes.\nAs an important part of the paper, it is concluded that the precession has a significant impact on the cyclicity of the productivity and thus on CO2 release. What would be interesting to include in figure showing atmospheric CO2 is if this actually also contains the precession component. If the reservoir in the north Pacific is so large, this should be present in the overall signal as well."}, "author_response": "Reply to the reviewers’ comments on “Ice sheet and precession control on biological\nexport production and upwelling in the subarctic Pacific over the last 550,000 years” by\nYao et al.\nReviewer #1\nComment:\nYao et al. have substantially revised their manuscript and addressed all my comments,\nespecially my questions on cyclicities. The manuscript is much more organized, readable and\nclear than the first version, and excellent in sections. With that said, I am accepting this version\nas ready for publication. I want to emphasize that this manuscript covers a lot of interesting\nand new grounds at the intersection of biological, physical, and climate changes of the recent\nglacial cycles in the Pacific, and I believe this should be published.\nReply:\nWe sincerely appreciate your thoughtful and constructive comments during the initial\nreview, which has significantly enhanced the quality of our manuscript. Your positive\ncomments and the confirmation of acceptance for publication greatly encourage us.\nReviewer #2\nComment:\nThe main complaint in my original review was about Yao et al.'s insistence that the\nstrongest upwelling in the subarctic North Pacific occurs during maxima in the northern\nsummer insolation. I identified work by Galbraith et al. (2008) that seemed to show the opposite\ntendency, i.e., stronger upwelling during minima in the northern summer insolation. I would\nagree that the upwelling in the subarctic North Pacific varies with the precessional forcing, but\nit is far from certain that the strongest upwelling occurs during the insolation maxima.\nMaxima and minima in the precessional forcing are separated by only 10-11,000 years.\nThe chronologies of the cores in question, LV76-16-1 in this paper and MD2416 in the\nGalbraith et al. paper, are constrained by spikes in the Ca or Ca/Ti contents that occur every\n100,000 years or so during the big deglaciations. The chronologies at all other times are based\non an assumption that the sedimentation between the spikes is constant over time. It would not\ntake much variability in the sedimentation between the deglaciations to obscure 10-11,000-yr\n1\ndifferences in the timing of events.\nA third record has since come to my attention in Yao et al.'s response to Reviewer 1. This\nrecord is from IODP site U1343 in the Bering Sea, as described in Worne et al. (2019). Worne\net al. see evidence of upwelling variations that are similar to the variations seen at the other\nsites. The sedimentation rate at site U1343 is higher than at the other locations, however, and\nthe sediments were recovered from a shallower depth. The sediments also contain benthic\nforams from which a d18O record was constructed. This is important because Worne et al.'s\nd18O record provides time markers for the intervals between the deglaciations.\nAccording to the upwelling index in the red curve in Worne et al.'s Figure S3, the strongest\nupwelling during stages 5 and 7 in the Bering Sea occurred around 70, 95, 190, and 210 ka. As\nI was inclined to think in my original view, the strongest upwelling t site U1343 occurred when\nthe when the northern summer insolation was relatively weak. Yao et al. use the same upwelling\nindex, yet their results seem completely different. The differences, it seems to me, arise from\nthe chronologies of the two cores. One is better resolved than the other.\nYao et al.'s entire paper is built around their certainty about the upwelling variations at\nsite LV76-16-1. I don't see any justification for this level of certainty. My recommendation at\nthis point is that the paper should be rejected.\nReply:\nWe would like to extend our gratitude for your patience and for undertaking another\nreview of our manuscript. Following your comments, we have implemented the following\nrevisions and clarifications.\nChronology framework: To address concerns regarding the age model of our study core\n(LV76-16-1), we have enriched our manuscript with additional details to reinforce the\ncredibility of our age model. The chronology for ODP Site 882 was established by matching\nthe high-resolution XRF scanning Ba/Al ratios with the millennial-suborbital variations of δD\nin the Antarctic ice core (Jaccard et al., 2005; Galbraith et al., 2008). For core MD2416, the\nage model was developed by aligning XRF scanning Ca counts with those from ODP Site 882.\nThis methodology underscores that the age model for core LV76-16-1 has a high enough\nresolution to allow for discussions of suborbital-scale variations. While we recognize that there\nis some degree of uncertainty in our chronology and that of others, the consistent trends in\nproductivity and other parameters across numerous sites in the North Pacific lend strong\nsupport to the robustness of our age model. This information has been added in the revised\nmanuscript (P4, Lines 91-96).\n“Specifically, the chronology for ODP Site 882 was established by correlating high-\n2\nresolution XRF scanning Ba/Al ratios with the millennial-suborbital variability of δD from\nAntarctica ice core12,22. The age model for MD2416 were developed by aligning XRF scanning\nCa counts with those from ODP Site 88222. This demonstrates that the chronology of core\nLV76-16-1 has a high enough resolution to allow for discussions of suborbital-scale variations.”\nConcerns on the relationship between upwelling intensity and insolation: In light of\nthe feedback regarding the upwelling index presented by Worne et al., 2019, Reviewer #2\nclaimed that the most significant upwelling events in the Bering Sea during stages 5 and 7\noccurred around 70, 95, 190, and 210 ka, coinciding with periods of relatively low northern\nsummer insolation. We believe this interpretation might have been an oversight. To clarify,\nwe\nhave produced a plot comparing Worne's upwelling index with Northern Hemisphere summer\ninsolation, which clearly shows that strong upwelling phases mostly align with periods of\nstrong northern summer insolation (refer to Supplementary Figure S1 shown below). This\npattern is in agreement with observations from our study location, LV76-16-1. While there are\nslight discrepancies in certain periods, likely attributable to chronological precision, these do\nnot detract from the fundamental relationship observed. In our manuscript, we have already\nillustrated the comparison between the upwelling index at LV76-16-1 and northern summer\ninsolation (as seen in Figure 2d), alongside the correlation between Worne's U1343 site and\nLV76-16-1 (presented in Figures 3c and d). To maintain the manuscript's conciseness and avoid\nredundancy, we have chosen not to include Supplementary Figure S1 in the main text.\nIn conclusion, our findings and the methodological approach to establishing a\nchronological framework are robust, providing a sound basis for discussing upwelling\nvariations in relation to precession cycles.\nU1343 with boreal summer insolation. (a) Reconstructed upwelling index of core LV76-16-\n1 (this study; black line) and summer insolation at 65°N (Berger and Loutre, 1991; red line).\n(b) Upwelling index of IODP Site U1343 from the Bering Sea (Worne et al., 2019; blue line)\nand summer insolation at 65°N (Berger and Loutre, 1991; red line). The vertical green dashed\nlines represent 70 ka, 95 ka, 190 ka and 210 ka, four dates which were explicitly mentioned by\nReviewer #2.\nComment:\nThis is a comment to the manuscript by Yao et al. and specifically the reviewer’s\ncomments, which show some opposite opinions after two rounds of review.\nReviewer 1 provides an extensive and constructive review, which is was answered in\ndetail by Yao et al. Reviewer 2 has a problem with the age model presented by Yao et al. and\nargues that the paper by Galbraith et al. shows opposite results because of the age model\nproblem. However, Yao et al. argue that the interpretation of the reviewer 2 was not correct and\nprovide an explanation that appears valid.\nWhat additionally argues for me that the age model, and therewith the story, is valid, is\nthat their age model and trends in productivity on themselves are not something new for the\nnorth Pacific. These productivity peaks have been identified at many locations following a\nsimilar mechanism, including reference to the Galbraith et al. paper. This is also where their\nage model is based on, i.e. comparing these peaks to the neighboring ODP Site 882 and align\nthem to construct the age model.\nObviously, and also determined by setting, the age model would not be valid for high-\nresolution conclusions on leads and lags during the termination. But Yao et al. are also staying\naway from this. One, as a minor side comment, thing that I noticed in the age model plot was\nthat the variations in sedimentation rate are not very large due to the limited number of tie-\npoints. Although difficult to constrain, these peaks in productivity probably increased\nsedimentation rates significantly during the time intervals that these conditions prevailed. But\nthis would not change the correlation of the peaks with those at other locations in the north\nPacific as these would have experienced the same changes.\nAs an important part of the paper, it is concluded that the precession has a significant\nimpact on the cyclicity of the productivity and thus on CO release. What would be interesting\n2\nto include in figure showing atmospheric CO is if this actually also contains the precession\n2\ncomponent. If the reservoir in the north Pacific is so large, this should be present in the overall\nsignal as well.\nReply:\nThank you very much for your thorough and objective evaluation of our manuscript. We\ndeeply appreciate the constructive comments and the time you have dedicated to our work. We\nagree with your perspective on the age model and are grateful for your support of our\nexplanations against the concerns raised by Reviewer 2. Your deep understanding of our\n5\napproach reinforces the validity of our findings and the robustness of the age model we\nemployed, particularly in the context of the North Pacific's productivity peaks and their\ncorrelation with climate changes.\nFurthermore, in response to your insightful suggestion about examining the presence of a\nprecessional cycle in atmospheric CO , we have added a figure (Extended Data Figure 4d\n2\nshown below) in the revised manuscript. This figure presents a wavelet analysis of CO\n2\nvariations, and our results indeed show that CO changes contain a significant precessional\n2\ncomponent. This finding further strengthens our argument regarding the significant impact of\nprecession on productivity cycles and consequent CO release, underscoring the critical role of\nWe believe that this additional analysis not only addresses your valuable feedback but also\nenriches the scientific discussion around our findings, providing a clearer picture of the\nintricate relationship between orbital forcing, productivity, and carbon cycles in the North\nPacific.\n2\n(a) Comparison between Ba/Ti (dark line), 20-kyr filtering (blue line), Precession (red line)\nand Eccentricity (green line)60. (b) Wavelet analysis of Ba/Ti from core LV76-16-1. (c)\nWavelet analysis of upwelling index from core LV76-16-1. (d) Wavelet analysis of pCO 1,66.\n2\nBerger, A. & Loutre, M-F. Insolation values for the climate of the last 10 million of years. Quat.\nSci. Rev. 10 (4), 297–317 (1991).\nGalbraith, E. D. et al. Consistent relationship between global climate and surface nitrate\nutilization in the western subarctic Pacific throughout the last 500 ka. Paleoceanography\n23, 1–11 (2008).\nJaccard, S. L. et al. Glacial/interglacial changes in subantarctic North Pacific stratification.\nScience308, 1003–1005 (2005).\nWorne, S., Kender, S., Swann, G. E. A., Leng, M. J. & Ravelo, A. C. Coupled climate and\nsubarctic Pacific nutrient upwelling over the last 850,000 years. Earth Planet. Sci. Lett.\n522, 87–97 (2019)."}]

Decision letter