Sampling, DNA extraction, partial uracil DNA glycosylase library preparation and sequencing
We obtained permission from the Kunstkamera, Peter the Great Museum of Anthropology and Ethnography in St Petersburg for the sampling and ancient DNA analysis of 7 tooth specimens, excavated between 1885 and 1892 from the medieval cemeteries of Kara-Djigach and Burana (Supplementary Information 2). No statistical methods were used to predetermine the number of samples used in this study. All laboratory procedures were carried out in the dedicated aDNA facilities of the Max Planck Institute for the Science of Human History and the Max Planck Institute for Evolutionary Anthropology. The detailed procedures used for tooth sampling can be found in ref. 52. In brief, teeth were sectioned in the dentin–enamel junction using an electric saw with a diamond blade. After tooth sectioning, approximately 50 mg of powder was removed from the surface of the pulp chamber of each tooth using rounded dental drill bits.
The recovered tooth powder was used for DNA extractions using a previously established protocol optimized for the recovery of short fragments of DNA53. The exact steps and modifications of the procedure used have been made available in ref. 54. In brief, the tooth powder was incubated overnight (12–16 h) at 37 °C in 1 ml of DNA lysis buffer containing EDTA (0.45 M, pH 8.0) and proteinase K (0.25 mg ml−1). After incubation, DNA binding and isolation was performed using a custom GuHCl-based binding buffer and purification using High Pure Viral Nucleic Acid Large Volume Kit (Roche). Finally, DNA was eluted in 100 μl of Tris-EDTA-Tween containing Tris-HCl (10 mM), EDTA (1 mM, pH 8.0) and Tween-20 (0.05%). For procedure monitoring, extraction blanks and positive extraction controls were included throughout the laboratory processing steps.
All DNA extracts were converted into one-to-two double-stranded DNA libraries for Illumina sequencing, using 25 μl of input extract per library with an initial partial uracil DNA glycosylase (UDG) and endonuclease VIII treatment (USER enzyme; New England Biolabs) according to established protocols55,56. The detailed library preparation procedure, including the blunt-end repair, adaptor ligation and adaptor fill-in reaction steps can be found in ref. 57. After library preparation, each library was quantified using a quantitative PCR system (LightCycler 96 Instrument) using the IS7 and IS8 primers55. For multiplex sequencing, we performed double indexing of all libraries using previously published procedures58, outlined in detail in ref. 59. A combination of unique index primers containing 8 base pair (bp) identifiers were assigned to each library. To aid amplification efficiency, libraries were then split into multiple PCR reactions for the indexing step based on their initial IS7/IS8 quantification. The number of indexing PCR reactions performed for each library was determined so that every reaction was assigned an input of no more than 1.5 × 1010 DNA copies. Each reaction was set up using the Pfu Turbo Cx Hotstart DNA Polymerase (Agilent Technologies) and was run for 10 cycles using the following conditions: initial denaturation at 95 °C for 2 min followed by a cycling of 95 °C for 30 s, 58 °C for 30 s and 72 °C for 1 min, as well as a final elongation step at 72 °C for 10 min. All PCR products were purified using the MinElute DNA Purification Kit (QIAGEN), with some modifications to the manufacturer’s protocol59. Finally, all indexing PCR products were qPCR-quantified (LightCycler 96 Instrument) using the IS5 and IS6 primer combination58,59. To avoid heteroduplex formation, indexed libraries were amplified to 1013 DNA copies per reaction with the Herculase II Fusion DNA Polymerase (Agilent Technologies) and quantified using a 4200 Agilent TapeStation Instrument using a D1000 ScreenTape system (Agilent Technologies). Libraries were diluted to 10 nM and pooled equimolarly for sequencing. We performed shotgun DNA sequencing on an Illumina HiSeq 4000 platform using a 76-cycle kit (1 × 76 + 8 + 8 cycles).
Shotgun next-generation sequencing read processing and metagenomic screening
After demultiplexing, raw shotgun sequenced reads were preprocessed in the EAGER pipeline v.1.92.58 using AdapterRemoval v.2.2.0 (ref. 60), which was used to remove Illumina adaptors (minimum overlap of 1 bp), as well as for read filtering according to sequencing quality (minimum base quality of 20) and length (retaining reads ≥30 bp). Subsequently, all datasets were screened for the presence of pathogen DNA traces using the metagenomic pipeline HOPS29. First, preprocessed reads were aligned against a custom RefSeq database61 (November 2017) containing all complete bacterial and viral genome assemblies, a subset of eukaryotic pathogen assemblies and the GRCh38 human reference genome. Genome assemblies that contained the word ‘unknown’ were removed from the database, retaining a total of 15,361 entries. The database retained a number of Yersinia species entries: Yersinia aldovae (n = 1), Yersinia aleksiciae (n = 1), Yersinia enterocolitica (n = 16), Yersinia entomophaga (n = 1), Yersinia frederiksenii (n = 3), Yersinia intermedia (n = 1), Yersinia kristensenii (n = 2), Y. pestis (n = 39), Yersinia phage (n = 17), Yersinia pseudotuberculosis (n = 13), Yersinia rohdei (n = 1), Yersinia ruckeri (n = 4), Yersinia similis (n = 1) and Yersinia sp. FDA-ARGOS (n = 1). MALT v0.462 was run using the following parameters: -id 90 -lcaID 90 -m BlastN -at SemiGlobal -topMalt 1 -sup 1 -mq 100 -verboseMalt 1 -memoryMode load -additionalMaltParameters. The resulting alignment files were post-processed with MALTExtract for a qualitative assessment against a predefined list of 356 target taxonomic entries (https://github.com/rhuebler/HOPS/blob/external/Resources/default_list.txt). Specifically, reads were assessed according to their edit distance against a specific pathogen sequence in the database and the potential occurrence of mismatches that could signify the presence of aDNA damage29. In cases in which both parameters were met, the corresponding pathogen alignment was considered a strong candidate. Preprocessed reads were mapped against the Y. pestis CO92 (NC_003143.1) and human (hg19) reference genomes with the Burrows–Wheeler Aligner (BWA). Mapping parameters were set to 0.01 for the edit distance (-n) and seed length was disabled (-l 9999). Subsequently, we used SAMtools v.1.3 (ref. 63) to remove reads with mapping quality lower than 37 (for CO92) or 30 (for hg19); PCR duplicates were removed with MarkDuplicates v1.140 (http://broadinstitute.github.io/picard/). Finally, patterns of aDNA damage were assessed with mapDamage v.2.0 (ref. 64).
Single-stranded DNA library preparation and hybridization capture
For specimens BSK001 and BSK003, extra single-stranded DNA libraries were constructed from an input DNA extract of 30 μl. We performed library preparation at the Max Planck Institute for Evolutionary Anthropology using an automated protocol that is publicly available65. Single-stranded and double-stranded libraries from individuals BSK001, BSK003 and BSK007 were enriched using DNA probes covering the whole Y. pestis genome, as well as 1.24 million genome-wide SNP sites of the human genome66,67. For capture preparation, all libraries were amplified for the necessary number of PCR cycles to achieve 1–2 μg of input DNA. PCR reactions were carried out using the Herculase II Fusion DNA Polymerase. They were then purified using the MinElute DNA Purification Kit and eluted in EB elution buffer containing 0.05% Tween 20. Finally, library concentrations (ng μl−1) were quantified using a NanoDrop spectrophotometer (Thermo Fisher Scientific). For the in-solution Y. pestis captures, the probe set design was based on a set of publicly available modern genomes, specifically the Y. pestis CO92 chromosome (NC_003143.1), CO92 plasmid pMT1 (NC_003134.1), CO92 plasmid pCD1 (NC_003131.1), KIM10 chromosome (NC_004088.1), Pestoides F chromosome (NC_009381.1) and the Y. pseudotuberculosis IP32953 chromosome (NC_006155.1). For the in-solution human DNA captures, the probe set design was created to target 1,237,207 variants across the genome that are informative for studying the genetic history of worldwide human populations28,67. Both human DNA and Y. pestis hybridization captures were carried out for two rounds as described previously28,69,68,67,66, in which partially UDG-treated libraries from the same individual were pooled in equimolar ratios for capture and single-stranded libraries were captured separately.
Post-capture Y. pestis data processing
After Y. pestis whole-genome capture, libraries were sequenced on a HiSeq 4000 platform (1 × 76 + 8 + 8 cycles or 2 × 76 + 8 + 8 cycles) at a depth of approximately 11–27 million raw reads. The preprocessing of raw demultiplexed reads was carried out as described in the ‘Shotgun next-generation sequencing read processing and metagenomic screening’ section. At this stage, the datasets produced from partially UDG-treated libraries from the same individual were pooled and terminal bases were trimmed using fastx_trimmer (FASTX Toolkit 0.0.14, http://hannonlab.cshl.edu/fastx_toolkit/) to avoid damaged site interference with SNP calling during further processing. The following steps for read mapping, PCR duplicate removal and aDNA damage calculation were carried out in the EAGER pipeline70. We performed read mapping with BWA v.0.7.12 against the Y. pestis CO92 reference genome (NC_003143.1). For the pooled and trimmed partial UDG-treated libraries, BWA parameters were set to 0.1 for the edit distance (-n) and seed length was disabled (-l 9999). Given that the single-stranded libraries constructed for this study retained aDNA-associated damage, the BWA parameters were set to 0.01 for the edit distance (-n) to allow for an increased number of mismatches that could derive from deamination; seed length was disabled (-l 9999). We performed read mapping against the plasmids using the same parameters against a concatenated reference of all three Y. pestis plasmids (pMT1: NC_003134.1; pPCP1: NC_003132.1; and pCD1: NC_003131.1), masking the problematic pPCP1 region between nucleotides 3000 and 4200 that was shown to have high similarity to expression vectors used in laboratory reagents71. SAMtools v.1.3 (ref. 63) was used to remove all reads with mapping quality lower than 37 (-q), whereas MarkDuplicates was used to remove PCR duplicates. Deamination patterns associated with aDNA damage were retrieved with mapDamage v.2.0 (ref. 64). We used MALT62 for a taxonomic classification of mapped reads, to attempt a retention of reads that are more likely to be endogenous Y. pestis. MALT was run against the same database as described in the section ‘Shotgun next-generation sequencing read processing and metagenomic screening’, using the following parameters: -m BlastN -at SemiGlobal -top 1 -sup 1 -mq 100 -memoryMode load -ssc -sps. The minimum percentage identity parameter was set to default (-id 0.0), as opposed to a 90% identity filter used for running HOPS29, to avoid any reference bias that might arise from the removal of endogenous reads with a higher number of mismatches. After run completion, to retain the maximum number of reads accounting for the naive lowest common ancestor algorithm, we extracted reads that were assigned to the Yersinia genus node or summarized under the Y. pseudotuberculosis complex node. Reads were extracted in FASTA format from MEGAN v.6.4.12 (ref. 72). Subsequently, FASTA files were converted into FASTQ format with the reformat.sh script in BBMap from the BBtools suite (version 38.86, https://sourceforge.net/projects/bbmap/). FASTQ files were then remapped against the CO92 reference genome using the same parameters as described previously in this section. For single-stranded libraries, mapDamage v.2.0 (ref. 64) was used to rescale quality scores in read positions at which potential deamination-associated mismatches to the reference were identified. Subsequently, BAM files corresponding to the same individual were concatenated after mapping quality filtering and PCR duplicate removal. We performed concatenation using the SAMtools ‘merge’ command and with the AddOrReplaceReadGroups tool in Picard (http://broadinstitute.github.io/picard/) for assigning a single read group to all reads in each new file.
SNP calling, heterozygosity estimates and SNP filtering
Variant calling was carried out for BSK001 and BSK003, both before and after MALT62 filtering using the UnifiedGenotyper in the Genome Analysis Toolkit (GATK) v.3.5 (ref. 73). GATK was run using the EMIT_ALL_SITES option, which produced a call for every position on the chromosomal CO92 reference genome. The resulting genomic profiles of BSK001 and BSK003 were compared against a set of 233 modern and 46 historical Y. pestis genomes, as well as against the Y. pseudotuberculosis reference genome IP32953 (NC_006155.1), using the Java tool MultiVCFAnalyzer v.0.85 (https://github.com/alexherbig/MultiVCFAnalyzer). MultiVCFAnalyzer v.0.85 was run with the following parameters. SNPs were called at a minimum coverage of threefold and in cases of heterozygous positions, calls were made at a 90% minimum support threshold. In addition, SNPs were called at a minimum genotyping quality of 30. Furthermore, previously defined non-core and repetitive regions, as well as regions containing homoplasies, ribosomal RNAs, transfer-messenger RNAs and transfer RNAs were excluded from comparative SNP calling16,32. A set of 6,567 total variant sites were identified in the present dataset.
To investigate the extent of possible exogenous contamination within the BSK001 and BSK003 datasets, we estimated the number of ambiguous heterozygous variants beyond the SNP calling threshold. For this, MultiVCFAnalyzer v.0.85 (ref. 74) was used to generate an SNP table of alternative allele frequencies ranging between 10 and 90%. The results were then used to create ‘heterozygosity’ histogram plots of the estimated frequencies in R v.3.6.1 (ref. 75). Heterozygosity plots were created both before and after MALT filtering (see ‘Post-capture Y. pestis data processing’) to investigate whether taxonomy-informed filtering could aid the elimination of contaminant sequences in the investigated datasets (Supplementary Fig. 7).
An SNP table created with MultiVCFAnalyzer v.0.85, containing all variant positions across the present dataset, was filtered to identify SNP differences between the BSK001 and BSK003 genomes. The identified differences (n = 20) were then evaluated with the Java tool SNP_Evaluation30 (build date 13 August 2018; https://github.com/andreasKroepelin/SNP_Evaluation). The variant table and the VCF files of each genome were used as input for SNP_Evaluation. Furthermore, each identified private variant was evaluated within a 50 bp window and was considered ‘true’ when fulfilling the following criteria established in studies published previously17,21,30,76: (1) no multi-allelic sites were permitted within the evaluated window unless they were consistent with aDNA deamination (signified as spurious C-to-T or G-to-A substitutions); (2) the evaluated SNP position itself was not consistent with aDNA damage (no bases overlapping the SNP were downscaled by mapDamage v.2.0 (ref. 64)); (3) no gaps in genomic coverage were identified in the evaluated window; (4) reads overlapping the SNP sites showed specificity to the Y. pseudotuberculosis complex when screened with BLASTn (https://blast.ncbi.nlm.nih.gov/Blast.cgi).
Finally, to gain phylogenetic resolution, the BSK001 and BSK003 Y. pestis datasets were concatenated. We performed concatenation of BAM files, MALT62 filtering and aDNA damage rescaling (with mapDamage v.2.0 (ref. 64)) as described in the section ‘Post-capture Y. pestis data processing’. Moreover, the dataset was included in the comparative SNP analysis using MultiVCFAnalyzer v.0.85 (ref. 74) as described above. Finally, unique SNPs were evaluated with SNP_Evaluation30 according to the four criteria listed above.
Phylogenetic reconstruction and diversity estimations
Phylogenetic analysis was used to explore 233 Y. pestis genomes as part of the modern comparative dataset. An SNP alignment produced by MultiVCFAnalyzer v.0.85 (ref. 74) was used to construct a phylogenetic tree in MEGA7, using the maximum parsimony approach with 95% partial deletion (6,032 SNPs). Of the 233 modern Y. pestis genomes in the current dataset, 30 displayed extensive private branch lengths (Supplementary Fig. 12). Such an effect in bacterial phylogenies could result either from true biological diversity or from technical artefacts associated with false SNP incorporation during computational genome reconstruction. Although we cannot exclude the presence of several strains with exceedingly higher mutation rates in the current dataset, previous studies showed that modern Y. pestis strains with ‘mutator’ profiles are uncommon16,36. In this study, 27 out of 30 genomes that showed disparities in their private SNP counts compared to the rest of the dataset, were derived from assemblies for which the quality of SNP calls could not be evaluated (raw data unavailable). Because potential mis-assemblies or false-positive SNP calls can affect evolutionary inferences and diversity estimations, these genomes were excluded from further analyses. Therefore, we performed phylogenetic analysis using a subset of 203 modern Y. pestis genomes (Supplementary Table 13). The list of excluded genomes is as follows: 2.MED1_139 (ref. 19), 2.MED1_A-1809 (ref. 18), 2.MED1_A-1825 (ref. 19), 2.MED1_A-1920 (ref. 19), 2.MED0_C-627 (ref. 19), 2.MED1_M-1484 (ref. 19), 2.MED1_M-519 (ref. 19), 0.ANT5_A-1691 (ref. 18), 0.ANT5_A-1836 (ref. 18), 0.PE2_C-678 (ref. 77), 0.PE2_C-370 (ref. 77), 0.PE2_C-700 (ref. 77), 0.PE2_C-746 (ref. 77), 0.PE2_C-535 (ref. 77), 0.PE2_C-824 (ref. 77), 0.PE2_C-712 (ref. 77), 0.PE2b_G8786 (ref. 16), 0.PE4_I-3446 (ref. 78), 0.PE4_I-3517 (ref. 78), 0.PE4t_A-1815 (ref. 18), 0.PE4_I-3447 (ref. 78), 0.PE4_I-3518 (ref. 78), 0.PE4_I-3443 (ref. 78), 0.PE4_I-3442 (ref. 78), 0.PE4_I-3519 (ref. 78), 0.PE4_I-3516 (ref. 78), 0.PE4_I-3515 (ref. 78), 0.PE4_Microtus91001 (ref. 79), 0.PE5_I-2238 (ref. 80) and 0.PE7b_620024 (ref. 16).
A genome-wide SNP alignment consisting of 203 modern-day and 48 historical Y. pestis genomes (Supplementary Table 14), as well as the Y. pseudotuberculosis IP32953 genome, was used as input to construct a maximum likelihood phylogeny including 2,960 SNPs and up to 4% missing data. We performed phylogenetic analysis with RAxML81 v.8.2.9 using the generalized time-reversible (GTR) substitution model with 4 gamma rate categories. Finally, 1,000 bootstrap replicates were used to estimate node support for the resulting tree topology. After run completion, the maximum likelihood phylogenies were visualized with FigTree v.1.4.4 (http://tree.bio.ed.ac.uk/software/figtree/) and GrapeTree (v1.5.0)50.
To estimate the proportion of modern Y. pestis diversity descending from BSK001/003, we used the R package picante v1.8.282 to compute the MPD and FPD83 from the reconstructed maximum likelihood substitution tree. Measures made on a subset of the tree corresponding to the subclade descending from BSK001/003 (branches 1–4) were compared to that of the complete Y. pestis phylogeny. In both cases, only modern strains were included in the calculation. We used a bootstrapping approach to assess the sensitivity of our results with regard to sampling and phylogenetic uncertainty84. For each of the 1,000 RAxML bootstrap trees, we randomly resampled modern strains with replacement and only kept branches of the tree corresponding to the sampled strains. Diversity measures were performed for each of the obtained resampled bootstrap trees, from which median estimates and 95% percentile intervals were derived.
To assess the potential impact of uneven sampling among branches (branches 1–4 contained 130 modern strains whereas branch 0 contained only 73), we repeated the same analysis but adding an initial step intended to equalize the number of genomes in both parts of the tree. We subsampled branches 1–4 to the same number of strains as in branch 0 using sequence clustering in branches 1–4 to obtain representative subsamples. We performed hierarchical clustering based on pairwise phylogenetic distances (derived from the maximum likelihood phylogenetic tree) and the resulting tree was cut to define 73 clusters (functions hclust85 and cutree in R v.4.0.3). For each bootstrap tree, clusters were randomly downsampled to one strain, resulting in an equal number of strains between branch 1–4 and branch 0. Resampling with replacement was then applied as previously to each of the downsampled trees before computing diversity measures.
Plasmid SNP analysis
To investigate possible genetic variation among the plasmids of historical genomes, we performed read mapping of BSK001, BSK003 and BSK001/003 with BWA as well as SNP calling with GATK v.3.5 as described in the above section ‘SNP calling, heterozygosity estimates and SNP filtering’ against each of the three Y. pestis plasmids (pMT1: NC_003134.1; pPCP1: NC_003132.1; and pCD1: NC_003131.1). We then performed comparative SNP calling using MultiVCFAnalyzer v0.85 (ref. 74) against a set of 46 historical Y. pestis genomes as well as the modern reference strains CO92, KIM5 and 0.PE4-Microtus91001. Variants were filtered in individual genomes using SNP_Evaluation according to previously defined criteria (see the ‘SNP calling, heterozygosity estimates and SNP filtering’ section). In the present dataset, we identified ten variants in pCD1, eight in pMT1 and two in pPCP1 (Supplementary Table 15).
Time-calibrated phylogenetic analysis
To estimate the timing for the divergence of Y. pestis branches 1–4 using the BSK001/003 genomes as a new calibration point, we used a dataset comprising all modern genomes from branches 1–4 used for phylogenetic analysis (n = 130), genomes of the ancestral branching lineage 0.ANT3 (n = 8) and all 29 historical (fourteenth–eighteenth century) genomes in our dataset representing unique genotypes. In cases of identical genomes, the highest coverage genome was chosen for this analysis. We applied a molecular clock test using a maximum likelihood method in MEGA7 (ref. 86), using a GTR substitution model in which differences in evolutionary rates among sites were estimated using a discrete gamma distribution with four rate categories. On the basis of this molecular clock test, the null hypothesis of equal evolutionary rates across tested phylogenetic branches was rejected, which is consistent with previous studies showing substitution rate variation across Y. pestis lineages16,17. Therefore, a log-normal relaxed clock model was used for all subsequent molecular dating analyses.
For the molecular dating analysis, we used the Bayesian statistical framework BEAST2 v.6.6 (ref. 87). The ages of all ancient isolates were used as calibration points to construct a time-calibrated phylogeny with their radiocarbon or archaeological context age ranges set as uniform priors (see Supplementary Table 19 for all used age ranges). The ages of all modern isolates were set to 0 years before the present. We tested a number of coalescent tree priors such as the coalescent constant size, Bayesian skyline88 and exponential population models, all of which have been used or tested in previous ancient pathogen genomic studies17,89,90,91. We also tested the birth–death skyline tree prior, which has gained traction in recent years91,92,93 because it can account for epidemiological variables and can model sampling disparities through time94. Moreover, we used jModelTest v.2.1.10 (ref. 95) to identify the substitution model of best fit for our dataset. The indicated transversion model was implemented in BEAUti by using a GTR model (4 gamma rate categories) and the AG substitution rate parameter fixed to 1.0 (as indicated previously93). All tree priors were used in combination with a log-normal relaxed clock rate with a uniform prior distribution ranging between 1 × 10−3 and 1 × 10−6 substituions per site per year for the SNP alignment (1,405 sites after a 95% partial deletion), corresponding to a range of 3 × 10−7 to 3 × 10−10 across the entire genome, which is within the range of previous estimates17. As part of the phylogenetic topology set-up, all branch 1–4 genomes (ancient plus modern) as well as the 0.ANT3 lineage were constrained to be independent monophyletic clades. For the constant population size and exponential population tree priors, all other parameters were set to default. For the coalescent skyline tree prior, a Jeffreys prior distribution (1/x) was used for the population sizes and a dimension of 5 was used to permit variations in the group and population sizes through time, with an upper bound of 380,000 for the effective population size (default). Moreover, for the birth–death skyline tree prior, we used a uniform prior for the rate to become non-infectious that ranged between 0.03 and 70, to account for possible infectious periods ranging from 30 years (lifelong infections in rodent reservoirs96,97) to 5 days (average infectious period for bubonic plague98). We used a prior beta distribution with mean = 0.1 (alpha = 10.0, beta = 90.0) for the sampling probability ρ at time 0 and a uniform distribution ranging between 0 and 0.1 for the sampling proportion s. For the latter, two shifts were allowed through time. Finally, the reproductive number R was allowed to vary between 0 and 4.0 using a long normal prior distribution of median = 1.0 and s.d. = 0.7, which is within the range of previous estimates for bubonic and pneumonic plague during medieval epidemics98.
The suitability of all tree priors was evaluated using path sampling as implemented in the model selection package of BEAST2 v.6.6. Path sampling was run in 50 steps, with 20 million states as the chain length for each step. The resulting log-marginal likelihoods favoured with ‘strong support’99 the coalescent skyline model for the present analysis (log Bayes factor= 8.35 when compared against the second best model) (Supplementary Table 20). Therefore, the coalescent skyline model was chosen for further analysis. To evaluate the temporal signal in the present dataset, we used TempEst v.1.5.3 to estimate the root-to-tip distance against specimen ages in a linear regression analysis100. For TempEst, we used a maximum parsimony tree computed in MEGA7 (ref. 86) in NEXUS format. Moreover, we used the midpoint of the archaeological or radiocarbon date ranges for all ancient genomes as tip dates. All modern genome ages were set to 0 years before the present. The resulting correlation coefficient r (0.39) and R2 (0.16) values supported the existence of a temporal signal in the present dataset. Furthermore, we used the BETS approach101 for a temporal signal assessment that takes into account all analysis parameters. BETS compares the (log)-marginal likelihood estimations produced from an isochronous model (all sampling dates set to 0 years before the present) against a heterochronous model (including real sampling times). As previously, path sampling was run in 50 steps with 20 million states as the chain length for each step. The estimated (log)-Bayes factor of 129.33 was in strong support of the heterochronous model; therefore, it indicated the presence of a temporal signal in the present dataset.
For the molecular dating analysis using a coalescent skyline model set-up, we performed Markov chain Monte Carlo sampling using 2 independent chains of 300–400 million states each. After completion, runs were combined using LogCombiner v.2.6.7 and convergence was evaluated using Tracer v.1.6 (http://tree.bio.ed.ac.uk/software/tracer/) ensuring that the effective sample sizes were greater than 200 for each estimated posterior distribution after a 10% burn-in. Maximum clade credibility trees were constructed using TreeAnnotator in the BEAST2 v.6.6 package87 with a 10% burn-in and were then visualized in FigTree v.1.4.4. In parallel with the molecular dating analysis, we performed a sampling from the prior analysis to test for possible overfitting of the prior to the data. We performed Markov chain Monte Carlo sampling for 2 independent chains of 600 million states each. After run completion, runs were combined and convergence was evaluated after a 30% burn-in. The results indicate that the posterior distributions of the uncorrelated log-normal relaxed clock and the time to the most recent common ancestor estimates are not concordant with those obtained when using a data-informed analysis (Supplementary Fig. 13).
Because most Bayesian phylogenetic frameworks (such as BEAST2) are based on bifurcating trees and hence are poor at resolving multifurcating nodes, we complemented our approach by using TreeTime v.0.8.4 (ref. 35) to infer a time-calibrated phylogeny using a maximum likelihood approach. TreeTime has been shown to resolve polytomies in a way that is consistent with specimen tip dates. We generated a rooted maximum likelihood phylogeny using RAxML (Supplementary Fig. 10) from the same SNP alignment as the one used for BEAST2 (95% partial deletion). The maximum likelihood tree was then used as input for TreeTime, which was run using all known sampling dates for modern genomes and the midpoint of the age range for the ancient genomes (Supplementary Table 22). TreeTime was run using the Kingman coalescent tree prior with the skyline setting. An appropriate substitution model was chosen for the data using the -gtr infer option. The time-scaled phylogeny was inferred using an uncorrelated relaxed clock and with the branch length optimization, keep-root and keep-polytomies options. Moreover, the divergence time intervals were estimated from the highest likelihood tree using the -confidence option. Analyses were run using a maximum number of 500 and 1,000 iterations (maximum number of iterations option) and produced consistent outputs. The resulting time tree can be found in Supplementary Fig. 11.
Further information on research design is available in the Nature Research Reporting Summary linked to this paper.