Multiplexing PCR allows the identification of within‐species genetic diversity in ancient eDNA

Sedimentary ancient DNA (sedaDNA) has rarely been used to obtain population‐level data due to either a lack of taxonomic resolution for the molecular method used, limitations in the reference material or inefficient methods. Here, we present the potential of multiplexing different PCR primers to retrieve population‐level genetic data from sedaDNA samples. Vaccinium uliginosum (Ericaceae) is a widespread species with a circumpolar distribution and three lineages in present‐day populations. We searched 18 plastid genomes for intraspecific variable regions and developed 61 primer sets to target these. Initial multiplex PCR testing resulted in a final set of 38 primer sets. These primer sets were used to analyse 20 lake sedaDNA samples (11,200 cal. yr BP to present) from five different localities in northern Norway, the Alps and the Polar Urals. All known V. uliginosum lineages in these regions and all primer sets could be recovered from the sedaDNA data. For each sample on average 28.1 primer sets, representing 34.15 sequence variants, were recovered. All sediment samples were dominated by a single lineage, except three Alpine samples which had co‐occurrence of two different lineages. Furthermore, lineage turnover was observed in the Alps and northern Norway, suggesting that present‐day phylogeographical studies may overlook past genetic patterns. Multiplexing primer is a promising tool for generating population‐level genetic information from sedaDNA. The relatively simple method, combined with high sensitivity, provides a scalable method which will allow researchers to track populations through time and space using environmental DNA.

An alternative source of aDNA is sedimentary ancient DNA (se-daDNA), DNA that can persist in cave, permafrost or lake sediments (Capo et al., 2021;Parducci et al., 2017;Pedersen et al., 2015) for extended periods of time (Kjaer et al., 2022).SedaDNA has primarily been used for the identification of taxa or even reconstructing palaeoecologies (Epp et al., 2012;Parducci et al., 2019;Rijal et al., 2021;Willerslev et al., 2014).However, a number of studies have explored the potential of sedaDNA for the retrieval of intraspecific variation through the usage of either shotgun metagenomic (Lammers et al., 2021;Meucci et al., 2021;Pedersen et al., 2021;Seersholm et al., 2016), qPCR melting curve assays (Nota et al., 2022) or hybridization capture techniques (Schulte et al., 2021;Zavala et al., 2021;Zhang et al., 2020).These methods, however, can be time intensive in terms of data generation and analysis, and in case of shotgun sequencing require a high concentration of the species of interest.
More scalable methods are thus required for large-scale reconstructions of past populations.
Amplicon-based methods such as DNA metabarcoding have proven to be scalable when working with sedaDNA, but as a singlebarcode approach, have been limited to species-level or higher identifications.Multiplexing PCR methods combine multiple primer sets into a single reaction to simultaneously obtain information from different groups of organisms, for example, mammals and plants (De Barba et al., 2014;Taberlet et al., 2018).However, different primer sets can also be included in a reaction that is designed to amplify intraspecific regions, thus obtaining population-level genomic information.So far, the usage of multiplexing PCRs to obtain populationlevel data is limited to contemporary DNA (Andres et al., 2021;De Barba et al., 2017;Skrbinšek et al., 2010), whereas multiplexing PCRs on ancient sediments have been limited to species-level identifications (Côté et al., 2016).Applying the multiplexing PCR method to sedaDNA comes with its own challenges given the damaged and highly fragmented nature of ancient DNA (Dabney et al., 2013;Orlando et al., 2021).
In this study, we investigate the potential of multiplexing PCRs for the retrieval of population-level genomic information from ancient sediments.We have selected the dwarf shrub Vaccinium uliginosum as our test species.V. uliginosum has well-known present-day intraspecific lineages (Alsos et al., 2005;Eidesen et al., 2007) and is frequently detected in sedaDNA studies (Alsos et al., 2022;Clarke et al., 2019), making it an ideal candidate for method development.
We designed multiple intraspecific primers for V. uliginosum and tested these on a range of sedaDNA samples to explore the palaeophylogeography of the species and evaluate the applicability of the method.

| Study species
The dwarf shrub V. uliginosum (Ericaceae) has a present-day circumpolar and circumboreal distribution but can also be found in more southern mountain ranges such as the Alps and Pyrenees (Hultén, 1986).The species is long-lived, bird-dispersed and insect-pollinated, though self-pollination has been observed (Jacquemart, 1996).Furthermore, the species has been reported in sedaDNA studies from Scandinavia (Alsos et al., 2022;Rijal et al., 2021), the Ural mountains (Clarke et al., 2019) and the Alps (Garcés-Pastor et al., 2022;van Vugt et al., 2022).A number of different subspecies or synonym species have been described, but these can be categorized into the following three lineages based on Amplified Fragment Length Polymorphism (AFLPs), ploidy and chloroplastic variation: Amphi-Atlantic, Arctic-Alpine and Beringian (Figure 1; Alsos et al., 2005;Eidesen et al., 2007).
The Amphi-Atlantic lineage is tetraploid and is widespread throughout Europe at lower elevations.Its distribution reaches from the northern Alps to Scandinavia, Iceland, southern Greenland and the northern Atlantic coast of North America (Alsos et al., 2005;Eidesen et al., 2007).
The Arctic-Alpine lineage is primarily diploid and can be split into two sublineages, namely the Arctic and Alpine sublineage.The Arctic sublineage has a circumpolar distribution and is found in northern Asia, Greenland, Svalbard and North America, with single populations in the Carpathians and northern Norway (Alsos et al., 2005).
The Alpine sublineage on the other hand has a more limited distribution and is found in the Alps and Pyrenees (Figure 1b).The split of the Arctic-Alpine lineage is most likely the result of isolation during the Quaternary cold periods, where the Alpine sublineage became isolated in southern Europe in relation to the Arctic sublineage (Alsos et al., 2005;Eidesen et al., 2007).
Finally, there is the Beringian lineage, which contains a mixture of di-, tetra-and hexaploids.This lineage has a geographical range from northern Japan and the Kamchatka Peninsula to Alaska, with smaller populations extending down the North American west coast (Alsos et al., 2005;Eidesen et al., 2007;Hirao et al., 2011).

| Vaccinium uliginosum variant discovery
The V. uliginosum variant discovery was based on the genome skim data for 18 modern individuals that were either part of the PhyloNorway project (n = 16; Alsos et al., 2020;Wang et al., 2021) or the PhyloAlps project (n = 2; Table S1).A reference chloroplast genome was assembled based on individual TROM_V_355013, an Alpine individual from the French Alps, following the methods described in Alsos et al., 2020.Intraspecific variable positions were identified by mapping each of the 18 individuals separately onto the assembled reference genome using bowtie2 v2.3.4.1 (Langmead & Salzberg, 2012) with default settings (Figure 2a).The resulting alignments were processed with SAMtools v1.7 (Li et al., 2009), in order to remove unaligned reads and duplicate reads with the view and markdup functions, respectively.Variants were called and filtered with BCFtools v1.9 (Li et al., 2009), using the mpileup, call and index functions, with the variants-only and multiallelic-caller options.The resulting variable positions were further filtered to only include those that had a mapping quality ≥30, coverage ≥20, and an alternative allele that was more abundant than the reference allele in case of haplotypic variation (Figure 2a).
The variable positions for the 18 individuals were combined into one set and processed to find candidate regions for primer design (Figure 2a).First, variants that were located within 100 bp of a contig end were discarded, as these variants could both be the result of incorrect mapping and limit the space available for primer design.
Second, variants located within 5 bp of an indel or homopolymer stretch of at least five bases were excluded to minimize the chance of misaligned reads.For the remaining variable positions, the haplotypes were called for each of the 18 individuals using the filtered bowtie2 alignments.Each haplotype was constructed using the bases with a base quality ≥13 for assemblies with a coverage ≥20, and alleles with a frequency ≥0.1 in case of haplotypic variation (Figure 2a).Individual variable positions were grouped into larger windows if two adjacent variants were within 50 bp of each other and treated as a single variable position onwards.To remove singleton haplotypes that carry limited population-level information, only positions that had at least two distinct haplotypes across the 18 individuals, with a minimum of two individuals per haplotype, were retained.
To ensure that each haplotype was specific to V. uliginosum, the variable positions were compared against all assembled chloroplast genomes in the PhyloNorway (n = 1969) and PhyloAlps (n = 3923) genome skim reference databases (Alsos et al., 2020;Wang et al., 2021).For each variable position, the haplotypes present, extended 20 bp up and downstream of the variant, were aligned against the reference databases with the blastn function from NCBI-BLAST+ v2.2.18 (Camacho et al., 2009).Any variable position that contained a haplotype that was not specific to V. uliginosum was removed (Figure 2a).

| Ancient sediment samples
Twenty ancient lake sediment samples were selected for V. uliginosum multiplex amplification.The sedaDNA samples originated from five different lake sediment cores that are known to contain V. uliginosum DNA based on previous trnL P6-loop metabarcoding of the sediments, and they cover the present-day extent of the different European V. uliginosum lineages (Data S1).Two lakes, Jøkelvatnet and Nordvivatnet, are located in northern Norway (Rijal et al., 2021) where the Amphi-Atlantic lineage is dominating today; two, Hopschusee and Krumschnabelsee, are from the Alps (Garcés-Pastor et al., n.d.) where the Alpine sublineage is most prevalent today, and one, Bolshoye Shchuchye, from the Polar Urals, Russia (Clarke et al., 2019), where the Arctic sublineage is dominant today.et al., 2013) for the Bolshoye Shchuchye, Jøkelvatnet and Nordvivatnet samples (Clarke et al., 2019;Rijal et al., 2021) and IntCal20 curves (Reimer et al., 2020).The age-depth models for Hopschusee and Krumschnabelsee were based on 25 and 12 radiocarbon dates, respectively, and followed the methods above (Garcés-Pastor et al., n.d.).
The material for Bolshoye Shchuchye, Jøkelvatnet and Hopschusee was previously extracted, along with three negative extraction controls (Garcés-Pastor et al., n.d.;Clarke et al., 2019;Rijal et al., 2021).The eight extracts used for Nordvivatnet and Krumschnabelsee were re-extracted for this study.All sampling and extractions were performed in a dedicated ancient DNA laboratory at the Arctic University Museum of Norway, UiT-The Arctic University of Norway, Tromsø, Norway.DNA was extracted from ~0.3 g sediment subsamples using a modified DNeasy PowerSoil kit (Qiagen, Germany).Modifications include a bead beating step and an initial lysis step with proteinase K, as described in (Alsos et al., 2020).
A negative extraction control was included for the re-extraction of the Nordvivatnet and Krumschnabelsee sediments.

| Multiplex PCR
Multiplex primers were developed for all variable positions that passed the filtering criteria.Primer sets were developed using a 100 bp flanking region around the variable positions and were designed such that they included both the V. uliginosum intraspecific variable positions and any positions required for V. uliginosum specificity.The primer sets were optimized to have a short amplicon length of an average of 79 bp (min 46 bp, max 116 bp; Table S2), to ensure amplification of highly fragmented ancient material (Orlando et al., 2021).Furthermore, the melting temperature was kept similar across all primer sets at an average of 55.7°C (min 55°C, max 57.5°C) to promote even primer annealing and amplification.In total, 61 primer sets were developed and tested (Data S1; Table S2), which resulted in a final set of 38 primers (Table S2).Each primer was modified by adding one of two 3 bp tags on the 5′ end to allow for pooling of up to four PCR products (Binladen et al., 2007).These tags had an edit distance of three to avoid misassignment based on sequencing errors (Table S3).The 38 primer sets were pooled together, with different concentrations per primer set to account for the different amplification efficiencies (Data S1; Table S2).S2) and 5 μL of DNA extract.Two negative PCR controls were included in the sample preparation.The following amplification protocol was used: enzyme activation step (2 min at 95°C), 45 PCR cycles of 30 s at 95°C, 90 s at 53°C and 30 s at 72°C, followed by a final elongation step of 10 min at 72°C.A annealing temperature of 2°C lower than the calculated 55°C temperature was used to increase the efficiency of the reaction (Taberlet et al., 2018).The amplification of the samples and controls was carried out using four multiplex PCR replicates.Positive amplification for a subset of the sample replicates was verified by electrophoresis on a 2% agarose gel.For each sample, an equal volume of PCR product from the four uniquely tagged replicates was pooled together and cleaned following (Voldstad et al., 2020).The pools were converted into DNA libraries with the Illumina TruSeq DNA PCR-Free protocol (Illumina Inc., CA, USA), with unique dual indexes and sequenced at 2 × 150 cycles on an Illumina NextSeq platform at the Genomics Support Centre Tromsø, UiT-The Arctic University of Norway.
For each library, the data were demultiplexed with the ngsfilter tool from OBITools (v1.2.12; Boyer et al., 2016) using default settings (Figure 2b).Identical reads were collapsed with obiuniq, and singleton sequences were removed with obigrep.PCR artefacts were identified with obiclean using a head: internal ratio of 0.05.Finally, the sequences were identified using ecotag and a curated reference database (Figure 2b).The reference database was created with ecoPCR v1.0.1 (Ficetola et al., 2010), using the multiplex primer sets and the reference genomes from the PhyloAlps (4604 specimens of 4437 taxa) and PhyloNorway (2051 specimens of 1899 taxa) projects (Alsos et al., 2020;Wang et al., 2021), allowing for a maximum of two mismatches between the primer sets and references.The identified data were exported with obitab for downstream analysis.
The data were further analysed with a custom Python script.Rare sequence detections, defined as those with less than three reads per PCR repeat, were removed from the dataset, as well as any sequences identified as "internal" by obiclean.For each sample, the four replicates were merged together, calculating the replication for each sequence in the process.Any sequence that was not present in at least two replicates was removed.The sequences that were identified as V. uliginosum were compared to the known haplotypes of the 18 reference individuals.For each sample, a reference individual, or individuals in case of identical haplotypes, was considered present if it was supported by at least one unique haplotype.The reads for each haplotype were split among the present reference individuals.Furthermore, for each sample, the V. uliginosum lineages present were determined as well.A lineage was considered present if there was at least one haplotype specific to the lineage detected in a sample (Figure 2b).

| Vaccinium uliginosum variant discovery
Phylogenetic analysis of the 18 Vaccinium uliginosum reference chloroplast genomes indicated that all four lineages and sublineages were covered.Seven individuals belonged to the Arctic sublineage, five to the Alpine sublineage, five to the Amphi-Atlantic lineage and finally one individual from the Bergingian lineage (Figure 1a).
A total of 1059 variable positions were found across the 18 Vaccinium uliginosum reference individuals.This was reduced to 704 variable positions after those close to contig edges, indels or homopolymers were removed.Grouping of adjacent variable positions into windows and removal of non-descriptive positions reduced the number to 292 variable positions.Out of these, 88 were fully specific to V. uliginosum and 81 were selected for multiplex primer design (Table S2).The remaining seven variable regions were excluded due to a high (>0.2) proportion of haplotypic variation among the individuals.

| Multiplex PCR performance
Multiplex-compatible primer sets were developed for 61 of the variable regions.Initial multiplex PCR tests indicated that 52 of these had both successful amplification and informative results (Data S1).
A final set of 38 primers that had the highest read counts in testing were used.Phylogenetic analysis of just the amplifiable regions revealed that they are capable of identifying the four Vaccinium uliginosum lineages, though five out of the seven Arctic, two out of the five Alpine and two out of the five Amphi-Atlantic reference individuals had identical haplotypes within their lineage and thus cannot be identified down to the exact reference individual by the primer sets (Figure S1).
We generated 98,947,000 paired-end reads for the sediment samples and controls.After merging and demultiplexing, we obtained 22,376,000 reads, out of which 19,050,000 reads remained after filtering and identification.Of the identified reads, 4,177,000 were identified as V. uliginosum and the remainder as off-target bycatch (Table S4).The most common off-target taxon was identified as Vaccinium, containing sequences that could either be identified as V. myrtillus or V. vitis-idaea.Vaccinium uliginosum was identified in all 20 sediment samples and was absent from the controls.The controls instead contained sequences identified to higher taxonomic levels (Spermatophyta, asterids, Pentapetalae) and common food contaminants (Solanum) (Table S4).Each sediment sample contained on average 28.1 (SD 11.07) and 34.15 (SD 15.4) on-target primer sets and variant sequences, respectively (Table S4).Seven sediment samples had on-target amplification for all 38 primer sets, while the lowest number of on-target primer sets was six, for a sediment sample from Nordvivatnet (Table S4).
All 38 primer sets could be observed in the sediment samples, where on average each primer could successfully amplify V. uliginosum in 14.79 (SD 3.11) sediment samples (Table S5).Furthermore, three primer sets had successful on-target amplification in all sediment samples (Table S5).
After filtering and identifying the reads, on average, each primer had 501,000 reads assigned to it, of which 110,000 reads were identified as V. uliginosum.The most abundant primer set was Vu_62, with 7,228,000 reads, though only 132,000 reads were on-target.
The most abundant on-target primer set was Vu_16, with 699,000 V.
In total, 89 different on-target variant sequences were obtained for the samples, with an average of 2.34 (SD 0.58) sequences per primer (Table S5).The average on-target replication for these sequences was 3.51 (SD 0.75), with 86 out of the 89 variant sequences being fully replicated for at least one sample (Table S7).
All four V. uliginosum lineage and sublineages were represented in our data.The Amphi-Atlantic lineage was the most common lineage with 15 detections, while the Beringian lineage was the rarest, with one detection.The Arctic and Alpine sublineages had 5 and 10 detections, respectively.On average, we detected 1.55 lineages per sample, with a minimum of one and a maximum of three lineages per sample (Tables S8-S10).All reference individuals within the lineages were present for at least one sediment sample and were represented by at least one fully replicated primer (Table S8).
When subsampling the data, the main lineage, as defined by containing more than >1% of the unique reads, is present in most samples from 500 reads per repeat onwards.Most additional lineages are also present from 500 reads per repeat onwards, though some lineages, which are only represented by one repeat, are never detected at stricter filtering regimes (Figure S2).The number of identified references on the other hand shows a declining trend for all filtering methods as the number of reads per subsample increases (Figure S3).

| Vaccinium uliginosum sedaDNA results per region
The Bolshoye Shchuchye sediment samples primarily contained reads assigned to the Arctic V. uliginosum sublineage.The two oldest samples exclusively contained sequences and repeats assigned to the Arctic sublineage, while the two youngest samples contained material assigned to both the Arctic sublineage and the Amphi-Atlantic lineage, where the Arctic sublineage was the more abundant lineage, with 99.9% and 92.2% of the unique reads and repeats, respectively (Table 1; Figure 3; Tables S8 and S9).
Both northern Norwegian lakes mainly contained material assigned to the Amphi-Atlantic lineage, the most dominant lineage in the region today.For Nordvivatnet, the Amphi-Atlantic lineage was the only lineage in the oldest three samples, but in the most recent sample 10% and 16% of the reads and repeats, respectively, were assigned to the Arctic sublineage.Jøkelvatnet, again mainly contained material identified to the Amphi-Atlantic lineage, but in addition contained material identified as the Alpine sublineage in all samples.
The Alpine sublineage, however, was present, on average, in only 0.6% and 7% of the reads and repeats, respectively, across the four samples (Table 1; Figure 3; Tables S8 and S9).
Similar to the Norwegian sites, the Alpine sites contained a mixture of V. uliginosum material.The Alpine sublineage was the main lineage in Hopschusee.It was the only lineage present from 10,500-3100 cal.yr BP and the dominant lineage in the youngest sample.The youngest sample also contained material assigned to the Amphi-Atlantic and Beringian lineages, at 32.9% and 0.7% of the reads, respectively.The sediment samples from Krumschnabelsee only cover the period 2900 cal.yr BP to present, and it contained primarily Amphi-Atlantic material, which was the only lineage in the oldest and youngest samples from this site.The two samples from 1800 to 600 cal.yr BP had an even split between the Alpine sublineage and Amphi-Atlantic lineage for both the reads and repeats (Table 1; Figure 3; Tables S8 and S9).

| Multiplex PCR results strengths and weaknesses
Our V. uliginosum results indicate that we can reliably amplify multiple informative target regions per sedaDNA sample.Despite the low coverage of some of our samples, we were still able to determine the main V. uliginosum lineages present.Narrowing down the exact reference genomes present is more challenging due to haplotype sharing between a number of reference individuals (Figure S1).As a result, we conclude that the method presented here is most robust for identification of the V. uliginosum lineages and that identifications down to an individual reference genome should be interpreted with caution.
The samples with the lowest number of on-target primer sets tended to have a lower on-target replication and relatively few reads assigned to V. uliginosum (Tables S4 and S7).The low V. uliginosum read coverage observed is mainly the result of off-target amplification.A number of primer sets, most notably Vu_62, are able to amplify taxa outside of V. uliginosum and Ericaceae.The amplification of off-target material will both reduce the relative amount of V. uliginosum material in these primer sets and reduce the overall amount of V. uliginosum within the multiplex reaction.That combined with the already low proportion of V. uliginosum template material for some samples, as estimated from previous trnL P6-loop metabarcoding results, can result in a limited amount of usable data (Table S4).Removal of less specific primer sets should increase the amount of assigned V. uliginosum sequences and result in improved identifications.
This study utilized 38 primer sets to identify the different V. uliginosum lineages.Given the experimental nature of this study, redundancy was built into this set of primers, to ensure that the main lineages could be identified even in the event of non-amplification of the target regions.For follow-up studies, this set of primers can be reduced to improve the read coverage per primer at equal sequencing depths or alternatively open up the possibility to include different primer sets targeting different regions or species.
Finally, the primer sets designed for this study primarily differentiate between the Amphi-Atlantic, Arctic and Alpine lineages, given the focus on European sedaDNA samples.Although a Beringian V. uliginosum sample was included in the set used for TA B L E 1 The proportion of unique reads assigned to each V. uliginosum lineage for the 20 lake seda DNA samples.primer design and as a reference, the Beringian region contains substantial variation that likely is not fully captured by a single individual (Alsos et al., 2005;Alsos, 2003;Brochmann et al., 2004;Eidesen et al., 2007).Applying the multiplex method and the primer sets developed here to sedaDNA samples from Beringia would require sequencing of additional reference individuals and inspecting the amplifiable variation present prior to analysis of sedaDNA samples.

| Vaccinium uliginosum phylogeography
Due to the proximity of the potential V. uliginosum Arctic sublineage refugia east of the Scandinavian ice sheet (Hughes et al., 2016) to northern Norway, we assumed that the first V.  1 and Table S9.The oldest sample from Krumschnabelsee ( 5) is from 2900 cal.yr BP and is thus empty for the periods >9 ka and 8.9-5 ka.The chart in period <1.9 ka for Krumschnabelsee is an average of three samples.
Cape (Alsos et al., 2005;Eidesen et al., 2007) is an ancient or more recent introduction in Norway, more sites and samples need to be analysed.
The Alps provide a more dynamic story.Although this dataset already provided some novel insights regarding V. uliginosum lineage migration, more work needs to be done to provide a complete overview of the species post-glacial migration.
This study is limited in its spatial and temporal resolution.Samples were selected to cover a geographically wide area to increase the chance of observing the three European V. uliginosum lineages and allow for the evaluation of the multiplex PCR's performance.For a full reconstruction of the species migration, more sedaDNA sampling locations are required along the migration paths in Scandinavia, and potential refuge locations such as Western Europe, the Urals and Pyrenees.Furthermore, the temporal resolution needs to be improved so that both the timing of migration events can be narrowed down and increase the ability to observe short-term population turnover.

| Future applications of multiplexing PCRs
The success of the multiplexing method for sedaDNA has implications for other taxa and studies.First, a set of multiplex PCR primers can be used to improve the taxonomic resolution of a group of interest (Côté et al., 2016) or resolve cryptic species identifications (Brosseau et al., 2019), compared to single-barcode methods.Second, combining multiple metabarcoding primer sets into a single multiplex PCR can greatly increase the taxonomic information obtained from a single amplification run (De Barba et al., 2014;Schuette et al., 2022;Weber et al., 2023).A traditional metabarcoding approach can achieve similar results by running multiple distinct metabarcoding runs using different primer sets (Epp et al., 2012;Garcés-Pastor et al., 2022;Willerslev et al., 2014); however, this comes at the cost of additional laboratory expenses and, especially in the case of sedaDNA, usage of available DNA extracts.
The multiplex PCR method on the other hand requires some careful planning.First, the method, similar to shotgun metagenomics and hybridization capture methods, requires the presence of an extensive reference database.This database is needed to both identify the intraspecific regions of interest and avoid false positives caused by regions shared among closely related species.Furthermore, the identification of extinct lineages, regardless of sedaDNA method, can be problematic if no appropriate reference material is available in the reference databases.Second, to ensure sufficient template material for amplification, multicopy regions are preferred, such as plastid or ribosomal DNA.In addition, taxa with a low biomass might not produce the template DNA required for reliable amplification across all primer sets.Third primer sets must be carefully designed to ensure that they have comparable annealing temperatures and do not form dimers. Fourth, primer sets need to be specific to the target taxa, which requires access to complete reference libraries for both design and identification of the amplicons.Finally, different primer sets and taxa will differ in available template DNA, which can lead to uneven amplicon concentrations and read counts.This can be solved by adjusting the primer concentrations to match the available template, but this can take time to fully optimize (Taberlet et al., 2018).
Furthermore, since multiplexing does not retain deamination profiles that are commonly used to authenticate ancient DNA (Dabney et al., 2013;Orlando et al., 2021), sufficient negative controls need to be incorporated to avoid detections based on contamination (Pedersen et al., 2015;Taberlet et al., 2018).
The multiplex PCR method presented here for generating population-level genomic data from sedaDNA samples has some advantages compared to the more commonly used shotgun metagenomic and hybridization capture methods.The main advantage is the scalability and time effectiveness of the method.A multiplex reaction, similar to a metabarcode reaction, is relatively simple to set up compared to more complex methods such as hybridization capture (Schulte et al., 2021).Furthermore, the multiplex method requires a lower sequencing depth compared to shotgun metagenomics (Kjaer et al., 2022;Wang et al., 2021), which results in lower sequencing costs.The bioinformatic analysis of the multiplex data is also relatively simple compared to shotgun metagenomics and hybridization capture, as the data generated can be processed to some extent by well-established metabarcoding pipelines, which reduces the analysis time considerably.All these factors combined allow the multiplex PCR method to generate population-level genomic data for a large number of samples and thus enable the potential to track populations through time and space, opening up the field of palaeo-phylogeography.

F
I G U R E 1 (a) Chloroplast maximum likelihood tree for the 18 V. uliginosum reference individuals.Node values indicate bootstrap support.(b) Map with the location of the European Arctic, Alpine and Amphi-Atlantic reference individuals indicated by the coloured boxes.Coloured circles indicate the location of previously genotyped individuals by Alsos et al. (2005).The black dots indicate the location of the five lake sedaDNA cores, (1) Bolshoye Shchuchye, (2) Nordvivatnet, (3) Jøkelvatnet, (4) Hopschusee and (5) Krumschnabelsee.All sites are close to present-day meeting points between the lineages (Figure1b).
Each site was represented by four sedaDNA samples, one sample for the oldest known V. uliginosum detection, one for the most recent detection, and the remaining two samples spaced in between.The lake sediment cores were previously dated by radiocarbon dating of plant macrofossils using Accelerator Mass Spectrometry (AMS) at the Poznań Radiocarbon Laboratory.Bayesian age-depth models were constructed with "Bacon" v2.3.4 (Blaauw & Christen 2011) using the calibrated radiocarbon dates based on the IntCal13 curves (Reimer F I G U R E 2 (a) Workflow utilized for the variant discovery.Starting with aligning N genomes or genome skims to a shared reference genome.The resulting variants are combined, and unsuitable positions (rare, near contig edges or indels) are removed.Each variant is compared to other reference taxa, and only those that are distinct are retained for primer design.(b) Workflow for the multiplex data analysis.Each replicate for a sample is demultiplexed separately, and for each primer present the variants are identified.After identification, the replicate data for a sample are combined.Off-target identifications and non-replicated variants are removed.Variants are retained if they are supported by at least one unique haplotype.

All
PCRs were prepared in a dedicated ancient DNA lab at the Arctic University Museum of Norway, UiT-The Arctic University of Norway, Tromsø, Norway.Each amplification was carried out in a 50 μL final volume, containing 25 μL of Platinum Multiplex PCR master mix (Life Technologies, Carlsbad, CA, USA), 15 μL of molecular grade water, 5 μL of the multiplex primers with a final concentration of 1.735 μM (Table uliginosum colonization of this region would originate from the Arctic sublineage, with later replacement of the boreal Amphi-Atlantic lineage observed in the region today.Our results, however, indicate that the Amphi-Atlantic lineage was first to arrive in northern Norway, as was detected in Jøkelvatnet and Nordvivatnet.In the case of the more western Jøkelvatnet, this was paired with a background presence of the Alpine sublineage, which could be the result of unknown Arctic variation not present in our reference set, suggesting potential longterm presence of the Arctic-Alpine lineage.On the other hand, at Nordvivatnet only the most recent sample contains a small component of Arctic sublineage reads and repeats.Thus, to determine if the presence of a modern Arctic V. uliginosum population at North F I G U R E 3 Multiplex results for the five lake sedaDNA cores: (1) Bolshoye Shchuchye, (2) Nordvivatnet, (3) Jøkelvatnet, (4) Hopschusee and (5) Krumschnabelsee.Each lake is represented by a pie chart indicating the proportion of unique reads supporting each lineage.Raw values are provided in Table Hopschusee was first colonized by the Alpine sublineage, which remained the only lineage present until the arrival of the Amphi-Atlantic lineage in the youngest sample, where both lineages are now present in roughly even amounts.The introduction of the Amphi-Atlantic lineage could represent the migration of the more boreal lineage to higher altitudes.The Beringian lineage is also present in the youngest sample, though at a low abundance.This detection is based on just two haplotype sequences and could be the result of unknown Alpine or Amphi-Atlantic variation that is overlapping with the variation found within the Beringian lineage.Vaccinium uliginosum is a relatively recent detection in Krumschnabelsee, where it was first detected at 3000 cal.yr BP although the core goes back to 8600 cal.yr BP.The Amphi-Atlantic lineage was the first to appear but was joined by the Alpine sublineage between 1800 and 600 cal.yr BP.The most recent sediment sample, representing the present, however, indicates that the Amphi-Atlantic lineage is currently the only lineage remaining.The Bolshoye Shchuchye samples, representing the Polar Urals, are the most stable in terms of V. uliginosum lineages present.The Arctic sublineage is the main lineage present, where only the two most recent samples have a minor representation of the Amphi-Atlantic lineages.The present-day extent of the Amphi-Atlantic lineages reaches the Urals (Eidesen et al., 2007), so the more recent detection of this lineage in the Polar Urals samples could indicate a recent background presence in the region.