Extracting and authentication of aDNA from ancient mastics

We performed DNA extraction on eight mastic samples from Huseby Klev (Fig. 1, Supplementary Note 1, Supplementary Data 3). First, we used Yang-Urea extraction following the modifications employed by Svensson et al.29 of the published protocol30. We named the successful samples ble004, ble007, and ble008. For the ble008 sample we tested one more extraction method using an extraction kit designed to process samples with high-inhibitor content (QIAamp PowerFecal DNA Kit, Qiagen), with slight modification to the protocol provided with the kit (described in the Methods chapter). Blunt-end repaired libraries were built on the concentrated DNA31, four of which were UDG treated (for ble004 sample). The libraries were amplified using polymerase chain reaction (PCR) and shotgun sequenced at the Science for Life laboratories (SciLife) in Stockholm, on Illumina Hiseq X platform. Initial tests showed that the individual libraries from Yang-Urea extracts contained authentic ancient human DNA, ranging from 2 to 8% (Supplementary Table 1). The library built on the PowerFecal DNA Kit extract contained 23% endogenous DNA. The tenfold increase in endogenous DNA content makes the PowerFecal DNA extraction kit a valid approach to process mastic samples. By repeating the process with successful extracts, we obtained genome-wide data from three of the mastic pieces, ranging from 0.1× to 0.49× coverage (Table 1). We merged individual libraries using SAMtools32 and calculated library statistics (before and after PMD filtering, Supplementary Tables 2 and 3) and produced damage plots (Supplementary Fig. 1).

Fig. 1 The studied material and its origin. a One of the chewing gums from Huseby Klev, (Fynd 2037), with two plastelina casts for each side. The cast to the left captures several teeth imprints from the left side of the maxilla, the one to the right is of the corresponding teeth from the mandible. The presence of the second molar and analysis of tooth wear suggest that the individual, who left these imprints was in the early teens (12–14 years old)11. Scale bar: 50 mm (photo by Verner Alexandersen). b The location of the sites, genomes from which were used in this study. 1—LaBrana; 2—Loschbour; 3—Bichon; 4—Hummervikholmen; 5—Huseby Klev; 6—Steigen; 7—Motala; 8—Stora Bjers; 9—Stora Förvar; 10—Yuzhnyy Oleni Ostrov; 11—Samara Full size image

Table 1 The library properties for ble samples after processing and pmd filtering Full size table

Contamination and pmd filtering

We estimated mitochondrial contamination rates using near-private consensus alleles as described by Green et al.33 To exclude the effects of sequencing errors, we used bases that have a base-call and mapping quality score of more than Q30. Also, we filtered the positions where we detected transition patterns to compensate the post-mortem damage (Supplementary Table 4). We used the PMDtools software to filter out the possible contaminant sequences34. This tool compares each aDNA sequence to its modern counterpart to calculate deamination specific nucleotide transitions and assigns a pmd score. This score is used to evaluate the authenticity of the sequence. We set the pmd threshold to 0 and removed contaminant DNA sequences below this threshold (Supplementary Figs. 2–6). After removing the potential contaminants from the merged libraries, we re-calculated the library statistics, deamination patterns, and MT contamination estimates to analyse the authenticity of the dataset. Deamination patterns (Supplementary Fig. 1) and MT contamination rates (Table 1, Supplementary Tables 4 and 5) present a strong case that the aligned data is authentic and represents the individuals that chewed the ancient mastic (Supplementary Tables 3 and 5).

Relationships between ancient individuals

We used READ35, to explore kinship between the individuals. READ compares the non-overlapping 1 Mb segments in the genome and calculates the nonidentical allele ratio between the samples (P0). Lower P0 values mean more shared chromosomal segments. We confirmed that none of the genomes are identical to each other (Supplementary Table 6). We also found that ble004 and ble007 have a possible second degree relationship. However, it should be noted that using three individuals for analysis is not recommended for this tool, and we refrain from using this result in further discussion. In summary, we can confirm that we sequenced DNA from three distinct individuals.

Mitochondrial DNA

We used samtools mpileup command to create mitochondrial consensus sequences using the nucleotides that have a base-call and mapping quality score of more than Q30. We assigned mitochondrial haplogroups with HaploFind36 and HaploGrep 237 (Table 1). We reviewed the results with PhyloTree (build 17)38. Mitochondrial genomes from all three individuals belong to the U5a2d haplogroup. ble008 was assigned to U5a2 by HaploGrep 2, but the same sequence is assigned to U5a2d haplogroup by HaploFind and we accept this result (Supplementary Excel Table 1). The mitochondrial U5a2d haplogroup is consistent with earlier published results for ancient individuals from Scandinavia, U5a being the most common within SHG. Of the 16 Mesolithic individuals from Scandinavia published prior to our study, seven belong to the U5a haplogroup, and nine share the U2 and U4 haplogroups5,6 (Table 1). The library properties for ble samples after processing and pmd filtering. Statistic for ble004 is for merged damage repair and blunt-end repair libraries, ble007 is for blunt-end-repair merged libraries, ble008 is for merged libraries made from extracts produced with different protocols (Supplementary Table 1, Methods)

Demographic history and population genomics

To explore the demographic history of the ancient individuals, we curated a set of ancient genomes (Supplementary Table 7) and compared with the publicly available Human Origins SNP reference set39,40. This set contains 594,924 single nucleotide polymorphisms (SNPs) from 2404 modern individuals and 203 different worldwide populations. We coded deamination patterns as missing data to compensate for possible biases caused by deamination patterns. We used principal component analysis to acquire an overview of the affinity of the ble individuals with selected ancient and modern populations41(Fig. 2). We merged ancient individuals with the Human Origins reference set, coded nucleotide transitions as missing data, and used Procrustes transformation to project ancient individuals on the principal component space (Supplementary Fig. 7, Supplementary Table 7). The ancient individuals are divided into three groups. The WHG, include individuals from the sites of LaBrana, Loschbour and Bichon; the EHG, include Yuzhnyy Oleni Ostrov and Samara; and the SHG, include Norwegian Steigen and Hummervikholmen and Swedish Motala, Stora Bjers and Stora Förvar. The projected ancient individuals (WHG, EHG and SHG cluster) show close affinity to modern day North, East and Western European populations, and BLE individuals from Huseby Klev (the earliest among the SHG) cluster with the ancient genomes originating from Scandinavia4. Our three samples are located between the two Hummervikholmen individuals (Norway), and the Stora Förvar SF9/12 (Sweden), all three with dates earlier than c. 9000 calBP. By reproducing EHG and WHG populations in this plot6, we confirm the close affinity of ancient individuals from Scandinavia to WHG and EHG (Fig. 2). To evaluate the relationship between the Huseby Klev individuals and other Mesolithic Scandinavians, we examined the relative shared allele frequencies and estimated shared drift among populations via performing F3 and F4 tests. These tests propose formal statistical frameworks to study the patterns of allele frequency correlation across populations40. We first tested the relative allele sharing of BLE individuals between EHG and WHG groups. Results show a high contribution of WHG ancestry to BLE individuals (the highest contribution observed in ble008 individual F4: 0.022, Z-score 1.8, Fig. 3). We compared the contribution between ancestry of EHG to SHG and WHG to SHG for BLE individuals. Results show that all tested individuals have relative high-allele sharing with the SHG group, with ble008 individual showing the significant value (F4: 0.03, Z-score: 2.95). We divided the SHG group into two groups: SHGa and SHGb, ancient individuals found in contemporary Norway (Steigen and Hummervikholmen) and Sweden (Motala, Stora Bjers and Stora Förvar), respectively. We based this on both the geographical distribution and the previous studies demonstrating the close relation of SHGa to EHG group and SHGb to WHG group6. It should be noted that the Z-score values for our samples are low for ble004 and ble007, due to low amount of SNPs for these individuals. Moreover, the F3 tests also support the affinity of BLE individuals to SHG group. To further explore the demography within the SHG group, we compared the ancestry of BLE individuals within SHGa and SHGb groups. Both F3 and F4 tests suggest a trend to a high relative shared drift between BLE individuals and the SHGb group (Supplementary Excel Table 2, for F4 and Supplementary Excel Table 3 for F3 tests); however, a symmetrical relationship between BLE and both SHGb and SHGa cannot be ruled out, since only one of the individuals has shown considerable deviation towards any of the groups.

Fig. 2 Principal component analysis of the Huseby Klev individuals within the diversity of Mesolithic individuals from Europe. The magnified section incaptures BLE individuals’ relation to Western hunter–gatherer (WHG), Eastern hunter–gatherer (EHG) and Scandinavian hunter–gatherer (SHG) individuals (Supplementary Table 7) Full size image

Fig. 3 Results of relative allele sharing (F4) test between Huseby Klev individuals and ancient population groups (the triangle marks the significant result in deviation from zero). WHG, Western hunter–gatherers; EHG, Eastern hunter–gatherers; SHG, Scandinavian hunter–gatherers Full size image

To further explore our findings from the D and F4 tests, we used the model-based clustering algorithm admixture42, which is helpful when exploring the genetic components in a given dataset. The essence of this algorithm is to calculate the admixture proportions as a parameter of a model. The results show that starting from K15, while there is a similarity in genetic elements between the BLE individuals and SHGa, we do not see the red and blue/green component appearing at the same time. Since our best sequenced individual is ble008, which lacks both red and green component, we interpret that the BLE individuals and SHGb population have more genetic similarity (Fig. 4, Supplementary Fig. PDF 1, Fig. PDF 2).

Fig. 4 Admixture analysis showing the major mode for K = 15. The figure represents 11 runs out of 20 replicates (Greedy algorithm ran with the Jaccard distance and a 0.97 similarity threshold). WHG, Western hunter–gatherers; EHG, Eastern hunter–gatherers; SHG, Scandinavian hunter–gatherers Full size image

Lithics results

The technological analysis conducted for the purpose of this study, shows that the lithic artefacts from the deep pit display clear affiliation with the eastern pressure blade technology, as documented from a large number of sites in northern and western Scandinavia, eastern Fennoscandia and the East European Plain8,9,10,43. No artefacts diagnostic to the preceding Early Mesolithic blade technology, that would indicate chronological or technological mixing, were observed (Supplementary Note 1). Based on the composition of lithic artefact types, the site appears to represent a production site to which lithic raw materials, in their more or less unworked condition, were transported, and where initial core preparation and exploitation was performed.

In addtion, some standardised blade production and re-tooling was performed on-site, visible in the presence of discarded tools (a relatively low number) and regular blade blanks. Blades were produced by the same overall concept: serial production from single-platform, sub-conical and conical cores with faceted and smooth platforms. No complete regular blade cores are present, but fragments of conical cores with visible scars deriving from the detachment of very regular thin blades, along with core rejuvenation flakes with small-flake faceting and a platform to front angle close to 90°, suggest that the eastern pressure blade concept was employed (Supplementary Note 1, Supplementary Fig. 14). Although the majority of the studied blades display features found in blades produced by direct and indirect percussion techniques, a selection of blades display diagnostic characteristics of the pressure technique, and the variation in knapping techniques is best explained as related to the different stages of the production process (Supplementary Note 1). Morphometric analysis shows the production of a consistent range of blade blanks, which in turn allowed the production of standardised tools, such as barbed points (hulling-type), slender lanceolate microliths, as well as blades with lateral retouch on one or both edges. The last mentioned were probably used as inserts in composite slotted tools, to which the inserts were attached using mastic made of birch bark pitch44 (Fig. 5). A bone point with remains of pitch retrieved from the deep pit shows that birch bark pitch mastic was part of tool production at the site, while fragments of slotted points, contextually dated to the same period as the finds from the deep pit, were found nearby45.