Introduction

Knowledge of the temperature variability of the past two millennia provides critical natural context and long-term perspective for the assessment of recent climate excursions1,2. Much research has leveraged the temperature-sensitive proxy records from tree-rings, ice cores, laminated sediments, corals, speleothems, and historical documentary archives to understand the past climate changes3,4,5,6,7,8. Over the last 2000 years, temperature reconstructions of the Northern Hemisphere (NH) exhibit a coherent long-term cooling from the Medieval Climate Anomaly (MCA; ca. 1000–1300 CE) to the Little Ice Age (LIA; ca. 1400–1850 CE), followed by abrupt industrial-era warming3,4,5,6. However, the majority of temperature proxies used to reconstruct this cooling trend are biased toward warm seasons9. Unfortunately, cold-season temperature reconstructions spanning the last two millennia are extremely rare, hampering our comprehensive understanding of the temperature evolution and its forcing mechanisms.

The LIA, which was used originally as a description of the most recent period of glacial advance10, has long been considered one of the coldest periods of the Holocene in the NH. Numerous proxy-based reconstructions and model simulations have attributed the LIA cooling at decadal to multidecadal timescales at least partly to decreased solar irradiance and/or increased volcanic activity7,8,11,12,13,14. On centennial timescales, the long-term cooling that lasted through the MCA and into the LIA was sustained by an orbitally-driven reduction in summer insolation15 and weakening of the North Atlantic subpolar gyre (SPG) circulation in response to increased Arctic freshwater export16,17,18,19. In cold winter-spring seasons, however, the NH’s insolation has steadily increased from the MCA to the LIA and the SPG strength is largely driven by the prevailing Atlantic atmospheric circulation via wind-stress forcing, such as winter storm conditions16,18,19. These differences in seasonal forcings create additional complexities for holistically understanding the temperature evolution of the last 2000 years. Remarkably, a recent reconstruction of North Atlantic winter-spring temperature from a freshwater lake sediment record in southwest Iceland has shown an overall warming from the MCA to the LIA20, in sharp contrast to the NH summer cooling21. However, the underlying forcing mechanism behind this seasonal warming trend remains unclear, and it is unknown whether it is spatially isolated or can also be documented in other regions.

The Tibetan Plateau (TP), the world’s highest and largest plateau, is highly susceptible to global climate changes, especially via atmospheric teleconnection associated with Atlantic climate variability22,23,24. In this study, we present a ~2000 year long, sub-decadal- to multidecadal-scale record of cold-season temperature from Lake Ngoring (Fig. 1) on the northeastern TP (NETP). Our temperature reconstruction is based on the unsaturation index (\({{{{{{\rm{U}}}}}}}_{37}^{{{{{{\rm{K}}}}}}}\)) of long-chain alkenones produced by the phylogenetically classified Group 1 Isochrysidales25,26,27,28, whose growth blooms during the spring transitional cold season in freshwater lakes27,29,30,31. We focus on temperature variability during the MCA and LIA, and compare our record with a suite of previously published temperature reconstructions, including regional- to NH-scale summer temperatures and North Atlantic sea surface temperatures. Our data and analysis provide insights into underlying forcing mechanisms of winter-spring temperatures through the MCA into the LIA.

Fig. 1: Map of study site.
figure 1

a Topographic map of the Tibetan Plateau and the ___location of study site. The map was created using BIGEMAP software and the geographic information data are from National Platform for Common Geospatial Information Services (https://www.tianditu.gov.cn). b Contour map showing the ___location of the sediment core (red star) collected from Lake Ngoring and the lake’s hydrological context. The contour map was created using Surfer software and the geographic information data are from National Platform for Common Geospatial Information Services (https://www.tianditu.gov.cn).

Full size image

Results and discussion

A Group 1 alkenone-based record of cold-season temperatures

To reconstruct cold-season temperature changes, we use the \({{{{{{\rm{U}}}}}}}_{37}^{{{{{{\rm{K}}}}}}}\) paleotemperature proxy (Eq. (1)), which is based on alkenone assemblages produced by Group 1 Isochrysidales. By conducting field sampling throughout the course of a seasonal cycle, numerous genomic and geochemical studies have demonstrated that Group 1 Isochrysidales bloom during the spring transitional season in freshwater lakes27,29,30,31, which may be triggered by lake ice break-up due to increased light intensity of lake waters. This indicates the generally consistent seasonality of Group 1 growth thriving in the relatively cold condition. Similarly, in Lake Ngoring, the timing of lake ice-off began in late April of spring transitional season32. The lake temperatures during the spring transitional season in the relatively cold regions are mainly dependent on lake ice thickness in the winter and early spring, as well as the timing and duration of spring ice melt27,28, which has been validated by a thermodynamic lake model to respond to winter-spring air temperature changes20,33. In Lake Ngoring, a lake model study has also shown that the lake temperatures at that time are dominantly influenced by the winter-spring air temperature changes32. Further, measurements of the Group 1 \({{{{{{\rm{U}}}}}}}_{37}^{{{{{{\rm{K}}}}}}}\) proxy from surface sediments in freshwater lakes of the northern mid- and high-latitudes reliably reflect the cold-season air temperature signal27,28, indicating the potential wide applicability of this proxy in recording this temperature signal. The proxy has also been successfully used to reconstruct past winter-spring temperatures in Arctic Alaska and Iceland20,33. Thus, the Group 1 \({{{{{{\rm{U}}}}}}}_{37}^{{{{{{\rm{K}}}}}}}\) proxy can be used as an indicator of cold-season temperature changes in our Lake Ngoring sediment record based on the above ecological and environmental reasons.

The Group 1 Isochrysidales can produce a distinctive C37:3b alkenone (tri-unsaturated isomer with Δ14,21,28 double bond positions) as a chemotaxonomic biomarker, in addition to the common C37:3a with Δ7,14,21 double bond that is also produced by Group 2 Isochrysidales thriving in saline lakes27,28,29. The isomer-based RIK37 index (Eq. (2)) can be used to infer the relative contributions of alkenones produced by Groups 1 and 2 Isochrysidales26,27,34. RIK37 values of 0.51–0.60 correspond to Group 1 Isochrysidales, whereas a value of 1 is diagnostic of Group 2, and an intermediate value of 0.60–0.99 to mixed Groups 1 and 226. Thus, higher RIK37 values reflect a greater input of Group 2 alkenones to lake sediments26,27,34.

In our sediment record from Lake Ngoring, the alkenones contain the C37:3b isomer (Supplementary Fig. 1) and RIK37 values vary between 0.58 and 0.78 (Supplementary Fig. 2a), indicating that alkenone contributions come from dominant Group 1 or mixed Groups 1 and 2. To reduce the potential impact of Group 2 mixing to the cold-season temperature signal derived from Group 1 alkenones, we removed the sediment samples with a relatively high input of Group 2 alkenones from the reconstruction by using an RIK37 cut-off of 0.7 (Supplementary Fig. 2b). When the samples with RIK37 > 0.7 are excluded, the filtered \({{{{{{\rm{U}}}}}}}_{37}^{{{{{{\rm{K}}}}}}}\) values are not correlated with RIK37 index (RIK37 ≤ 0.7; R2 = 0.06; Supplementary Fig. 3), indicating that a small amount of Group 2 alkenone mixing does not significantly affect the \({{{{{{\rm{U}}}}}}}_{37}^{{{{{{\rm{K}}}}}}}\) temperature record. Moreover, alkenone-producing Isochrysidales in lacustrine environments are only composed of two phylogenetically distinct Group 1 and Group 225,35, our record is thus not affected by other alkenone sources. This same approach (filtering the samples with RIK37 > 0.7) was applied to a 2000-year sediment record from southwest Iceland, which was demonstrated to be resistant to the phylotype mixing effects while still reliably recording changes in winter-spring temperatures from Group 1 alkenones20. We thus interpret the corrected \({{{{{{\rm{U}}}}}}}_{37}^{{{{{{\rm{K}}}}}}}\) values in our sediment core as a record of winter-spring temperature variability from the dominant Group 1 alkenones.

Winter-spring temperature variability of the Tibetan Plateau and a climate teleconnection

We employ three existing in-situ Group 1 \({{{{{{\rm{U}}}}}}}_{37}^{{{{{{\rm{K}}}}}}}\)-temperature calibrations based on suspended particulate matter and sediment trap samples from three freshwater lakes26,29,30 (Supplementary Fig. 4a), and calculate the dominant Group 1 \({{{{{{\rm{U}}}}}}}_{37}^{{{{{{\rm{K}}}}}}}\)-inferred temperatures from our sediment record (Supplementary Fig. 4b). Due to differences in y-intercept and slope of the three calibrations from different lakes, the calculated temperatures show obvious differences in our sediment record (Supplementary Fig. 4b). We acknowledge three potential factors influencing the calculated absolute temperatures: (1) differences in temperature responses of Group 1 \({{{{{{\rm{U}}}}}}}_{37}^{{{{{{\rm{K}}}}}}}\) for different lakes (Supplementary Fig. 4a); (2) the \({{{{{{\rm{U}}}}}}}_{37}^{{{{{{\rm{K}}}}}}}\) values in our sediment record (~0.1–0.3) beyond the range of the existing calibrations (\({{{{{{\rm{U}}}}}}}_{37}^{{{{{{\rm{K}}}}}}}\): −0.7–−0.3; Supplementary Fig. 4a); (3) a small input of Group 2 alkenones as inferred by the RIK37 proxy (Supplementary Fig. 2). Thus, we focus here on the changes in \({{{{{{\rm{U}}}}}}}_{37}^{{{{{{\rm{K}}}}}}}\) values in our sediment record.

The \({{{{{{\rm{U}}}}}}}_{37}^{{{{{{\rm{K}}}}}}}\) values in our sediment record abruptly rise before the MCA, followed by an overall increase but at a slower pace across the MCA to the LIA (Fig. 2a). During the LIA, a decrease between 1560–1650 CE is superimposed on the increasing trend. Overall, our 2000-year record shows a long-term warming (Fig. 2a), which is also reflected by the Ensemble Empirical Mode Decomposition analysis (Fig. 2d). The long-term trend may mainly respond, or in part, to increasing winter-spring insolation (Fig. 2g). The Ensemble Empirical Mode Decomposition analysis also decomposes our \({{{{{{\rm{U}}}}}}}_{37}^{{{{{{\rm{K}}}}}}}\) time series into two relatively main intrinsic mode functions (IMFs) IMF1 and IMF2 (Fig. 2e, f; note that the contributions of other IMFs to the total variance are very small; <10%). The IMF1 component appears to display ~100-year cycles of temperature variability, which approximately corresponds to the 80- to 90-year Gleissberg solar cycles36. The IMF2 contains only a full ~800-year cycle. However, the contributions of the IMF1 and IMF2 to the total variance are relatively small (14.5 and 13.9%; Fig. 2e, f). Thus, higher-resolution reconstructions of cold-season temperatures in the future are required to further confirm the observed centennial to multi-centennial signals.

Fig. 2: Cold-season temperature variations of the past two millennia and ensemble empirical mode decomposition analysis.
figure 2

a \({{{{{{\rm{U}}}}}}}_{37}^{{{{{{\rm{K}}}}}}}\) values from Lake Ngoring on NETP (this study; red line is a 5-point running mean). b \({{{{{{\rm{U}}}}}}}_{37}^{{{{{{\rm{K}}}}}}}\) values from Vestra Gíslholtsvatn in southwest Iceland20 (light blue line is a 3-point running mean). c \({{{{{{\rm{U}}}}}}}_{37}^{{{{{{\rm{K{{\hbox{'}}}}}}}}}}\) values from Lake Sihailongwan in northeastern China37 (blue-green line is a 3-point running mean). d–f Intrinsic temporal components (including intrinsic mode functions (IMFs)) of the Lake Ngoring \({{{{{{\rm{U}}}}}}}_{37}^{{{{{{\rm{K}}}}}}}\) record by ensemble empirical mode decomposition. Before the analysis, the raw data was interpolated at 5 years time slices using a spline method. g Winter-spring (December–May) and summer (June–August) insolation at 35°N. Black dashed lines depict long-term linear trends in temperature changes (as derived from linear regressions of \({{{{{{\rm{U}}}}}}}_{37}^{{{{{{\rm{K}}}}}}}\) or \({{{{{{\rm{U}}}}}}}_{37}^{{{{{{\rm{K{{\hbox{'}}}}}}}}}}\) values as a function of time). Gray and blue shading note the timing of the MCA and LIA, respectively.

Full size image

The observed warming trend in our record is broadly similar that of the Group 1 alkenone-based North Atlantic winter-spring temperature reconstruction recently reported from Vestra Gíslholtsvatn in southwest Iceland20 (Fig. 2b). We note that an alkenone-based long-term warming is also documented in a 1600-year sediment record from a freshwater volcanic lake (Lake Sihailongwan) of northeastern China37 (Fig. 2c). Our recent study has demonstrated the widespread occurrence of Group 1 alkenones in a series of volcanic lakes in this region28, and we thus interpret that the Lake Sihailongwan alkenone record may also reflect winter-spring temperature changes. The similarity of these long-term trends suggests a potential climate teleconnection between the North Atlantic region and the NETP (and northeastern China) on the timescale of our study.

Studies of model simulations, instrumental data, and paleoclimatic reconstructions have revealed the key role of Atlantic multidecadal variability (AMV) in modulating temperature changes over the TP as well as other geographic regions22,23,24,38, such as eastern Asia22,23,24 and the western tropical Pacific38. Our atmospheric general circulation model simulation clearly displays a zonal dipole pattern of sea level pressure (SLP) anomalies between the warm and cold AMV phases across the Atlantic-Eurasia region, with negative anomalies over the North Atlantic and positive anomalies over the TP (Fig. 3a). In response to the warm AMV anomaly, the negative SLP anomalies over the North Atlantic result from the strong ascending atmospheric motion in this region, whereas the positive (anticyclonic) anomalies over the TP are caused by the compensating subsidence24 (Fig. 3b). Such remote responses of atmospheric circulation may be the main physical mechanisms linking the AMV to climate impacts over the TP and other geographic region24,38.

Fig. 3: Climate teleconnection between the North Atlantic region and Tibetan Plateau.
figure 3

a Composite differences in the sea level pressure (SLP) between the warm (1995–2018) and cold AMV phases (1963–1994). b Schematic diagram of the climate teleconnection between the North Atlantic region and Tibetan Plateau. Warm North Atlantic temperature anomalies induce an east-west dipole pattern of SLP across the Atlantic-Eurasia region, with low-pressure anomalies over the North Atlantic and high-pressure anomalies over the Tibetan Plateau. This pattern generates strong upward motion and upper-level divergence over the North Atlantic, resulting in upper-level outflows that converge eastward to the Tibetan Plateau, and subsequently subside over the Tibetan Plateau. The base map in the schematic diagram was obtained from ETOPO1, a 1 arc-minute global relief model of Earth’s surface that integrates land topography and ocean bathymetry (https://ngdc.noaa.gov/mgg/global/relief/ETOPO1/tiled/).

Full size image

Opposing seasonal temperature trends across the MCA to LIA and potential driving mechanisms

The MCA and LIA are two notable climate features of the past two millennia1,2. Our winter-spring temperature record from the NETP, together with those from southwest Iceland20 (Fig. 2b) and northeastern China37 (Fig. 2c), exhibit an overall warming trend through the MCA into the LIA (Fig. 4a). In contrast, a tree-ring-based annual temperature record of the regional NETP39, a composite of summer temperature reconstructions from the broader TP40 (including the previously published lake records41,42,43), two composites of Chinese temperature reconstructions44,45, and a large-scale NH summer temperature reconstruction21 all consistently show a long-term cooling over this same time period (Fig. 4b–e). Similarly, in the geographically remote North Atlantic region, this cooling trend has also been observed in a regional sea surface temperature reconstruction5 (for which there is a warm season bias for the majority of the proxy records; Fig. 4f). The contrasting seasonal patterns of temperature trends suggest that seasonal biases should be considered to fully understand the temperature evolution of the past two millennia.

Fig. 4: Regional- to hemisphere-scale temperature variations and climate forcings through the MCA into the LIA.
figure 4

a \({{{{{{\rm{U}}}}}}}_{37}^{{{{{{\rm{K}}}}}}}\) values from Lake Ngoring on NETP (this study; red line is a 3-point running mean). b Composite of Group 2 alkenone-based summer temperature anomalies on the Tibetan Plateau (TP)40. c Tree ring-based annual temperature anomalies on NETP39. d Composite temperature anomalies from China based on multiple proxy-based studies44,45. e Tree ring-based NH summer temperature anomalies21. f Sea surface temperature (SST) anomalies in the Atlantic Multidecadal Oscillation (AMO) region5. g Relative abundance of hematite-stained grains (%HSGs) in the Labrador Sea marine sediment record18. h Marine-sourced sea salt sodium (Na+) concentrations in the Greenland ice core record49. i North Atlantic Oscillation (NAO) index51,52. j Winter-spring (December–May) and summer (June–August) insolation at 35°N. Black dashed lines depict long-term linear trends (as derived from linear regressions of each record as a function of time). Gray and blue shading note the timing of the MCA and LIA, respectively.

Full size image

Summer (or warm season) cooling through the MCA into the LIA could be partly attributed to an orbitally-driven reduction of summer insolation at the top of the atmosphere15 (Fig. 4j; 35°N as an example). In addition, increasing evidence suggests that the LIA cooling was caused by a weakening of the North Atlantic SPG circulation16,17,18,19, which reduced northward ocean heat transport from low to high latitudes19. Observational and modeling studies have shown that the freshening of the upper Labrador Sea is associated with an overall reduction in deep water formation due to a decrease in vertical ocean density gradients in the region, which plays a major role in weakening the SPG19,46. Such a freshwater anomaly can be caused by the southward export of Arctic freshwater from melting sea ice and/or glaciers, mainly via the East Greenland Current (Supplementary Fig. 5). From the MCA to the LIA, a marine sediment record in the Labrador Sea shows an overall increase in the relative abundance of hematite-stained grains18 (Fig. 4g). This indicates an intensification in the contribution of sea ice exported from the Arctic region, eventually leading to SPG slowdowns18 (Supplementary Fig. 5).

However, winter-spring temperatures in both the North Atlantic and NETP exhibit overall warming through the MCA into the LIA (Figs. 4a and 2b). In contrast to the decrease in summer insolation of the extratropical NH, winter-spring insolation steadily increases through the MCA into the LIA (Fig. 4j; 35°N as an example). The seasonality of insolation trends at the middle latitudes is caused by the precession of the equinoxes depending on Earth’s wobbling of the axis and turning of its elliptical orbit. We propose another potential mechanism for the observed warming trend. During the cold season, the strength of the SPG is largely controlled by wind stress forcing via heat losses from the surface of the Labrador Sea47,48. An increase in winter storm conditions over the North Atlantic can largely remove its sea surface heat, enhancing deep water formation of the Labrador Sea and intensifying the SPG16 (Supplementary Fig. 5b). From the MCA to the LIA, the intensity of winter storms over the North Atlantic increased, as indicated by a Greenland ice core record of marine-sourced sea salt sodium (Na+) concentrations49 (Fig. 4h). Note that the frequency of winter storms associated with the phase of the North Atlantic Oscillation (persistent positive phase of North Atlantic Oscillation could enhance winter storminess)50 has not undergone a long-lasting shift from the MCA to the LIA51,52 (Fig. 4i). We propose, therefore, that the increase in winter storm intensity leading to enhanced SPG may counterbalance or overwhelm the freshwater-induced SPG slowdowns (Supplementary Fig. 5). This mechanism, together with increasing winter-spring insolation, may explain our observed winter-spring warming trend during the MCA and LIA, in contrast to the summer cooling.

Over the past millennia, solar irradiance and volcanic eruptions have been considered two main forcings affecting mean annual and warm-season temperature variability on decadal to multidecadal timescales8,11,12. We perform two cross-wavelet analysis to assess the role of volcanic and total solar irradiance forcings on our \({{{{{{\rm{U}}}}}}}_{37}^{{{{{{\rm{K}}}}}}}\) record on different timescales (Supplementary Fig. 6). Our result shows that the \({{{{{{\rm{U}}}}}}}_{37}^{{{{{{\rm{K}}}}}}}\) values have a positive correlation with volcanic forcing during the period of ~1200–1350 CE on multidecadal timescales (Supplementary Fig. 6a), but have a negative correlation with total solar irradiance during some periods of ~600–1600 CE (Supplementary Fig. 6b). The negative correlation is in contradiction with the impacts of solar activity on temperature changes, suggesting that total solar irradiance cannot account for the multidecadal-scale temperature variability in our record. During ~1200–1350 CE, our observed volcanic-induced cooling is associated with the largest volcanic signal occurring in 1257 CE (Samalas eruption, Indonesia53) over the past two millennia (Supplementary Fig. 7). We therefore suggest that the cooling response of our regional multidecadal-scale cold-season temperature variations may be only sensitive to very large global eruptive episodes.

We note that winter-spring, summer, and annual temperatures on the NETP show a relatively cold episode between 1560–1650 CE during the LIA (Fig. 4a–c). This cooling event was also recorded by a record of sea surface temperature variability in the North Atlantic (Fig. 4f). The cold anomalies cannot be explained by the solar or volcanic radiative forcing, as solar irradiance did not decrease and volcanic activities did not substantially increase at that time (Supplementary Fig. 7). This event is generally consistent with the massive Arctic freshwater discharges from melting sea ice between ~1500–1600 CE (Fig. 4g), although there is a slight mismatch of the chronologies, possibly due to dating uncertainties of the sediment records. We thus infer that the cooling event of 1560–1650 CE may be associated with the Arctic meltwater pulses. The meltwater pulses may have led to substantial weakening of the SPG throughout the course of a year, eventually causing the cooling teleconnection during this period of the LIA.

Overall, our winter-spring temperature reconstruction of the past two millennia indicates a long-term and pronounced warming on the NETP, likely reflecting high-elevation amplification. The TP, as the Asian Water Tower, provides abundant freshwater resources for downstream populations of more than 1.4 billion, due to its vast areas of glaciers and snow54. If the warming amplification continues in the future, this could accelerate the regional retreat of glaciers as well as spring snow melt, which could affect the downstream water availability and food security.

Materials and methods

Study site and sampling

Lake Ngoring (34°46′–35°50′N, 97°32′–97°54′E; 4,272 m a.s.l; Fig. 1), located in the alpine environment of the NETP, is the largest freshwater lake (salinity measured on August 2015 CE was ~448 mg l–1)55 in the source region of the Yellow River (the second largest river in China). The Yellow River flows through the lake from the southwest to the northeast. The lake has a surface area of ~610 km2 and a catchment area of ~18,188 km2. The mean water depth is ~18 m, and the maximum water depth is ~31 m. In recent decades, the lake has generally maintained hydrological balance, with annual inflow of ~14.43 × 108 m3, outflow of ~6.36 × 108 m3, and evaporation of ~8.07 × 108 m3 56. The regional mean annual air temperature was –3.7 °C, and the mean annual precipitation was 318 mm (1961–2010; from nearby Madoi station), with a maximum in summer and minimum in winter.

A 85 cm sediment core (35°3′N, 97°43′E; ELLC14-1) was retrieved from Lake Ngoring at a water depth of ~16 m in September 2014 CE (Fig. 1b). The core was subsampled at 0.5 cm intervals, and all samples were frozen at –20 °C in the laboratory prior to analysis. The chronology of the core has been well established by 137Cs-210Pb dating and a 14C model57.

Alkenone analysis

All sediment samples were freeze-dried and extracted by sonication (3×) with dichloromethane (DCM): methanol (MeOH) (9:1, v/v). Total lipid extracts were purified by column chromatography with silica gel using the following sequence of eluents: n-hexane, DCM, and MeOH. Following the method reported by Zheng et al.58, alkenones in the DCM fractions were analyzed by an Agilent 8890N GC system equipped with a flame ionization detector (FID) and a Restek Rtx-200 GC column (105 m × 250 µm × 0.25 µm) at Xi’an Jiaotong University, China. The following GC-FID oven program was used: initial temperature of 50 °C (hold 2 min), ramp 20 °C min–1 to 255 °C, ramp 3 °C min–1 to 312 °C (hold 35 min). Helium was used as the carrier gas and was held at a constant flow rate of 1.3 ml min–1. The well-defined alkenone profile from our previously reported Lake Wudalianchi samples28 was used as a standard for the comparison of alkenone peaks. Samples with compounds co-eluting with alkenones were further purified using column chromatography with sliver thiolate silica59 with n-hexane: DCM (1:1, v/v) and acetone.

The alkenone unsaturation index \({{{{{{\rm{U}}}}}}}_{37}^{{{{{{\rm{K}}}}}}}\)60 and alkenone isomer-based RIK3726 were calculated as follows:

$${{{{{{\rm{U}}}}}}}_{37}^{{{{{{\rm{K}}}}}}}=\frac{{{{{{{\rm{C}}}}}}}_{37:2}-{{{{{{\rm{C}}}}}}}_{37:4}}{{{{{{{\rm{C}}}}}}}_{37:2}+{{{{{{\rm{C}}}}}}}_{37:3}+{{{{{{\rm{C}}}}}}}_{37:4}}$$
(1)
$${{{{{{\rm{RIK}}}}}}}_{37}=\,\frac{{{{{{{\rm{C}}}}}}}_{37:3{{{{{\rm{a}}}}}}}}{{{{{{{\rm{C}}}}}}}_{37:3{{{{{\rm{a}}}}}}}+{{{{{{\rm{C}}}}}}}_{37:3{{{{{\rm{b}}}}}}}}$$
(2)

where the “a” and “b” subscripts refer to the Δ7,14,21 and Δ14,21,28 tri-unsaturated alkenones, respectively.

Statistical analyses and model experiments

We used an Ensemble Empirical Mode Decomposition analysis61 to extract the IMFs of our \({{{{{{\rm{U}}}}}}}_{37}^{{{{{{\rm{K}}}}}}}\) record. The statistical analyses were performed by R hht package62.

Following the method from Shi et al.24, we performed an atmospheric general circulation model simulation using the observed SLP composite analysis between the warm and cold AMV phases (1995–2018 minus 1963–1994) in the MATLAB software. The interannual high-frequency variability (10 year) was denoised and the global warming trend of sea surface temperature was removed before we performed the composite analysis. The SLP during 1950–2018 CE was obtained from NCEP/NCAR R1 dataset (https://www.psl.noaa.gov/data/gridded/data.ncep.reanalysis.html). Sea surface temperature was obtained from ERSST v5 dataset (https://www.ncei.noaa.gov/products/extended-reconstructed-sst).