Sample preparation and DNA extraction

We sampled c. 250 mg from the specimen for DNA analysis. Briefly, the sample was washed in 5% bleach solution to remove any surface contamination, rinsed in molecular biology grade water and left to dry. We tested three different extraction methods using between 20–50 mg of starting material: For method (1), 1 ml of lysis buffer containing 0.45 M EDTA (pH 8.0) and 0.25 mg/ml Proteinase K was added to the sample and left to incubate on a rotor at 56 °C. After 12 h the supernatant was removed and concentrated down to ~150 µl using Amicon Ultra centrifugal filters (MWCO 30 kDa), mixed 1:10 with a PB-based binding buffer51, and purified using MinElute columns, eluting in 30 µl EB. For method (2) the sample was digested and purified as above, but with the addition of a phenol-chloroform clean-up step. Briefly, 1 ml phenol (pH 8.0) was added to the lysis mix, followed by 1 ml chloroform:isoamyl alcohol. The supernatant was concentrated and purified, as described above. For method (3) the sample was dissolved in 1 ml chloroform:isoamylalcohol. The dissolved sample was then resuspended in 1 ml molecular grade water and purified as described above. DNA extracts prepared using a Proteinase K-based lysis buffer followed by a phenol-chloroform based purification step produced the best results in terms of the endogenous human DNA content (see Supplementary Table 1); however, following metagenomic profiling the extracts were found to be contaminated with Delftia spp., a known laboratory contaminant52. The contaminated libraries were excluded from metagenomic profiling.

Negative controls

We included no template controls (NTC) during the DNA extraction and library preparation steps. The NTCs prepared with the additional phenol-chloroform step were also found to be contaminated with Delftia spp., suggesting that the contaminants were introduced during this step. In addition, we included two soil samples from the site, weighing c. 2 g each, as negative controls. DNA was extracted as described above using 3 ml EDTA-based lysis buffer followed by 9 ml 25:24:1 phenol:chloroform:isoamyl alcohol mixture to account for the larger amount of starting material. The sequencing results are reported in Supplementary Table 1.

Library preparation and sequencing

16 µl of each DNA extract were built into double-stranded libraries using a recently published protocol that was specifically designed for ancient DNA53. One extraction NTC was included, as well as a single library NTC. 10 µl of each library were amplified in 50 µl reactions for between 15 and 28 cycles, using a dual indexing approach54. The optimal number of PCR cycles was determined by qPCR (MxPro 3000, Agilent Technologies). The amplified libraries were purified using SPRI-beads and quantified on a 2200 TapeStation (Agilent Technologies) using High Sensitivity tapes. The amplified and indexed libraries were then pooled in equimolar amounts and sequenced on 1/8 of a lane of an Illumina HiSeq 2500 run in SR mode. Following initial screening, additional reads were obtained by pooling libraries #2, #3, and #4 in molar fractions of 0.2, 0.4, and 0.4, respectively and sequencing them on one full lane of an Illumina HiSeq 2500 run in SR mode.

Data processing

Base calling was performed using Illumina’s bcl2fastq2 conversion software v2.20.0. Only sequences with correct indexes were retained. FastQ files were processed using PALEOMIX v1.2.1255. Adapters and low quality reads (Q < 20) were removed using AdapterRemoval v2.2.056, only retaining reads >25 bp. Trimmed and filtered reads were then mapped to hg19 (build 37.1) using BWA57 with seed disabled to allow for better sensitivity58, as well as filtering out unmapped reads. Only reads with a mapping quality ≥30 were kept and PCR duplicates were removed. MapDamage 2.0.959 was used to evaluate the authenticity of the retained reads as part of the PALEOMIX pipeline55, using a subsample of 100k reads per sample (Supplementary Fig. 6). For the population genomic analyses, we merged the ancient sample with individuals from the Human Origin dataset11 and >100 previously published ancient genomes (Supplementary Data 1). At each SNP in the Human Origin dataset, we sampled the allele with more reads in the ancient sample, resolving ties randomly, resulting in a pseudohaploid ancient sample.

MtDNA analysis and contamination estimates

We used Schmutzi60 to determine the endogenous consensus mtDNA sequence and to estimate present-day human contamination. Reads were mapped to the Cambridge reference sequence (rCRS) and filtered for MAPQ ≥ 30. Haploid variants were called using the endoCaller program implemented in Schmutzi60 and only variants with a posterior probability exceeding 50 on the PHRED scale (probability of error: 1/100,000) were retained. We then used Haplogrep v2.261 to determine the mtDNA haplogroup, specifying PhyloTree (build 17) as the reference phylogeny62. Contamination estimates were obtained using Schmutzi’s mtCont program and a database of putative modern contaminant mitochondrial DNA sequences.

Genotype imputation

We used ANGSD63 to compute genotype likelihoods in 5 Mb windows around 43 SNPs associated with skin, eye, and hair colour10 and lactase persistence into adulthood (Supplementary Data 2). Missing genotypes were imputed using impute264 and the pre-phased 1000 Genome reference panel65, provided as part of the impute2 reference datasets. We used multiple posterior probability thresholds, ranging from 0.95 to 0.50, to filter the imputed genotypes. The imputed genotypes were uploaded to the HIrisPlex-S website10 to obtain the predicted outcomes for the pigmentation phenotypes (Supplementary Data 3).

Principal component analysis

Principal component analysis was performed using smartPCA66 by projecting the ancient individuals onto a reference panel including >1000 present-day Eurasian individuals from the HO dataset11 using the option lsq project. Prior to performing the PCA the data set was filtered for a minimum allele frequency of at least 5% and a missingness per marker of at most 50%. To mitigate the effect of linkage disequilibrium, the data were pruned in a 50-SNP sliding window, advanced by 10 SNPs, and removing sites with an R2 larger than 0, resulting in a final data set of 593,102 SNPs.

D- and f-statistics

D- and f-statistics were computed using AdmixTools67. To estimate the amount of shared drift between the Syltholm genome and WHG versus EHG and Neolithic farmers, respectively, we computed two sets of f 4 -statistics of the form f 4 (Yoruba, X; EHG/Barcın, WHG) where “X” stands for the test sample. Standard errors were calculated using a weighted block jackknife. To confirm the absence of EHG and Neolithic farmer gene flow in the Syltholm genome and to contrast this result with those obtained for other Mesolithic and Neolithic individuals from Scandinavia, we computed two sets of D-statistics of the form D(Yoruba, EHG/Barcın; X, WHG) testing whether “X” forms a clade to the exclusion of EHG and Neolithic farmers (represented by Barcın), respectively.

qpAdm

Admixture proportions were modeled using qpAdm12, specifying Mesolithic Western European hunter-gatherers (WHG), Eastern hunter-gatherers (EHG) and early Neolithic Anatolian farmers (Barcın), as possible ancestral source populations. We present the model with the lowest number of source populations that fits the data, as well as the model with all three admixture components (see Supplementary Table 6). When estimating the admixture proportions for WHGs and EHGs, the test sample was excluded from their respective reference populations.

MetaPhlan

We used MetaPhlan213 to create a metagenomic profile based on the non-human reads (Supplementary Data 4). The reads were first aligned to the MetaPhlan2 database13 using Bowtie2 v2.2.9 aligner68. PCR duplicates were removed using PALEOMIX filteruniquebam58. For cross-tissue comparisons 689 human microbiome profiles published in the Human Microbiome Project Consortium14 were initially used, comprising samples from the mouth (N = 382), skin (N = 26), gastrointestinal tract (N = 138), urogenital tract (N = 56), airways and nose (N = 87). The oral HMP samples consist of attached/keratinised gingiva (N = 6), buccal mucosa (N = 107), palatine tonsils (N = 6), tongue dorsum (N = 128), throat (N = 7), supragingival plaque (N = 118), and subgingival plaque (N = 7). Pairwise ecological distances among the profiles were computed at genus and species level using taxon relative abundances and the vegdist function from the vegan package in R69. These were used for principal coordinate analysis (PCoA) of Bray–Curtis distances in R using the pcoa function included in the APE package70. Subsequently, we calculated the average relative abundance of each genus for each of the body sites present in the Human Microbiome Project and visualised the abundance of microbial orders of our sample and the HMP body sites with MEGAN615.

MALT

To further characterise the metagenomic reads we performed microbial species identification using MALT v. 0.4.1 (Megan ALignment Tool)16, a rapid sequence-alignment tool specifically designed for the analysis of metagenomic data. All complete bacterial (n = 12,426) and viral (n = 8094) genomes were downloaded from NCBI RefSeq on 13 November 2018, and all complete archaeal (n = 280) genomes were downloaded from NCBI RefSeq on 17 November 2018 to create a custom database. In an effort to exclude genomes that may consist of composite sequences from multiple organisms, the following entries were excluded:

GCF_000922395.1 uncultured crAssphage

GCF_000954235.1 uncultured phage WW-nAnB

GCF_000146025.2 uncultured Termite group 1 bacterium phylotype Rs-D17

The final MALT reference database contained 33,223 genomes and was created using default parameters in malt-build (v. 0.4.1). The sequencing data for the ancient pitch sample, two soil control samples and associated extraction and library blanks were de-enriched for human reads by mapping to the human genome (hg19) using BWA aln and excluding all mapping reads. Duplicates were removed with seqkit v.0.7.171 using the ‘rmdup’ function with the ‘–by-seq’ flag. The remaining reads were processed with malt-run (v. 0.4.1) where BlastN mode and SemiGlobal alignment were used. The minimum percent identity (–minPercentIdentity) was set to 95, the minimum support (–minSupport) parameter was set to 10 and the top percent value (–topPercent) was set as 1. Remaining parameters were set to default. MEGAN615 was used to visualise the output ‘.rma6’ files and to extract the reads assigned to taxonomic nodes of interest for our sample. A taxon table of the raw MALT output for all samples and blanks, as well as species level read assignments to bacteria, archaea and DNA viruses for the ancient pitch sample are shown in Supplementary Data 5, where reads listed are the sum of all reads assigned to the species node, including reads assigned to specific strains within the species. Reads assigned to RNA viruses were not considered for further analyses, since our dataset consisted of DNA sequences only. Due to the limited number of reads assigned to archaeal species (Supplementary Data 5), we did not consider Archaea in downstream analyses of species identification. To validate the microbial taxa, we aligned the assigned reads to their respective reference genomes and examined the edit distances, coverage distributions, and post-mortem DNA damage patterns (see Supplementary Note 5).

Pneumococcus analysis

We reconstructed a S. pneumoniae consensus genome (Fig. 4) by mapping all reads assigned to S. pneumoniae by MALT16 to the S. pneumoniae TIGR4 reference genome (NC_003028.3). To investigate the presence of multiple strains we estimated the number of heterozygous sites using samtools57 mpileup function, filtering out transitions, indels, and sites with a depth of coverage below 10. Coverage statistics of the individual alignments (MQ ≥ 30) were obtained using Bedtools72 and plotted using Circos73 in 100 bp windows. Mappability was estimated using GEM274 using a k-mer size of 50 and a read length of 42, which is comparable to the average length of the trimmed and mapped reads in the ancient pitch. Virulence genes were identified by assembling the ancient S. pneumoniae MALT extracts into contigs using megahit75. The contigs were aligned against known S. pneumoniae TIGR4 virulence genes in the Virulence Factor Database17 (downloaded 22/11–2018) using BLASTn76. Only unique hits with a bitscore >200, >20% coverage, and an identity >80% were considered as shared genes (Supplementary Data 6).

To identify all streptococcus virulence factors in the ancient pitch, we aligned the contigs against the full Virulence Factor Database17 (downloaded 22/11–2018) using BLASTn76 and the same filtering criteria as described above (Supplementary Data 6). To validate the approach we repeated the analysis with five modern oral microbiome samples (SRS014468; SRS019120; SRS013942; SRS015055; SRS014692) from the Human Microbiome Project (HMP)14 using only the forward read (R1) (Supplementary Data 6). We find that the number of virulence genes we recovered directly correlates with sequencing depth (Supplementary Fig. 16).

Holi

For a robust taxonomic assignment of reads aligning to Metazoa (animals) and Viridiplantae (plants), all non-human reads were parsed through the ‘Holi’ pipeline18, which was specifically developed for the taxonomic profiling of ancient metagenomic shotgun reads. Each read was aligned against the NCBI’s full Nucleotide and Refseq databases (downloaded November 25th 2018), including a newly sequenced full genome of European hazelnut (Corylus avellana, downloaded April 10th 2019)77. The alignments were then parsed through a naive lowest common ancestor algorithm (ngsLCA) based on the NCBI taxonomic tree. Only taxonomically classified reads for taxa comprising ≥1% of all the reads within the two kingdoms and a declining edit distance distribution after edit distance 0 were parsed for taxonomic profiling and further validation. To validate the assignments, we aligned the assigned reads to their respective reference genomes and examined the edit distances, coverage distributions, and post-mortem DNA damage patterns (see Supplementary Note 5; Supplementary Data 7).

Reporting summary

Further information on research design is available in the Nature Research Reporting Summary linked to this article.