Research ArticleCLIMATOLOGY

Evolution of vegetation and climate variability on the Tibetan Plateau over the past 1.74 million years

See allHide authors and affiliations

Science Advances  06 May 2020:
Vol. 6, no. 19, eaay6193
DOI: 10.1126/sciadv.aay6193


The Tibetan Plateau exerts a major influence on Asian climate, but its long-term environmental history remains largely unknown. We present a detailed record of vegetation and climate changes over the past 1.74 million years in a lake sediment core from the Zoige Basin, eastern Tibetan Plateau. Results show three intervals with different orbital- and millennial-scale features superimposed on a stepwise long-term cooling trend. The interval of 1.74–1.54 million years ago is characterized by an insolation-dominated mode with strong ~20,000-year cyclicity and quasi-absent millennial-scale signal. The interval of 1.54–0.62 million years ago represents a transitional insolation-ice mode marked by ~20,000- and ~40,000-year cycles, with superimposed millennial-scale oscillations. The past 620,000 years are characterized by an ice-driven mode with 100,000-year cyclicity and less frequent millennial-scale variability. A pronounced transition occurred 620,000 years ago, as glacial cycles intensified. These new findings reveal how the interaction of low-latitude insolation and high-latitude ice-volume forcing shaped the evolution of the Tibetan Plateau climate.


The Tibetan Plateau has long been a focus of geoscientific studies because of its importance in global tectonics and Asian and global climate change across a wide range of time scales (1). However, with only few available paleoarchives of coarse resolution [>8 thousand years (ka)] (2, 3), little is known about its environmental history through the Quaternary ice ages. To understand the mode and tempo of changes and, ultimately, the underlying drivers during this period, we need long-term high-resolution records from the elevated plateau with well-constrained chronologies.

The Zoige Basin, occupied by a huge lake until the latest Pleistocene (3) and located on the eastern Tibetan Plateau within the South Asian monsoon zone (fig. S1), represents a potential site to fill this gap. Mean annual precipitation (MAP) at Zoige is ~600 to 650 mm, and the basin is primarily covered by alpine meadows, and the surrounding mountains have scattered forests up to ~4000 m above sea level (a.s.l.) (fig. S1). A sediment core extending over the past 0.9 million years (Ma) was previously recovered, but its analytical resolution and chronological reliability were insufficient to resolve orbital- and suborbital-scale changes (3).


New drilling (33°58.163′N, 102°19.855′E, 3434 m a.s.l.) was undertaken in 2013 in the central basin guided by a seismic survey. A 573.39-m core (ZB13-C2) was obtained with 96% recovery, mostly consisting of fine-grained freshwater lacustrine sediments. Only the upper 50 m contains two episodic fluvial sandy layers, 10.11 and 10.4 m thick, respectively.

Independent age control derived from magnetostratigraphy in combination with radiocarbon [accelerator mass spectrometry (AMS)] and luminescence [optically stimulated luminescence (OSL)] dating provides an initial chronological framework (Materials and Methods; figs. S2 and S3), according to which the ZB13-C2 core extends back to 1.74 Ma before present (BP) (Materials and Methods; fig. S3). Fluctuations in arboreal pollen abundances (AP%) based on an initial age model using a combination of 14C, OSL, and paleomagnetic control points (table S3) show clear ~100-, 40-, and 20-ka cyclicities (fig. S3), suggesting possible eccentricity (E), tilt (T), and precession (P) powers. The presence of astronomical frequencies in the Zoige Basin record is further supported by spectral analyses of AP% in the depth domain, which indicate the occurrence of ~34-, 15-, and 7.5-m cycles, whose ratios are close to those of 100:40:21; the ~7.5-m cycle appears to be stronger in the lower ~75 m, while the ~34-m cycle is stronger in the top ~200 m (Materials and Methods and fig. S3). On this basis, a more detailed age model was constructed by tuning the AP% record to an ETP record that is generated by normalizing and averaging variations in eccentricity, tilt, and reversed precession (4). As this approach may artificially introduce astronomical frequencies in our record, we compared the ETP age model against an age model constructed by aligning the Zoige AP% to the Chinese speleothem δ18Ocalcite record, an independently dated (U-Th) archive of changes in Asian monsoon intensity over the past 640 ka (5). Comparison of the ETP- and speleothem-based age models of Zoige reveals a close correspondence (fig. S3).


Pollen and sediment analyses (2787 and 3274 samples, respectively, with a mean sampling resolution of ~530 to 620 years), as well as x-ray fluorescence (XRF) scanning (~6-year resolution) data, yield a detailed multiproxy record back to 1.74 Ma BP.

The vegetation in the eastern Tibetan Plateau is strongly influenced by the Asian summer monsoon. A stronger monsoon with warmer and moister climate would cause an expansion of tree populations. However, in the alpine Zoige region, the density and elevational limits of forests are primarily controlled by temperature, as moisture availability is relatively plentiful (6, 7). This is confirmed by the distinct elevational distribution of modern vegetation (fig. S1), the good match of AP% and axis 1 of principal components analysis (PCA) on the pollen data (fig. S4), and the close relationship between summer temperature and PCA axis 1 in the significance test of the quantitative temperature reconstructions for core ZB13-C2 (fig. S4). Thus, on long time scales, variations in AP% in the core are mostly a reflection of changes in temperature, particularly summer temperature, but drought stress, for example, during weak monsoon or glacial intervals, would also have an impact on tree populations.

Pollen-based quantitative reconstruction of past climate variables was undertaken (Materials and Methods), providing the first independent paleotemperature history for the Tibetan Plateau because of the dominant control of temperature on local vegetation (figs. S1 and S4) (6, 7). The mean temperature of the coldest month (MTCM) and precipitation reconstructions failed the relevant statistical significance test and are thus less reliable (Materials and Methods).

Rubidium/strontium (Rb/Sr) ratios and carbonate content (Carb%) are used as supplementary proxies. Rb/Sr primarily reflects the chemical weathering intensity of the catchment or strength of summer monsoon and associated run-off (8). When weathering/run-off is stronger, there is greater Sr input into the lake, leading to a lower Rb/Sr ratio. However, Sr contained in carbonate can influence Rb/Sr ratios from bulk scanning data, rendering climate interpretations less reliable. We therefore measured Rb/Sr ratios on bulk samples after removing the carbonate content (Materials and Methods). This showed good agreement with the high-resolution XRF scanning data (fig. S5), supporting the view that the ZB13-C2 Rb/Sr signal derived from XRF scanning is largely independent of changes in carbonate content. Grain size changes could also distort the climate signal of Rb/Sr, but examination of core ZB13-C2 shows a weak correlation between Rb/Sr and grain size changes (fig. S5). The Zoige lake sediments are generally fine and have no large variations, except the two sandy layers near the top. The high correlation between Rb/Sr and the chemical index of alteration (CIA) at Zoige further support the weathering interpretation (fig. S5) that the proxies are sensitive to both summer precipitation and temperature conditions.

Previous studies from core RH nearby ZB13-C2 suggest that carbonate content mainly represents authigenic chemical precipitation, as both detrital carbonate content and shell carbonate are in trace amounts (3). The measured Carb% in core ZB-C2 should mainly reflect processes of chemical precipitation in the lake, which are largely related to summer temperature and precipitation. High temperature could increase the precipitation of carbonate through changing the precipitation-dissolution equilibrium and photosynthesis process. Precipitation could also enhance carbonate content by washing more Ca2+ and HCO3 into the lake through chemical weathering. Carb% therefore indicates warm and wet climate. Carbonate content from core RH shows good positive correspondence with hydrogen index, a proxy of the effect of lake water depth (3). Our loss-on-ignition (LOI) measurements of surface mud samples, taken along water-depth transects from four lakes near ZB13-C2, also reveal that Carb% generally increases with water depth in each lake (fig. S5). High Carb% from core ZB13-C2 likely agrees with high lake level, which depends on the balance between precipitation and evaporation.

The coherent variations in AP%, Rb/Sr, and Carb% (fig. S5), with low Rb/Sr and high Carb% corresponding to high AP%, therefore suggest that coupled changes of temperature and precipitation occurred over the past 1.74 Ma.


Long-term vegetation and climate changes

The Zoige Basin record reveals three major intervals with transitions at ~1.54 and ~0.62 Ma BP (Fig. 1). The interval of 1.74–1.54 Ma BP is marked by conifer and deciduous mixed forest (mostly Picea, along with Pinus, Betula, and Quercus) for interglacial times and steppe (dominated by Artemisia) for glacial times (Fig. 2). Interglacial biomes shift to conifer forest, and glacial biomes shift to meadow (dominated by Cyperaceae/steppe) in the 1.54–0.62 Ma BP interval with increasing steppe (mainly Poaceae and Ranunculaceae) at ~1.03 Ma BP. Over the past 0.62 Ma, forest expanded less frequently. A less prominent change is observed at ~0.41 Ma BP, after which Taraxacum-type and Poaceae percentages declined (Fig. 2).

Fig. 1 Comparison of multiproxy records from the Zoige Basin core ZB13-C2 with regional insolation of the past 1.74 Ma.

(A) Mean June insolation at 30°N (red) (4). (B) Chinese speleothem δ18Ocalcite (gray) (5). (C) Arboreal pollen percentages AP% (~600-year resolution) from ZB13-C2 (this study), a drill core from the Tibetan Plateau consisting of lacustrine sediments. ZB, Zoige Basin. (D) Reconstructed mean temperature of the warmest month (MTWM) and (E) mean annual temperature (MAT) based on pollen data from ZB13-C2. The horizontal lines indicate the mean values of MTWM and MAT for the intervals of 0–0.62, 0.62–1.03, 1.03–1.54, and 1.54–1.74 Ma BP. Bootstrap sample-specific estimates of uncertainties are shown in fig. S4. (F) Sixty-year smoothed Rb/Sr ratio based on XRF scanning from ZB13-C2. (G) Carbonate percentage (~600-year resolution) determined by loss on ignition from ZB13-C2. The gray bars mark the two fluvial sandy layers.

Fig. 2 Pollen percentage diagram of major taxa and biome reconstruction from the Zoige Basin core ZB13-C2.

(A) Pollen percentage diagram. Tree taxa mainly include Picea, Pinus, Betula, and deciduous Quercus. Major meadow/steppe taxa consist of Cyperaceae, Artemisia, other Asteraceae, Poaceae, and Ranunculaceae. The consistent presence of aquatic pollen of Myriophyllum suggests that the study site was occupied by a shallow-intermediate lake over the past 1.74 Ma. The pollen zonation is based on CONISS results aided by a multivariate regression tree analysis, which shows the biggest changes at 1.54 and 0.62 Ma BP. (B) Combined biome. The biomes were combined into five types. The results indicate a less frequent occurrence of forest from 0.62 Ma BP onward.

Mean annual temperature (MAT) and mean temperature of the warmest month (MTWM) pollen-based reconstructions (Fig. 1) also show a stepwise cooling trend with the major transitions occurring at ~1.54, ~1.03, and ~0.62 Ma BP. Because winters at Zoige are extremely cold and dry with seasonal temperatures below 0°C, these signals would mainly reflect summer conditions and are, hence, mostly related to the summer monsoon circulation. The similarity between MAT and MTWM supports this interpretation.

Rb/Sr shows a generally increasing trend with shifts at ~1.54 and 0.62 Ma BP (Fig. 1), suggesting the reduced precipitation associated with weakening summer monsoon. Carb% also indicates a reduced summer temperature and a decreasing lake-level trend. Although gradual filling of the basin by sediment may have contributed to the long-term changes in the lake environment, the similarity between the Carb%, Rb/Sr, and AP% records implies the presence of a long-term cooling and drying trend.

Orbital cyclicity

On orbital time scales, AP%, Rb/Sr, and Carb% all reveal ~100-, ~40-, and ~20-ka cycles with varying intensities through time (Fig. 3 and fig. S6), pointing to a link with astronomical frequencies (4). The 1.74–1.54 Ma BP interval is marked by an exceptionally strong ~20-ka cycle and a relatively weaker ~40-ka cycle. Weakening but still strong, ~20-ka along with ~40-ka signals are observed for the 1.54–0.62 Ma BP interval. The ~100-ka cycle appears at ~1.03 Ma BP and becomes established in the past ~0.62 Ma, along with an intermittent ~20-ka cycle and a weak 40-ka cycle. Sensitivity analyses, using the ETP age model and the age model based only on 14C and OSL dates, and magnetostratigraphy produce similar evolution in periodicities (fig. S3), underlining the robustness of these features at Zoige.

Fig. 3 Continuous wavelet and spectral results of arboreal pollen (AP%) from the Zoige Basin core ZB13-C2.

(A) The continuous wavelet and (B) the multitaper method spectral results for the intervals of 0–1.74, 0–0.62, 0.62–1.54, and 1.54–1.74 Ma BP. The data were resampled at equally spaced 0.6-ka intervals and detrended before analysis. AP% generally shows powers of ~40, 20, 10.6 to 9.6, 8, 6.6, 5.2, and 4.4 to 4.1 ka. Black lines denote the 90% confidence level. The ~100-ka cycle becomes established in the past ~0.62 Ma, along with an intermittent ~20-ka cycle and a weak ~40-ka cycle.

Low-latitude summer insolation, dominated by ~20-ka periodicity, is commonly considered a main forcing of the Asian summer monsoon. High summer insolation enhances the summer monsoon through stronger continent-ocean thermal gradient that increases the atmospheric humidity and wind circulation intensity, or latitudinal oscillations of the intertropical convergence zone (9) via the associated meridional flow. A stronger summer monsoon brings more heat and moisture to the Tibetan Plateau, leading to warmer and moister summers. The ~20-ka cyclicity throughout the Zoige Basin record implies the strong influence of summer insolation via the summer monsoon on the climate and vegetation of the Tibetan Plateau.

Because Early Pleistocene ice-volume changes were dominated by the ~40-ka cycle (10), the particularly strong and continuous ~20-ka cycle in the 1.74–1.54 Ma BP interval indicates the dominant influence of insolation versus ice-volume forcing. During the 1.54–0.62 Ma BP interval, the ~20-ka insolation cycle was still pervasive, but with the 40-ka glacial cycle becoming stronger. The first appearance (~1.03 Ma BP) and establishment (~0.62 Ma BP) of the ~100-ka cycle at Zoige are consistent with the rise of the ~100-ka glacial cycle (11).

Studies have also highlighted the interaction of the westerly jet and the Asian monsoon system (12). Seasonally, the role of the latitude of the westerly jet relative to the Tibetan Plateau is critical in determining the stepwise transitions in East Asian rainfall seasons. At the glacial-interglacial scale, climate model studies suggest that changes to the westerly jet stream intensity and its position relative to the Tibetan Plateau might be closely coupled with the large-scale Asian monsoon changes (12). The presence of large ice sheets could lead to the strengthening of the westerly jet and its displacement south of Tibet, which, in turn, could prevent the penetration of monsoonal precipitation into East Asia, especially over the past ~0.62 Ma.

The evolution in orbital spectral pattern at Zoige differs from that of other long-term records in Asia, including the loess-soil sequence in northern China (1315) and the lower-latitude and low-elevation Heqing lacustrine record (8). Most of the loess records show only one major transition at ~0.9 to 0.7 Ma BP from ~40-ka (13, 14) to ~100-ka cycles and less prominent ~20-ka periodicity, perhaps a reflection that the loess records contain climate signals related to both Asian summer and winter monsoon circulations, while the Zoige climate signals are mainly related to the summer monsoon circulation. The pollen record from Heqing Basin shows a shift in the amplitude of Tsuga (hemlock) abundances after ~0.9 Ma BP, with larger interglacial expansions and major glacial contractions (8). Spectral analyses of the Tsuga record show strong 40-ka and 20-ka cycles during the interval of 1.82–0.9 Ma BP, with a rise in the 100-ka power after that. By comparison, the Heqing Indian summer monsoon (ISM) index (derived from combining total organic content and Rb/Sr records) shows a shift toward reduced amplitude changes after ~0.9 Ma BP, while spectral analyses show a shift from 40- and 20-ka cyclicity to 100- and 20-ka cyclicity at 0.9 Ma BP. Unlike the Tsuga record, the ISM index does not have a 40-ka spectral peak during the interval of 0.9–0.13 Ma BP but contains peaks at 73-, 55-, and 29-ka periods. An et al. (8) attributed the ISM changes to the interplay of the southern Mascarene high (influenced by Antarctic temperature) and the northern Indian low (influenced by ice volume). Zoige Basin located to the north of Heqing is less subject to cross-equatorial transfers of heat, which may account for the differences between the two records. The clear ~100-ka cycle in the past 0.62 Ma at Zoige also differs from the Chinese composite speleothem δ18O record (5) of the Asian Monsoon changes over the past 0.64 Ma; the cave records are obtained from South China and have an almost pure ~20-ka cyclicity. In summary, these divergences may arise from differences in the dominant mode of climate variability at the geographical locations of these records.

Millennial-scale variability

Superimposed on the long-term and orbital-scale changes are millennial-scale oscillations with frequencies centered at ~10 ka (10.6, 10.2, 9.6, and 8 ka) and ~5.5 ka (6.6, 5.2, and 4.4 ka) (Fig. 3). These signals are weak in the 1.74–1.54 Ma BP interval, but increase in frequency and amplitude in the 1.54–0.62 Ma BP interval, and decrease in frequency in the past 0.62 Ma. The amplitude of millennial-scale changes in the pollen record over the past ~0.62 Ma is subdued, but changes in the Rb/Sr and Carb% values (Fig. 1 and figs. S6 and S7) indicate that millennial-scale variability has been a permanent feature of the Tibetan Plateau climate since 1.54 Ma BP. Such temporal evolution of millennial features is not obvious in other long Asian records (8, 13, 14), perhaps as a result of geographical location and/or lower analytical resolution (~1.5 to 2 ka).

The origin of suborbital climate variability in paleoclimate archives remains controversial and is linked either to millennial-scale freshwater discharges and disruptions of the Atlantic meridional overturning circulation (AMOC) (11, 16) or to tropical insolation changes at half- (~10-ka) and quarter-precession (~5.5-ka) bands (17, 18). Although the Zoige Basin millennial signals exhibit cyclicities around these bands, their quasi-absence in the 1.74–1.54 Ma BP interval and the frequency/amplitude changes in the past 1.54 Ma do not support the link with tropical insolation, because insolation values have not substantially changed around 1.54 and 0.62 Ma BP. Half- and quarter-precession cyclicities are also detectable in Atlantic ice-rafting proxies (11), further suggesting that linking these signals in the geological record with tropical forcing is not always straightforward. A North Atlantic origin of millennial signals at Zoige is also supported by the strong similarity between the Rb/Sr and the North Greenland Ice Core Project (NGRIP) ice δ18O (19) and Sanbao Cave δ18O records (5) during the last glacial (fig. S8).


The Zoige Basin temperature reconstruction shows a long-term cooling trend over the past 1.74 Ma, starting with warm glacial and interglacial conditions and followed by cooling steps at ~1.54, ~1.03, and ~0.62 Ma BP. The overall cooling is in line with sea surface temperature (SST) reconstructions from the North Atlantic, equatorial Pacific, and Arabian Sea (Fig. 4 and fig. S9) (20, 21). Some of the SST records show declines after ~1.6 to 1.5 Ma BP, most show a prominent cooling during marine isotope stages (MIS) 24–22 (866 to 936 ka BP) (22), and none indicate a particularly distinctive change during MIS 16 (~676 to 621 ka BP). In contrast to most SST reconstructions, the Zoige shows much higher glacial temperatures before 1.54 Ma BP, perhaps a reflection of smaller ice sheets during that interval. In terms of ice volume changes, marine records (11, 23) indicate (i) a transition toward larger ice sheets and increased ice rafting from ~1.54 Ma BP; (ii) a prominent glaciation ~0.9 Ma BP; (iii) and the establishment of the 100-ka cycle, with growth of large ice sheets and the initiation of large iceberg discharges through the Hudson Strait from 0.64 Ma BP, in broad agreement with the long-term changes observed at Zoige Basin.

Fig. 4 Correlation of vegetation and temperature records from the Zoige Basin core ZB13-C2 with global and regional climate records over the past 1.74 Ma.

(A) LR04 global benthic δ18O stack (blue) (10). The numbers denote the MIS. (B) Seawater δ18O (23) at ODP 1123 (gray), east of New Zealand, which is a direct indicator of ice volume. (C) CO2 records are shown as follows: orange, δ11B-based data from the Caribbean (27) and eastern tropical Atlantic (31); blue squares, low-resolution δ11B record from the equatorial Atlantic (all with 2σ error bars/envelopes) (29); purple circle, ice core CO2 measurements from blue ice (30); purple line, ice core compilation (36). (D) Southern Ocean dust mass accumulation rate (MAR) (33). (E) Bulk carbonate δ18O at site U1308, north Atlantic (11). The dashed lines indicate the mean values for three intervals. (F) Estimates of SST at site DSDP 607 (20), north Atlantic, and tropical stack (21). (G) Reconstructed MAT and 9-point running mean (black line) from ZB13-C2 (this study). The dashed lines indicate the mean values of MAT for three intervals. (H) Reconstructed MTWM and 9-point running mean (black line) from ZB13-C2 (this study). The dashed lines indicate the mean values of MTWM for three intervals. (I) Arboreal pollen percentages (~600-year resolution) from ZB13-C2. Black line shows the 9-point running mean.

A secular decline in atmospheric CO2 concentrations may have contributed to the observed long-term global cooling (24). This would have gradually raised the solar energy threshold required for a complete deglaciation starting ~1.5 Ma BP, leading to skipped insolation peaks after 1 Ma BP and longer glacials with larger ice sheets (25), although progressive regolith erosion and attendant changes in ice sheet dynamics may have also contributed to this change (2628). Boron-based CO2 reconstructions (Fig. 4) do not support a long-term decrease in pCO2, but a more abrupt change around 0.9 Ma BP, from consistently greater than 200 to ~180 parts per million (27, 2931), although existing data are generally of low resolution. The 0.9 Ma BP decline in atmospheric CO2 may be related to a weakening of AMOC during MIS 24–22 and increased deep-ocean carbon storage (32). On the other hand, a gradual increase in dust and iron flux from ~1.5 Ma BP in the subantarctic Southern Ocean (Fig. 4) suggests increased iron fertilization and possible CO2 drawdown during over several glacial intervals (33). In addition, vertical carbon isotope gradients between the intermediate and deep Atlantic increased ~1.5 Ma BP, implying greater carbon storage in the deep ocean during glacials and attendant lowering of atmospheric CO2 (11, 34). After 0.9 Ma BP, glacial atmospheric CO2 concentrations reached relatively similar minima (29, 35, 36).

With respect to millennial-scale variability, comparison of the Zoige Basin record with North Atlantic ice rafting history reveals parallel trajectories (Fig. 4 and fig. S7). Before ~1.54 Ma BP, Atlantic millennial ice rafting only occurred episodically with much limited intensity and spatial extent due to the smaller ice sheets (11). This is in line with the quasi-absence of millennial signals at Zoige during the 1.74–1.54 Ma BP interval. The 1.54 Ma BP boundary marks the emergence of millennial variability at Zoige with larger amplitudes and higher frequencies. It is correlative with a climate transition in the North Atlantic characterized by more active and permanently established millennial variability associated with more extensive ice sheets (11, 15). Another transition occurred during MIS 16 in the North Atlantic when Heinrich events were initiated (37). Since then, the intensity of millennial variability increased, but its frequency decreased (37), which can account for the less frequent millennial-scale changes at Zoige in the past 0.62 Ma BP.

Last, it is worth noting that a major shift in the Zoige pollen record occurred at ~0.62 Ma BP, near the end of MIS 16, one of the largest Pleistocene glaciations, which included the most extensive ice advance in the southern Russian Plain (38). In addition to the change in the frequency of millennial-scale changes, the pollen record shows a shift to lower arboreal values during both glacials and interglacials, compared with pre–0.62 Ma BP levels. The glacial declines are in line with our understanding of increased glacial ice volumes from MIS 16 onward, but the reduced interglacial tree population expansions (mainly in Picea) in the Zoige Basin do not fit global trends (39). Uplift of the Tibetan Plateau at ~0.62 Ma BP could provide an alternative explanation for both the glacial and interglacial vegetation shifts by permanently reducing temperatures (40). However, the eastern part of the Tibetan Plateau is dominated by strike-slip faulting with small rates of vertical motion (41), and plate convergence considerations and paleoaltimetry reconstructions indicate decreased and relatively constant uplift rates in the Quaternary (4244). It is possible that the emergence of longer and more extreme glacials in the past ~0.62 Ma led to tree populations seeking refugial locations at increasingly lower elevations. Longer migrational distances could then account for their subdued expansions, with alpine meadows (as shown by the Cyperaceae curve) becoming more important around Zoige, not only during glacials but also during interglacials. The shift toward lower arboreal pollen percentages at ~0.62 Ma BP is paralleled by changes in carbonate content and Rb/Sr records (Fig. 1), although mainly in glacial values. Together, the data indicate a major environmental transition on the Tibetan Plateau toward cooler and drier summers, as global ice volumes increased.


Changes in ocean circulation, carbon cycle, climate, and the cryosphere during the Early-Middle Pleistocene have led to the prolongation and intensification of glacial cycles. Pollen and sediment records from the Zoige Basin reveal the manifestation of these global changes in the eastern Tibetan Plateau as a stepwise cooling trend over the course of the past 1.74 Ma. More specifically, the evolution of the Tibetan Plateau climate and vegetation can be described by three successive modes with different combinations of low-latitude insolation forcing and high-latitude ice forcing (Fig. 5). The 1.74–1.54 Ma BP interval corresponds to an insolation-dominated mode, while the past 0.62 Ma correspond to an ice-driven mode. The 1.54–0.62 Ma BP interval represents a transitional insolation-ice mode. The relative impact of ice forcing versus insolation forcing increased in a stepwise manner, with a prominent transition occurring at ~0.62 Ma BP. The influence of North Atlantic millennial-scale variability was established from ~1.54 Ma BP. The interaction of these forcings shaped the expression of orbital- and millennial-scale climate variability on the Tibetan Plateau during the course of the Quaternary. While other long terrestrial records from China and marine SST records show a major transition in the mode and tempo of climate variability centered near 0.9 Ma BP, the Zoige Basin records indicate a major shift toward cooler conditions at ~0.62 Ma BP, coeval with the intensification of glacial cycles. Our records advance understanding of the spatial complexity of long-term environmental responses in monsoonal Asia and the pervasive influence of millennial-scale climate variability over the past 1.54 Ma.

Fig. 5 Schematic diagram of the forcings of Zoige Basin vegetation and climate change modes at orbital and millennial scales over the past 1.74 Ma.

The dark colors in the bars represent the trends of the strength of the forcings (low-latitude insolation and ice volume). The triangles demonstrate the frequency and strength of northern Atlantic ice-rated debris (IRD) variabilities. The curves show the Zoige-filtered AP% values centered around 20-, 40-, 100-ka bands, and the millennial variabilities of AP% (green) and Rb/Sr ratio (brown) at the powers between 15 and 3 ka.


Study site description

The Zoige Basin is a low-relief tectonic basin (32°10′-34°10′N, 101°45′-103°25′ E, ca. 3350 to 3450 m a.s.l.) on the eastern Tibetan Plateau. The surrounding mountains exceed 4000 m in elevation. At present, the Yellow River flows into the western part of the basin. However, a huge lake occupied the basin during the Pleistocene, until it was finally drained about 30 to 40 ka ago when the Yellow River cut through the mountain barrier to the east (3). A new AMS 14C date on terrestrial plant macrofossils at 11.28 m (near the bottom of the upper fluvial sandy layer) from core ZB13-C2 constrains the timing of the eventual drying to younger than 28 ka ago at the coring site.

MAT at the nearby Zoige meteorological station at the elevation of 3439 m a.s.l. (based on data from 1981 to 2010) is ca. 1.1°C, with July temperature of ca. 10.8°C and January temperature of ca. −10.2°C. MAP is about 650 mm. Most precipitation falls as rain during the summer months, owing to the influence of the Asian summer monsoon.

The basin is primarily covered by alpine meadows dominated by Kobresia spp. and other Cyperaceae taxa. The surrounding mountains are covered by scattered forests up to ca. 4000 m a.s.l. (fig. S1), mainly composed of Picea asperata, Picea wilsonii, Picea purpurea, Abies faxoniana, Pinus densata, Betula platyphylla, and Quercus liaotungensis and by scrub dominated by Salix and Rhododendron (6, 7). The zone between ca. 4000 and 4400 m is occupied by alpine shrub meadow. Alpine periglacial desert occurs from ca. 4400 to 4600 m. The peaks above ca. 4600 m are covered by ice and snow. Temperature is the dominant climate control of the major vegetation types in this forest-alpine meadow ecotonal region (7).

Seismic survey, drilling operation, and sediment lithology

The seismic survey was conducted in the sedimentary center of the paleolake using Tromino 3G seismometers. Two 1-km transects along the north-south and west-east directions were surveyed, with 40 points along each transect. Another six locations within 3 km of this region were also investigated. The passive seismic data were analyzed using the spectral ratio of horizontal and vertical component data. The results indicate ca. 790 m (with 20% error range) of continuous sediments in this region.

Scientific deep drilling in the sedimentary center of the paleolake of the Zoige Basin was performed in August to October 2013. A 573.39-m core (ZB13-C2: 33°58.163′N, 102°19.855′E, 3434 ± 4 m) was retrieved, with the recovery up to 96%. The depth was calibrated from the initial drilling after considering errors caused by sediment expansion in some sections, loss of samples, and drilling deviation. The polyvinyl chloride (PVC) tubes were cut into two halves, one-half for XRF scanning and laboratory analysis, and the other for archiving.

The lithology of the sediment core mainly consists of clay, silt-clay, and clay-silt, except for two intervals of sand layers (fig. S2). In general, the sediments are dominated by fine bluish gray clay at depths of 573.39 to 490.52 m, and grayish green clays and silty clays with thin-bedded silt at 490.52 to 50.4 m and 40.29 to 11.28 m. Two sand layers occur at 50.4 to 40.29 m and 11.28 to 1.88 m. The top 1.88 m is dominated by brown-black silty clay, with soil at the surface that is disturbed by grazing. Mollusk shells indicating shallow water are found intermittently, except at the depths below 490.52 m and in the two sand layers.

Dating methods

AMS radiocarbon dating. Eleven samples of terrestrial macrofossils and one sample of bulk sediment from the top 20 m were picked for AMS 14C dating, following standard methods. Calibration was calculated using the database associated with the 2013 INTCAL program. The results are shown in table S1.

OSL dating. OSL measurements of 11 levels of samples from the top 80 m were carried out using an automated Risø TL/OSL-DA-15 reader. Burial age was estimated by dividing the equivalent dose (De) by the environmental dose rate.

All raw samples were first treated with 10% HCl and 20% H2O2 to remove carbonate and organic matter. The samples were then sieved in water to pick out the grain size fraction of 63 to 90 μm. This grain fraction was separated by heavy liquids to obtain quartz using densities between 2.62 and 2.75 g/cm3 and K-feldspar using densities lower than 2.58 g/cm3. After drying, the quartz grains were treated with 40% hydrofluoric acid (HF) for 40 min to remove the outer layer irradiated by alpha particles and also to remove any remaining feldspars. The K-feldspar grains were treated with 10% HF for 40 min. All grains were then treated with 1 M HCl for 10 min to remove fluorides created during the HF etching.

Laboratory irradiation was carried out using 90Sr/90Y sources mounted within the reader with a dose rate of 0.084 billion years (Ga)/s. De determination was done adhering to both the post-IR (infrared) IRSL (IR-stimulated luminescence) (pIRIR) protocol (45) and the multiple elevated temperature (MET) stimulation (MET-pIRIR) protocol (46).

K-feldspar is known to have two advantages over quartz for optical dating: The IRSL signal (per unit absorbed dose) is usually much brighter than the OSL signal from quartz, which leads to lower measurement errors, and the IRSL traps saturate at a much higher dose than the OSL does, making it possible to date older samples using K-feldspars. However, the routine dating of K-feldspars using the IRSL signal has been hampered by the phenomenon of “anomalous fading” (47), which gives rise to substantial underestimates of age. Recent progress in understanding anomalous fading of the IRSL signals in K-feldspar has raised the prospect of isolating a nonfading IRSL component for dating Quaternary deposits. A high-temperature (ca. 290°C) pIRIR signal, which is stimulated after a low-temperature (ca. 50°C) IRSL stimulation, exhibits a negligible anomalous fading rate. This has led to the development of pIRIR protocols (45) and also a MET stimulation (MET-pIRIR) protocol (46). The two protocols obtain almost the same results for samples younger than 80 ka, and the MET-pIRIR protocol can be applied to old Quaternary deposit events (>100 ka).

In this study, the quartz OSL traps are expected to be saturated, owing to the old ages of some samples (>100 ka). We applied the pIRIR protocol to the samples younger than 80 ka to save laboratory time and applied the MET-pIRIR protocol to older samples. We also allowed for any residual dose at the time of sediment deposition, considering that pIRIR traps are less easily bleached than the “fast”-component OSL traps in quartz. The resulting MET-pIRIR ages should therefore be reliable estimates of the time of sediment deposition.

To confirm that there is a negligible residual dose for our samples, we bleached five aliquots of sample #IOSL 2014-25 for a total of ~24 hours (~4 hours per day for six sunny days in June). The residual dose was then measured using the same MET-pIRIR procedure. The calculated residual dose for MET-pIRIR signals at 250°C is 17.40 ± 0.40 Ga, which is negligible compared to the De values of 401.11 ± 0.49 Ga in this study. We therefore ignore the residual dose for all of the samples. At the simulation temperatures of 200° and 250°C, the De values form a plateau, which indicates that the fading for the 200° and 250°C signals in the MET-pIRIR method has a negligible effect on De estimates for our samples. Because the MET-pIRIR signals from K-feldspar have a much slower bleaching rate compared with the quartz OSL signal (46), all the K-feldspar De values are consistent with a central value within 2σ on a radial plot, indicating that the measured K-feldspar grains were well bleached before burial.

The radioactive elements of the U, Th, and K contents in sediments were determined by means of neutron activation analysis. All measurements were converted to beta and gamma dose rates according to the conversion factors of Adamiec and Aitken (48). The dose rate from cosmic rays was calculated on the basis of burial depth, latitude, longitude, and altitude of the samples. For the calculations of the internal dose rate of K-feldspar separates, a value of 13 ± 1% was chosen as the K content (49).

All the dating results are shown in table S2. The dates at the depths of 35.5, 41.96, and 61.22 m are associated with minor reversals. However, Bayesian model shows that all the OSL dates, except the date at the depth of 58.06 m, which was dated as 39.25 ka BP and is obviously too young, are in order within their error ranges (fig. S3). The seemingly out-of-order dates would not cause a problem for the precision and accuracy of the age model for such a long record, although with some acceptable uncertainties as shown in the error range. There are two OSL dates (108.7 and 136.1 ka BP) to constrain the position of the Last Interglacial.

Paleomagnetic measurements. A total of 3018 samples of 2-cm cubes were used for paleomagnetic measurements. Anisotropy of magnetic susceptibility (AMS) measurements was first performed using a KLY-4s Kappabridge (Agico Ltd., Brno). The susceptibility tensor for each sample was calculated from measurements in 15 positions using Jelínek’s (50) method. Most minimum susceptibility axes (Kmin) of the AMS ellipsoid are close to the vertical, perpendicular to the bedding plane, whereas most maximum axes (Kmax) are close to the horizontal, parallel to the bedding plane (fig. S2). The AMS results are typical for a primary sedimentary magnetic fabric, indicating that the Zoige sedimentary sequence has not been disturbed since deposition (51) and suitable for paleomagnetic analysis.

Second, high-temperature magnetic susceptibilities were measured to identify the mineral type. The results show a major decrease in magnetic susceptibility at about 580°C, the Curie point of magnetite, suggesting that magnetite is the major contributor to the susceptibility (52). Hysteresis parameters, such as saturation magnetization (Ms), saturation remanence (Mrs), coercivity field (Bc), and the coercivity of remanence (Bcr), were also determined using a MicroMag 2900 Alternating Gradient Magnetometer to decide on the demagnetization approach. The results show that majority of the samples obtained >95% of saturation remanence at ca. 300 mT, and the coercivity of remanence is generally low (ranging from 39.78 to 65.52 mT), suggesting that magnetite with low coercivity is the major component. The magnetic hysteresis loops are closed in the shape of the coarse waist type, also indicating that the mineral with low coercive force controls hysteresis behavior (fig. S2). The hysteresis ratios (Mrs/Ms versus Bcr/Bc) show that the average magnetic grain size falls in the single-domain range. All these features reveal that alternating field (AF) demagnetization approach is more suitable for the Zoige samples.

Remanence measurements were then made using a 2G-760 cryogenic superconducting magnetometer in the magnetic shielded space (<300 nT). All 3018 samples were subjected to stepwise AF demagnetization up to 70 mT at 5-mT intervals. A total of 2350 samples gave reliable characteristic remanence (ChRM) directions. For the 668 samples whose AF results were not satisfactory, thermal demagnetization was further measured using a TD-4 thermal demagnetizer. They were stepwise heated to 620°C at 20° to 50°C increments. In general, both AF and thermal demagnetizations can remove viscous components of magnetization after 15 to 20 mT or 300°C, respectively. A total of 86% samples give stable and reliable ChRM directions (fig. S2).

The magnetism stratigraphy results were then correlated with the geomagnetic polarity time scale (table S3 and fig. S2). A clear geomagnetic reversal (the Brunhes/Matuyama boundary) was identified at 272.26 m. The Jaramillo (343.28 to 364.52 m) and Cobb Mountain (390.43 to 398.01 m) subchrons were identified in the Matuyama chron.

In addition, there are four possible excursions at the depths of 72.72 to 73.90, 176.19 to 176.49, 378.20 to 378.90, and 521.35 to 522.17 m, probably corresponding to the Blake, Emperor, Panaruu, and Gilsa events, respectively. These excursions could provide helpful information for our age model establishment discussed below, although they are not used as control points.

Age-depth model establishment

First-order tie points are provided by the AMS 14C, OSL, and magnetostratigraphic results (table S3), which were used in constructing an initial age model. The AP% record was then tuned using linear interpolation to the ETP record that is generated by normalizing and averaging variations in orbital eccentricity, tilt, and reversed precession from Laskar et al. (4).

Three pieces of evidence support the robustness of ETP tuning chronology. First, wavelet analysis on the AP% using the untuned age model based on 14C, OSL, and magnetic reversal stratigraphy shows ~100-, 40-, and 20-ka cyclicities and their evolution over the past 1.74 Ma (fig. S3). Second, spectral analyses of the AP% in the depth domain show the presence of ~34-, 15-, and 7.5-m cycles (above 90% significance level) (fig. S3). The ratios of 34:15:7.5 (=4.5:2:1) for depth cycles are close to those of astronomical frequencies of 100:40:21 (=4.7:1.9:1). Examination of changes over three sections of the sequence (493 to 573.39, 216 to 493, and 0 to 216 m) reveals that the ~7.5-m cycle appears to be stronger in the lower ~75 m, while the ~34-m cycle appears to be stronger in the top ~200 m. Third, we aligned AP% to the Chinese speleothem record of the Sanbao, Hulu, and Dongge Caves for the past 640 ka, which has independent U-Th dating. Comparison of the ETP- and speleothem-based age models of Zoige reveals a close correspondence (Fig. 3).

On the basis of the ETP tuning age, the four possible excursions of the Blake, Emperor, Panaruu, and Gilsa events are centered at 0.11, 0.49, 1.15, and 1.62 Ma, respectively, further supporting the robustness of the chronology.

Laboratory analyses

Rb/Sr ratio and CIA estimates. The Rb/Sr ratios were determined on core halves of ZB13-C2 using an XRF core scanner (ITRAX, Cox Ltd., Sweden), equipped with an Mo tube that was set to 30 kV and 30 mA. Scanning was performed at 2-mm resolution using an integration time of 10 s per measurement.

To test the climate indication of Rb/Sr ratios retrieved from XRF scanning and the relationship between Rb/Sr and the CIA indicating weathering, we measured the elements (Al, K, Na, Ca, Rb, and Sr) of 481 bulk sediment samples on an Axios advanced wavelength-dispersive XRF spectrometer (WD-XRF; PANalytical, Ea Almelo, The Netherlands). Samples were first air dried and then ground and homogenized in an agate mortar. The carbonate fraction of the homogenized samples was removed using 2 M HCl and then rinsed thoroughly four times in deionized water and dried. Each 5-g sample was compacted into a disc of 32-mm diameter under 30 metric ton force for measurements. The CIA defined as the molecular ratio was calculated using the formula CIA = 100 × Al2O3/(Al2O3 + Na2O + K2O + CaO).

LOI measurements. Volumetric subsamples of ~2 cm3 of 3273 samples from core ZB13-C2 were used for LOI analysis. The samples were weighed and oven dried at 105°C to estimate moisture content. LOI at 500°C was then used to estimate organic matter content. The weight lost at 950°C was multiplied by 1.36 to derive the carbonate content (53).

LOI measurements were also performed on the samples from the transects in the four lakes near the coring region to investigate the relationship between carbonate content and water depth on fig. S5. In the Zoige region, there are few deep lakes in the present day suitable for modern process study. The gradients of water depths of the selected shallow water lakes may be not large enough compared with the Zoige paleolake; however, they can provide parsimonious information about the general carbonate–water depth relationship.

Pollen analysis. A total of 2787 pollen subsamples of ~2 to 8 cm3 in volume were taken at ca. 20-cm intervals. The subsamples were processed following standard procedures (54), including HCl, KOH, HF, heavy liquid flotation, and acetolysis treatments, as well as fine sieving to remove clay-sized particles. Known quantities of Lycopodium spores were added at the beginning of processing to estimate pollen concentrations. At least 300 to 500 terrestrial pollen grains were counted for each sample.

The pollen percentages are calculated on the basis of the sum of the terrestrial pollen taxa. Aquatic/spore pollen percentages are calculated on the basis of the sum of terrestrial plus aquatic spores and pollen. TGViewer software was used for the calculation of percentages and for drawing initial diagrams. The zonation division was conducted using CONISS (55) with the aid of multivariate regression tree analysis.

PCA, a multivariate numerical method for dimensionality reduction (ordination analysis), was performed on the pollen data to identify the major pattern of palynological changes. The dominant pollen taxa (with >2% in at least one sample) were used in the PCA using the CANOCO program.

Biome reconstruction

Square root transformation was applied to the pollen percentage values to increase the importance of the minor pollen taxa before biomization. A modified scheme for China by Ni et al. (56) was used to assign taxa into plant functional types (PFTs) and assemblages to principal vegetation types (biomes) on the basis of the modern ecology, bioclimatic tolerance, and spatial distribution of pollen-producing plants. The biome score calculation was performed following the standard equation and reconstruction procedure by using the PPPBase software (57). The biome with the highest affinity score or the one defined by a smaller number of PFTs (when scores of several biomes are equal) was assigned for each pollen assemblage. The biomes were then combined into five types: forest, meadow, steppe, shrubland, and desert as shown in Fig. 2.

Quantitative climate reconstruction

A total of 5405 modern pollen samples from across China, where the percentages of all the taxa cover the range in the fossil pollen samples, were used as a reference or training dataset. Locally weighted weighted-average partial least squares (LWWAPLS) (58), a combination of the modern analogue technique and the transfer function approach, was used to derive quantitative climate reconstructions. This approach helps to minimize the problem of modern sample selection for traditional transfer function methods [e.g., weighted-average (WA) and weighted-average partial least squares (WAPLS)] and analog failure (e.g., too few or too many analogues) for the modern analogue technique. It is a novel and ideal approach for climate reconstructions at the broad glacial-interglacial scale when changes can be large (58).

First, for each fossil sample, a “local” training set of analogs in the modern sample dataset was generated on the basis of the smallest squared-chord distance using the modern analogue technique. The results of leave-one-out cross-validation analysis show that the root mean square error of prediction (RMSEP) value of selecting 20 analogs is significantly lower than for 25, 30, or 50 analogues (table S4), while the represented climatic gradient does not decrease. Therefore, 20 modern analogs for each fossil sample are the optimal choice. In total, there are 597 samples in the analogue training set for all the fossil samples. The analog evaluation (fig. S4) suggests that all fossil samples have good analogs within the modern training set.

To calculate the modern climate at the sampling sites of the analogue training set, we applied a thin plate spline regression approach to interpolate climate data from 756 stations across China from the China Meteorological Data website ( using the open-source R software package Fields. Second, on the basis of the local analogue training set, WAPLS was used to reconstruct the climate variables: MAT, MTWM, MTCM, MAP, mean precipitation of the warmest month, and mean precipitation of the coldest month.

Performance statistics for the WAPLS climate reconstructions are listed in table S5; however, the RMSEP can be misleading for a large heterogeneous dataset. Therefore, a significance test, which is increasingly being used, was performed to evaluate the statistical significance of the quantitative paleoclimate reconstructions. In this test, a number of reconstructions (e.g., 999) based on random variables are reconstructed using the pollen data from the training set. The proportion of variance explained by these random reconstructions is estimated using redundancy analysis. The results show that MTWM of ZB13-C2 explains >95% of the variance of the random reconstructions, suggesting that it is statistically significant (fig. S4). In addition, MAT is the variable that is the closest to the test line. In this significance test, the summer temperature reconstructions are near the PCA axis 1 that agrees with AP%, indicating that they are highly correlated. All these results suggest that summer temperature reconstructions from ZB13-C2 are statistically significant. MTCM and all the precipitation reconstructions, however, fail to pass the test, suggesting that they are not reliable.

Bootstrap sample-specific estimates of uncertainties for each reconstructed temperature value are shown in fig. S4. Note that the possible muted tree response due to migrational distances during the interglacials in the past ~0.62 Ma might have caused the underestimation of temperature range to some degree.

Time series analysis

Time series analysis was conducted using kSpectra software. Before the analysis, the data of AP%, MTWM, Carb%, and Rb/Sr ratio were resampled at equally spaced 0.6-ka intervals using AnalySeries. They were then detrended using singular spectrum analysis (SSA) in kSpectra. A window length of 300 was selected following the suggestion that the window length should be less than about one-fifth of the point number in the time series (59). Burg covariance estimation and heuristic significant tests were used for the SSA decomposition. The multitaper method was then applied to the detrended data, and spectral peaks were evaluated against a red noise background.

To track the time-varying amplitude of orbital and suborbital periods, we calculated the continuous wavelet transform using the MATLAB code of Grinsted et al. (60). The statistical significance of wavelet power was tested relative to a red noise background power spectrum.


Supplementary material for this article is available at

This is an open-access article distributed under the terms of the Creative Commons Attribution-NonCommercial license, which permits use, distribution, and reproduction in any medium, so long as the resultant use is not for commercial advantage and provided the original work is properly cited.


Acknowledgments: We gratefully acknowledge D. Za for logistical support and the villages (Nenwa and Maiqi) that hosted our fieldwork. We thank X. Sun and Y. Tang (field and laboratory assistance); Z. Zheng, X. Xiao, C. Ma, Q. Xu, and U. Herzschuh (access to modern pollen data); B. M. Benito, A. W. R. Seddon, and R. J. Telford (advice on numerical analysis); and D. B. Jiang (plotting the regional airstream map). A. Copley, G. Ramstein, F. Chen, S. Wang, and L. Qin are acknowledged for the discussions. We thank D. Lea and two journal reviewers for their constructive comments and useful suggestions. Funding: Financial support of this research was provided by the National Natural Science Foundation of China (nos. 41690113, 41888101, and 41330105), the National Key Research and Development Program of China (no. 2016YFA0600501), and the Strategic Priority Research Program of the Chinese Academy of Sciences (no. XDA20070101). P.C.T. acknowledges funding from the Natural Environment Research Council (NE/R000204/1), H.J.B.B. and V.A.F. from the European Research Council Advanced grant 741413 Humans on Planet Earth (HOPE), and V.A.F. from the VISTA (Norwegian Academy of Science and Letters and Equinor) project IGNEX-eco (6166). Author contributions: Y.Z. initiated and designed the study. Y.Z., Q.L, F.Q., W.R., M.C., and H.L. undertook the fieldwork. Y.L., F.Q., C.L., Q.C., Q.L., H.L., Z.Z., H.W., Y.Z., H.Y. and J.Z. conducted the proxy analyses. J.G. and C.D. supervised the paleomagnetism analysis. H.Z. carried out the OSL dating. C.L. and Q.C. performed the quantitative climate reconstructions and numerical analysis with guidance/advice from H.J.B.B., V.A.F., and Y.Z. Y.Z. performed the spectral and wavelet analyses on the multiproxy data. Y.Z. and P.C.T. constructed the age model. Y.Z. led the writing of the paper in collaboration with P.C.T., while Z.G. and H.J.B.B. provided additional input. All co-authors contributed to the ideas in the paper. Competing interests: The authors declare that they have no competing interests. Data and materials availability: All data needed to evaluate the conclusions in the paper are present in the paper and/or the Supplementary Materials. The R code for pollen-based climate reconstructions as well as clustering of pollen-based and reconstructed time series are available from the corresponding author upon request.

Stay Connected to Science Advances

Navigate This Article