Postglacial species arrival and diversity buildup of northern ecosystems took millennia.

What drives ecosystem buildup, diversity, and stability? We assess species arrival and ecosystem changes across 16 millennia by combining regional-scale plant sedimentary ancient DNA from Fennoscandia with near-complete DNA and trait databases. We show that postglacial arrival time varies within and between plant growth forms. Further, arrival times were mainly predicted by adaptation to temperature, disturbance, and light. Major break points in ecological trait diversity were seen between 13.9 and 10.8 calibrated thousand years before the present (cal ka BP), as well as break point in functional diversity at 12.0 cal ka BP, shifting from a state of ecosystem buildup to a state where most habitat types and biotic ecosystem components were in place. Trait and functional diversity stabilized around 8 cal ka BP, after which both remained stable, although changes in climate took place and species inflow continued. Our ecosystem reconstruction indicates a millennial-scale time phase of formation to reach stable and resilient levels of diversity and functioning.

INTRODUCTION

Understanding the establishment of contemporary ecosystems over millennia is essential to anticipate how ecosystems may respond to climate changes (, ). While climate defines the broad-scale distribution of species and is the major driver of past range changes (, ), it interacts with other range-limiting processes such as dispersal, disturbance, and biotic interactions, which are not predictable from a climate-only perspective (). Reconstructing past environments in terms of edaphic and biotic factors () may facilitate understanding of the interplay among environmental factors, biodiversity, and ecosystem functioning that determines the time scale of ecosystem buildup (). Species traits reflect habitat characteristics () and, combined with pollen records, have been used to reconstruct past environments yet with limited power due to taxonomic limitations (). The higher taxonomic resolution, local provenance, and more persistent detection of sedimentary ancient DNA (sedaDNA) (, ) provide more power to understand the influence of past environmental change on ecosystem buildup and functionality.

At the present rate of climate change, the northern flora will have to move hundreds of meters per year to stay within their climate niches (). Short-term observations show elevational and latitudinal movements (), and species are predicted to move northward (). Yet, addressing how the timings of these movements translate into effective ecosystem buildup and diversity requires long-term data. However, there has been a vigorous debate about the role of postglacial dispersal lags based on pollen studies, with some arguing for tree species tracking climate changes (), while others argue for time lags (). Many pollen records show no increase in diversity over time, which is inferred as no migration lag (, ), while others show an increase that may be due to dispersal limitation (). A recent study using sedaDNA suggests an almost continuous immigration of species to northern Fennoscandia, suggesting a legacy of the last glaciation (), that could be substantial for ecosystem buildup and functionality.

Here, we reconstructed the establishment of plant taxa in arctic ecosystems using sedaDNA across 16 millennia. We focus on their traits related to colonization, providing unique information about the formation of ecosystems. Fennoscandia was almost completely glaciated by the Scandinavian Ice Sheet during the Last Glacial Maximum (LGM) 26 to 19 thousand years (ka) ago (), and we study 10 sites representing the increasing availability of colonizable land as the ice subsequently retreated (Fig. 1 and figs. S1 to S3). We expanded plant sedaDNA data from () and constructed a nearly complete regional DNA reference library (PhyloNorway, 1822 vascular plant taxa) that we used for identification. We combined the taxonomic data with a large trait dataset () to test which edaphic and biotic factors were important in postglacial ecosystem buildup. Because dwarf shrubs play a key role in the functioning of Fennoscandian ecosystems (), we further investigate the dynamics of their arrival in relation to other groups.

Fig. 1.

Appearance of 315 taxa in northern Fennoscandia in relation to deglaciation and climate.

(A) The 10 study sites (black dots) in relation to the reconstructed ice sheet (white) and mean daily 2-m air temperature for July (other colors) (). Dark gray lines indicate reconstructed surface elevation at 1000-m steps above sea level (a.s.l.) starting at 0 m a.s.l. () with modification for the Baltic region (). For additional time slices, see fig. S1. (B) Accumulated richness of all growth forms aggregated into 500-year time slices. The gray lines represent numbers of species appearing in nine samples per time slice (based on 100 resampling of all lakes within 500-year time slices; lines are not visible when they overlap with the full data), and dotted lines represent estimated first arrival date (taking into account variation in detectability; see Materials and Methods), whereas solid lines are the observed first appearance. (C) Box plot of observed first appearances (Ob, filled, vertical line = median) and estimated first arrival (Es, unfilled, vertical line = median) for each growth form. The estimated first arrival takes into account uncertainties in detection of taxa. (D) The relative abundance of each growth form is estimated as the proportion of all PCR repeats for each 500-year time slice across all lakes.

RESULTS AND DISCUSSION

The dataset

By extracting and quality checking the P6 loop of the trnL (UAA) intron from a dataset of recent low-coverage genome skims of 2051 specimens from 1899 taxa [PhyloNorway ()], we obtained a P6 loop reference library of 1969 specimens from 1822 taxa (data S1). This DNA reference library not only covers most of the species in Norway but also includes taxa from other polar regions. Our reanalyses of the sedaDNA data of () coupled with additional samples (see Materials and Methods) resulted in a total dataset of 355 samples containing 42,093,315 reads and 64,338 polymerase chain reaction (PCR) repeats of 315 taxa (data S2). For up to 227 taxa, data on species traits were available (data S3). As we ran eight PCR repeats per sample, we used these replicates as a proxy for taxon abundance (see Materials and Methods). Sedimentary ancient DNA is a powerful source of information about past floras but, just as with detecting the disappearance of a species, the first arrival is likely to go undetected because of rare and scattered occurrences. Therefore, we estimated likely arrival times on the basis of detectability of each species and used these to analyze the possible effect of species traits on arrival times (data S4).

Ecosystem buildup

Our data show that the buildup of the vegetation composed of different growth forms was characterized by both gradual increase and phases of more rapid accumulation. The earliest vegetation consisted of arctic bryophytes, forbs, and graminoids, of which the first two dominated (Fig. 1 and figs. S4 to S6). These represent the earliest post-LGM plant records known from the region (). The first dwarf shrubs were detected at ca. 14.0 calibrated thousand years before the present (cal ka BP), but their abundance increased only later from 11.0 to 10.0 cal ka BP. Aquatic macrophytes were the last to appear, from around 12.9 cal ka BP (Fig. 1), but rapidly accumulated and made up >30% of the abundance at 12.0 to 10.0 cal ka BP before they dropped to <20%. Birch, the most dominant woody taxon in the region today, appeared from 12.9 cal ka BP, with the first record probably representing shrub-birch Betula nana and later records also representing the tree-birch Betula pubescens (). In the period 13.0 to 12.0 cal ka BP, there was a rapid increase of taxa, especially forbs, graminoids, aquatic macrophytes, and vascular cryptogams. While there were an increasing number of sites and samples from around this time period (fig. S3), our resampling analyses indicate that the increase in number of taxa is not an artifact of sampling (Fig. 1B) but represents an increase in regional diversity as temperatures rise and more land becomes available (Fig. 1A). Order of appearance of taxa was similar to what has been detected in pollen and macrofossils studies (table S1), strengthening the interpretation that the observed order of appearance in our DNA data reflects the order in which they became common in the ecosystem.

The Scandinavian Ice Sheet finally disappeared around 10.3 to 9.9 cal ka BP (). The major change in vegetation composition was around 7.8 cal ka BP, when the last tree species, Alnus incana, appeared (figs. S4 and S5). Around this time, pollen studies indicate a spread of pine and birch forest in the region (). During the past 8000 years, the relative importance of different growth forms showed smaller fluctuations (Fig. 1D and fig. S6), although changes in species composition were seen (figs. S5), as also seen in pollen-based studies at many sites in the region ().

Both observed first appearance dates and estimated arrival times showed differences among growth forms (Fig. 1C and fig. S5). The median first appearance date was earliest for aquatic macrophytes, although this was the last growth form to appear. It was intermediate for forbs, trees/shrubs, and bryophytes and was latest for dwarf shrubs (9.2 cal ka BP) and graminoids (8.8 cal ka BP; Fig. 1C). We regard the data for dwarf shrubs as especially robust, as all common and even some rare dwarf shrubs in the region were detected (fig. S5). The average estimated arrival times were 1500 to 3000 years later for dwarf shrubs than other growth forms, although confidence intervals overlapped except for bryophytes (Fig. 1C, fig. S7, and data S4). Most of the dwarf shrubs in the region are Ericaceae. The late arrival of some Ericaceae species may have gone undetected in pollen studies because of low pollen production and because they are often taxonomically lumped in pollen analyses (, ). Thus, little is known about the postglacial arrival of individual species. Most of the Ericaceae have adaptations to wind or animal dispersal, but even some of these arrived late. For the arctic wind-dispersed Cassiope tetragona, the late arrival (4.6 detected, 6.0 estimated cal ka BP) may be caused by long migration distance from LGM refugia in Beringia, as indicated by contemporary genetic patterns (). The reason for late arrival of the relatively cold-tolerant wind-dispersed Phyllodoce caerulea (detected 5.4, estimated 6.3 cal ka BP) is less clear, but the discontinuous global distribution () suggests that factors other than climate may have caused restrictions in its dispersal speed. Thus, while dwarf shrubs on average show later arrival than other growth forms, the cause of these time lags may differ among the species.

Overall, our data expand knowledge of arrival due to the high number of taxa identified (315), many of which cannot be identified to the same taxonomic level from pollen () and were never detected as macrofossils in the region (). The two growth forms best studied in the past, trees and aquatic macrophytes (, , ), had early median appearance dates. Here, we see a continuous appearance of newly arrived taxa, in contrast to what has been inferred by some pollen studies (, ). Our study suggests that the ability of species to track climate changes may have been underestimated by pollen analyses, which is highly relevant for prediction of species response to ongoing climate changes.

To better understand the variation in arrival time, we focused on nine plant traits important for dispersal and establishment including microenvironmental preference (light, phosphorus, and soil disturbance) and biotic requirements (seed dispersal, mycorrhiza, and pollinator dependence). We further looked at the macroenvironmental requirement of the species including temperature optimum, continentality, and moisture requirements (fig. S8). We determined the relative importance of these traits, together with growth form, for the estimated arrival time using boosted regression tree (BRT) analysis. Temperature, as estimated from CHELSA-TraCE21k V1.0 dataset (see Materials and Methods), had the highest relative influence on the estimated arrival time and taxa with a lower temperature optimum generally arrived earlier (Fig. 2). This result was robust to different weightings that took into account uncertainty in estimated arrival time and to different complexities of the models (tree length; fig. S9). Adaptation to soil disturbance and light further determined arrival time. Species adapted to open areas with free soil patches and high light conditions tended to arrive earlier (Fig. 2). This was in accordance with univariate relationships (weighted linear regressions) that showed a positive correlation between earlier estimated arrival times and higher light requirement as well as increasing tolerance to soil disturbance (fig. S10), suggesting that factors related to establishment were determining the arrival of taxa. The relative importance of the remaining plant traits were each estimated as less than 10%. Our analyses show that, in addition to temperature, a suite of traits, rather than single traits, influenced the estimated arrival time, supporting the idea that models of future distribution of species and ecosystems may need to include additional factors (), with requirements for light and the ability to grow in already established vegetation being the most influential.

Fig. 2.

Impact of trait variables on the estimated arrival time.

(A) The relative importance of different trait variables in explaining the estimated arrival time of 227 vascular plant taxa to northern Fennoscandia based on BRT analysis. (B to D) Relationship between the trait value and establishment time predicted by the BRT analysis for the three variables with the highest relative importance: temperature optimum (B), soil disturbance (C), and light requirement (D).

The pulses of species arrivals had consequences for ecosystems’ trait diversity and hence to functional diversity. We used the abundance of taxa in each trait category as an estimate of niche space presence and occupation through time (Fig. 3 and fig. S11). There were clear shifts in all traits at the onset of the Holocene. The vegetation before 12.0 cal ka BP was dominated by cold-tolerant taxa adapted to soil disturbance and full light, as also observed during deglaciation elsewhere (). The largest fluctuations in trait composition were seen during the Late Glacial and Early Holocene (Fig. 3), with break points in trait diversity earliest in phosphorus requirement (13.9 cal ka BP) and latest in mycorrhiza (10.8 cal ka BP; Fig. 4, fig. S12, and tables S2 and S3). Pollination dependence deviated, with the highest diversity in early communities and a low-diversity break point at 10.8 cal ka BP when insect-pollinated species were at their most dominant. We also found the highest abundance of taxa independent of mycorrhiza before 10.8 cal ka BP, and variation in adaptation to mycorrhiza had a low impact on the estimated arrival time (Fig. 3). Similarly, Felde and Birks () found highly fluctuating mean trait values at 14.0 to 9.0 cal ka BP based on analyses of 58 taxa recorded as macrofossils in southwestern Norway. This is as expected, given that this was a period with large climate fluctuations through the Late Glacial Interstadial-Stadial Complex (Figs. 1A and 5A) and transition from tundra to birch forest (, ). An important question is if there was a possible migration mismatch between plants and their pollinators. We found the highest proportion of taxa independent of pollinators at >12 cal ka BP. This may be explained by the high proportion of arctic plant species adapted for asexual reproduction (, ). However, exclusively insect-pollinated taxa dominated both in abundance and in number of taxa during the period of rapid colonization 12 to 10 ka cal BP (Fig. 3 and fig. S11), and only a low proportion of the variation in arrival time was explained by variation in adaptation to pollination (Fig. 2). Similarly, variation in adaptation to mycorrhiza had a low impact on the estimated arrival time. This suggests that mismatch in colonization between plants and their pollinators or fungi symbionts is of low importance for the long-term ability of plants to track climate changes. Thus, although short-term studies raise concerns about mismatch in colonization between symbiont taxa (, ), our data show that this is of low importance for northern plants’ long-term ability to track climate changes.

Fig. 3.

Northern ecosystem changes through 16 millennia based on abundance of plant taxa with different trait values.

(A to I) Dashed vertical lines indicate the break point in diversity of each trait (see fig. S8). The break points in cal ka BP are given in brackets. Note that the break point in pollinator dependence is negative, indicating low diversity. The plots include both terrestrial and aquatic taxa except for moisture, which is for terrestrial taxa only. Taxa are quantified as detection in proportion of PCR repeats. Gray areas represent undefined traits for some taxa. For full description of traits and categories, see the Supplementary Materials. AM, arbuscular mycorrhiza; EM, ectomycorrhiza; ErM, ericoid mycorrhiza.

Fig. 4.

Breakpoints estimated by segmented regression on Shannon diversity of traits based on the proportions of PCR repeats.

Estimated break points follow facet names, and are indicated by dashed line. The blue lines represent fitted values of break point analysis, the solid circles represents the Shannon diversity, and pink shadings represent approximate 95% confidence intervals (Estimates ± 2 SE) of the break points.

Fig. 5.

New appearance of taxa in relation to climate, niche saturation, and ecosystem diversity.

(A) Mean July temperature (°C) for the 10 sites (lines) and observed number of new taxa per 500 years according to growth form (bars). (B) Accumulated trait values for each of the nine traits. (C) Shannon diversity of trait values. (D) Functional diversity measured as functional dispersion (pink fill for 95% confidence intervals of break point).

Although new taxa appeared over the entire period, taxa representing most of the trait values were found within 12.0 cal ka BP (Fig. 5, A and B). This finding was robust even when either only looking at the three lakes with the earliest records or using resampling analyses, suggesting that the observed accumulation of traits was not due to increased sampling intensity (fig. S13). The most competitive and shade-tolerant taxa appeared at 11.0 cal ka BP and the most continental at 10.5 cal ka BP. Taxa of all temperature optima appeared by 12.0 cal ka B.P. except Huperzia that has intermediate temperature requirement and appeared at 10.2 cal ka BP. Thus, most of the habitats were likely available by 12.0 cal ka BP, as also seen by the overall break point in functional diversity at this time (Fig. 4, C and D). However, 59% of observed first appearances and 41% of estimated first arrival time of the taxa were after this date, and there were large fluctuations in trait diversity (Fig. 5C) also after 12.0 cal ka BP. This suggests that, although the habitats were present, their availability may have fluctuated in time and space, which may have limited the chance for any arriving propagule to actually land in the right place for establishment. A similar inference has also been suggested for dispersal to the high arctic archipelago of Svalbard ().

Trait diversity became relatively stable only after the maximum temperature in the region was reached at about 8.0 cal ka BP (Figs. 3 and 5C). When quantifying functional diversity from all nine traits, we see similar stabilization in the past 8000 years (Fig. 5D). This suggests that it took several thousand years from when all habitats were available to actually reach stable ecosystem function. This has consequences for the evaluation of ecosystem equilibrium state (), as our data indicate that it may take millennia from a major change in climate to reach a stable ecosystem. Thus, the first change observed after an environmental change may not represent the new equilibrium. For example, the changes observed due to ongoing climate warming may not have reached a new ecosystem equilibrium and may continue to change even if no further climate changes take place.

During the past 8000 years, only minor changes were seen in ecosystem traits, although many new taxa appeared and some changes in temperature were seen (Fig. 4). Thus, the ecosystem might already have sufficient diversity for resistance to minor climate changes and further introductions, as the newcomers did not bring in new traits nor change the relative abundance of trait categories, which is also reflected by stable functional redundancy (Figs. 3 and 5 and fig. S14). Notably, most late appearing taxa were of small life forms like forbs and bryophytes (Fig. 4 and fig. S5) and may not have had the same influence as the earlier arriving trees and ferns, which cause shade and alter the structure of the ecosystems. The stability of the Late Holocene ecosystem is in strong contrast to the fluctuating ecosystem in Late Glacial and Early Holocene and supports the hypothesis that high diversity ensures ecosystem resilience (, ).

While the ability of species and ecosystems to adapt to environmental conditions in the future is expected to differ from the past because of fragmentation of habitats, disturbance, and species introduced by humans (), the polar latitudes remain the most pristine regions on Earth. Ecosystem processes similar to those that have been shown here may be expected for areas of large ice sheet melt, like Greenland, whereas areas with smaller ice caps may be closer to seed source areas that may speed up the ecosystem buildup process. The long time needed for ecosystem buildup documented for the past millennia suggests that other processes beyond dispersal will delay how ecosystems will track ongoing climate change geographically, and functioning will take time to rebuild. Larger environmental changes may cause large fluctuations in species and ecosystem traits, and if an ecosystem is severely disrupted, then it may take thousands of years to rebuild.

MATERIALS AND METHODS

Regional flora and vegetation

Northern Fennoscandia represents the most remote area that plants migrated to from unglaciated areas in Europe after the last glaciation (Fig. 1 and fig. S1). Palaeoecological reconstruction based on pollen analyses shows a sparse tundra vegetation from 14.6 cal ka BP, with the first birch woodland developing 12.0 to 10.0 cal ka BP, and low human impact (, , ). The region now has the lowest vascular plant diversity in Europe (), although it is richer than expected for its northern latitude (250 to 550 species per 75 km by 75 km), possibly because of contribution from northern refugia ().

The flora of the northernmost counties of Norway (Troms and Finnmark) is well known, although a modern, comprehensive flora is only available for Troms (areal extent 25,863 km2). A past (1934) survey of the flora of Finnmark (48,631 km2) is largely based on material collected more than a hundred years ago and is outdated. A database maintained by T.A. covering Finnmark includes 598 taxa of vascular plants (660 if subspecies are included, excluding apomictic taxa in Ranunculus auricomus coll., Hieracium and Taraxacum, each counted as one species) and a further 256 introduced taxa, in which some are well established (data S5). The dominant vegetation types in the region are birch forest in the western lowlands (up to 300 to 400 m above sea level) with an oceanic climate, pine forest up to 350 m in the interior with a more continental climate, and alpine tundra above and north of the tree line ().

Climatic data

Climate data ranging from current conditions (1950 to 1990) to the LGM (21 cal ka BP) were represented by monthly mean daily near-surface (2-m) air temperatures for the month July from the CHELSA-TraCE21k V1.0 dataset (, ). The monthly mean daily near-surface (2-m) air temperatures for July are shown for three selected time slices in Fig. 1 and for the entire time period from 22 to 0 cal ka BP in fig. S1.

The 10 lake cores and the time covered

Site and lake core characteristics are given in (). For each of the 10 lake cores, Bayesian age-depth models were previously established including for depths older than the Holocene (). We therefore reused these models, which allowed for the dating of all but six samples, which fell below basal-most radiocarbon dates. Five of these samples were from the undated Unit A of EG17 [Jøkelvatnet; core JØP-112; ()], and so we explored extrapolation of the existing age-depth model to assign sample ages. All age-depth modeling was performed following (), except as otherwise noted. First, we reran the existing Jøkelvatnet age-depth model by extending the d.max parameter to 258 to enable extrapolation, which yielded an age of 12,018 cal BP for the basal sample (257 cm; fig. S2A). However, we noted that the extrapolated curve was dependent on the accumulation rate prior (accMean), with a reduced accMean of 20 giving an age of 10,453 cal BP for this basal sample (table S4 and fig. S2B). The impact of prior choice is likely driven by two nonconcordant basal radiocarbon dates for Jøkelvatnet (). To inform which curve was more reasonable, we constructed age-depth models for Lake Storvatnet (STP-112), located 5 km from Jøkelvatnet in the same valley, and using the same parameters listed above and radiocarbon dates from (). Crucially, the transition between Units A and B in both JØP-112 and STP-112 shows the same magnetic susceptibility and loss-on-ignition profiles (), allowing for proxy curve anchoring between the lake core profiles. The accMean prior did not affect the STP-112 age-depth model, which gave anchor-inferred ages of 10,151 and 10,125 cal ka BP for the 257-cm sample from JØP-112 (table S5 and fig. S2, C and D). We therefore selected the extrapolated JØP-112 age-depth model with an accMean of 20 for Jøkelvatnet. Revised ages for all samples from the Jøkelvatnet core are presented in table S4.

The core records begin at the following: Langfjordvatnet, 16.9 cal ka BP (with the oldest sample in the dataset at 16,147 cal BP); Nordvivatnet, 12.9 cal ka BP; Sandfjorddalen, 12.4 cal ka BP; Eaštorjávri South, 11.4 cal ka BP; and Horntjernet, 11.1 cal ka BP. From around 10.5 cal ka BP, there are samples from 8 to 10 cores in each millennium (data S2 and fig. S3). Two cores cover the Late Holocene period only (Nesservatnet and Sierravannet), which could cause a late detection of taxa restricted to habitats found in these two catchments. However, as only two taxa (Potentilla egedei/anserina and Alopecurus geniculatus) were unique to these two cores, this has a negligible impact on the overall first detection of taxa.

The sedaDNA dataset

Our sedaDNA dataset consists of a supplemented version of that presented in (). All samples had been analyzed at The Arctic University Museum of Norway, where the dedicated ancient DNA laboratory is in a separate building to post-PCR facilities. DNA was extracted from 0.25 to 0.35 g of sediments following () and amplified using the gh primer targeting the vascular plant trnL p6-loop locus of the chloroplast genome () and sequenced on Illumina NextSeq platform at either Fasteris, SA (Switzerland) or the Genomics Support Centre Tromsø at the UiT The Arctic University of Norway [see full details in (, )]. We used all 371 samples from (), and in contrast to (), we included all extraction duplicates (n = 17). We supplemented this dataset with 53 samples from five of the same lake cores that either predated the Holocene or fell outside of the age-depth models used in (table S6) (). We also included 76 samples from three cores taken from Lake Øvre Æråsvatnet lakes on Andøya (table S7) (), which were later excluded because of poor DNA quality (see below). All additional samples were processed using the methods in () with core depth and other sample metadata given in table S6 and library, sequencing details, and tag-sample lookup presented in table S8. This gave a total pre–quality control (QC)–filtered dataset of 517 samples, together with 76 negative and 28 PCR laboratory (positive and negative) controls (table S9).

Construction of reference library

A comprehensive reference library is essential for correct identification to lower taxonomic levels, especially when working with short barcodes as sequence-sharing among taxa is common, and if the true species is missing from the reference library, then a match to a closely related species is likely. For the P6 loop used here, the taxonomic resolution to species level is 33% at a circumpolar scale due to the many congeneric species in Beringia and the Amphi-Atlantic region (), but resolution increases to 77 to 93% when applied to local catchments (). Thus, from our recent low-coverage genome skims of 2051 specimens from 1899 taxa [PhyloNorway ()], we extracted the P6 loop of the trnL (UUA) intron from the assembled genomes by running the ecoPCR program (). We proofread the sequences, and erroneous identifications were removed. In addition, we had access to an early version (and at time of use not yet proofread) PhyloAlps P6 loop reference library (), which we used for cross-checking our identifications of the sedaDNA sequences.

Identification and final filtering of sequences

We matched against four DNA reference libraries: (i) our new PhyloNorway library (provided in data S1), (ii) a combination of 815 arctic () and 835 boreal () vascular plant taxa and 455 arctic-boreal bryophytes () from the circumpolar region (ArcBorBryo, n = 2280 sequences of which 1053 are unique), (iii) PhyloAlps [n = 4604 specimens of 4437 taxa collected in the Alps and Carpathians (); http://phyloalps.org/], and (iv) EMBL (release 143, n = 159,748 sequences of 74,936 taxa). Matches to PhyloNorway were given priority. As 85% of the PCR repeats and 91% of the reads matched PhyloNorway and/or ArcBorBryo using a 100% identity criterion, we disregarded sequences not matching any of these two reference libraries with the following exceptions: (i) sequences matching Ericaceae as we especially focus on the immigration of dwarf shrubs, (ii) sequences with >100 repeats, or (iii) >10,000 reads. The sequences were categorized as follows: 1, if biogeographically plausible; 2, if biogeographically marginal and/or assumed sequence variant of a more common sequence; and 3, if regarded as not plausible. The remaining low frequency sequences were assumed to mainly represent PCR and/or sequencing errors, although it may also include rare taxa (table S10)—these were given category 4 and not further checked. Further, we removed sequences that were only found in positive controls. As a last step, we checked average reads and repeats in samples compared to negative controls and removed sequences with higher average number of reads in negative controls than in samples. We used this approach rather than removing all sequences detected in the negative controls because the most abundant sequences in samples are regularly also detected in negative controls (). For sequences with a 100% match to several taxa, taxa were listed if they were recorded in the region according to a regional checklist for Northern Norway (), the Norwegian Biodiversity Information Centre (www.biodiversity.no/), circumpolar distribution maps (), and/or the PanArctic Flora Checklist (); the taxonomy follows the latter. While it is possible that different species were present in the region during the past 12,000-year time period as indicated by, for example, the past occurrence of Larix sibirica, indicated by a fossil record in northern Sweden (), we are not aware of any non–Northwestern European plant taxa found in pollen or macrofossil records from the region. Thus, in cases of sequence sharing among closely related taxa, we assume that the taxa in the region are most likely correct matches but also list all matches in table S10. The P6 loop contains homopolymers that make identification difficult within some families. Our almost complete reference library made it worth checking the sequences in more detail, and we aligned sequences within the largest families using Geneious Prime 2.2 to check the nature of variation. Other taxa that regularly occur as somewhat unexpected exotics in sedaDNA studies [reviewed in ()] were Anacardiaceae (cashew nut family) and Juglandaceae (walnut family), which we now consider to be due to species of these families being used in plastic production, and Hevea brasiliensis, which is used in latex gloves. Further comments to inclusion or exclusion of taxa are given in table S10.

In total, 612 sequences had a 100% match to at least one taxon in PhyloNorway, PhyloAlps, ArcBorBryo, and/or EMBL (tables S10 and S11). Of these, 256 had a 100% match to PhyloNorway (57,821 repeats; 45,938,251 reads). A further 115 sequences had a 100% match to ArcBorBryo (10,746 repeats; 2,033,869 reads), 79 of them to bryophytes (4939 repeats; 844,277 reads). An additional 34 sequences had a 100% match to PhyloAlps (3066 repeats; 1,066,760 reads). The remaining 208 sequences had a 100% match only to EMBL (9263 repeats; 3,597,452 reads). We manually cross-checked all sequences with a 100% match to PhyloNorway and/or ArcBorBryo as well as 26 sequences with a 100% match to Ericaceae, 19 sequences with >100 repeats, and 35 sequences with >10,000 reads (some of them with both high reads and high repeats and/or match to Ericaceae). In total, 432 sequences representing 99.7% of the reads (52,481,639 reads) and 97.2% of the PCR repeats (78,568 repeats) were manually checked for consistency across libraries, biogeographical plausibility, and potential contamination and were classified in categories 1 to 4 (see above). We kept only sequences of category 1. After collapsing sequences matching the same taxon, our dataset consisted of 342 taxa (48,783,645 reads; 68,937 repeats).

We next removed low-quality samples by applying the QC approach described in (), which uses minimum metabarcoding technical and analytical quality (MTQ/MAQ) score thresholds as a filter. All negative controls associated with the additional samples were reported in (), and the Andøya negative controls did not affect the thresholds (table S7). We therefore applied the same MTQ/MAQ score filtering thresholds (MTQ ≥ 0.75; MAQ ≥ 0.1). We retained the 316 samples reported in (), five extraction duplicates, and 34 of the additional samples from four cores (tables S6 and S9). As only 16 of the 76 samples from Lake Øvre Æråsvatnet passed the QC thresholds (table S7), we omitted this site from further analysis. From the other 10 sites, a total of 355 of the 441 samples passed the QC thresholds (table S9 and fig. S3). These 355 samples contained 319 taxa (43,220,432 reads; 65,664 repeats). As the reference library for algae is poor, we removed four taxa of algae (but these data are included in data S2).

Quantification of sedaDNA data

Regarding abundance, sedaDNA studies may be taxonomically biased for various reasons, for example, unequal preservation, taphonomy, primer binding, marker length, or bias in the DNA reference library used (, ). The PCR step in the metabarcoding approach may especially cause biases in quantities and can underestimate graminoids (, ), although the close similarity of past circum-arctic vegetation reconstructed using metabarcoding () and shotgun metagenomic sequencing () suggests that this might be less of a problem than previously assumed. Possible ways to quantify metabarcoding data are using raw read counts, proportion of reads, or number of PCR repeats where species occur (any of these can also be transformed or weighted). Raw reads have the advantage that they show the data exactly as they are, but read numbers vary strongly with sequencing depth of individual samples and, while useful for evaluation of sample quality and strength of signal (, ), it may give a strongly biased estimate of abundance. Proportion of reads are widely used to quantify sedaDNA (–). This approach uses all available quantitative information in the data and works well for most purposes. However, if some taxa dominate the record, as is commonly found for Salix and aquatic macrophytes in northern ecosystems [including in shotgun metagenomics data suggesting that this may be due to other reasons than PCR bias, such as taphonomy ()], then the dominant sequence may mask patterns of other taxa. If this is the case, then using proportion of PCR repeats is a more conservative alternative (). This requires that a sufficient number of PCR repeats have been analyzed () and that they have not been merged before sequencing (). Thus, while the actual contribution of the dominant taxa may be underestimated when using proportion of PCR repeats, this approach highlights the full diversity of the taxa detected. In the Results and Discussion section, we therefore chose to quantify the sedaDNA data using our eight PCR repeats. Similar patterns in growth form and traits were observed when using the proportion of reads (data S2) or only the qualitative information, which is the presence of taxa (figs. S6 and S11). The total number of repeats across multiple samples and lakes within a time slice was used while calculating the proportion of repeats for each 500-year time slice. A similar approach was used after taking the sum of repeats of all species within a growth form while calculating repeat proportion for growth forms within a 500-year time slice.

Number of samples and sites over time

As our selection of sites represented the regional variation in deglaciation times, the earliest sample varied among sites. Both the number of sites and the number of samples increased over time (fig. S3). To assess the impact of increasing sampling intensity on species accumulation curves, we resampled the data to a maximum of nine samples for each time slice of 500 years representing all the sampled lakes. Thus, if a time slice contained more than nine samples, then we randomly selected only nine samples by choosing at least one sample from each lake that was sampled during this time. The resampled curves fitted very well with the curves generated from the full dataset (Fig. 1B). Similarly, trait accumulation based on resampled data (see below) showed no effect of sampling (fig. S13). Thus, we used the full dataset for all analyses.

Observed temporal patterns of taxa occurrence

The temporal patterns of occurrences are plotted as the number of PCR repeats (1 to 8) using the rioja (v0.9.26) and ggplot2 (v3.3.2) packages in R. Note that sequences assigned to Salix are included in “Trees & shrubs” although this genus also includes some dwarf shrubs. We performed multivariate regression tree (MRT) analysis (fig. S4) to find optimal partition in the plant sedaDNA data as it does not make any assumptions about species-environment relations and provides robust results (). While performing MRT in the mvpart (v1.6.2) package in R, we considered species assemblage as the response and sample age as the predictor variables. A major change in the vegetation was found at 7.8 cal ka BP, with secondary changes at 13.0, 9.7, 3.6, and 0.8 cal ka BP (fig. S4).

Plant trait data

We used species traits from the Flora of Sweden (). If a sequence matched to several taxa or was identified only to genus level, then we included the relevant taxa if they were found in the region today (). In this case, we used mean trait values if there was a <3 category difference within the genus/taxon group, except for moisture where we used <4 and climatic variables where we used <5. If the trait variation within a genus/taxon group was larger than these values, then we did not average but considered the trait value for this taxon as missing, as averaging may reduce model fit (). In a few instances, the sequence matched to several taxa, but only one of them was present in (). In this case, we used the trait values of the taxon present in the database for all the taxa represented by the sequence (see notes in data S3). For species lacking traits in (), we used floras to add traits where possible (data S3). Trait data were each available for up to 227 taxa. We also grouped the taxa into the following growth forms: trees and shrubs (including Salicaceae that could be either tree, shrub, or dwarf shrub), dwarf shrubs, forbs, graminoids, vascular cryptogams, bryophytes, and aquatic macrophytes. The aquatics were defined as vascular plants with a moisture index of 10 (temporarily inundated) to 12 (deep permanent water) of ().

We initially selected 14 of the 30 traits available in () that we assumed relevant for arrival time and ecosystem buildup (see the Supplementary Materials): heat requirement, cold requirement, temperature optimum, continentality, light requirement, moisture requirement, soil reaction (pH), nitrogen, phosphorus, salinity, soil disturbance, pollinator dependence, seed dispersal, and mycorrhiza. To minimize redundancy, we carried out a multivariate analysis with mixed quantitative variables and factors [fig. S8; using the function dudi.mix in the package ade4 in R () to identify correlated variables and keep a simplified set of variables]. For this analysis, missing values for continuous traits were replaced by the mean and taxa with missing values for categorical traits were excluded (two taxa). A single taxon with orchid mycorrhiza was also excluded.

Estimated first arrival time based on observed date

To estimate the likely arrival time, we used the approach developed by Wang et al. () [see also Wang and Marshall () for a review]. The method was developed to estimate extinction dates, but it can be used to estimate arrival time by reversing the time parameter. Compared to other approaches such as that in (), it allows for flexible modeling of recovery potential. The model has two components, a function f(x,λ) describing the recovery potential as a function of date x and parameter λ (given that the species is present) and a distribution for the true arrival time, θ. The recovery function f in the study by Wang et al. () was a modified Beta distribution called the reflected Beta distribution. Given a set of dates x = (x1, x2, xn), the likelihood L(θ,λ|x) is the product for i = 1 to n of f(xi|θ,λ). Posterior distributions of arrival time θ were obtained using a Bayesian approach, that is, assuming a prior distribution for θ and λ given in the study by Wang et al. ()

We used the R script provided by Wang et al. () to estimate arrival time. This method allows the estimation of arrival time and the associated uncertainty (90% credible interval—the interval is bounded on the earliest observed date). To account for the uncertainty in the estimated arrival times in the analyses relating traits and arrival times, we considered the width of credible intervals as an approximation of the precision of each estimate.

For some species, estimated arrival times predated the oldest analyzed sample. However, because our analysis focuses on the studied lakes and the time frame covered by the retrieved cores, these dates were adjusted to the date of the oldest analyzed sample at 16,147 cal BP. Credible intervals were also given an upper bound at this date for consistency. For species detected in the oldest sample, this resulted in a credible interval of 0, which we replaced by 200 years, corresponding roughly to the precision of the age-depth models.

The difference between detected and estimated first arrival was lowest for aquatic macrophytes (data S4), which are generally well detected as their sedaDNA is readily deposited in the lake sediments (). Bryophytes [which are regarded as bycatch using this method ()] and graminoids [which are often underrepresented in metabarcoding studies (, )] showed larger uncertainties, whereas the other growth forms showed intermediate uncertainties.

Effect of species traits on the estimated first arrival time

The effect of species traits on the estimated arrival time was assessed with a BRT analysis. BRTs provide a flexible method for ecological predictions based on multiple independent variables that can handle nonlinear relationships, interactions, and missing values and do not require the choice of a single best model (). We used estimated arrival time as the response variable and the 10 selected plant traits as predictor variables for BRTs. To account for differences in the precision of estimated arrival times, we weighted the observations by the inverse of the 90% credible interval. The BRT model was fitted using the gbm.step function of the dismo package in R () according to the recommendations in (). We adjusted the learning rate to obtain an optimal number of trees of >1000 (learning rate = 0.001) and ran the selection process for tree complexities between 2 and 4. As a final result, we chose the model with a tree complexity of 4 as it resulted in the lowest holdout deviance (Fig. 2), but the results from the other analyses are shown in fig. S9. To test for sensitivity to the chosen weighting of observations, we repeated the BRT analysis with weaker weighting (inverse of the square root of credible intervals) and stronger weighting (inverse of the square of credible interval).

In addition, we estimated univariate relationships of estimated arrival time with growth form and the nine variables selected as representing the likely most important traits using univariate weighted linear models, where the weights accounted for the uncertainty in estimated arrival times. An approximation of the estimated variance for each arrival time was obtained by considering that the length of 90% credible intervals was roughly 4 SE, so that the variance was given by ((credible interval)/4)2. This variance corresponds only to the uncertainty in arrival times (aka measurement error), thus using those values as weights would exaggerate the importance of arrival times with high precision and lead to high leverage of these samples. Hence, we modeled arrival times by considering that an additional process (residual) error determined the overall residual variance of the regression. The magnitude of the process error was chosen such as to result in a total variance of 1 for the weighted regression model, i.e., that the overall variance was given by the sum of the measurement error variance and the process variance. Results of the linear models were displayed as parameter estimates with 95% confidence intervals.

Reconstruction of past environments based on species traits

To reconstruct past environments, we grouped the data into 500-year time bins and calculated the proportion of PCR repeats of each trait category for each of the traits (temperature optimum, continentality, moisture, light, nitrogen, soil disturbance, pollinator dependence, mycorrhiza, and seed dispersal). We chose the same traits as used for analyzing factors important for arrival, as we focus mainly on factors important for ecosystem buildup. Because of the potential bias in sedaDNA data (see above), in the Results and Discussion section, we semiquantify the sedaDNA data using our eight PCR repeats. Similar patterns in traits were observed when using only the presence of taxa (fig. S13) or proportion of reads (data S3).

Trait and functional diversity through time

To identify the major changes in plant community traits, we performed break point analysis for trait diversity through time (Fig. 4 and fig. S11). We first calculated Shannon diversity using the proportion of PCR repeats for each trait value rounded to the whole number within each 500-year time bin using the library “vegan” (). We then fitted linear models with Shannon diversity as the response and time bins as the predictor. We repeated the entire process stated above for trait diversity calculated on the basis of the proportion of species within each 500-year time bin. We considered linear models as the base models to perform single break point analysis of Shannon diversity of individual traits.

To quantify the aspect of changes in functional diversity, we computed abundance-weighted functional dispersion that is the average trait distance of individual species to the trait centroid of all species in the community (). The metric is very similar to other commonly used metrics such as Rao’s functional diversity (, ). We used the mFD package () to calculate functional dispersion. We then used break point analysis to explore whether a joint measure of functional diversity, i.e., functional dispersion, provides comparable break points to the break points of individual traits. As trait similarity of the assemblage may enhance temporal stability (), we also calculated functional redundancy to assess stability of northern ecosystem using adiv package () and performed break point analysis. The break point analysis was performed using the library “segmented” (). All the analyses in this work were performed in R v4.1.3 ().