Central Asian modulation of Northern Hemisphere moisture transfer over the Late Cenozoic

Earth’s climatic evolution over the last 5 million years is primarily understood from the perspective of marine mechanisms, however, the role of terrestrial feedbacks remains largely unexplored. Here we reconstruct the last 5 million years of soil moisture variability in Central Asia using paleomagnetism data and isotope geochemistry of an 80 m-thick sedimentary succession at Charyn Canyon, Kazakhstan. We identify a long-term trend of increasing aridification throughout the period, along with shorter-term variability related to the interaction between mid-latitude westerlies and the Siberian high-pressure system. This record highlights the long-term contribution of mid-latitude Eurasian terrestrial systems to the modulation of moisture transfer into the Northern Hemisphere oceans and back onto land via westerly air flow. The response of Earth-surface dynamics to Plio-Pleistocene climatic change in Central Asia likely generated terrestrial feedbacks affecting ocean and atmospheric circulation. This missing terrestrial link elucidates the significance of land-water feedbacks for long-term global climate. Late Cenozoic variation in Central Asian hydroclimate resulted from the interaction between mid-latitude westerlies and the Siberian high-pressure system and may have driven terrestrial feedbacks, according to analyses of sediments from Charyn Canyon, Kazakhstan.

O ver the last 5 million years (Ma), Earth has evolved from a largely ice-free, warmer world sustained by high atmospheric carbon dioxide concentration, to a cooler planet oscillating between ice ages and interglacials on hundred thousand years (100 ka) timescales 1 . The onset of glaciations, initiating in the Northern Hemisphere 3.6 million years ago (Ma) and intensifying between 3.3 and 2.7 Ma 2,3 , required cold temperatures at high latitudes combined with sufficient relocation of moisture onto land for snowfall and ice accumulation. In the early phases of this process, changes in Earth's obliquity triggered highlatitude cooling but were insufficient to sustain ice cover 2 . Only after 3.3 Ma 3,4 did polar summer insolation and atmospheric CO 2 reduce adequately for widespread ice sheet expansion, indicating the contribution of additional feedback mechanisms beyond orbital variability.
This sequence of events is presently best understood through the lens of marine mechanisms. The shoaling at 4.6 Ma of the Panama gateway 5,6 altered ocean dynamics by reduced mixing between Atlantic and Pacific waters, thereby strengthening Atlantic meridional overturning circulation 7 . This process transported warmer near-surface waters to the North Atlantic via the Gulf Stream 5 and increased moisture supply to the high latitudes. High-latitude glaciation is proposed to have become widespread following reduced global oceanic warm water exchange, due both to constriction of equatorial Indonesian throughflow from 3.5 Ma 8 and establishment of the Panama Isthmus by 2.8 Ma 9 .
It is unlikely that continental regions were passive during this process; understanding the contribution of land−ocean-feedbacks is imperative in order to elucidate global climate feedback loops over the late Cenozoic. Once hypothesis invokes moisture, transported by westerlies onto Eurasia, increasing the outflow of the large Siberian rivers 10 . The subsequent injection of freshwater to the oceans is proposed to have expedited sea ice formation at high latitudes 10,11 . Another potential feedback involves wholesale removal of surface sediment from the continents, which may have changed ice sheet response to orbital forcing, enhanced global dust flux, and reduced albedo at the mid-latitudes 12 .
Here, we interrogate the continental contribution to Northern Hemisphere climate feedbacks by reconstructing the last 5 Ma of moisture variability and landscape response in Arid Central Asia (ACA). Our record derives from an 80-m-thick sedimentary succession at Charyn Canyon located in the Ili Basin in southeast Kazakhstan ( Fig. 1, Supplementary Fig. 1, Supplementary Note 1, and Supplementary Table 1). This region is uniquely situated to investigate long-term land-climate feedbacks since its semiarid climate is driven by the interplay between the mid-latitude westerlies 13 and high-latitude polar front 14 (Fig. 1a). The East Asian Summer Monsoon is unlikely to have penetrated this far north beyond the mountains over these timescales 15,16 , and the site is situated too far east for the incursion of rainfall from the Indian monsoon system (Fig. 1a).
The Neogene story of ACA is one of aridification. Uplift of the Himalaya and Tian Shan to the south, and westward retreat of the Paratethys Sea 17 , mitigated monsoon penetration and drove longterm drying 18 . The region's longest and most continuous post-Paratethyan archives are predominantly aeolian and have accumulated as thick blankets of loess along the mountain piedmonts. This aeolian dust has been transported by air masses associated with the major climate subsystems of the region, chiefly the westerlies 15 , with additional contributions from northerly transport or local katabatic winds 13,19 . Despite the potential of its widespread aeolian deposits 20 to record palaeoclimatic information [21][22][23][24][25] , palaeoclimate datasets generated for ACA, thus are either coarse in timescale 15 or focus on recent, short-term intervals insufficient to inform us about long-term feedbacks 14,26 . Continuous archives of terrestrial response to atmospheric drivers in ACA over Plio−Pleistocene timescales so far remain largely unexplored. The Charyn Canyon sedimentary succession (43°16′ N, 78°59′ E) addresses this imbalance by providing a quasi-continuous record of late Cenozoic change for central Eurasia (Fig. 1b, Supplementary Fig. 1, and Supplementary Note 1). The sequence enables detailed interrogation of the continental contribution to the global climate system over these timescales. We explore how ACA landscape response to moisture availability contributed to global dust flux, and may represent a proxy for freshwater inflow to the Arctic. These two terrestrial feedbacks might have influenced sea ice formation and modulated the dynamics of the westerly air masses.

Results
The Charyn Canyon sequence records long-term aridification in ACA (Supplementary Note 1). It comprises 13 m of distal alluvial fan, composed of gravels and silts which accumulated during the early Pliocene ( Supplementary Fig. 2e), superposed by 24 m of alternating aeolian silts and alluvial sheetwash gravel deposits indicative of drier climates with short-lived wet pulses corresponding to the late Pliocene (Fig. 1b, Supplementary Fig. 2d). The uppermost 43 m consists of almost purely aeolian sediment which accumulated from 2.6 Ma (Fig. 1, Supplementary Fig. 2a-c). Detailed stratigraphic and facies description is reported in Supplementary Note 1 and grain size distributions ( Supplementary  Fig. 1c) are available in Supplementary Data 1. Low-field magnetic susceptibility is lithologically driven, with sharp increases corresponding to mafic gravel sheetwash layers (Fig. 1b, Methods, and Supplementary Data 2).
Our chronology for the Charyn Canyon deposits (Fig. 1b) was determined using palaeomagnetic analyses (Methods, Supplementary Data 3, Supplementary Fig. 3), combined with absolute 238 U− 206 Pb dating of carbonate nodules and a bovid tooth (Methods, Supplementary Data 4 and 5, and Supplementary  Fig. 4). We derived an age model using linear interpolation between the absolute age control points (Methods). The combined chronology and stratigraphy indicate a quasi-continuous record spanning ca. 5.0−0.5 Ma (Fig. 1b).
The timing of changes in sedimentation rates based on the age model ( Fig. 2a, Methods) appears to pace global climatic events. Sedimentation at Charyn Canyon increased ca. 3.3 Ma, coeval with increased dust input to the oceans 27 , and again around the onset of Northern Hemisphere glaciation at ca. 2.7 Ma. Aeolian dust sedimentation rates also increased from 0.8 Ma, coeval with the onset of hundred-thousand-year glacial-interglacial cycles during the Pleistocene (Fig. 2a). We assume on the basis of these results, combined with our observations of horizontal stratification and nature of river incision in the Ili Basin 28 , that climate, rather than tectonics, regulated sedimentation at Charyn Canyon.
We evaluated the origin of the calcium within pedogenic carbonates down the sequence using strontium isotope ratios, and strontium and calcium concentrations (Methods, Supplementary  Supplementary Fig. 5).
The carbon and oxygen isotope composition of pedogenic carbonates (Methods) provides complementary hydroclimate information. The carbon isotopic composition of the pedogenic carbonates (δ 13 C pc ) is primarily sensitive to the concentration of atmospheric and soil-respirated CO 2 15,18,29,30 . The 1-D production−diffusion model developed by refs. 18,30,31 indicates that as rainfall decreases, so too does soil respiration flux and δ 13 C values of soil-respirated CO 2 32,33 ; the depth of carbonate formation correspondingly increases within the stratigraphical layer 34 , resulting in greater δ 13 C pc values. Consequently, δ 13 C pc values provide a proxy for long-term moisture variability. Mean δ 13 C pc values in Charyn Canyon increase upward through the sequence overall, indicating a general aridification trend over the last 5 Ma. However, it is clear that the process of aridification was not linear and developed over multiple phases (Fig. 2b, Supplementary Data 6, and Supplementary Fig. 5). The earliest drying phase, from 3.3 to 2.0 Ma, is evidenced by increased sedimentation rates and δ 13 C pc values, and oversaw the establishment of a purely aeolian depositional regime ca. 2.7 Ma. Between 2.0 and 1.5 Ma, precipitation increased concomitant with a decrease in aeolian sediment flux. Arid conditions intensified from 1.5 Ma, as indicated by an increase of δ 13 C pc values (Fig. 2b).
By contrast, δ 18 O values of pedogenic carbonates (δ 18 O pc ; Methods) vary as a function of both soil temperature and soil water isotopic composition (δ 18 O sw ), which is assumed to be in isotopic equilibrium with the meteoric water; i.e. precipitation (δ 18  Clumped isotope (Δ 47 ) analysis directly measures the temperature of carbonate formation within the soil, which is independent of the δ 18 O composition of the meteoric water 38 Table 3, and Supplementary Data 7). These temperature estimates for the Pliocene and early Pleistocene fall within the range typical of modern near-surface temperatures for this region-mean annual temperature (MAT) is 8.4°C and summer temperature is 19.3°C 40 -which suggests that T Δ47 from pedogenic carbonates represent summer soil temperatures rather than MAT 36 . We calculated meteoric water δ 18 O mw values using paired T Δ47 and δ 18 O pc data (Supplementary Data 7). These vary by 3‰ and may suggest a change in moisture source over the Plio−Pleistocene.
A bovid molar found west of the main profile is stratigraphically correlated at −6.25 m and corresponds to a time slice at ca. 2.6 Ma (Methods, Supplementary Fig. 6, and Note 3). A direct dating of the enamel with U−Pb methods yields a weighted mean age of 2.65 ± 0.15 Ma, which supports our previous stratigraphically correlated age (Methods, Fig. 1b, Supplementary  Fig. 4). Limited stable isotopic variability (<2%; Supplementary Table 4) throughout the lifespan of the animal suggests a period of reduced seasonal variation combined with semiarid, relatively stable conditions. The climate regime at this time was likely to be less continental than the present day 41 .

Discussion
The Charyn Canyon sequence records late Cenozoic variability in moisture and sediment transport over ACA, as modulated by the interplay between westerly and northerly air masses. Presently, most rainfall at Charyn Canyon derives from the westerlies and their interaction with the spring and autumn latitudinal migrations of the polar jet 42,43 . Winter precipitation is limited by the expansion of the Siberian high-pressure system against the Tian Shan 44 (Fig. 1a). Over the last 5 Ma, an intensified Siberian high, triggered by reduced Arctic Oscillation 45 and/or weakened westerly circulation, would further reduce rainfall in ACA. In this way, Charyn Canyon could serve as a litmus test for the interaction of the Siberian high with the rain-bearing westerlies. We hypothesize that moisture variability in our record may be used as a proxy for Siberian river activity to the north of our site, since the hydroclimate regime at Charyn Canyon echoes that of the steppes to the north where a number of Siberian river headwaters, such as the Irtysh and Ob, are situated. These are likewise influenced by the dynamics of the Siberian high and by the westerly air masses 46,47 . The Charyn Canyon sediments preserve landscape response to changes in the relative influences of the westerlies and Siberian high, and provides an indication for the emergence of land-based feedbacks into the atmospheric circulatory systems.
Relatively warm and wet conditions prevailed in the early Pliocene (Fig. 2b, c). Rainfall most likely derived from westerly air masses bearing moisture boosted by an amplified North Atlantic current following the 4.6 Ma shoaling of the Panama seaway 6 . During the warm Pliocene period, the westerlies lay northward of their present trajectory 48 . Westerly dominance across ACA at this time was strengthened by a relatively weak Siberian high which facilitated persistent warmer temperatures at high latitudes 49 (Fig. 3, Step 1). Relatively high rainfall during this phase promoted fluvial activity at Charyn Canyon and was likely echoed by the Siberian rivers further north. Increased Siberian river outflow facilitated the freshening of Arctic ocean waters 11 and stimulated sea ice formation in the Atlantic. This terrestrial response to increased rainfall in Eurasia may have exceeded a tipping point for the largest Pliocene cooling event (marine isotope stage (MIS) M2: 3.3−3.2 Ma), and very likely set the scene for sustained Northern Hemisphere glaciation at 2.7 Ma 2,10 .
Cool, arid conditions, coeval with increased sedimentation rates, initiated at Charyn Canyon during the Mid-Pliocene (ca. 3.3 Ma) and persisted until ca. 2.0 Ma (Fig. 2b, c). These conditions were widespread across central Eurasia, having been recorded in sediment records to the northeast at Lake Baikal 50 and Arctic Siberia 49 , as well as to the south in the Tarim Basin 51,52 and to the west in Tajikistan 53 and northern Iran 54 . The distribution and nature of this phase likely reflects the intensification of the Siberian high, which brought cooler temperatures and effectively reduced westerly moisture transport over Eurasia (Fig. 3, Step 2). It is feasible that weakened more southward trajectories of westerly air-flow was linked to the M2 (3.3 Ma) and 2.7 Ma Northern Hemisphere glaciations 48 . These early glaciations are attributed to the coevolution of high latitude feedbacks, obliquity-driven climate shifts 55 , and wholesale sediment transport from the continents 12,56 . Once established, these glaciations oversaw decreased sea-surface temperatures (SST) in the North Atlantic (Fig. 2d, e), further reducing evaporation and weakening moisture transport onto continental Eurasia. Model simulations attribute cooler continental temperatures to increasingly active Arctic Oscillation and periodic penetration of cold Arctic air to the midlatitudes 57 . We propose that the redistribution of global water vapor linked to North Atlantic Ocean-atmosphere feedbacks, coupled with an intensified Siberian high, contributed to the conditions observed in ACA at this time (Fig. 3, Step 2). Reduced vegetation cover in ACA, in response to drier climates, would increase land-surface albedo and further decrease soil moisture through direct evaporation, creating an additional possible feedback. Drier conditions additionally increase aeolian entrainment and transport, increasing sedimentation rates across Eurasia 56 .
The early Pleistocene (ca. 2.0−1.5 Ma) oversaw an abrupt reversal into warmer, wetter conditions in ACA (Fig. 2b, c). These conditions are also recorded in Northern Iran by the development of thick palaeosol horizons which indicate milder climates 54 . Sedimentation rates remained high until ca. 1.8 Ma, followed by an abrupt decrease, which may reflect reduced sediment availability for transport and suggests a change in climate dynamics (Fig. 2). This amelioration phase corresponds to the global establishment of 41 kyr, orbitally forced glacial−interglacial cycles. The influence of astronomical forcing at this time is evident worldwide, dictating a change in North Atlantic thermohaline circulation and driving Arctic glacial melting during interglacials 55,58 . Coeval reduction in polar jet activity resulted in particularly warm interglacial temperatures in the Arctic 59 . This, along with a more northward path of strengthened westerlies 48 associated with increased evaporation over slightly warmer North Atlantic waters 55 (Fig. 2d, e), would likely have weakened the Siberian high, promoting warmer conditions and increasing precipitation at Charyn Canyon (Fig. 3 Step 3).
Aridification in ACA intensified from ca. 1.5 Ma and persisted through the Mid-Pleistocene transition (MPT; ca. 1.25−0.7 Ma) (Fig. 2b, c). This stepwise aridification is also reported in sediments of the Tarim basin 51 . The MPT was characterized by a shift in the periodicity of orbitally paced glacial−interglacial climate cycles, from 41 to 100 kyr. Its onset was accompanied by gradual global SST cooling 60 and substantially increased continental ice volume 61 . Models suggest that the transition from 41 to 100 kyr cycles was controlled by the combined effect of long-term declining atmospheric CO 2 and large-scale regolith removal from the continents, which altered surface albedo and ice-sheet behavior 12 . We propose that the increased sedimentation rates observed at Charyn Canyon are emblematic of the modeled terrestrial feedbacks driving this cooler world. Colder conditions and widespread glaciation furthermore led to a decrease in the evaporation potential of moist air masses originating over the North Atlantic, resulting in a more southward trajectory of the westerlies 48 and a reduced moisture transport over Eurasia. Concomitant intensification of the Siberian high further exacerbated ACA aridity (Fig. 3 Step 4).
To conclude, the Charyn Canyon sequence provides the first clear evidence for long-term interactions between mid-and highlatitude climate systems over the Eurasian continent. At key periods over the last 5 Ma, the response of ACA land surfaces to climate dynamics likely generated terrestrial feedbacks affecting atmospheric and ocean circulation, leading to tipping points in the regulation of moisture transfer and sediment flux. Since the Pliocene is our best analog for Anthropocene climatic conditions 62 , the missing terrestrial mechanisms provided by the Charyn Canyon record are of essential significance for understanding future climate thresholds.

Methods
Stratigraphy. Sediment description and sampling were undertaken by abseil with fixed ropes to ensure continuous collection and observation of the sequence. The 80 m-thick composite section and position of samples along the profiles were located in three-dimensional topographic space using a Leica total station (Supplementary Note 1, Supplementary Fig. 1, and Supplementary Table 1). Grain size analysis of Charyn sediments was performed following recommendations by ref. 63 . A total of 55 samples, ranging in weight from 33 to 36 g and spaced 20−25 cm through the whole composite section, were treated with 10 ml of 30% H 2 O 2 in order to remove any organic matter content. After 24 h, samples were dried in an oven at 70°C and treated for the other 24 h with 10 ml of 10% HCl to remove carbonate content. Dried again at 70°C, samples were gently disaggregated in an agate mortar. Sample weight was measured after each step. Samples were then dispersed in 0.05% (NaPO 3 ) 6 solution and passed through a Microtrac MRB Bluewave analyzer, coupled with a Sample Dispersion Controller, at the CORE Laboratory of Instituto Oceanografico (USP, São Paulo) (Supplementary Data 1).
Low field-magnetic susceptibility was measured in situ using a Bartington Instruments MS3 magnetic susceptibility meter with MS2K sensor (Supplementary Data 2), and discrete samples have been replicated in the IAG laboratory at the Universidade de São Paulo (USP, Brazil) with an AGICO Kappabridge MFK-FA.
Palaeomagnetism. Fine-grained and cohesive sediments for palaeomagnetic analyses were sampled in the field using a non-magnetic cordless drill yielding 2.5 cm diameter cores. The sampling positions were all placed within threedimensional space using the total station and palaeomagnetic samples were oriented in the field with a magnetic compass (Supplementary Data 3). Palaeomagnetic analyses were performed at the Alpine Laboratory of Palaeomagnetism (ALP, Peveragno, Italy) on a total of 121 cores collected at 50 cm intervals along the sequence. The structure of the natural remanent magnetization (NRM) was investigated for both thermal and alternating fields (AF) progressive demagnetization in 10−15 steps, using a 2G-Enterprises 600 Degausser and an ASC TD48 furnace, respectively. The NRM was measured with a 2G-Enterprises 755 cryogenic magnetometer located in a magnetically shielded room (ambient field < 500 nT). Intensity generally ranged from 5 to 15 mA/m, with a few samples yielding up to 30 mA/m. Orthogonal projections of demagnetization data 64 indicate the existence in most samples of a low-unblocking temperature (AF soft) component (erased at 150-200°C or at 5-15 mT). This bears only normal polarity parallel to the present field, superimposed on a moderate-high unblocking temperature (AF hard) component, and can be completely removed at 500-580°C (50-60 mT) or 625-650°C ( Supplementary Fig. 3a). This second component bears both normal and reverse magnetic polarity and is interpreted as the characteristic magnetization (ChRM). The range of maximum unblocking temperatures indicates the main magnetic carriers as magnetite and hematite. Magnetization components were resolved from the demagnetization data using the least-squares method 65 applied to at least four linear points on the orthogonal projections. Reliable ChRM directions were isolated in 102 samples (84%) and show antipodal distributions (Supplementary Fig. 3b) that deviate from being coaxial by only 3.5°, giving a positive Class B reversal test (95% confidence angle = 8.3°6 6 ). The ChRM is interpreted to be the primary magnetization of the Charyn Canyon sediments, acquired at or soon after deposition (DRM). The virtual geomagnetic polar (VGP) latitudes (i.e. the latitude of the magnetic pole which accounts for the direction of the ChRM at the site) were used to build the magnetostratigraphy of Charyn Canyon (Supplementary Data 3). 238 U/ 206 Pb dating. The absolute constraint of the palaeomagnetic age model was undertaken using 238 U− 206 Pb dating of pedogenic carbonates and a bovid tooth. Four large pedogenic carbonate nodules from key layers within the magnetostratigraphy and a bovid tooth stratigraphically correlated to −6.25 m and corresponds to a time slice at ca. 2.5 Ma, were prepared for measurement on polished thick sections. Isotope measurements were made using laser ablation-sector fieldinductively coupled plasma-mass spectrometry (LA-SF-ICPMS) 67 Fig. 4). Uranium, thorium, and lead isotopes were analyzed using a ThermoScientific Neptune+ MC-ICP-MS coupled to a RESO-Lution S-155 (Resonetics) 193 nm ArF Excimer laser (CompexPro 102, Coherent), equipped with a two-volume ablation cell (Laurin Technic, Australia). Results are summarized in Supplementary Data 4 for the samples and for the reference materials and Supplementary Data 5 for the data report. Samples were ablated in a helium atmosphere (0.3 l/min) and mixed in the ablation funnel with 1.0 l/min argon and 0.08 l/min nitrogen. Static ablation used a spot size of 104 µm in diameter for carbonate samples and 50 µm for phosphate samples, a crater depth of 15 µm and a fluence of ca. 2 J/cm 2 at 8 Hz for carbonates samples and 2 J/cm 2 at 6 Hz for phosphate samples. The detection limits for 206 Pb and 238 U were~0.1 and 0.03 ppb, respectively. Prior to analysis, each spot was pre-ablated (1 µm) to remove surface contamination. For each spot, the analysis comprised 20 s of background acquisition followed by 18 s of sample ablation. Soda-lime glass SRM-NIST614 and -NIST612 were used as reference glasses together with four carbonate and two apatite standards, to bracket sample analysis (Supplementary Data 4). Raw data were corrected offline using an in-house MS Excel © spreadsheet program 67,68 . Following background correction, outliers (±2 s) were rejected based on the timeresolved 207 Pb/ 206 Pb and 206 Pb/ 238 U ratios together with the Pb and U signal. SRM-NIST 612 was used as the primary standard for carbonate and phosphate samples to correct mass bias for the 207 Pb/ 206 Pb ratio (0.4%) and the final driftcorrected 206 Pb/ 238 U ratio for an inter-element fraction (ca. 30%). The WC-1 carbonate 69 (254 Ma) and Durango apatite 70 (31.44 Ma) reference materials were used to determine the effect of the inter-element fractionation when injecting a carbonate and phosphate matrix, respectively, relative to a synthetic glass (SRM-NIST) in the ICP. In the case of carbonate samples, the calculated offset of the 206 Pb/ 238 U between those two matrices were 2 and 4%, respectively. For the first sequence, ASH-15D 71  Age model and calculation of sedimentation rates. Our age model was built on the magnetostratigraphy of the Charyn Canyon sequence, tied to absolute ages produced by 238 U− 206 Pb laser ablation dating. We undertook a linear interpolation between nine successive age control points using the software AnalySeries 73 (Fig. 1b and Supplementary Fig 4).
The lower section correlates with the Gilbert chron. We identified two age control points for this part of the sequence, the upper part of Nunivak (C3n.2n 74 , 4.5 Ma) and the lower part of Cochiti (C3n.1n 74 , 4.3 Ma). This correlation is supported by independent age control provided by a 238 U− 206 Pb laser ablation age of 5.0 ± 1.2 Ma.
The middle section correlated with the Gauss and Matuyama chrons. We determined three age control points to build the age model: the upper part of the chron C2An.3n 74 (3.3 Ma), and the normal polarity interval of C2An.1n 74 (3.0 and 2.6 Ma, Gauss−Matuyama transition). Independent age control for this part of the section was provided by a U/Pb laser ablation age of 2.6 ± 0.9 Ma on pedogenic carbonate and by two U/Pb laser ablation ages of 2.61 ± 0.2 Ma and 2.7 ± 0.2 Ma from the tooth enamel, which constrains the normal polarity interval C2An.1n.
The upper section correlates with the normal polarity Olduvai chron (C2n 74 , 2−1.78 Ma), extending upward to the Matuyama−Bruhnes boundary (C1n 74 , 0.75 Ma). This interpretation is supported by two 238 U− 206 Pb laser ablation ages (0.78 ± 0.78 Ma and 0.9 ± 0.7 Ma). The uppermost exposed surface of this section is characterized by a horizontal surface strewn with gravels and may represent an unconformity. Moreover, we observed a typical alternating loess-palaeosol succession approximately 1 km to the southeast which likely stratigraphically overlies our upper section, but could not be included in our analyses due to irregular topography which complicated accurate surveying of key horizons. Both the gravel surface and apparently overlying loess deposits suggest that the Charyn Canyon sequence analyzed in this study extends to ca. 0.5 Ma.
Assuming a linear interpolation of our age model, we calculated sediment accumulation rates ranging between 10 and 20 mm ka −1 for water-lain deposits, with a mean of 13 mm ka −1 (−51 to −39 m of depth, Fig. S5), and between 15 and 40 mm ka −1 (−33 to −10 m of depth, Fig. S5) for interbedded alluvial−aeolian sediments, yielding a mean of 20 mm ka −1 . The upper aeolian deposits yield sedimentation rates of around 36.3 mm ka −1 for the depth interval −10 to +22 m and between 10 and 20 mm ka −1 for the depth interval +22 to +34 m ( Supplementary Fig. 5). The sedimentation rates are consistent within the three different depositional environments at Charyn Canyon. We infer this result to imply that variations in accumulation rates are likely to have been triggered by climate change over late Cenozoic timescales. 87 Sr/ 86 Sr isotope ratio, calcium and strontium concentrations. Twelve representative samples were selected for strontium isotope analysis and the measurement of elemental concentrations of calcium and strontium. The sample suite included two soil matrix and ten pedogenic carbonate samples (Supplementary Table 2, Supplementary Fig. 5). Of the latter, three samples derived from fluvial deposits, one from alluvial sediment, three from caliche, and three from carbonate nodules. Sequential leaching, modified from ref. 75 was performed in the clean laboratory of the Max Planck Institute of Chemistry, Mainz (MPIC), to isolate the carbonate (labile) phase from the silicate (residue) phase. Fifty milligrams of powdered nodule and soil matrix were leached with 0.5 M acetic acid (HOAc) in centrifuge tubes. Samples were placed on a wrist action shaker for 2 h. The residual fraction was separated from the carbonate soluble fraction by a 5-mincentrifugation at 6000 rpm and repeated twice. The two purified carbonate phases of each sample were then merged, evaporated to dryness, and taken up in 6 M HCl. The residues were rinsed with distilled water and digested in a 27 M HF-16 M HNO 3 mixture overnight at 90°C. After evaporation, residues were taken up in 1 mL 1 M HCl. After final evaporation, both carbonate and residual solutions were dissolved in 3 M HNO 3 . Twenty-four aliquots (12 from carbonates and 12 from residual phases) were analyzed to determine the radiogenic 87 Sr/ 86 Sr ratio. Strontium was separated from the other elements using a column in shrink Teflon filled with 100 µm of Sr-Spec resin (50−100 mesh). The isotope measurements were performed using thermal ionization mass spectrometry ThermoFisher Triton at the Max Planck Institute of Chemistry, Mainz. Measurement of the standard NBS 987 was performed regularly to control the accuracy of the machine. The NBS 987 yields an 87 Sr/ 86 Sr of 0.710257 ± 0.000005 (2σ), in agreement with the certified value. The calcium and strontium concentrations were measured in carbonate and residual solutions. Solutions were analyzed using inductively coupled plasma mass spectrometry (Agilent 7900 ICP-MS) at the Max Planck Institute of Chemistry, Mainz. Analytical uncertainties were below 5% based on the repeated standard analyses of certified BCR-2 and BHVO-2 rock standards 76 .
Stable isotope analysis of pedogenic carbonates. Stable isotope analysis was undertaken on 201 pedogenic carbonates, collected with a sampling resolution ranging between 20 and 80 cm (Supplementary Data 6, Supplementary Fig. 5). Pedogenic nodule samples were cut transversally with a circular saw and three replicates per sample were analyzed. Samples were collected using a handheld dentist drill with a diamond-coated drill bit, which was carefully cleaned between samples. Between 0.3 and 0.5 mg of carbonate powder was digested with phosphoric acid at 70°C in an exetainer vial and flushed with helium gas. The liberated CO 2 was analyzed using a Thermo Finnigan GasBench II preparation device interfaced with a Thermo Delta V advantage mass spectrometer at the Max Plank Institute of Chemistry in Mainz. All δ 13 C and δ 18 O values are reported relative to the Vienna-Peedee Belemnite (V-PDB) standard. The long-term standard deviation of a routinely analyzed in-house CaCO 3 standard was <0.1‰ (1σ) for both carbon and oxygen isotope ratios. The CaCO 3 standards were calibrated using NBS19 and IAEA 603. Internal variability was tested by the analysis of twelve samples along with a single pedogenic carbonate nodule (CHA17-CC14400) (Supplementary Table 5); isotopic composition intra-variability is substantially smaller than the stable isotope variability records down the profile.
Bovid tooth stable isotope record. An upper M1 or M2 molar of the bovid (Supplementary Fig. 6 and Note 3) was analyzed for stable carbon and oxygen isotopes in a cross-sectional transect perpendicular to growth. Each sample was split into two specimens. The first set of specimens was analyzed without pretreatment and the second set was leached in order to remove soluble carbonates. Samples were leached by soaking in 1 M of ammonium acetate-buffered acetic acid (Hac NH4Ac) and intermittently mixed in a vortex shaker for 4 h, then washed 3 times with milliQ water and dried in a 40°C oven. The isotopic composition of the non-leached and leached samples from the bovid tooth are presented in Supplementary Table 4. The measured δ 13 C and δ 18 O values were analyzed using the same stable isotope procedure as the pedogenic carbonates, in the same laboratory. The isotopic composition of the leached and non-leached samples yields the same general trends, suggesting a low diagenesis effect. The stable carbon isotopic composition indicates a preferred diet of C3 plants (−23 to −24‰), which likely grew in an open and dry environment. We observed minimal variability in oxygen isotopic composition across the transect (around 2‰), suggesting limited seasonal variability and annual temperature change over the lifespan of the animal. Conditions at ca. 2.6 Ma appears to have been semiarid, relatively stable, and less continental than at present.
Clumped isotope analyses. Carbonate clumped isotope (Δ 47 ) analyses were performed at Goethe University Frankfurt using a fully automated gas extraction and purification line connected to a Thermo Scientific 253 Plus gas source isotope ratio mass spectrometer, following published protocols 77 . We selected three pedogenic carbonate samples representing maximum, minimum, and mean δ 18 O values respectively. Five replicates from each sample were measured. Between 18 and 20 mg of homogenized carbonate powder per sample was digested in >106% phosphoric acid at 90 ± 0.1°C in a common acid bath. The resulting CO 2 was first purified and then analyzed. CO 2 gases equilibrated at 1000°C and 25°C, respectively, were analyzed along with the samples to test the linearity of the mass spectrometer and to determine the empirical transfer function (ETF). Each sample gas was analyzed in 13 acquisitions consisting of ten cycles with an ion integration time of 20 s each. Background with no gas in the ion source was not determined. The negative background below m/z 47 was continuously monitored on the m/z 47.5 half-cup, and subtracted from the corresponding m/z 47 raw-intensities of sample and reference gas measurements using a scaling factor of −0.988. In δ 47 vs Δ 47 space, background-corrected at 1000°C and 25°C gases displayed-within errors-no further residual slopes anymore (Supplementary Data 7). All data were processed using IUPAC parameters 78 . Intercepts displayed by backgroundcorrected equilibrated gases in δ 47 vs Δ 47 space and Wang/Dennis values 78 were used to constrain the ETF. An acid fractionation factor of 0.088‰ 78 was finally applied to report Δ 47 on the CDES 25°C scale. The resulting Δ 47 values were then converted to clumped isotope temperatures using the regression equation Δ47 − RFAC(Br,WD,newAFF) of ref. 78 (1): Δ47 À RFACðBr; WD; newAFFÞ ¼ ð0:0387 ± 1:7E À6 Þ Ã ð10 6 =T 2 Þ þ ð0:257 ± 1:7E À5 Þ ð1Þ The measured δ 47 and Δ 47 values of CO 2 equilibrium gases (HG: CO 2 gas equilibrated at 1000°C and 25G: CO 2 gas equilibrated at 25°C) and ETF applied during this study were reported in Supplementary Data 7. Each sample was analyzed in five replicates, corresponding to a shot-noise limit of 0.003‰. Whenever the external reproducibility (1 sdev) was better than 0.003‰, the shotnoise limit was adopted instead of the external reproducibility to discuss temperature uncertainties. The measured isotopic compositions of the three carbonate samples δ 47 and Δ 47 values and the calculated temperature (T Δ47 ) were reported in Supplementary Table 3. δ 18 O mw was calculated using the oxygen isotope equilibrium fractionation equation 35,79 using the pairs of T Δ47 and δ 18 O pc values (Supplementary Table 3).

Data availability
All data needed to evaluate the conclusions in the paper are present in the paper and/or the Supplementary Material. And the Supplementary Data are also available here: https:// doi.org/10.17632/v9s3bhn27k.1.