Late Quaternary dynamics of Arctic biota from ancient environmental genomics

Abstract

During the last glacial-interglacial cycle, Arctic biotas experienced substantial climatic changes, yet the nature, extent and rate of their responses are not fully understood^1-8. Here we report a large-scale environmental DNA metagenomic study of ancient plant and mammal communities, analysing 535 permafrost and lake sediment samples from across the Arctic spanning the past 50,000 years. Furthermore, we present 1,541 contemporary plant genome assemblies that were generated as reference sequences. Our study provides several insights into the long-term dynamics of the Arctic biota at the circumpolar and regional scales. Our key findings include: (1) a relatively homogeneous steppe-tundra flora dominated the Arctic during the Last Glacial Maximum, followed by regional divergence of vegetation during the Holocene epoch; (2) certain grazing animals consistently co-occurred in space and time; (3) humans appear to have been a minor factor in driving animal distributions; (4) higher effective precipitation, as well as an increase in the proportion of wetland plants, show negative effects on animal diversity; (5) the persistence of the steppe-tundra vegetation in northern Siberia enabled the late survival of several now-extinct megafauna species, including the woolly mammoth until 3.9 ± 0.2 thousand years ago (ka) and the woolly rhinoceros until 9.8 ± 0.2 ka; and (6) phylogenetic analysis of mammoth environmental DNA reveals a previously unsampled mitochondrial lineage. Our findings highlight the power of ancient environmental metagenomics analyses to advance understanding of population histories and long-term ecological dynamics.

Fig. 1. Site distribution (North Pole-centred view).

Samples (n = 535) from a total of 74 circumpolar sites were grouped into four geographical regions (Supplementary Information 2). The grey dashed circle indicates the Arctic Circle (66.5° N). Site IDs are labelled on the map. The corresponding information is provided in Supplementary Data 1.

Fig. 2. Climate and vegetation changes over the past 50 kyr.

a, Pan-Arctic climate changes and vegetation variations. LGM (26.5−19 ka) and Younger Dryas (YD) (12.9−11.7 ka) are indicated by grey bars. The six time intervals are indicated by light blue bars (Supplementary Information 2). The error bands denote s.e. From top to bottom (see Methods for detailed calculations): the Greenlandic ice-core δ¹⁸O ratio and snow accumulation rate; the plant Shannon diversity and the Greenlandic ice-core calcium concentration; the average modelled annual temperature and precipitation for all eDNA sampling sites; the proportion of plant growth forms; the proportion of the herbaceous plant growth forms; and the vegetation turnover rates. b, The number of observed genera in different regions. c, Regional vegetation turnovers. d, Regional vegetation morphological compositions. The sample sizes for each region and time interval are provided in Supplementary Information 2. Calculations are supplied in the Methods.

Fig. 3. Spatiotemporal models to retrodict the explanatory factors for animal distribution.

The values indicate posterior parameter estimates of covariate effects for the models explaining the presence–absence of each animal’s eDNA. Only covariates included in the model with lowest Watanabe–Aikaike information criterion are shown (Methods). The dots represent the posterior means, and the whiskers represent the posterior 2.5% and 97.5% quantiles. Covariate effects of which the 2.5% and 97.5% quantiles are both negative (red), and effects of which the 2.5% and 97.5% quantiles are both positive (blue) are indicated.

Fig. 4. Mammoth distribution and mitochondrial haplotypes.

A total of 78 mammoth mitochondrial genomes and 159 eDNA-identified mammoths (79 among them were assigned to mitochondrial haplotypes) are shown. Records of dated mammoth fossils are also plotted. All samples older than 26.5 ka were combined into the pre-LGM interval.

Extended Data Fig. 1

Circum-Arctic plant abundance variations and vegetation similarity clustering. a, Pan-Arctic plant abundance heatmap. b,  Spatial vegetation dissimilarities. Pairwise spatial beta-diversities (dissimilarities between every two plant communities) against the geographical distances between the two communities. c, Non-metric Multidimensional Scaling (NMDS, k=3) on vegetation communities.

Extended Data Fig. 2

Regional vegetation differences and climate changes. a, Vegetation similarities between each two regions. All identified plant genera across sites in a region during a time interval were merged as a plant assemblage. Spatial beta-diversity between every two assemblages were calculated and illustrated. NAt, North Atlantic; WcS, Northwest and central Siberia; ES, Northeast Siberia; Nam, North America. b and c, Modelled annual temperature and precipitation in different regions. Means of the modelled annual temperature and precipitation values (Methods) at all eDNA sampling sites within a region at each 1,000-year time step were calculated. The changing trends are illustrated.

Extended Data Fig. 3

Regional plant abundance heatmaps. Heatmaps show the relative abundances of the 40 abundant plant genera in each region.

Extended Data Fig. 4. Environmental explanatory factors for animal distribution, and plant NMDS components.

a, Posterior parameter estimates of covariate effects for the models explaining the presence/absence of each animal’s eDNA using climate, human presence and plant NMDS as explanatory variables. The dots represent the posterior means, and the whiskers represent the posterior 2.5% and 97.5% quantiles. The colour red denotes covariate effects whose 2.5% and 97.5% quantiles are both negative, while the colour blue denotes covariate effects 2.5% and 97.5% quantiles are both positive. b, The plant genera and morphological forms correlated to the 3 components of plant NMDS. Plant genera (morphological forms) are ranked by the p-value of t-test, and only the top 20 Pearson correlations are shown. The colour red denotes negative correlations while the colour blue denotes positive correlations.

Extended Data Fig. 5. Distribution chronologies for woolly rhinoceros, bison, horse, caribou, hare, wolf, and vole.

We combined our DNA results and the fossil records (available for woolly rhinoceros, bison, and caribou). Samples older than 26.5 ka were combined into Pre-LGM; samples younger than 4.2 ka were combined into the Late Holocene.

Extended Data Fig. 6. Mammoth mitochondrial phylogenetic tree.

For placed eDNA samples the number of supporting single-nucleotide polymorphisms is given in braces (Methods). IDs for the Wrangel Island population are underlined.