Repeated plague infections across six generations of Neolithic Farmers

Sampling, library building and sequencing

We sampled a total of 133 teeth, 38 petrous bones and three femur samples from the collections at the University of Gothenburg (Gothenburg, Sweden), Statens Historiska Museum (Stockholm, Sweden) and Laboratory of Biological Anthropology (University of Copenhagen, Denmark). Laboratory processing of samples was carried out at the clean laboratory facilities at Globe Institute, University of Copenhagen and at the Human Evolution ancient DNA laboratory at Uppsala University, following strict ancient DNA guidelines. Subsampling and mechanical cleaning of the samples were carried out by removing the surface layer using a drill bit, followed by sampling of the dense bone material of the outer cementum layer for teeth, or the cochlea for petrous bones. DNA extraction was carried out on approximately 100–150 mg of bone material as previously described (ref. ⁶, modified after ref. ³⁸), followed by a uracil-DNA glycosylase treatment step to remove deaminated cytosines for the majority of samples (Supplementary Table 3). Subsequently, double-stranded DNA-sequencing libraries were built following ref. ³⁹. DNA sequencing was performed on a NovaSeq 6000 instrument (Illumina) in paired-end configuration at the GeoGenetics Sequencing Core (Copenhagen) using the S4 200 cycle kit (v.1.5). Furthermore, for selected samples with promising results but low content of endogenous DNA, single-stranded libraries were built as described in ref. ⁴⁰ followed by sequencing on the Illumina NovaSeq 6000 platform (Supplementary Table 3).

The Uppsala samples were decontaminated using ultraviolet irradiation, 0.5% sodium hypochlorite and deionized Millipore water. The outer surface layers were removed as described above, and 50–150 mg of either bone powder or bone pieces was collected using drill bits or diamond cutting wheels. DNA was extracted using silica-based protocols^41,42,43. Double-stranded DNA libraries were prepared as described above, with some being uracil-DNA glycosylase treated (Supplementary Table 3)⁴⁴. Some libraries were enriched using European MYbaits from MYcroarray following the manufacturer’s instructions. Sequencing was performed at SciLifeLab SNP&Seq in Uppsala using Illumina HiSeq 2500, HiSeq X Ten or NovaSeq platforms and 125 or 150 base pair paired-end chemistry.

Capture enrichment

For libraries in which shotgun-sequencing results suggested relatively high concentrations of plague DNA, in-solution capture enrichment of Y. pestis DNA was carried out. We used the myBaits kit from Arbor Sciences described in ref. ⁴⁵ and carried out hybridization capture following the manufacturer’s High Sensitivity protocol, but using only a single round of enrichment and 1.1–4.4 μl of baits per reaction. We pooled between one and 11 samples per capture reaction and reamplified the post-capture pools for 16 cycles (Supplementary Table 6).

Basic bioinformatic processing and PCAs

Demultiplexed fastq files were trimmed for adapter sequences, and overlapping paired-end reads were collapsed using AdapterRemoval (v.2.0)⁴⁶. Next, trimmed reads were mapped to the human reference genome (hg19) using bwa aln (v.0.7.17). The resulting alignments were converted to BAM files and merged, sorted and filtered using Samtools (v.1.12). The resulting BAM files were merged at sample level (Supplementary Table 2), and duplicates were identified using MarkDuplicates (v.2.27.4) from Picard Tools. Next, duplicate reads and reads with mapping quality less than 30 were removed from the sample-level BAM files using Samtools. Lastly, samples deriving from either the same individual or a monozygotic twin were merged based on both an initial run of ngsRelate and results from CrosscheckFingerprints (Picard; Supplementary Table 1). Because it is impossible, based on genetics, to distinguish whether two identical samples are derived from the same individual or from monozygotic twins, we merge all genetically identical samples. Nevertheless, in the case of FRA029, both samples belonging to this genetic individual are actually derived from the upper right first molar. This observation demonstrates that we have at least one pair of monozygotic twins in our dataset.

Average depth of coverage values was calculated using BEDTools genomecov (v.2.30.0), and mitochondrial haplogroups were assigned using mutserve (v.1.3.0) and haplogrep (v.2.1.25). Chromosome Y haplogroups were assigned using a tool developed in house and described in Supplementary Note 2.

For PCA we called pseudohaploid genotypes by random selection of one of the variants from mpileup (Samtools) that passed mapping and base quality filters (over 30). The resulting panel of pseudohaploid genotypes were merged with a panel of 1,430 previously published ancient shotgun genomes. We used smartpca from eigensoft to carry out PCA on 2,086,279 transversion-only SNPs with minor allele frequency of over 0.1%. All samples were projected onto the variation of the reference panel.

IBD analysis

For identity-by-descent analysis we imputed the genotypes of 82 individuals with sufficient coverage (over 0.1×) using Glimpse⁴⁷ (v.1.1.1; https://github.com/odelaneau/GLIMPSE/tree/glimpse1), as described in ref. ²³. Imputed genomes were then merged with a panel of imputed shotgun-sequenced genomes (Supplementary Table 7). Genomic segments shared via identity-by-descent were identified in all pairs of ancient samples using IBDseq⁴⁸ (v. r1206). Next we excluded the lowest-coverage individual in all pairs of close relatives (shared IBD segments over 1,500 centimorgans (cM), resulting in a total of 54 individuals from this study remaining for clustering (Supplementary Table 8). Genetic clustering was carried out based on IBD-sharing patterns in all pairs of samples, as described in ref. ²³. For analysis of ROHs we filtered the imputed panel described above for genotype posterior values over 0.95 and identified segments of homozygosity-by-descent using IBDseq (v.r1206). These segment lengths were converted to cM using the HapMap Phase II genetic map (https://github.com/johnbowes/CRAFT-GP/tree/master/source_data/genetic_map_HapMapII_GRCh37) and summarized in six bins based on segment length.

Kinship reconstruction

Initially we tested five different approaches to reconstruct kinship based on ancient DNA data (KING⁴⁹, ngsRelate^50,51, READ⁵², KIN⁵³ and the relationship inference algorithm based on pairwise identity-by-state sharing presented in ref. ⁵⁴). We found that all five approaches produced very similar results, but in some edge cases we observed that ngsRelate outperformed the four other approaches. Accordingly we decided to use ngsRelate for our analysis. We ran ngsRelate on a subset of the imputed panel described above, focusing on ancient Northern Europeans (816 individuals) and on 2,086,279 transversion-only SNPs with minor allele frequency over 0.1% in the panel. Next we classified all pairs of individuals as either zero-, first-, second- or third-degree relatives or as unrelated, based on pairwise relatedness (R_AB). We used the midpoint between the theoretical estimates of pairwise relatedness to distinguish between each degree of relatedness. For example, for zero-, first- and second-degree relatives the theoretical pairwise relatedness values are 1.00, 0.50 and 0.25, respectively. Accordingly we define first-degree relatives as pairs of individuals with R_AB values between 0.75 (midpoint between 1.0 and 0.5) and 0.375 (midpoint between 0.50 and 0.25; Supplementary Table 9). Lastly, to distinguish relationship type for second-degree relatives we used results from the likelihood files produced by KIN. To run KIN we first ran KINgaroo on filtered BAM files, focusing on 2,086,279 transversion-only SNPs with minor allele frequency over 0.1%. KIN was then run using default settings.

For pedigree reconstruction we used a combination of the automated tool PRIMUS (pedigree reconstruction and identification of a maximum unrelated set⁵⁵) and the manual triangulation approach described in ref. ³⁶ (Supplementary Note 2).

Pathogen screening

We carried out a screening for hits to all known human pathogens on all samples sequenced. The detailed description of the pathogen-screening pipeline can be found elsewhere⁵⁶. In brief, we classified all reads using KrakenUniq⁵⁷ (v.0.5.8) on a custom database of microbial reference genomes, with a focus on human pathogens and common environmental microbes. For genera represented by more than 50 k-mer in a sample, traditional mapping was carried out. For a given genus/sample combination we used minimap2 (v.2.17-r941)⁵⁸ to map all reads classified to that genus to corresponding reference genomes of all species within that genus. Next we calculated standard mapping statistics for each sample/species combination, such as depth of coverage, ancient DNA damage, average nucleotide identity, mean edit distance and mean read length. We also calculated an estimate of coverage evenness (actual breadth of coverage/expected breadth of coverage given the depth of coverage) and ranked each species within a genus based on its number of unique k-mers identified by krakenUniq. We define a positive microbe detection as species represented by more than 50 reads, that are ranked as having the highest number of unique k-mers within their genus, have an average number of soft clippings under eight, have an average nucleotide identity of over 0.97 and have a coverage evenness score of over 0.8.

Plague phylogenetic analysis

For all samples in which Y. pestis was detected in the pathogen screening, we mapped trimmed reads from all libraries of that sample to the reference plague genome (CO92; GCA_000009065.1) using bowtie2 (v.2.3.2) with the following parameters: -D 20 -R 3 -N 1 -L 20 -i S,1,0.50 –end-to-end –no-unal. Next we merged and sorted BAM files, and marked duplicates using MarkDuplicates from Picard Tools. Lastly, duplicate reads and reads of mapping quality under 30 were removed from the BAM files using Samtools.

The reference plague panel was generated by downloading raw read data from 400 modern plague strains, 30 ancient plague strains and data from the Y. pseudotuberculosis reference genome (GCF_000834295.1) at the European Nucleotide Archive. Raw fastq read files were trimmed with AdapterRemoval and mapped to the reference plague genome as described above, followed by filtering and calculation of the average depth of coverage using BEDTools genomecov. We called genotypes for plague samples with an average depth of coverage over ten for modern samples and over four for ancient samples, with the exception of Gok2, which was included despite having lower coverage. Genotype calling was carried out using Genome Analysis Toolkit⁵⁹ (v.4.1.9.0). First, sample-wise genomic variant call formats (GVCFs) were generated using HaplotypeCaller with a record at every single site (-ERC BP_RESOLUTION). Next, GVCFs were merged with CombineGVCFs and genotypes were called jointly with GenotypeGVCFs. Using VariantFiltration, genotypes were set to no-call if genotype quality was below 50 or if depth was over 1,000 or under four or ten, for ancient and modern samples, respectively. Furthermore, genotypes were set to no-call if genotype was supported by less than 90% of the reads at a given position. Lastly we filtered variants using VariantFiltration based on variant confidence/quality by depth, Phred-scaled P value, root mean squared mapping quality, z-score of Alt versus Ref read mapping qualities, z-score of Alt versus Ref read position bias, as follows: –filter-expression “QD < 2.0||FS > 60.0||MQ < 40.0 || MQRankSum < −12.5 || ReadPosRankSum < −8.0”.

To generate the phylogenetic tree we converted the vcf file to a multifasta using bcftools consensus, replacing missing data with N (flag: -M N) and filtering out regions with high proportions of reads of mapping quality zero (flag: -m [bedfile]; Supplementary Note 2). We used RAxML-NG (v.0.9.0)⁶⁰ to generate the phylogenetic tree, applying the GTR + G model and using the Y. pseudotuberculosis reference genome (GCF_000834295.1) as an outgroup. Bootstrap support values were calculated with raxml using both the transfer bootstrap expectation and Felsenstein’s bootstrap proportions metric. Both sets of bootstrap values produced similar support for the main nodes of the phylogeny. Lastly we plotted the best tree generated by RAxML-NG using the ggtree^61,62 (v.2.0.4) package in R (v.4.2.2).

For lower-coverage plague genomes (0.01–1.0×) for which traditional genotype calling was not feasible, we used EPA-ng (v.0.3.8)⁶³ to place each sample onto the phylogeny described above. For this approach we used major allele genotypes called by selecting the most common allele passing filters (mapping quality (MQ) over 30, base quality (BQ) over 30) from Samtools mpileup, and converted the vcf to fasta format as described above. Next we subset the fasta to the same SNPs as in the reference panel using BEDTools getfasta. We ran the evolutionary placement algorithm on the resulting fasta file, with default parameters supplying the best tree and the associated model from RAxML-NG as reference. Placements were plotted using the ggtree^61,62 (v.2.0.4) package in R (v.4.2.2).

Pangenome analysis

We used a pangenome variation graph framework to investigate variation in pangenome content among ancient Y. pestis lineages. We selected 81 complete reference assemblies of the Y. pseudotuberculosis complex (56 Y. pestis, 24 Y. pseudotuberculosis and one Y. similis) from NCBI (assembly level ‘chromosome’ or ‘complete’) and built a pangenome variation graph using pggb⁶⁴ with the following parameters:

• segment length: 5 kb (-s 5,000)

• pairwise identity: 0.9 (-p 0.9)

• minimum match length: 47 (-k 47)

• target sequence length for POA: 700, 900, 1,110 (-G 700,900,1,100)

• score parameters for POA: asm5 (-P asm5)

We built graphs separately for the chromosome and plasmids and merged them into a single-reference graph using the ‘vg combine’ function of the vg toolkit⁶⁵. The resulting graph was then indexed for short-read mapping with the Giraffe⁶⁶ mapper, with default parameters except k-mer length (-k 15) and window length (-w 5) for the minimizer. Shotgun-sequencing reads for each sample were then mapped to this index using ‘vg giraffe’. We filtered the resulting graph alignment files (GAM format) for mapping quality of 30 or above (-q 30) and normalized alignment score of 0.9 or above (-r 0.9) and removed reads with tandem repeat motifs at the read ends (-E 3). To obtain coverage statistics we then computed read support for each node in the reference graph using ‘vg pack’ with the ‘-d’ option. We considered a reference graph node covered in an ancient sample if the proportion of the node length covered by sequencing reads was 0.2 or above. To classify nodes into groups based on their presence or absence among the three Y. pseudotuberculosis complex species in the reference graph, we obtained a path coverage haplotype matrix for the reference contigs across all nodes in the graph using ‘odgi paths’⁶⁷ with the ‘-H’ option. Each node was then labelled according to the combination of Y. pseudotuberculosis complex species with coverage observed. For each ancient sample we aggregated coverage across nodes and classes to obtain the total amount of pangenome sequence contained in each class.

To identify the gene content of the Y. pseudotuberculosis complex pangenome in ancient Y. pestis lineages absent from modern Y. pestis strains, we extracted all nodes of the reference graph from classes not containing modern Y. pestis (that is, shared by Y. similis and Y. pseudotuberculosis or observed only within one of the two) in which at least one ancient sample had coverage. The reference graph start/end coordinates for those nodes were then translated into reference path positions using ‘odgi position’, and nodes were collapsed into regions if their translated position was within 200 base pairs of each other on the same reference contig. From this list of initial candidate regions we selected all of length greater than 500 base pairs, extended them by 10 kb on each side and merged overlapping regions to yield a final list of genomic regions harbouring putative pangenome content retained in ancient Y. pestis. To further investigate coverage across genes in those regions we split the reference graph into region-specific graphs using ‘vg chunk’ and produced single-reference alignments (BAM format) for each ancient sample and region using ‘vg surject’. The resulting BAM files were filtered by removal of duplicates (picard) and retaining only reads with mapping quality over 30, minimum number of aligned bases over 30 and, at most, five soft-clipped bases. Final per-gene coverage statistics were obtained from the filtered BAM files using ‘bedtools coverage’⁶⁸.

Radiocarbon dating and modelling

For Frälsegården a total of 91 radiocarbon dates are available, of which 46 are produced within the present project. The new datings targeted individuals in the pedigree, those showing plague and those with Steppe ancestry or HG (hunter-gatherer) admixture. Some samples were also redatings of individuals with anomalous previous dates. Seven new datings were made on samples from the Landbogården passage grave and six from the Nästegården site (Supplementary Table 4).

New datings were performed at the Keck carbon cycle AMS facility, University of California, Irving. The samples were decalcified in 1 N HCl, gelatinized at 60 °C and pH 2.0 and ultrafiltered to select a high molecular weight fraction (over 30 kDa). δ¹³C and δ¹⁵N values were measured to a precision level of below 0.1‰ and below 0.2‰, respectively, on aliquots of ultrafiltered collagen, using a Fisons NA1500NC elemental analyser/Finnigan Delta Plus isotope ratio mass spectrometer. Further details can be found at https://sites.ps.uci.edu/kccams/education/protocols/.

All samples were dated successfully, and showed C:N ratios well within the accepted range for well-preserved collagen, 2.9–3.6 (ref. ⁶⁹). δ¹³C and δ¹⁵N values indicate a generally terrestrial diet with little or no contribution from marine or freshwater protein, suggesting that reservoir effects are negligible (Supplementary Table 4).

Datings were calibrated in Oxcal 4.4.4 using the Intcal20 calibration curve^70,71. Because the time period of interest shows a large plateau in the calibration curve, the calibrations result in broad ranges, usually approximately 5250–4850 cal. bp (2σ). Chronological modelling of individuals within the pedigrees was made, allowing considerably higher precision. The model makes the simple assumption that genetically identified generations can be regarded as sequential phases. Tooth dates were offset by 20 ± 5 years to make them comparable to those from bone, based on formation time. The assumption of a generational interval of 29 ± 19 years, as suggested in ref. ²⁷, was tested but did not improve the result. Four dates were identified as outliers by both low-agreement indices and outlier analysis (HJE012, FRA023, FRA002 and FRA009). These were marked as outliers in the Oxcal codes but not modelled. Detailed results and Oxcal code can be found in Supplementary Note 3 and Supplementary Tables 11–13.

Strontium isotope analysis

Ten new Sr isotope values were produced within the present project, all on samples from Frälsegården (samples F11500–F11510; Supplementary Table 5). In addition, 24 samples from Frälsegården, Firse sten, Hunnebostrand and Rössberga were previously published^11,72,73. A further 34 samples were previously analysed and are reported here for the first time (Supplementary Table 5). The measurements are bulk values taken on small pieces of tooth enamel, representing an average over the time of enamel formation. All measurements of ⁸⁷Sr:⁸⁶Sr ratios were conducted at the Geochronology and Isotope Geochemistry Laboratory, Department of Geological Sciences, University of North Carolina, Chapel Hill. Samples were dissolved in nitric acid and the strontium fraction purified by ion selective chromatography (Eichrom Sr resin) before analysis by thermal ionization mass spectrometry (TIMS) on a VG Sector 54 mass spectrometer run in dynamic mode. Internal precision in the laboratory is consistently around 0.0007% standard error (or 1σ = 0.00006 in the ratio of a particular sample). Long-term, repeated measurements of SRM-987 are around 0.710260—an acceptable difference from the recognized value of 0.710250—and raw sample values from individual runs are standardized to the recognized value of SRM-987. Baseline values from Falbygden and surrounding areas in western Sweden were published in ref. ⁷⁴.

Reporting summary

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