Glacial legacies on interglacial vegetation at the Pliocene-Pleistocene transition in NE Asia

Abstract

Broad-scale climate control of vegetation is widely assumed. Vegetation-climate lags are generally thought to have lasted no more than a few centuries. Here our palaeoecological study challenges this concept over glacial-interglacial timescales. Through multivariate analyses of pollen assemblages from Lake El'gygytgyn, Russian Far East and other data we show that interglacial vegetation during the Plio-Pleistocene transition mainly reflects conditions of the preceding glacial instead of contemporary interglacial climate. Vegetation-climate disequilibrium may persist for several millennia, related to the combined effects of permafrost persistence, distant glacial refugia and fire. In contrast, no effects from the preceding interglacial on glacial vegetation are detected. We propose that disequilibrium was stronger during the Plio-Pleistocene transition than during the Mid-Pliocene Warm Period when, in addition to climate, herbivory was important. By analogy to the past, we suggest today's widespread larch ecosystem on permafrost is not in climate equilibrium. Vegetation-based reconstructions of interglacial climates used to assess atmospheric CO₂ -temperature relationships may thus yield misleading simulations of past global climate sensitivity.

ARTICLE

Received 28 May 2015 | Accepted 17 May 2016 | Published 24 Jun 2016

Broad-scale climate control of vegetation is widely assumed. Vegetation-climate lags are generally thought to have lasted no more than a few centuries. Here our palaeoecological study challenges this concept over glacialinterglacial timescales. Through multivariate analyses of pollen assemblages from Lake Elgygytgyn, Russian Far East and other data we show that interglacial vegetation during the Plio-Pleistocene transition mainly reects conditions of the preceding glacial instead of contemporary interglacial climate. Vegetationclimate disequilibrium may persist for several millennia, related to the combined effects of permafrost persistence, distant glacial refugia and re. In contrast, no effects from the preceding interglacial on glacial vegetation are detected. We propose that disequilibrium was stronger during the Plio-Pleistocene transition than during the Mid-Pliocene Warm Period when, in addition to climate, herbivory was important. By analogy to the past, we suggest todays widespread larch ecosystem on permafrost is not in climate equilibrium. Vegetation-based reconstructions of interglacial climates used to assess atmospheric CO2temperature relationships may thus yield misleading simulations of past global climate sensitivity.

1 Periglacial Research Section, Alfred Wegener Institute, Helmholtz Centre for Polar and Marine Research, Telegraphenberg A43, 14473 Potsdam, Germany.

2 Institute of Earth and Environmental Sciences, Faculty of Sciences, University Potsdam, 14479 Potsdam-Golm, Germany. 3 Department of Biology, University of Bergen, Postboks 7803, 5020 Bergen, Norway. 4 Bjerknes Centre for Climate Research, Bergen, Norway. 5 Environmental Change Research Centre, University College London, London WC1E 6BT, UK. 6 Institute of Geology and Mineralogy, Faculty of Mathematics and Natural Sciences, University of Cologne, 50674 Kln (Cologne), Germany. 7 Institute of Geology and Petroleum Technologies, Faculty of Natural Sciences, Kazan Federal University, Kazan 420008, Russia. 8 Department of Geosciences, University of Massachusetts Amherst, Amherst, Massachusetts 01003-9297, USA. Correspondence and requests for materials should be addressed to U.H. (email: mailto:[email protected]

Web End [email protected] ).

NATURE COMMUNICATIONS | 7:11967 | DOI: 10.1038/ncomms11967 | http://www.nature.com/naturecommunications

Web End =www.nature.com/naturecommunications 1

DOI: 10.1038/ncomms11967 OPEN

Glacial legacies on interglacial vegetation at the Pliocene-Pleistocene transition in NE Asia

Ulrike Herzschuh1,2, H. John B. Birks3,4,5, Thomas Laepple1, Andrei Andreev6,7, Martin Melles6

& Julie Brigham-Grette8

ARTICLE NATURE COMMUNICATIONS | DOI: 10.1038/ncomms11967

Natural vegetation patterns in space and time are assumed to reect mainly climate variability on all environmentally relevant temporal and spatial scales1,2. For example, it

has been shown that vegetation turnover in space can substitute for turnover in time on millennial timescales3, which ts with the ndings that the climatic niches of individual plant taxa are generally stable over even longer timescales4. Although lag times between vegetation and climate changes have been extensively discussed with respect to tree migration patterns in Europe and North America following the last glacial stage57, most proxy studies do not suggest major vegetationclimate disequilibria8. Relevant model evidence is lacking. Patterns of Late Pleistocene interglacial vegetation histories (marine isotope stage (MIS) 11, 9, 7 and 5) as inferred from long-term pollen records from Greece9 and the Massif Central10 are, in addition to their established relationship with interglacial climate, hypothesized to reect the different locations of glacial refugia, creating unique interglacial migrational patterns. This accords with the suggestion that in addition to climate, post-glacial migration limitations might also explain the current distribution of many European plant species11. Furthermore, an analysis integrating modern plant community composition, their climatic niches and climate change since the Last Glacial Maximum implies that North and South American forests are still responding to last glacial climate change12. However, systematic investigations of changes in vegetationclimate relationships over glacial-interglacial timescales are lacking and thus the concept of glacial legacies on interglacial vegetation has not been demonstrated using proxy data.

High-northern latitudes are ideally suited to detect former non-equilibrium patterns as climate change is particularly strong due to polar amplication13. Furthermore, the southern expansion of permafrost soils during glacial climate states14 may have increased the distance between glacial refugia and potential interglacial refugial areas. However, non-climatic vegetation drivers such as herbivory15,16, re17 or soil disturbance18 should also be considered. Finally, tectonic impacts and their climatic effects may represent a hidden driver of long-term vegetation change in the Arctic19.

The only high-latitude pollen record continuously spanning several glacial-interglacial cycles is from Lake Elgygytgyn in Chukotka, Russian Far East ranging from 3,530 to 2,150 kyr ago BP (refs 20,21), covering both the Mid-Pliocene Warm Period (MPWP: 3,5302,900 kyr ago) and the Plio-Pleistocene transition (PPT: 2,9002,150 kyr ago). The identication of vegetationclimate disequilibrium requires proxy records that reliably reect regional climate change and are independent of the Lake Elgygytgyn pollen data and of sufcient temporal length and resolution. Unfortunately, no independent climate record from Lake Elgygytgyn or any other arctic terrestrial site of sufcient data quality exists. Accordingly, only marine records can be used. Here we assume and later conrm that the globally integrating benthic isotope stack (LR04)22 and nearest sea-surface temperature (SST)23 records from the Pacic can be used to approximate the climate conditions at our study site. While all climate records indicate that the relative intensities of glacials and interglacials vary over time, a behaviour likely related to differences in orbital conguration24 and internal climate dynamics, these intensities appear to be largely globally coherent.

To assess the impact of herbivory, re and soil disturbance on vegetation, common non-pollen palynomorphs in lake sediments can be used including Sporormiella25 (a coprophilous fungal spore), Gelasinospora26 (a fungal spore characteristic of burned soils) and Glomus27,28 (an endomycorrhizal fungal spore characteristic of disturbed soils), respectively.

Here through multivariate analyses we show that variation among interglacial pollen assemblages at Lake Elgygytgyn through the PPT is best explained by long-term drivers (that is, proxies for climate and the vegetation condition of the preceding glacial), which we interpret to indicate major vegetationclimate disequilibrium. In contrast, short-term drivers (that is, proxies for contemporaneous climate, soil disturbance/erosion and/or herbivory) best explain statistically the variations in the pollen data of the PPT glacial and MPWP interglacials and glacials. From these results we conclude that the severe glacial climate and related extensive permafrost during some PPT led to distant glacial refugia, particularly of evergreen tree taxa, that hindered their subsequent interglacial establishment and the development of forests attaining equilibrium with climate. Our results imply that the current widespread larch ecosystem actually represents a transitional vegetation type reecting severe last glacial (MIS2) conditions rather than contemporaneous Holocene interglacial (MIS1) climate.

ResultsVegetation change on glacialinterglacial timescales. A principal curve of the Lake Elgygytgyn pollen record, which summarizes complex multivariate data in one dimension, shows strong temporal variations (Fig. 1 and Supplementary Figs 1 and 2, 47% variation explained). The most obvious turnover at the glacialinterglacial timescale has already been described and related to 41 kyr obliquity insolation cyclicity20. Furthermore, the principal curve is consistent with a landscape openness curve based on the pollen-based biome-reconstruction technique of the same pollen record29. This suggests that millennial-scale vegetation variability generally reects global climate shifts between glacial and interglacial states known from many proxy-based reconstructions22,23,30. The sequence of pollen types along the principal curve (Fig. 2) generally reects their present-day relative occurrences in vegetation ranging from arctic steppe-tundra to dark needle-leaf forests31. As we focus on variability at glacialinterglacial timescales, similar principal curve analyses were performed with pollen spectra averaged for each marine isotope stage (see Methods). They show a very similar trend (Fig. 1, 75% variation explained: Supplementary Fig. 3). Inferred glacialinterglacial vegetation variability generally conrms SST reconstructions23, suggesting that the glacialinterglacial amplitude during the MPWP was rather low compared with that of the PPT. Furthermore, the principal curve suggests that differences in pollen composition among glacials and among interglacials are in the same range as glacialinterglacial variability, particularly during the PPT. As expected, the pollen spectra-specic contribution to compositional turnover ( Local-Contribution-to-Beta-

Diversity32), here used as a measure of the uniqueness of the pollen composition in a marine isotope stage, is generally high for stages with extreme principal curve values.

A globally coherent pattern of stage intensities. To test our assumption that we can approximate the climate conditions at our site with remote marine temperature and isotope records, we rst investigate the global pattern of glacial and interglacial intensities. We compare the time series of glacial peak temperatures of a global set of SST records (Supplementary Table 1 and Supplementary Fig. 4) with the glacial intensities recorded in the LR04 stack (see Methods). The same analysis is also performed on the interglacial intensities. To ensure that most of the records overlap, we analyse the correlation during the PPT (2,1502,900 kyr ago). The glacial intensities in all records cores show a signicant (Po0.05) positive correlation to the inverted

2 NATURE COMMUNICATIONS | 7:11967 | DOI: 10.1038/ncomms11967 | http://www.nature.com/naturecommunications

Web End =www.nature.com/naturecommunications

NATURE COMMUNICATIONS | DOI: 10.1038/ncomms11967 ARTICLE

Vegetation

Age Large herbivores Fire Global ice volume Pacic SSTDense Common Unique Rare Seldom Large Cold

Abundant Frequent Small Warm

Open

Soil erosion

Weak Strong

2,100 2,200 2,300 2,400 2,500 2,600 2,700 2,800 2,900 3,000 3,100 3,200 3,300 3,400

MIS 82

MIS 82

MIS 84c MIS 83

MIS 84c MIS 83

MIS 84a MIS 84b

MIS 84a MIS 84b

Mid-Pliocene Warm Period Plio-Pleistocene-transition

MIS 86 MIS 85

MIS 85

MIS 87

MIS 86

MIS 88 MIS 89

MIS 88 MIS 87

MIS 89

MIS 90 MIS 92

MIS 94

MIS 91

MIS 90

MIS 92 MIS 94

MIS 91

MIS 93

MIS 93

MIS 95

MIS 99 MIS 101

MIS 103 G 1 G 3

G 7 G 9

MIS 96

MIS 95

MIS 98 MIS 97

MIS 98 MIS 97

MIS 96

MIS 100 MIS 99

MIS 100

MIS 102 MIS 101

MIS 102

MIS 104

MIS 103

MIS 104

G 2 G 1

G 2 G 4

G 4

G 6 G 5

G 3

G 6 G 5

G 8

G 7

ka BP

G 8

G 10 G 9

G 10

G 12

G 11

G 12

G 11

G 14 G 13

G 16 G 15

G 14 G 13

G 16 G 15

G 17

G 17

G 18

G 20

G 19

G 18

G 19

G 21

G 20

G 21

G 22 K 2 K 1

G 22 K 2 K 1

KM 2 KM 1

KM 3

KM 2 KM 1

KM 4 KM 5

KM 4 KM 3

KM 5

KM 6

M 1

M 2

KM 6 M 1

MG 2

M 2 MG 2

MG 1

MG 7

MG 3

MG 4

MG 4 MG 3

MG 5

MG 5

MG 6

MG 6

3,500

MG 8

MG 1

MG 7

MG 8

1.5

2.0 0.0

0.5

1.0

1.5

2.0 0 0.02 0.04

Local-contribution-to-diversity (LCBD) based on pollen data from Lake Elgygytgyn

1.0

0.5

0.0

4.0

O of marine cores (Lisiecki & Raymo, 2005)22

3.6

3.2

2.8

C

Principal curve (PC) values based on Lake Elgygytgyn pollen data (single sample data-set)

Principal curve (PC) values based on Lake Elgygytgyn pollen data (averaged stage data-set)

Sporormiella (sqrt%) in Lake Elgygytgyn

0 1 2 0.0 0.4 0.8 1.2

Gelasinospora (sqrt%) in Lake Elgygytgyn

0.0 0.5 1.0 1.5

Glomus (sqrt%) in Lake Elgygytgyn

20 25

Alkanone-based SST left: ODP 1208 (Venti et al., 2013)

right: ODP 846 (Lawrence et al., 2006)

Figure 1 | Proxy records of vegetation change and of changes of potential vegetation drivers during the MPWP and the PPT. Stratigraphic plot including the principal curve of the original pollen data set and of the averaged stage data set, and the Local-Contribution-to-Beta-Diversity (LCBD) of each stage, all based on analyses of the Lake Elgygytgyn pollen data. It also includes environmental data of inverted globally averaged oxygen stable isotope measurements of benthic foraminifera (d18Ostack)22, SST of ODP1208 (mid-latitude Pacic)30 and ODP846 (Tropical Eastern Pacic)23, and Sporormiella sqrt% (herbivory indicator), Gelasinospora sqrt% (re indicator) and Glomus sqrt% (soil disturbance/erosion indicator) relative to the vascular-plant pollen sum from Lake Elgygytgyn21. The names of the glacials and interglacials are indicated. Red and blue labels indicate PPT interglacials with preceding weak and strong glacials, respectively, that were included in the comparison of inter-taxa relationships (see text and Fig. 5). Dotted horizontal lines are drawn to aid visual comparisons. sqrt, square root.

LR04 stack (Fig. 3) with a mean correlation of r 0.81.

Interglacial intensities also show a positive correlation at all core sites with the inverted LR04 stack (mean r 0.69). The weakest

correlations are obtained for the tropical cores from the Indian Ocean (ODP722) and Atlantic (ODP662). As the interglacial temperatures at these cores are also at the temperature limit of the specic (UK33) proxy (Supplementary Table 1), it is unclear if this is a climatic signal or a proxy artefact.

While a systematic Plio-Pleistocene analysis, omitting pollen-based records to avoid circular reasoning, is limited to marine temperature proxies, analysing the Late Quaternary period also allows the inclusion of other climate proxy records. Here Antarctic ice cores provide independent evidence for the coherency of glacialinterglacial intensities34. Our comparison of stage extremes in isotope, CO2 and methane records (Supplementary Fig. 5) show a signicant correlation of the glacial intensities (mean correlation of ice-core records with the LR04 stack r 0.52) and a very strong correlation of the

interglacial intensities (r 0.91). This result strongly suggests that

the synchronizing factors (including greenhouse gases) are dominant over the local imprint of seasonal insolation24 in shaping the strength of the interglacials.

Drivers of vegetation change. To separate the effects of long-term and short-term drivers on interglacial vegetation, redundancy analysis (RDA) was applied to the data set of PPT averaged interglacial pollen spectra and proxy-based environmental variables (Fig. 1). We nd that variables reecting contemporaneous interglacial climate and re, that is, short-term

ecological drivers, explain only a small portion of the variation in the pollen data when included separately in RDA (15% Po0.1, 18% Po0.05, respectively, Fig. 4 and Supplementary Table 2) compared with variables reecting climate and vegetation type from the preceding glacial stage (42% Po0.001, 29% Po0.01, respectively). Though a large part of pollen variation is shared by all or is a combination of several variables, each variable explains a unique portion of variation as inferred from variation partitioning using a specic ensemble of RDA runs (Fig. 4 and Supplementary Figs 6 and 7). Preceding glacial climate explains the highest unique portion (18%) of pollen variation.

To test whether this result is sensitive on the choice of the climate record, RDAs of PPT interglacial pollen assemblages were repeated for each of the globally available climate records (Supplementary Table 1) with reasonable data quality. Results (Supplementary Table 3) reveal that the preceding glacial climate in all tested cases explains a signicant unique portion while interglacial climate does not explain a unique portion or only a small portion. Consistent with the global coherency of the glacialinterglacial intensities, this indicates that neither the selection of a specic proxy type nor the selection of a specic site would impact the major outcome of this study.

In addition, similar RDAs were applied to the other averaged data sets. We nd that the variation in the PPT glacial pollen spectra is most parsimoniously explained by the contemporary climate variable and to a lesser extent by the soil disturbance/ erosion proxy, with both explaining a unique and statistically signicant portion of the variation (Fig. 4). MPWP interglacial pollen are signicantly correlated to climate, vegetation of

NATURE COMMUNICATIONS | 7:11967 | DOI: 10.1038/ncomms11967 | http://www.nature.com/naturecommunications

Web End =www.nature.com/naturecommunications 3

ARTICLE NATURE COMMUNICATIONS | DOI: 10.1038/ncomms11967

0.6

0.6

0.6

0.8

Tsuga

Abies

Picea

Pinaceae_undif

Larix

Pinus_Haplox

0.20

0.4

0.4

0.4

0.4

0.10

0.4

0.2

0.2

0.2

0.2

0.00

0.0

0.0

0.0

0.0

0.0

0.0 1.0 2.0

0.0 1.0 2.0

0.0 1.0 2.0

0.0 1.0 2.0

0.15

0.20

0.20

0.00

Populus

Myrica

Corylus

Alnus_sp

Betula_undif

Lycopodium_clav

0.20

0.10

0.08

0.0 1.0 2.0

0.10

0.10

0.10

0.05

0.10

0.04

0.00

0.00

0.00

0.00

0.00

0.0 1.0 2.0

0.0 1.0 2.0

0.0 1.0 2.0

0.0 1.0 2.0

0.0 1.0 2.0

0.0 1.0 2.0

0.30

0.20

0.20

0.20

0.30

0.6

Lycopodium_sp

Equisetum

Osmunda

Betula_alba

Alnus_frut

Betula_nana

0.6

0.0 1.0 2.0

0.20

0.10

0.10

0.4

0.4

0.10

0.10

0.2

0.2

0.00

0.00

0.00

0.00

0.0

0.0

0.0 1.0 2.0

0.0 1.0 2.0

0.0 1.0 2.0

Salix

0.12

0.20

Ericales

Lamiaceae

0.20

Polygonum_viv

Saxifraga

Cyperaceae

0.20

0.0 1.0 2.0

0.4

0.4

0.08

0.10

0.10

0.10

0.2

0.04

0.2

0.00

0.0

0.00

0.00

0.00

0.0

0.0 1.0 2.0

0.0 1.0 2.0

0.0 1.0 2.0

0.0 1.0 2.0

0.0 1.0 2.0

0.20

0.20

0.15

0.08

Valeriana

Rosaceae

Papaveraceae

Lycopodium_ann

Pteridium

Anthoceros

0.4

0.10

0.10

0.10

0.10

0.04

0.2

0.05

0.00

0.00

0.00

0.0

0.00

0.00

0.0 1.0 2.0

0.0 1.0 2.0

0.0 1.0 2.0

0.12

0.0 1.0 2.0

Botrychium

Encalypta

0.4

0.8

0.30

0.0 1.0 2.0

Polypodiaceae

Sphagnum

Riccia

Selaginella_rupestris

0.4

0.08

0.4

0.10

0.15

0.2

0.04

0.2

0.00

0.00

0.0

0.0

0.00

0.0

0.0 1.0 2.0

0.0 1.0 2.0

0.0 1.0 2.0

0.0 1.0 2.0

0.0 1.0 2.0

0.0 1.0 2.0

0.30

Trilets_undif

0.15

Selaginella_sang

Huperzia

0.20

0.0 1.0 2.0

Indeterminata

0.20

0.0 1.0 2.0

Ranunculaceae

0.30

0.0 1.0 2.0

Caryophyllaceae

0.20

0.10

0.20

0.10

0.10

0.10

0.10

0.15

0.05

0.00

0.00

0.00

0.00

0.00

0.00

0.0 1.0 2.0

0.0 1.0 2.0

0.0 1.0 2.0

0.0 1.0 2.0

0.0 1.0 2.0

0.15

Brassicaceae

Asteraceae

Poaceae

0.6

Artemisia

Thalictrum

0.20

0.0 1.0 2.0

0.6

0.0 1.0 2.0

0.20

0.10

0.4

0.4

0.10

0.05

0.2

0.10

0.2

0.00

0.00

0.0

0.0

0.00

0.0 1.0 2.0

0.0 1.0 2.0

0.0 1.0 2.0

0.0 1.0 2.0

0.0 1.0 2.0

Figure 2 | Fitted response curves of each pollen or spore taxon along the principal curve. Open circles are the observed Hellinger-transformed abundance along their principal curve locations, and the solid red line is the optimized smoother from the nal iteration of the principal curve.

The taxa are arranged according to visual inspection of their response curves optima (from right to left and from top to bottom). The vertical axis is the Hellinger-transformed pollen taxon abundance (%) and the horizontal axis is distance along the principal curve.

the preceding glacial, herbivory and re when each variable (or variable set) is considered separately in RDA. However, variation partitioning shows that the major part of the explained

variation is shared by all or a combination of several variables. MPWP glacial pollen spectra are most strongly related to contemporaneous glacial climate and herbivory (Fig. 4).

4 NATURE COMMUNICATIONS | 7:11967 | DOI: 10.1038/ncomms11967 | http://www.nature.com/naturecommunications

Web End =www.nature.com/naturecommunications

NATURE COMMUNICATIONS | DOI: 10.1038/ncomms11967 ARTICLE

Correlation glacial intensities

180 90 W 0 90 E 180

Correlation interglacials intensities

180 90 W 0 90 E 180

60 N 60 N

ODP1239

ODP1208

ODP1208

ODP1143

ODP1143

ODP662

ODP722

ODP1239

ODP662

ODP722

EQ EQ

ODP846

ODP982

DSDP607

ODP846

ODP982

DSDP607

60 S 60 S

1 0.7 0.3 0 0.3 0.7 1

Correlation R

Figure 3 | Correlation of glacial and interglacial intensities of marine UK33 temperature records with the inverted LR04 stack. All correlations are calculated during the core period (2,1502,900 kyr ago), except ODP-982, which is evaluated for 2,5003,500 kyr ago as the data resolution between 2150 and 2500 kyr ago does not allow determination of the minima and maxima. Black circles indicate signicance at P 0.05; green circles indicate

records with interglacial peaks near the limit of UK33 thermometry (427 C). (See Supplementary Fig.4 and Supplementary Table 1 for details and references of records.)

Differences among PPT interglacial vegetation histories. On the basis of the RDA biplots (Supplementary Fig. 6a,b) where variables for the preceding glacial are constrained to the rst RDA axes and all the remaining known environmental variation is partialled out, we conclude that interglacials with preceding weak glacials and prevailing tree-tundra, such as G9, G7, G3, MIS101 and MIS93, are characterized by Pinus (pine) and Picea (spruce) forests (see also interglacial sequences in Supplementary Figs 8 and 9). Furthermore, these interglacials are closely located in the principal curve versus Local-Contribution-to-Beta-Diversity plot (Fig. 5a). Interglacials with preceding warm glacials are also characterized by similar inter-taxa relationships as inferred from the statistically signicant t of all interglacial pairs of ordination-derived taxa scores shown by Procrustes analyses and associated Protests (see Methods, Fig. 5b, see also taxa-specic residuals in Supplementary Fig. 10). This suggests that they reect similar interglacial vegetation histories. In summary, during the PPT, interglacials with preceding warm glacial stages share similar evergreen forest vegetation, which we assume is driven by and is in equilibrium with the contemporaneous climate.

In contrast, Pinus and Picea are absent or rarely occur during interglacials with a preceding cold glacial-stage climate and extensive arctic steppe such as G1, MIS103, MIS99, MIS95 and MIS87 while deciduous shrubs (in particular Alnus (alder)) and wetlands (indicated by high Cyperaceae (sedge) and Sphagnum (bog-moss) values) are abundant (Supplementary Figs 2, 8 and 9; see also taxa-specic residuals in Supplementary Fig. 10). Some of these interglacial spectra strongly contribute to the Local-Contribution-to-Beta-Diversity when compared with spectra from other stages with similar principal curve values (Fig. 5a). Also, the inter-taxa relationships differ strongly among single interglacials as indicated by high Protest P values (Fig. 5b and Supplementary Figs 8 and 9). In summary, interglacials with preceding strong (cold) glacial stages lack Pinus and Picea and their vegetation characteristics are non-analogous among each other.

DiscussionWith respect to PPT interglacials, we infer from our analyses of the Lake Elgygytgyn pollen record that the environmental state of the preceding glacial was more relevant for subsequent interglacial vegetation development, in particular for the presence or absence of evergreen conifers, than the contemporaneous interglacial climate, contradicting the hypothesis of a close vegetationclimate equilibrium. One alternative explanation for the observed weak relationship between proxy series for interglacial climate and interglacial pollen data would be that the analysed d18Ostack and mid-latitude Pacic SST do not represent the regional interglacial climate around Lake Elgygytgyn due to proxy-specic or regional-specic differences. However, our analysis shows that the strength of an interglacial is largely coherent in the PPT between all available high-quality SST records from the Northern Hemisphere (Fig. 3), and is coherent between different proxy records, including Antarctic temperature and greenhouse gases in the last 800 kyr (Supplementary Fig. 4). This is consistent with theoretical and model considerations that climate variations on long timescales are typically associated with broad spatial scales35,36 and with evidence from climate-model simulations using orbital forcing as well as greenhouse-gas forcing37.

Such a local climate hypothesis would also be inconsistent with our nding that the climate of the preceding glacial as recorded by the d18Ostack and mid-latitude Pacic SST should inuence the local vegetation of the interglacial. Hence, we conclude that the most parsimonious explanation for our results is that glacial climate and vegetation conditions are the main drivers of interglacial vegetation.

Accordingly, we interpret the absence of Pinus and Picea in Lake Elgygytgyn pollen spectra from interglacials following strong glacial stages to be related to three main factors: permafrost extent, distant glacial refugia, and extensive res, plus possible interactions (Fig. 6). Today the vast continuous permafrost zone in NE Asia is covered by mono-specic larch forests with a deciduous-shrub understorey9. Occurrences of

NATURE COMMUNICATIONS | 7:11967 | DOI: 10.1038/ncomms11967 | http://www.nature.com/naturecommunications

Web End =www.nature.com/naturecommunications 5

ARTICLE NATURE COMMUNICATIONS | DOI: 10.1038/ncomms11967

PPT interglacials

V:13**

PPT glacials

F:4*

CFPV:3

FPV:9

CPV:15

C:10*

C:23**

S:5* SC: 3

FP:3

***

**

*

40

20

40

20

**

.

P:18**

*

Total: 61***

Total: 46***

Total: 31**

C H F V P S

MPWP interglacials

MPWP glacials

F:1

HF:7

CF:13

H:3

CHFV:5

H:12*

C:18*

HV:9

C:1

40

20

CV:15

40

20

*

*

* * *

V:9.

CH:9

*

Total: 30**

C H F P

V

Climate Herbivory Fire

Preceeding Stage Climate

Preceeding Stage Vegetation

Soil Disturbance/

S Erosion

*** P < 0.001 ** P < 0.01* P < 0.05. P < 0.1

Figure 4 | Drivers of vegetation change of interglacials and glacials of the MPWP and the PPT. Bar charts indicating the variation of the pollen spectra explained by single environmental variables (climate (C), herbivory (H), re (F), preceding-stage climate (P), preceding-stage vegetation (V) and soil disturbance erosion (S)) for each of the analysed data sets for glacials and interglacials of the PPT and of the MPWP. Venn diagrams indicate the results of variation partitioning with those variables included that separately explain a statistically signicant part of the variation (Po0.1). Total variation explained (in %) is indicated in the lower left part of each sub-plot (see Supplementary Table 2 and Supplementary Fig. 6/7 for further detailed RDA results).

Results imply that vegetation conditions of interglacials of the PPT are best explained by long-term drivers such as the climate and vegetation condition of the preceding glacial, instead of contemporaneous conditions such as climate or herbivory. In contrast, vegetation conditions of the PPTglacials and MPWP interglacials can be best explained by short-term drivers. Accordingly, we assume that vegetation conditions of interglacials with preceding cold glacials were in strong vegetationclimate disequilibria.

evergreen taxa (for example, Picea obovata and Pinus sibirica: for Pinus pumila see Supplementary Note 1) are restricted to the southern and western permafrost limit characterized by a summer thaw depth 41.5 m and azonal vegetation, such as on deep soils in river beds9. Permafrost with a shallow active layer is also assumed to limit the northward extent of Picea in Alaska33. In contrast, Larix gmelinii and Larix cajanderi can grow on soils with an active-layer depth of o40 cm (ref. 38). Thus, the question arises whether glacial-stage permafrost conditions are able to persist through strong interglacial warming and hence hinder or slow down the invasion of evergreen taxa and are

therefore the key environmental link causing vegetationclimate disequilibrium. It can be assumed that permafrost during strong glacials was particularly extensive and deep because widespread continental-scale ice sheets were absent here and would not insulate soils from the very cold and harsh glacial climate39. Furthermore, during strong glacials, permafrost was likely more widespread due to enhanced continentality related to the absence of in-land moisture transport (that is, no snow in winter), perhaps as a result of lower sea level (6070 m)40 and more persistent sea-ice over the Arctic ocean41. Though permafrost was inferred to persist even to a few metres below the surface in

6 NATURE COMMUNICATIONS | 7:11967 | DOI: 10.1038/ncomms11967 | http://www.nature.com/naturecommunications

Web End =www.nature.com/naturecommunications

NATURE COMMUNICATIONS | DOI: 10.1038/ncomms11967 ARTICLE

a b

0.004

mid-Pliocene warm period Plio-Pleistocene-Transition

MIS95

MIS99

MIS85

G03

G07

MIS103

MIS93

MIS101

G09

G01

0.5

0.003

0.1

Protest -values

0.05

LCBD

0.002

0.01

0.005

0.001

0.001

0.0 0.5 1.0 1.5 Principal curve value

Interglacials with

preceding warm glacials

Interglacials with

preceding cold glacials

Figure 5 | Comparison of pollen-based vegetation characteristics of PPT interglacials with preceding weak glacial conditions to PPT with preceding strong glacial conditions. (a) Principal curve values versus Local-Contribution-to-Beta-Diversity (LCBD) values for all MPWP and PPT interglacials and glacials. Red labels denote interglacials with a preceding warm glacial (G09, G07, G03, MIS101 and MIS93), and blue labels denote interglacials with preceding cold glacials (G01, MIS103, MIS99, MIS95 and MIS85). It is clear that averaged pollen spectra of interglacials with preceding cold glacials strongly contribute to the LCBD when compared with spectra from other stages with similar principal curve values, in contrast to averaged spectra from interglacials with preceding warm glacials. (b) Boxplots of P values obtained for pairwise Procrustes comparison of ordination-derived variable scores of 13 taxa among 5 selected PPT interglacials with preceding weak (warm) glacials and among 5 selected PPT interglacials with preceding strong (cold) glacials. Whiskers represent the 1.5 time length of interquartile range. Results indicate that inter-taxa relationships among interglacials with a preceding cold glacial differ strongly among each other (that is, indicated by high Protest P values) while they show a signicant t among interglacials with preceding warm glacials. We assume that the vegetation histories of interglacials with preceding cold glacials were unique as a result of distant glacial refugia causing broad-scale vegetation-climate disequilibria.

Cold glacials with sparse vegetation and remote refugia of pine and spruce

Mild glacials with dense vegetation and near-by refugia of pine and spruce

Interglacials with evergreen forests on permafrost-free soils

Interglacials with larch forests and deciduous shrub on permafrost soils

Fire

Spruce Pine

Larch

Shrub

Active layer < 2m

Perma frost

Figure 6 | Visualization of the conceptual model of the Glacial legacy on interglacial vegetation. We assume that interglacial vegetation during the PPT mainly reects conditions of the preceding glacial instead of contemporary interglacial climate. During interglacials following mild glacials, NE Asia was colonized by evergreen trees from nearby refugia. In contrast, during interglacials following cold glacials larch and deciduous shrubs expand in the area related to the combined effects of permafrost persistence and distant glacial refugia of evergreen trees, and re. This implies that vegetationclimate disequilibrium can last for many millennia.

warmer climates during the Eemian interglacial in areas that are today located in the discontinuous permafrost zone42, modelling results indicate that the active-layer depth responds rather directly to climate43 while only deep permafrost may have millennial-scale lags to insolation maxima44. However, model studies have not been performed over long timescales and most

do not include all relevant Earth-system components and their interactions. Our knowledge of permafrost response to climate change during the PPT is thus very limited.

The more southerly extensive permafrost during strong glacial stages possibly shifted the glacial refugia of Pinus and Picea farther south in addition to the severe glacial climate.

NATURE COMMUNICATIONS | 7:11967 | DOI: 10.1038/ncomms11967 | http://www.nature.com/naturecommunications

Web End =www.nature.com/naturecommunications 7

ARTICLE NATURE COMMUNICATIONS | DOI: 10.1038/ncomms11967

For example, the vast steppes recorded at Lake Elgygytgyn during MIS96 stretched as far south as the Lake Baikal region, which accompanied the widespread extinction of evergreen needle-leaf taxa45. In contrast, Larix (larch) probably survived in NE Russia during all glacials, by analogy to its last glacial maximum (LGM) history46. The pollen compositional turnover of the PPT interglacials indicates the establishment of dense larch forests with shrub Alnus (Supplementary Fig. 9). These forests in particular, when occurring in combination with a thick Sphagnum cover47, insulate the soil from climate, thereby stabilizing the permafrost38 and likely hinder the germination and establishment of Picea and Pinus. Such an ecosystem with a high-standing biomass in a warm, continental early interglacial climate is likely to experience frequent res. This is indicated by the maximum Gelasinospora abundance coincident with the LarixAlnus-type optimum in the interglacial sequences (Supplementary Fig. 9). By a study of succession in present forests from central Siberia it has been shown that the ratio of deciduous (Larix and short-lived shrubs) to evergreen taxa cover (Picea and Pinus) is higher at sites characterized by frequent res48. Hence it is reasonable to assume that no single driver caused the expansion of deciduous shrubs and larch forests in warm interglacials but the combination of persistent permafrost, distant refugia and high re frequency may have slowed the invasion of Picea and Pinus and thus caused long-term vegetationclimate disequilibrium (Fig. 6).

The prevailing transient vegetation types during interglacials following strong glacials have unique characteristics, which probably originated from taxon reshufing. In contrast, during interglacials following weak glacial stages, Pinus and Picea remained in the area or in nearby refugia and invaded northern Siberia and the northern Russian Far East before shrub-rich larch forests became stabilized by feedbacks with permafrost and re. Accordingly, vegetationclimate equilibrium, characterized by a xed vegetation composition, was reached in a consistent manner. Hence, we conclude that disequilibrium versus equilibrium conditions differentiate interglacials with strong or weak preceding glacial states rather than a bifurcation with two stable states.

Our reasoning that permafrost presence likely contributed to the slow invasion of evergreen taxa is further supported by results from numerical analyses with the remaining averaged data sets. We nd that the effect of preceding glacial drivers is much smaller during MPWP interglacials when widespread permafrost was probably absent. In contrast, contemporaneous climate during the MPWP was a main driver. For example, the warm climate supported the expansion of dark needle-leaf forests during MPWP interglacials, in particular during the two oldest recorded interglacials (MG7 and MG5) that were also among the warmest with the smallest ice-sheet extent22,23. Furthermore, MPWP and PPT glacial vegetation is not related to interglacial climate or vegetation conditions but is best explained statistically by contemporaneous climate. In other words, severe climate change during a glacial-to-interglacial transition causes a leading-edge disequilibrium due to slow invasion8, while a slow die-back of long-lived trees (trailing-edge disequilibrium) is unlikely given the severe climate deterioration during an interglacialglacial transition; the Arctic is known for its strong feedback mechanisms13.

In addition to long-term drivers, our study provides evidence that short-term drivers can disturb the vegetationclimate relationship on glacialinterglacial timescales. Our analyses suggest that MPWP glacial and interglacial vegetation is signicantly related to changes in herbivory (Fig. 4). It is assumed that large herbivores, particularly mammoths and mastodons but also deer and bison, have the potential to act as

ecological keystones in high-latitude vegetation15,16,49. The herbivore fauna of the MPWP was potentially even richer than that of MIS2 and potentially exerted its impact far into the forested areas. This is suggested by, for example, fossil records of camel ancestors on Ellesmere Island (Canadian High Arctic) formerly covered by boreal larch forests50. The varying herbivore abundance, as inferred from the coprophilous Sporormiella spore record during the MPWP, may have led to marked vegetation variation. Our results suggest that herbivory during MPWP glacials turned a deciduous shrub-tundra (less Betula nana type (dwarf birch) in the pollen spectra, Supplementary Fig. 6f) into an open steppe forest (more Larix and Poaceae (grasses)) comparable to the mammoth steppe proposed for the LGM in Beringia51. A parallel decline of large herbivores and deciduous-shrub encroachment also occurs at the end of the late-glacial in North America52,53 and this pattern is supported by modern exclosure experiments53. Accordingly, our results indicate that herbivory may represent a signicant vegetation driver on glacialinterglacial timescales acting independently of climate, thereby weakening the vegetationclimate link. However, no vegetationherbivore relationship is indicated for the PPT, which contradicts the exploitationecosystem hypothesis54,55. According to this hypothesis unproductive ecosystems can only support food-limited grazers that exert a strong control on vegetation. More productive ecosystems such as occurred during the MPWP, in contrast, should be characterized by community-level trophic cascades, which should, according to this hypothesis, result in a weaker vegetationherbivore relationship. However, the weak vegetationherbivore ties during the PPT can potentially be explained by a local functional extinction of the megafauna, which was previously assumed when Sporormiella declined to o2% of the pollen sum53. However, it remains unknown if our observation of generally lower herbivore densities during the PPT compared with the MPWP represents a regional peculiarity or an Arctic-wide picture. Interestingly, the two major Sporormiella declines, that is, at the M2/M1 and G18/G17 transitions, occur when the magnitude of late-glacial warming was particularly strong. By analogy with the climate-driven extinction of some megaherbivores during the MIS2/Holocene transition56 we can only speculate that random processes such as extinctions during severe warming, at least partly, explain the observed herbivore density pattern around Lake Elgygytgyn. Furthermore, the potential refugia for cold-restricted plants and animals in the Asian Arctic are, compared with the North American Arctic, particularly small because the arctic coast is located further to the south and because the Bering Land Bridge was submerged during interglacials40, which may have enhanced random effects in interglacial colonization processes.

Our results suggest that vegetationclimate disequilibrium may have lasted for millennia. With the background of our results and the severity of the LGM (MIS2) we hypothesize that todays widely distributed larch forests lacking Picea and Pinus may not be in equilibrium with contemporary climate. Support for this comes from the fact that compared with Larix, which colonized Yakutia during the late-glacial, evergreen taxa expanded only during the early- and mid-Holocene57. However, they did not signicantly extend their ranges beyond the present-day limit during the Holocene thermal maximum, in contrast to Larix46. Furthermore, evergreen taxa in southern and western Siberia probably represent late-successional vegetation in present-day larch-forest regions but they only become established when re frequency is low48.

Vegetation records are commonly used to infer past temperature quantitatively, which is then compared with modelled or reconstructed atmospheric CO2 concentrations to

8 NATURE COMMUNICATIONS | 7:11967 | DOI: 10.1038/ncomms11967 | http://www.nature.com/naturecommunications

Web End =www.nature.com/naturecommunications

NATURE COMMUNICATIONS | DOI: 10.1038/ncomms11967 ARTICLE

assess Earth-system sensitivity58. When either the modern analogues or fossil records originate from vegetation in disequilibrium with climate, the inferred climate may be misleading. As expected from our results, the PPT temperature reconstruction based on Egygytgyn pollen data30 poorly reects the global SST23 and global ice-volume22 trends. This suggests that pollen-based quantitative climate reconstructions are less reliable than previously assumed, particularly for interglacials. However, independent regional climate records are necessary to quantify the uncertainties of pollen-based climate reconstructions. Likewise, only transient Earth-system simulations that include realistic long-term tree-population and permafrost dynamics and run over several glacialinterglacial cycles can resolve the timing and effectiveness of vegetation-related feedback mechanisms and thus provide realistic model-based estimates for Earth-system sensitivity.

With the appearance of permafrost and the functional disappearance of herbivory, the PPT at northern high latitudes represents a valuable natural experiment in a palaeoenviron-mental setting with the potential to understand better the dynamics of Earths past states when new system components are included.

Methods

Numerical analyses of original Lake Elgygytgyn pollen data. This paper uses a pollen data set from Lake Elgygytgyn21. Of the 1,166 samples that were analysed, 1,043 contained pollen. Of these, 844 samples with at least 150 terrestrial pollen grains and spores were included in the numerical analyses here. The pollen sum is 4300 in most samples.

Principal curves59 were derived to best summarize, in a mathematical sense, the overall pollen compositional turnover and thus to maximize the variation represented in one dimension. To increase the signal-to-noise ratio, we a priori decided to include only half of the original 92 reliably recorded taxa (that is, taxa recorded at least twice in the entire data set); thus, all taxa with occurrences in 474 samples were retained for the principal curve analyses. Sample scores on the rst correspondence axis were used as a starting curve. The principal curve was tted using cubic smoothing splines allowing the complexity of the individual smoothers used in the local averaging step to vary between pollen taxa, using generalized cross-validation to select the optimal degree of smoothing for each taxon. A penalty term of 1.4 was used to increase the cost of degrees of freedom in the generalized cross-validation calculations60. Analyses were implemented in R using the pcurve package.

Analyses of averaged interglacial and glacial pollen data. To obtain averaged interglacial and glacial data sets, principal curve minima and maxima were rst assigned to marine isotope stages22 based on the original age model20. MIS84 was split into two glacial stages (MIS84a/c) and an interglacial stage (MIS84b) because a full glacialinterglacial turnover is clearly reected in the pollen data. Boundaries between glacial and interglacials were set at the point where the mean principal curve value was crossed, calculated from the minimum principal curve values for interglacials and the maximum principal curve values for glacials. Original counts of those samples that are below/above the median PC values for each glacial/interglacial were then summed for all taxa. This step was done to ensure that the averaged data sets represent the characteristic pollen assemblages of each stage rather than the transitional assemblages between stages. To reduce the noise in the data set, percentages were calculated including vascular-plant pollen taxa that occur in at least 12 stages with a value of at least 2% once. Finally, a data set of 32 averaged interglacial pollen spectra and another of 32 averaged glacial pollen spectra were obtained.

Principal curve analysis of the combined glacialinterglacial averaged data set was run with a similar conguration as for the original data set, only the penalty term was set to 2 to account for the lower number of samples. Local-Contributions-To-Beta-Diversity32, that is, representing an estimation of the compositional uniqueness of a stage was likewise calculated using the combined averaged glacialinterglacial data set.

Correlation analyses of global climate data sets. We analysed the global LR04 benthic inverted d18Ostack (ref. 22) as well as a global collection of high-resolution

SST records. We required that the full PPT (2,1502,900 kyr ago) is covered with a mean resolution of at least 5 kyr per sample. Ten published UK33 records and the LR04 benthic stack meet these criteria (Supplementary Table 1).

The irregular raw data of every core were low-pass ltered using a Gaussian Kernel Smoother with an effective cutoff frequency of 1/15 kyr. The assignment of the glacial minima and interglacial maxima was performed in two steps. (1) Local maxima and minima were automatically identied from the smoothed time series.

(2) All records were visually compared using their published age model, and the corresponding glacial minima and interglacial maxima were manually chosen from the minima and maxima of step 1. We only interpreted time periods where a 100-kyr window did not contain a data gap 415 kyr to reduce the possibility that we missed major peaks or troughs. In some cases, where no clear assignment was possible, no assignment was made; thus most records only contain a subset of all extremes. For ODP1082 and ODP1012, only a small number of maxima and minima could be assigned, probably due to the poor data quality of the records in the time period of interest. We thus excluded both these cores from further analysis. As the ODP982 SST record does not have a high enough resolution until 2,600 kyr ago, we correlated ODP982 and LR04 on an extended time window (2,6003,500 kyr ago) to obtain a meaningful overlap with the LR04 stack. The results are summarized in Supplementary Fig. 4.

We also investigated whether the coherency between glacial and interglacial intensities was restricted to the marine SST records. For the last 800 kyr, we analysed several proxies derived from Antarctic ice cores34: water isotopes related to local temperature; a methane record potentially related to low-latitude wetland dynamics; and the CO2 record representing the global carbon cycle. We compared these records with the global benthic d18Ostack (ref. 22) using the same method as

in the previous section (Gaussian smoothing, automatic local extremes identication and manual assignment according to the marine isotope stages22) (Supplementary Fig. 5). We used the same technique of comparing local extremes instead of analysing the same points in time to be consistent with our main analyses of the Plio-Pleistocene time series.

Redundancy analyses. Constrained ordination analyses were run to relate statistically the variation in the averaged glacial and interglacial pollen data setsto the proxy environmental variables. RDA was chosen over canonical correspondence analysis because of the compositional gradient lengths of the pollen data being o2 s.d. as revealed by initial detrended correspondence analyses, indicating that linear-based ordination methods are appropriate with these data sets. Because the boundary conditions strongly changed, we a priori separated the data of the MPWP (3,5102,900 kyr ago, 12 glacials, 12 interglacials) from the PPT (2,900 2,150 kyr ago, 20 glacials, 20 interglacials) following denitions of the PLIOMAX group61. This boundary also accounts for the nding that the major decrease in atmospheric CO2 concentrations occurred around 2,900 kyr ago (ref. 62), that the phase of interglacials having higher than present sea-level ended40 and that ice-rafted debris markedly increased in glacial stages from MIS16 onwards63. Setting the boundary to that time also agrees with the results of the biome reconstruction using Lake Elgygytgyn pollen data, which indicates that MIS17 was the last interglacial that was at least temporarily dominated by cool coniferous forests29.

The climate proxy set includes the d18Ostack record22 and two alkenone-based

SST records from the Pacic as no other suitable regional climate proxy data are available. Total organic carbon (TOC) values or Si/Ti ratios from Lake Elgygytgyn were not used as proxy environmental variables because terrestrial vegetation changes could have strongly impacted lake chemistry64, thereby risking circular reasoning in our analyses. We selected the LR04 d18Ostack record22, which, in

addition to the temperature-sensitive ice-sheet extent, reects the sea-level position and thus variability in continentality in northern Asia. Furthermore, the proxy set includes the high-resolution alkenone-based SST record from the northern mid-latitude ODP site 1208 in the Pacic30. This SST record represents the closest record to Lake Elgygytgyn and is located near to the potential glacial refugia in northeastern Asia and has, furthermore, a high data quality. Because the SST record from ODP site 1208 does not fully cover the MPWP the alkenone-based SST from the eastern Pacic ODP846 site23, which is well correlated to ODP site 1208 during the PPT65, is used for analysing the MPWP. The d18O stack and SST records were interpolated to 2.5 kyr ago time slices and the maximum/minimum value for each interglacial/glacial was used. The same proxy records were used as variables to represent the preceding stage climate.

Square-root-transformed values of Sporormiella spore percentages (relative to the vascular-plant pollen sum), a coprophilous fungus commonly occurring in the palynomorph record, were used as a herbivore intensity proxy31. Likewise, percentages of Gelasinospora spores, which are often found at times of massive charcoal input66,67, were used as a re-frequency proxy. Support for this comes from the observation that the occurrence of Gelasinospora is signicantly related to burnt soils26, which is explained by its rapid growth rate and the resistance of its spores to high temperatures compared to other fungi. Preceding stage PC values of the combined glacial/interglacial averaged data set were used to explore long-term vegetation drivers. Furthermore, percentages of Glomus spores were used to indicate changes in the disturbance and erosion of soils28.

The following sets of RDAs were run on each of the Hellinger-transformed glacial and interglacial averaged pollen data sets. Those variables (or variable sets), which when solely included in an initial RDA yielded a statistically signicant relationship with the data set (Po0.1), were included together in variation partitioning to extract the uniquely explained variation by each single environmental variable (or variable set) and the shared explained variation. Explained variations represent adjusted r2 values68.

Analyses of PPT interglacial vegetation histories. We aimed to investigate whether pollen inter-taxa relationships among PPT interglacials with preceding

NATURE COMMUNICATIONS | 7:11967 | DOI: 10.1038/ncomms11967 | http://www.nature.com/naturecommunications

Web End =www.nature.com/naturecommunications 9

ARTICLE NATURE COMMUNICATIONS | DOI: 10.1038/ncomms11967

warm glacials (that is, warm interglacials) are different from those among inter-glacials with preceding cold glacials (that is, cold interglacials). For that purpose, warm interglacials (that is, having a combination of low d18Ostack and high SST

values compared with all interglacials during the PPT, and having a minimum of 12 samples: G9, G7, G3, MIS101 and MIS93) and cold glacials (that is, having a combination of high d18Ostack and low SST, and having a minimum of 12 samples:

G1, MIS103, MIS99, MIS95 and MIS83) were selected based on the climate proxy information and the minimum sample number. Inter-taxa relationships of a single interglacial are reected by the species scores in an ordination. Accordingly, principal component analyses (with scaling focused on variable correlations) were performed on the original pollen samples of the respective warm or cold interglacials. Analyses were based on Hellinger-transformed pollen percentages of taxa that occur at least three times in each stage. To compare the ordination results of two interglacials the correlation of variable scores on the rst four principal component analysis (PCA) axes was assessed using Procrustes analysis69 and its associated Protest test70. Finally, the Protest P values gained from pair-wise comparisons among all warm interglacials and among all cold interglacials were presented in boxplots.

Statistical programmes used and data availability. All analyses were implemented in R version 3.03 using functions in packages analogue 0.12-0, calibrate 1.7.2, pcurve 0.6-5, vegan 2.0-10 and the function beta.div. Averaged data sets and environmental proxy data used for redundancy analyses are available at http://www.pangaea.de

Web End =www.pangaea.de .

References

1. Shugart, H. H. & Woodward, F. I. Global Change and the Terrestrial Biosphere (Wiley-Blackwell, 2011).

2. Vazquez-Rivera, H. & Currie, D. J. Contemporaneous climate directly controls broad-scale patterns of woody plant diversity: a test by a natural experiment over 14 000 years. Global Ecol. Biogeogr. 24, 97106 (2014).

3. Blois, J. L., Williams, J. W., Fritzpatrick, M. C., Jackson, S. T. & Ferrier, S. Space can substitute for time in predicting climate-change effects on biodiversity. Proc. Natl. Acad. Sci. USA 110, 93749379 (2013).

4. Svenning, J.-C. Deterministic Plio-Pleistocene extinctions in the European cool-temperate tree ora. Ecol. Lett. 6, 646653 (2003).

5. Smith, A. G. Problems of inertia and threshold related to post-glacial habitat changes. Proc. R. Soc. London B 161, 331342 (1965).

6. Prentice, I. C., Bartlein, P. J. & Webb, T. Vegetation and climate change in eastern North-America since the Last Glacial Maximum. Ecology 72, 20382056 (1991).

7. Giesecke, T. et al. The pace of Holocene vegetation changetesting for synchronous developments. Quat. Sci. Rev. 30, 28052814 (2011).

8. Svenning, J.-C. & Sandel, B. Disequilibrium vegetation dynamics under future climate change. Am. J. Bot. 100, 12661286 (2014).

9. Tzedakis, P. C. & Bennett, K. D. Interglacial vegetation succession: a view from southern Europe. Quat. Sci. Rev. 14, 967982 (1995).

10. Cheddadi, R. et al. Similarity of vegetation dynamics during interglacial periods. Proc. Natl. Acad. Sci. USA 102, 1393913943 (2005).

11. Normand, S. et al. Postglacial migration supplements climate in determining plant species range in Europe. Proc. R. Soc. B 278, 36443653 (2011).

12. Blonder, B. et al. Linking environmental ltering and disequilibrium to biogeography with a community climate framework. Ecology 96, 972985 (2015).

13. Miller, G. H. et al. Arctic amplication: can the past constrain the future? Quat. Sci. Rev. 29, 17791790 (2010).

14. Vandenberghe, J. et al. The Last Permafrost Maximum (LPM) map of the Northern Hemisphere: permafrost extent and mean annual air temperatures, 25-17 ka BP. Boreas 43, 652666 (2014).

15. Owen-Smith, R. N. Megaherbivores: The Inuence of Very Large Body Size on Ecology (Cambridge Univ. Press, 1988).

16. Gill, J. L. Ecological impacts of the late Quaternary megaherbivore extinctions. New Phytol. 201, 11631169 (2014).

17. Bond, W. J. & Keeley, J. E. Fire as a global herbivore: the ecology and evolution of ammable ecosystems. Trends Ecol. Evol. 20, 387394 (2005).

18. Komrkov, V. Vegetation type hierarchies and landform disturbance in arctic Alaska and alpine Colorado with emphasis on snowpatches. Vegetatio 106, 155181 (1993).

19. Kienast, F. et al. Paleontological records indicate the occurrence of open woodlands in a dry inland climate at the present-day Arctic coast in western Beringia during the Last Interglacial. Quat. Sci. Rev. 30, 21342159 (2011).

20. Brigham-Grette, J. et al. Pliocene warmth, polar amplication and stepped Pleistocene cooling recorded in NE Arctic Russia. Science 340, 14211427 (2013).

21. Andreev, A. A. et al. Late Pliocene and Early Pleistocene environments of the north-eastern Russian Arctic inferred from the Lake Elgygytgyn pollen record. Clim. Past 10, 10171039 (2014).

22. Lisiecki, L. E. & Raymo, M. E. A Pliocene-Pleistocene stack of 57 globally distributed benthic d18O records. Paleoceanography 20, PA1003 (2005).

23. Lawrence, K. T., Liu, Z. & Herbert, T. D. Evolution of the Eastern Tropical Pacic through Plio-Pleistocene Glaciation. Science 312, 7983 (2006).

24. Prescott, C. L. et al. Assessing orbitally-forced interglacial climate variability during the mid-Pliocene Warm Period. Earth Planet. Sci. Lett. 400, 261271 (2014).

25. Commandini, O. & Rinaldi, A. C. Tracing megafaunal extinctions with dung fungal spores. Mycologist 18, 140142 (2004).

26. Widden, P. & Parkinson, D. The effects of a forest re on soil microfungi. Soil Biol. Biochem. 7, 125138 (1975).

27. Anderson, R. S., Homola, R. L., Davis, R. B. & Jacobson, G. L. Fossil remains of the mycorrhizal fungal Glomus fasciculatum complex in postglacial lake sediments from Maine. Can. J. Bot. 62, 23252328 (1984).

28. Van Geel, B. in Tracking Environmental Change Using Lake Sediments, vol. 3: Terrestrial, Algal and Silicaceous Indicators. (eds Smol, J. P. Birks, H. J. B. & Last, W. M.) 99119 (Kluwer, 2001).

29. Tarasov, P. E. et al. A pollen-based biome reconstruction over the last 3.562 million years in the Far East Russian Arctinew insights into climate vegetation relationships at the regional scale. Clim. Past 9, 27592775 (2013).

30. Venti, N. L., Billups, K. & Herbert, T. D. Increased sensitivity of the Plio-Pleistocene northwest Pacic to obliquity forcing. Earth Planet. Sci. Lett. 384, 121131 (2013).

31. Isaev, A. P. et al. in The Far North: Plant Biodiversity and Ecology of Yakutia. (eds Troeva, E. I. et al.) 143260 (Springer, 2010).

32. Legendre, P. & DeCaceres, M. Beta diversity as the variance of community data: dissimilarity coefcients and partitioning. Ecol. Lett. 16, 951963 (2013).

33. Lloyd, A. Ecological histories from Alaskan tree lines provide insight into future change. Ecology 86, 16871695 (2005).

34. Masson-Delmotte, V. et al. EPICA Dome C record of glacial and interglacial intensities. Quat. Sci. Rev. 29, 113128 (2010).

35. Mitchell, J. M. et al. An overview of climatic variability and its causal mechanisms. Quat. Res. 6, 481493 (1976).

36. Jones, P. D., Osborn, T. J. & Briffa, K. R. Estimating sampling errors in large-scale temperature averages. J. Clim. 10, 25482568 (1997).

37. Willeit, M., Ganopolski, A., Calov, R., Robinson, A. & Maslin, M. The role of CO2 decline for the onset of Northern Hemisphere glaciation. Quat. Sci. Rev.

119, 2234 (2015).38. Osawa, A., Zyryanova, O. A., Matsuura, Y., Kajimoto, T. & Wein, R. W. (eds) Permafrost Ecosystems: Siberian Larch Forests (Springer 502, 2010).

39. French, H. The Periglacial Environment (Wiley 478 p, 2007).40. Miller, K. G., Mountain, G. S., Wright, J. D. & Browning, J. V. A 180-millon-year record of sea level and ice volume variations from continental margin and deep-sea isotopic records. Oceanography 24, 4053 (2011).

41. Polyak, L. et al. History of sea ice in the Arctic. Quat. Sci. Rev. 29, 17571778 (2010).

42. Reyes, A. V., Froese, D. G. & Jensen, B. J. L. Permafrost response to last interglacial warming: eld evidence from non-glaciated Yukon and Alaska. Quat. Sci. Rev. 29, 32563274 (2010).

43. Jiang, Y., Zhuang, Q. & ODonnell, J. A. Modeling thermal dynamics of active layer soils and near-surface permafrost using a fully coupled water and heat transport model. J. Geophys. Res. 117, D11110 (2012).

44. Galushkin, Y. Numerical simulation of permafrost evolution as a part of sedimentary basin modelling: permafrost in the Pliocene-Holocene climate history of the Urengoy eld in the West Siberian basin. Can. J. Earth Sci. 34, 935948 (1997).

45. Demske, D., Mohr, B. & Oberhansli, H. Late Pliocene vegetation and climate of Lake Baikal region, southern East Siberia, reconstructed from palynological data. Palaeogeogr. Palaeoclim. Palaeoecol. 184, 107129 (2002).

46. Polezhaeva, M., Laxcoux, M. & Semerikov, V. L. Cytoplasmic DNA variation and biogeography of Larix Mill. in Northeast Asia. Mol. Ecol. 19, 12391252 (2010).

47. Harden, J. W., Manies, K. L., Neff, J. C. & Turetsky, M. R. Effects of wildre and permafrost on soil organic matter and soil climate in interior Alaska. Global Change Biol. 12, 113 (2006).

48. Schulze, E.-D. et al. Factors promoting larch dominance in central Siberia: re versus growth performance and implications for carbon dynamics at the boundary of evergreen and deciduous conifers. Biogeosciences 9, 14051421 (2012).

49. Gill, J. L., Williams, J. W., Jackson, S. T., Donnelly, J. P. & Schellinger, G. C. Climatic and megaherbivory controls on late-glacial vegetation dynamics: a new high-resolution, multi-proxy record from Silver Lake, Ohio. Quat. Sci. Rev. 34, 6680 (2012).

50. Rybczynski, N. et al. Mid-Pliocene warm-period deposits in the High Arctic yield insight into camel evolution. Nat. Commun. 4, 1550 (2013).

51. Blinnikov, M. S., Gaglioti, B., Walker, D. A., Wooller, M. J. & Zazula, G. D. Pleistocene graminoid-dominated ecosystems in the Arctic. Quat. Sci. Rev. 30, 29062929 (2011).

52. Gill, J. L., Williams, J. W., Jackson, S. T., Lininger, K. B. & Robinson, G. S. Pleistocene megafaunal collapse, novel plant communities, and enhanced re regimes in North America. Science 326, 11001103 (2009).

10 NATURE COMMUNICATIONS | 7:11967 | DOI: 10.1038/ncomms11967 | http://www.nature.com/naturecommunications

Web End =www.nature.com/naturecommunications

NATURE COMMUNICATIONS | DOI: 10.1038/ncomms11967 ARTICLE

53. Bakker, E. S. et al. Combining paleo-data and modern exclosure experiments to assess the impact of megafauna extinctions on woody vegetation. Proc. Natl. Acad. Sci. USA 113, 847855 (2016).

54. Oksanen, L. & Oksanen, T. The logic and realism of the hypothesis of exploitation ecosystems. Am. Nat. 155, 703723 (2000).

55. Anapuu, M. et al. Spatial pattern and dynamic responses of Arctic foodwebs corroborate the exploitation ecosystems hypothesis (EEH). Am. Nat. 171, 249262 (2008).

56. Lorenzen, E. et al. Species-specic responses of Late Quaternary megafauna to climate and humans. Nature 479, 359364 (2011).

57. Binney, H. A. et al. The distribution of late-Quaternary woody taxa in northern Eurasia: evidence from a new macrofossil database. Quat. Sci. Rev. 28, 24452464 (2009).

58. Schmittner, A. et al. Climate sensitivity estimated from temperature reconstructions of the Last Glacial Maximum. Science 334, 13851388 (2011).

59. Death, G. Principal curves: a new technique for indirect and direct gradient analysis. Ecology 80, 22372253 (1999).

60. Simpson, G. L. & Birks, H. J. B. in Tracking Environmental Change Using Lake Sediments, vol. 5: Data Handling and Numerical Techniques. (eds Birks, H. J. B., Lotter, A. F., Juggins, S. & Smol, J. P.) 249327 (Springer, 2012).

61. Pliomax Group http://pliomax.org/pliowiki/index.php/Welcome_to_PLIOMAX<mac_aq id=

Web End =http://pliomax.org/pliowiki/index.php/Welcome_to_ http://pliomax.org/pliowiki/index.php/Welcome_to_PLIOMAX<mac_aq id=

Web End =PLIOMAX .

62. Seki, O., Foster, G. L., Schmidt, D. N. & Mackensen, A. Alkenone and boron-based Pliocene pCO2 records. Earth Planet. Sci. Lett. 292, 201211 (2010).

63. Flesche Kleiven, H., Jansen, E., Fronval, T. & Smith, T. M. Intensication of Northern Hemisphere glaciations in the circum-Atlantic region (3.52.4 Ma) ice-rafted detritus evidence. Palaeogeogr. Palaeoclim. Palaeoecol. 184, 213223 (2002).

64. Herzschuh, U. et al. Siberian larch forests and the ion content of thaw-lakes form a geochemically functional entity. Nat. Commun. 4, 2408 (2013).

65. Nie, J. et al. Surface-water freshening: a cause for the onset of the North Pacic stratication from 2.75 Ma onward? Global Planet. Change 64, 4952 (2008).

66. Anderson, R. S., Ejarque, A., Brown, P. M. & Hallett, D. J. Holocene and historical vegetation change and re history on the north-central coast of California, USA. Holocene 23, 17971810 (2013).

67. Froyd, C. A. et al. The ecological consequences of megafaunal loss: giant tortoises and wetland biodiversity. Ecol. Lett. 17, 144154 (2014).

68. Peres-Neto, P., Legendre, P., Dray, S. & Borcard, D. Variation partitioning of species data matrices: estimation and comparison of fractions. Ecology 87, 26142625 (2006).

69. Peres-Neto, P. & Jackson, D. A. How well do multivariate sets match? The advantages of a Procrustean superimposition approach over the Mantel test. Oecologia 129, 169178 (2001).

70. Jackson, D. A. PROTEST: a Procrustes randomisation test of community environment concordance. Ecoscience 2, 297303 (1995).

Acknowledgements

Arne Ramisch and Stefan Kruse contributed ideas and discussion to a very early version of the data-set analyses. We acknowledge Chris Brierley, Petra Dekens, Zhonghui Liu as well as the other data contributors for making their SST data sets publicly available and David Naafs for discussion of SST records. Cathy Jenks is thanked for her meticulous editing. This study was funded by the German Federal Ministry of Education and Research (grant no.03G0642A) and the German Research Foundation (grant nos. HE 3622/20 and ME 1169/24).

Author contributions

U.H. designed and performed the analyses with advice from H.J.B.B.; U.H. wrote the rst version of the manuscript that was revised by H.J.B.B.; T.L. performed the analyses of global climate proxy data. All authors contributed to discussion and revision of the nal manuscript.

Additional information

Supplementary Information accompanies this paper at http://www.nature.com/naturecommunications

Web End =http://www.nature.com/ http://www.nature.com/naturecommunications

Web End =naturecommunications

Competing nancial interests: The authors declare no competing nancial interests.

Reprints and permission information is available online at http://npg.nature.com/reprintsandpermissions/

Web End =http://npg.nature.com/ http://npg.nature.com/reprintsandpermissions/

Web End =reprintsandpermissions/

How to cite this article: Herzschuh, U. et al. Glacial legacies on interglacial vegetation at the Pliocene-Pleistocene transition in NE asia. Nat. Commun. 7:11967 doi: 10.1038/ ncomms11967 (2016).

This work is licensed under a Creative Commons Attribution 4.0 International License. The images or other third party material in this article are included in the articles Creative Commons license, unless indicated otherwise in the credit line; if the material is not included under the Creative Commons license, users will need to obtain permission from the license holder to reproduce the material. To view a copy of this license, visit http://creativecommons.org/licenses/by/4.0/

Web End =http://creativecommons.org/licenses/by/4.0/

r The Author(s) 2016

NATURE COMMUNICATIONS | 7:11967 | DOI: 10.1038/ncomms11967 | http://www.nature.com/naturecommunications

Web End =www.nature.com/naturecommunications 11