Pierre Clavel,Lexane Louis,Clio Der Sarkissian,Catherine Thèves,Claudia Gillet,Lorelei Chauvey,Gaétan Tressières,Stéphanie Schiavinato,Laure Calvière-Tonasso,Norbert Telmon,Benoît Clavel,Richard Jonvel,Stéfan Tzortzis,Laetitia Bouniol,Jean-Marc Fémolant,Jennifer Klunk,Hendrik Poinar,Michel Signoli,Caroline Costedoat,Maria A. Spyrou,Andaine Seguin-Orlando,Ludovic Orlando
Improving the extraction of ancient Yersinia pestis genomes from the dental pulp
Archeology,Biochemistry,Molecular biology
Summary Ancient DNA preserved in the dental pulp offers the opportunity to characterize the genome of some of the deadliest pathogens in human history. However, while DNA capture technologies help, focus sequencing efforts, and therefore, reduce experimental costs, the recovery of ancient pathogen DNA remains challenging. Here, we tracked the kinetics of ancient Yersinia pestis DNA release in solution during a pre-digestion of the dental pulp. We found that most of the ancient Y. pestis DNA is released within 60 min at 37°C in our experimental conditions. We recommend a simple pre-digestion as an economical procedure to obtain extracts enriched in ancient pathogen DNA, as longer digestion times release other types of templates, including host DNA. Combining this procedure with DNA capture, we characterized the genome sequences of 12 ancient Y. pestis bacteria from France dating to the second pandemic outbreaks of the 17th and 18th centuries Common Era.
Improving the extraction of ancient Yersinia pestis genomes from the dental pulp

Ancient DNA preserved in the dental pulp offers the opportunity to characterize the genome of some of the deadliest pathogens in human history. However, while DNA capture technologies help, focus sequencing efforts, and therefore, reduce experimental costs, the recovery of ancient pathogen DNA remains challenging. Here, we tracked the kinetics of ancient Yersinia pestis DNA release in solution during a pre-digestion of the dental pulp. We found that most of the ancient Y. pestis DNA is released within 60 min at 37°C in our experimental conditions. We recommend a simple pre-digestion as an economical procedure to obtain extracts enriched in ancient pathogen DNA, as longer digestion times release other types of templates, including host DNA. Combining this procedure with DNA capture, we characterized the genome sequences of 12 ancient Y. pestis bacteria from France dating to the second pandemic outbreaks of the 17th and 18th centuries Common Era.

Ancient DNA preserved within archaeological and paleontological remains provides genetic information about our evolutionary past. The increasing throughput of next-generation sequencing instruments has considerably improved our capacity to characterize ancient genomes, despite environmental microbes often producing the dominant fraction of the DNA extracts. Specific types of remains, including petrosal bones, ossicles, and the tooth cementum, are known to be generally less prone to environmental DNA contamination. Together with improved decontamination and DNA extraction procedures, involving various pre-digestion treatments of sample powder through bleaching, and/or washing buffers,,, these biological remains have facilitated the characterization of ancient genomes. The development of DNA capture methodologies targeting million-scale SNP (Single Nucleotide Polymorphism) panels,, whole chromosomes or even the entire genome, have provided a further opportunity to focus sequencing efforts on the DNA fraction of interest, therefore, increasing sensitivity and reducing sequencing costs. As a result, there has been an acceleration since the first ancient human genomes have been sequenced in 2010,, to the thousands of ancient human genomes characterized, with recent human paleogenomic studies typically including hundreds of individuals., Besides host and environmental microbial material, ancient DNA extracts can also preserve molecules originating from pathogens that circulated in the blood system when the host died. In cases where the pathogenic load is high and environmental contamination is low, shotgun sequencing has allowed for the successful characterization of the entire genome sequence from ancient pathogens, including viruses (e.g., smallpox,, hepatitis B,), bacteria (e.g., Mycobacterium leprae), and protozoan parasites (e.g., Plasmodium vivax). Reconstructing ancient pathogen genomes, however, most commonly requires DNA capture technologies, both due to the molecular complexity of ancient DNA extracts and the generally limited genome size of pathogens, relative to their hosts., Additionally, the type of material considered for ancient pathogen genome characterization is critical, as tissues favoring host DNA preservation, such as petrosal bones, show virtually no bacterial pathogen DNA., Despite these limitations, more than five hundred ancient pathogen genomes have been sequenced, revealing important insights into the epidemiological outbreaks with substantive morbidity and mortality. Plague represents the best-studied infectious disease from an ancient DNA perspective. DNA analyses of human individuals buried during the Justinian plague and the Black Death have confirmed Y. pestis as the etiological agent responsible for massive human loss in the sixth-eighth century Common Era (c. CE), and the 14th c. CE, respectively. Such analyses have, however, revealed that plague already circulated during the late Neolithic and early Bronze Age,,,,,, and have started reconstructing the evolutionary trajectory of key genetic changes underlying virulence and transmission., Genetic evidence has also tracked the expansion of the second plague pandemic from a Central Asian homeland and the following diversification of several lineages through space and time. Despite such advances, the characterization of ancient pathogen genomes remains difficult. Relatively low detection rates and minimal target ancient DNA contents are commonly observed in mass graves, even when the association with known pandemic, such as the second plague pandemic, is documented historically. While the performance of various DNA capture systems is well established,,, no single study has attempted to optimize the molecular toolkit for pathogen DNA extraction and capture. In this study, we investigated the temporal dynamics of plague and host DNA release from dental pulp, using a total of 120 human teeth dated to the 17th and 18th c. CE. These remains originated from five French archaeological sites associated with the second plague pandemic. We found that the vast majority of the pathogen DNA can be recovered following a short pre-digestion wash of the dental pulp, while the recovery of the host DNA benefits from further digestion. Therefore, focusing on the pre-digestion wash improves the ratio of Y. pestis-to-host DNA. Whether such an approach could also enhance the recovery of other ancient pathogens in different preservation contexts and/or tissue types remains to be tested. Combined with targeted capture, pre-digestion washes helped us characterize 12 ancient Y. pestis genomes, including from material showing limited DNA preservation. Our work contributes to ongoing efforts aimed at mapping the genetic diversity of second pandemic plague strains.

We collected 120 teeth from 89 human individuals that were excavated from French mass graves and cemeteries related to the second plague pandemic (Figure 1, Data S1). More specifically, a total of 56 teeth originated from the sites of La Major (1720 CE, Marseille) and Les Rayettes (1721 CE, Martigues), both associated with the last major plague outbreak in France. This epidemic started after the ship Grand Saint Antoine departed from the Levant and landed on the Marseille harbor in 1720 CE. It rapidly decimated half of the population of the city, before it spread to the Provence region, hitting Marseille a second time (causing considerably lower mortality), and finally vanished in 1722 CE. Overall, it is estimated that this plague outbreak killed more than a quarter of the population of the Provence region. The remaining 64 teeth were collected from 58 individuals buried during the 17th c. CE, including from the Lariey Puy-Saint-Pierre cemetery (1629–1630 CE), as well as the Maladrerie Saint-Lazare (1630–1670 CE, Beauvais) and the Hôtel-Dieu Saint-Jean-Baptiste (1630–1670 CE, Amiens). While the Lariey cemetery was previously analyzed for Y. pestis and human ancient DNA, the latter two sites were genetically investigated here for the first time. These sites correspond to plague mass graves formed at the time of the Thirty Years’ War. A total of 33 teeth were first subjected to the standard procedure for ancient DNA extraction routinely applied at the Centre for Anthropobiology and Genomics of Toulouse (CAGT) ancient DNA facilities. Briefly, the dental pulp powder is collected by cutting the tooth at the cementum-enamel junction in order to drill the pulp chamber from inside with a diamond ball bit. The powder is first pre-digested in a DNA extraction buffer for 1 h at 37°C. The undigested tooth pellets are then collected and digested until completion in a fresh DNA extraction buffer, following overnight incubation with agitation at 42°C. This procedure was previously optimized for recovering ancient DNA from the host, but nonetheless led to the successful characterization of two ancient human and Y. pestis genomes from the Lariey cemetery, following shotgun sequencing on Illumina platforms. Stringent alignment of shallow sequencing data against the Y. pestis reference genome, revealed limited, if any, Y. pestis DNA preservation across all the samples investigated (< 0.93%), regardless of the archaeological site considered (Table S1). DNA capture procedures were, thus, implemented to attempt to enrich for Y. pestis and human mtDNA templates, using the synthetic RNA 80-mers probes designed by Wagner and colleagues. In this experiment, a total of 5 and 55 amplified DNA libraries were selected from individuals excavated at Lariey and La Major, respectively, as representatives of the entire range of DNA preservation conditions. Bioinformatic demultiplexing and collapsing of capture DNA reads provided between 14,386 and 2,250,565 sequences per individual, of which only those from two individuals from Lariey (LAR8 and LAR11) provided sufficient coverage of the Y. pestis genome (Table S1). Further metagenomics analyses in MetaPhlAn 4 confirmed a greater Y. pestis content in capture versus shotgun collapsed reads for 12 libraries prepared from six individual teeth (Figures 2A and 2B). This supported numbers of classified microbial species dropping with capture, as expected for a procedure aimed at focusing sequencing efforts on Y. pestis only. Accordingly, the fraction of environmental and oral bacteria identified following capture was reduced relative to shotgun experiments (Figure 2B). While targeted DNA capture improved the Y. pestis DNA recovery rate, the vast majority of the sequence data for other individuals consisted of PCR duplicates (37.72–98.03%), which limited the cost-effectiveness as an approach for characterizing the plague genome. The high proportion of PCR duplicates indicated that the total amount of plague DNA templates present in the original extracts was extremely limited. We, thus, decided to identify DNA extraction procedures that could enhance the recovery of ancient plague DNA molecules.

Previous work reported increased proportions of host DNA in libraries constructed following complete digestion, relative to those obtained pre-digestion., This suggested that pre-digestion could wash away a fraction of the environmental microbial DNA in contrast to the host DNA, which would be preserved deeper into the calcified matrix, and thus would require longer digestion time to be released. We reasoned that the Y. pestis DNA would be mostly restricted to the blood circulation and would not penetrate the deeper extracellular matrix. Under this assumption the pathogen DNA would follow release kinetics similar to the environmental microbial DNA. To test this theory, we purified the DNA from 17 1-h pre-digested fractions (1h) to prepare libraries for the capture of Y. pestis and human mtDNA genomes. The DNA contents of the pre-digested samples were compared to those obtained following capture of the DNA libraries prepared after complete pellet digestion (1h+O/N; Figures 3A and 3B). Normalizing for sequencing efforts after Y. pestis whole-genome capture, we found that most of the libraries constructed following 1 h pre-digestion (1h, without overnight digestion) contained larger proportions of Y. pestis sequences. We also observed inverse trends with clonality, which is reduced in 1-h pre-digested samples (Figure S1). Combined, these results support the presence of more Y. pestis DNA molecules in the pre-digested fraction. Human mtDNA content followed no particular trend, representing lower, similar or greater proportions of the sequence data depending on the specimen considered. This indicated that only some of the remains analyzed still contained a substantial fraction of human DNA in the undigested pellets and that DNA preservation was considerably more limited in others. Encouraged by the increased proportions of Y. pestis DNA present in the pre-digested fraction, we further investigated the kinetics of DNA release during the 60 min (min) pre-digestion wash. To achieve this, we generated fresh dental pulp powder for four additional teeth from individuals that showed positive plague identification in our previous experiments. The powder was then subjected to six consecutive partial digestion steps at 37°C of 10 min each (10 min, 20 min, 30 min, 40 min, 50 min, and 60 min; Figures 3C and 3D). The remaining undigested pellets were collected after each step, and incubated further in the same volume of a fresh extraction buffer. The resulting six pre-digested fractions were purified and used for library preparation, amplification, DNA capture, and sequencing (Figures 3C and 3D). Normalizing for sequencing efforts after Y. pestis whole-genome capture, we found that most of the Y. pestis DNA was released in the extraction buffer during the first 10 min pre-digestion. The proportion of Y. pestis DNA released during further pre-digestion steps was consistently lower than that released during the previous incubations. This indicated rapid Y. pestis DNA release in the extraction buffer, and diminishing returns with increasing pre-digestion times. Human mtDNA release followed a similar trend for one tooth (MAJ46t9), but exhibited no consistent trend across the four teeth investigated. Interestingly, MetaPhlAn 4 taxonomic profiling of microbial diversity for the single of the four teeth for which sufficient sequencing data were generated (≥ 100,000 collapsed reads; AMIENS268t2), revealed increasingly abundant off-target microbial hits with pre-digestion times (Figures 2C and 2D). This translated into an increasing number of classified microbial species, consistent with capture probes not being saturated by Y. pestis DNA templates only allowing aspecific annealing, including to oral bacterial DNA (e.g., Desulfobulbus oralis; Figure 2D). We repeated the same experiment for four additional teeth belonging to the same four individuals, except that shorter pre-digestion times were considered (0.5 min, 1 min, 2.5 min, 5 min, 7.5 min, and 10 min; Figures 3E and 3F). The proportions of Y. pestis and human mtDNA present in each resulting pre-digested fraction were compared to those found in the absence of pre-digestion (0 min), which were estimated after suspending the dental pulp powder in the pre-heated extraction buffer, and immediately centrifugating the pellets to collect the supernatant. Two of the teeth examined were found to contain almost no Y. pestis DNA considering pre-digestion times up to 10 min (MAJ46t11 and SQ8723t4; Figures 3E and 3F). Other teeth from these individuals, however, were tested positive for Y. pestis in previous experimental conditions. This highlights that the success of Y. pestis detection strongly depends on the tooth material considered, possibly due to differences in DNA preservation across tooth microenvironments and/or tooth physical integrity conditioning colonization success from the oral and environmental communities. The other two of the four teeth examined showed inconsistent kinetics of DNA release, with one sample (AMIENS268t4) releasing roughly twice as much Y. pestis DNA between 2.5 and 5 min than during the following 2.5 min (between 5 and 7.5 min), but the same proportion again between 7.5 and 10 min. For the last tooth (LAR27t4), the proportion of Y. pestis DNA collected in the last experimental condition was 6.43-fold–10.10-fold larger than during the two previous steps, despite lasting for 2.5 min each. It was also equivalent to the proportion of Y. pestis DNA obtained following a 1 min pre-digestion (Figure 3E). The AMIENS268t4 tooth provided sufficient numbers of collapsed reads for metagenomic profiling with MetaPhlAn 4, which revealed the presence of other microbes than Y. pestis in some but not all pre-digestion conditions (Figures 2C and 2D). Combined, these results indicated stochastic effects dominating DNA release of both Y. pestis and other microbes during the first 10 min of pre-digestion. No specific trend could again be observed for the release of human mtDNA. As the pre-digestion conditions between 10 and 60 min provided consistent results across the four teeth investigated, we decided to further characterize the underlying kinetics of DNA release in the pre-digestion buffer. We found that the proportion of unique stringent alignments characterized decreased exponentially through time (Figures 4A and 4B). This was true for both Y. pestis and human mtDNA, which followed trends that were not significantly different (ANCOVA test, p values ≥ 0.19). As a result, the cumulative proportion of Y. pestis and human mtDNA templates sequenced rapidly reached saturation, leaving only 5% (1%) of unreleased templates after 28.0–38.5 min (43.1–59.2 min) of pre-digestion (Figures 4C, 4D, and S2). This indicated that pre-digesting the dental pulp powder for 30 to 60 min in the experimental conditions tested would be sufficient to release most of the Y. pestis DNA present. We, thus, performed the DNA analysis of the remaining teeth following a conservative pre-digestion for 60 min at 37°C to extract the DNA fraction of the pathogen almost entirely. The supernatant was purified and concentrated for preparing a total of 176 DNA libraries for Y. pestis and human mtDNA enrichment.

The improved DNA extraction procedures developed in this study allowed us to characterize the genome sequence of Y. pestis strains present in 12 ancient individuals. Patterns of nucleotide mis-incorporations and base compositional profiles at the genomic positions preceding and following read termini were characteristic of authentic ancient DNA data generated following USER-treatment of extracts (Figure S3). Specifically, USER enzymes cleave those cytosine residues that are deaminated postmortem, leaving only slightly inflated C→T (G→A) mis-incorporation rates at read starts (ends), as well as an excess of ancient DNA templates starting immediately after cytosine residues. Read-to-reference edit distance distributions confirmed a closer genetic proximity to Y. pestis relative to the closest outgroup, Yersinia pseudotuberculosis (Figure S4). Phylogenetic reconstructions revealed that most of the second pandemic Y. pestis genomes clustered together within a monophyletic group, including strains from the Black Death to the 18th c. CE (Figures 5 and S5, Table S2). Within this cluster, the 12 ancient Y. pestis genomes characterized in this study grouped together with all the other post-Black Death ancient Y. pestis genomes hitherto sequenced. This was true genome-wide (Figure S6), or when considering only a subset of 29,609 parsimony informative sites located outside non-core, repetitive or highly conserved regions, which are reported to be less prone to environmental contaminant mis-mapping (Figures 5 and S5). The latter analysis confirmed the 1338–1339 CE genome from KaraDjigach, Kyrgystan, as basal to all subsequent strains sequenced, belonging to the second and third pandemic. Phylogenetic reconstructions also supported all modern branch 1 genomes responsible for the third pandemic descending from a subset of 14th c. CE strains identified in the Netherlands (BergenOpZoom), the United Kingdom (London), and Russia (Bolgar), as previously reported.,, Additionally, the LAR27 genome from Lariey characterized in this study (1629–1630 CE) appeared closely related to the two genomes previously reported from the same site (LAR8 and LAR11), confirming their genetic proximity with contemporary Y. pestis genomes from San Procolo a Naturno (Italy, 1636 CE). The Y. pestis genomes from Beauvais and Amiens (1630–1670 CE) also clustered with this group, suggesting common epidemiological origins. This lineage further spread during or in the aftermath of the Thirty Years’ War (Figures 5 and S5), which is recognized as one of the most destructive conflicts in European history, it affected Central Europe between 1618 and 1648 CE, a time contemporary with the ancient Y. pestis genome from Brandenburg, which is located basally within this group. The Y. pestis genome sequences retrieved from the samples of La Major (Marseille, 1720 CE) and Les Rayettes (Martigues, 1721 CE) formed a monophyletic group together with the previously reported genomes from Observance (Marseille, 1722 CE). This supports historical sources that describe an initial introduction of the plague in Marseille by the ship Grand Saint Antoine on May 1720 CE, and a subsequent rapid spread to the entire Provence region. Interestingly, our phylogenetic reconstructions indicated that the ancient Y. pestis genomes from Observance do not directly descend from those entering Marseille two years earlier and present at La Major. Taken at face value, this may have provided support for the introduction of two closely related plague strains. However, we caution that due to extreme DNA degradation, the genomes from La Major and Les Rayettes could only be characterized with limited coverage (between 1.19-fold and 5.25-fold), despite the improved methodology described in this study. Furthermore, the majority of ancient strains presented in our phylogenetic analysis appear to have excessively long terminal branches, even when compared to modern genetic diversity (Figures S5 and S6). The combination of low-genome coverage, remnant sequencing errors due to DNA degradation, and residual contamination from environmental microbes, can result in excessive terminal branch lengths in ancient genomes. The presence of such long branches can introduce significant biases in phylogenetic placement,, possibly impacting fine-scale evolutionary predictions. Additionally, we note that the Observance genomes were characterized using a different DNA capture system, introducing subtle technical batch effects, also potentially impacting phylogenetic reconstruction. This can be illustrated by the single Lariey genome characterized in this study, following DNA capture, which appears phylogenetic close, but does not form a monophyletic clade with those previously characterized by shotgun sequencing in Seguin-Orlando and colleagues. In light of these caveats, we consider our phylogenetic reconstructions consistent with the abundance of historical sources supporting a unique introduction of the plague in Marseille in 1720 CE. Determining whether the underlying strain first emerged in the Levant, or derived from a European source re-entering Europe through the Grand Saint-Antoine, requires further work. For example, mapping the genetic structure of Y. pestis strains at a finer-grained resolution during the 17th and the 18th c. CE may clarify the geographic sources and spread of related plague outbreaks. Integrating historical evidence will however be essential, given that the mutation rate of the Y. pestis genome provides genetic resolution at the 10-year timescale, which is considerably slower than the pace of human circulation and trade.

We further tested whether the ancient Y. pestis strains sequenced in this study showed specific genetic characteristics compared to their closest phylogenetic relatives. The average depth-of-coverage of the 12 individuals sequenced here ranged between 0.16 and 56.16-fold for the circular chromosome CO92 (median = 3.97-fold, average = 11.67-fold), while it ranged between 0.79 and to221.45-fold, 0.20– to73.46-fold and 5.71– to1092.67-fold for the pCD1, pMT1, and pPCP1 plasmids, respectively. Interestingly, the average proportion of high-quality alignments against the CO92 chromosome was lower-than-expected on the basis of the fraction of CO92 probes used for capture. The average proportions of high-quality alignment against the three different plasmids were, conversely, larger-than-expected on the basis of the probe distribution (Figures 6A–6D). This is in line with bacterial cells containing multiple copies of plasmids. The excess of plasmid DNA relative to probe proportions was, however, significantly more limited in the specimens dated from the 17th c. CE than in those from the 18th c. CE (Figures 6B–6D). This suggests that the latter contained a greater number of plasmid copies per cell than the former. Interestingly, the plasmid copy number can have an impact on virulence, in particular the pCD1 plasmid, which carries the T3SS genes responsible for delivering Yersinia outer membrane proteins into the host cells. The number of plasmids carried by the 17th and 18th c. CE strains may, thus, have driven differences in their virulence. Whether this genetic feature played a role in the excessive mortality observed during the Great Plague of Marseille, which decimated half the city’s population in only a few months, remains to be tested. Whether it also had equally dramatic consequences on the populations of rodent natural hosts, thus, precipitating disease disappearance is also unknown. Furthermore, a 49-kb region from the CO92 chromosome was found to be largely absent in the strains from La Major and Les Rayettes (Figures 7A–7E). The cspE gene (cspE_1883773–1883985), which is involved in the survival and the multiplication at low temperature, returned spurious non-null coverage, due to sporadic alignments resulting from environmental Yersinia bacteria. The 49-kb region spans the mgtA (mgtA_1888533–1891232, previously named mgtB) and mgtC (YPORS09280_1887434–1888132) genes involved in the survival of the bacteria under condition with low concentrations of Mg2+. It was previously reported to also be missing from the genomes of their closest phylogenetic relatives, including 18th c. CE strains from France (Observance, Marseille), Russia (Maist), and Sweden (Pestbacken), as well as mid-16th/mid-17th c. CE strains from the United Kingdom (London). Moreover, a 47-kb region overlapping 11 virulence and pathogenicity factors (inv_2040296–2042350, fyuA_2140840–2142861, ybtE_2142992–2144569, ybtT_2144573–2145376, ybtU_2145373–2146473, irp1_2146470–2155961, irp2_2156049–2162156, ybtA_2162347–2163306, ybtP_2163563–2165275, ybtQ_2165262–2167064, ybtX_2167057–2168337) was largely under-covered in all the genomes from La Major, Beauvais, Lariey and Amiens reported here (Figure 7A). They were, however, covered in two genomes from Lariey that were previously characterized using shotgun sequencing. We, thus, consider that the limited coverage observed most likely reflects under-performing probe-template annealing in those regions than genuine genomic rearrangements. We next used snpToolkit to characterize the single nucleotide mutations found among the new 17th and 18th c. CE strains sequenced here, and their closest phylogenetic relatives. This analysis revealed a total of 375 SNPs, comprising 148 non-synonymous, 97 synonymous, and 130 intergenic mutations (Figure S7, Table S3). Interestingly, the iron ABC transporter encoded by the YbtP gene, which is essential for the bacterial physiology and can deplete Y. pestis avirulence in mice, was found to be polymorphic in both strains from the 17th and 18th c. CE. No mutations were otherwise found in the 235 virulence and pathogenicity loci considered in the analyses above. The genomes from Amiens and Beauvais Y. pestis strains (1630–1670 CE) showed remarkable genetic similarity with those from Lariey (1629–1630 CE), as they shared only eight non-synonymous variants not found in the latter (Table S3). One such variant affected the autotransporter YapC gene, which mediates adhesion to human cells and, thus, may affect host innate immune response. In contrast, the number of non-synonymous differences found between the 1720–1722 CE genomes from La Major, Les Rayettes, and Observance was larger. Whether this is indicative of currently unidentified minor mis-mapping from other bacteria than Y. pestis, or faster mutation rates in this specific lineage requires further clarification through additional sequencing of complete genomes from Y. pestis strains responsible for the plague epidemic that spreads to the Provence region in the early 18th c. CE.

In this study, we report that ancient Y. pestis DNA is almost entirely released during a 60 min pre-digestion wash of the dental pulp. This finding might have important methodological implications for the workflow of ancient pathogen studies, as DNA extracts obtained following pre-digestion appear to maximize the ratio of ancient pathogen DNA to human DNA. This contrasts to the extracts obtained following complete pellet digestion, which are generally used for the study of ancient plagues, despite containing a much greater ratio of human to pathogen DNA. In the experimental conditions described herein, a pre-digestion of the dental pulp powder for 60 min at 37°C appears sufficient to release at least 99% of the ancient Y. pestis DNA generally preserved. Applying targeted capture to these DNA extracts allowed us to characterize 12 ancient Y. pestis genomes from five French archaeological sites dated to the 17th and 18th c. CE. We find that the number of plasmid copies within Y. pestis bacteria from the 17th c. CE was lower than in those from the 18th c. CE. The Thirty Years’ War likely facilitated the spread of a plague lineage throughout Germany, Italy, and France, all of which were eventually dead-end emergences. The last plague outbreak in France appears to follow a single emergence in Marseille, before spreading to the Provence region.

In this study, genome variation among ancient Y. pestis strains was only investigated from sequence alignments against one single reference genome. The full extent of structural variation and their possible consequences on transmission, pathogenicity and virulence remains to be characterized. Moreover, the quality of the data generated for most of the ancient genomes presented was limited, which precluded the examination of fine-scale strain evolution and divergence dates. Only a small number of samples were subjected to pre-digestion tests, they were also confined to a single geographic area and, only encompassed preservation conditions from a relatively recent period. Even if our results highlight a strong pattern for the DNA release of Y. pestis, experiments should be extended to various areas, periods and epidemics. Finally, whether the recommended methodology shows similar performance on other types of non-dental tissues, and will equally apply to the entire range of pathogenic bacteria that can circulate in the blood, remains untested.



Further information and requests for resources and materials may be directed to lead contact, Prof. Ludovic Orlando (ludovic.orlando@univ-tlse3.fr).

This study did not generate new unique reagents or materials.

The 120 teeth analyzed in this study originate from 89 individual skeletons buried in five archaeological sites (Figure 1, Table S1; see Data S1 for further description of the archaeological contexts). Prior to destructive sub-sampling, precise virtual 3D models were constructed for the majority of teeth investigated (94/120 teeth), using the ArtecMicro© surface scanner (Artec 3D, 1 mm voxel size). Samples and laboratory procedures were registered in the Laboratory Information Management System CASCADE. The plague cemetery of Lariey Puy-Saint-Pierre dating to 1629-1630 CE is located in the French Alps, and was excavated in 2002 CE (Minimal number of individuals, MNI = 34). Two out of 10 individuals (LAR8 and LAR11) investigated in this study were previously sequenced by Seguin-Orlando and colleagues in 2021, who demonstrated that they were infected with a Y. pestis strain genetically close to those infecting 1636 CE individuals from San Procolo a Naturno, Italy (1636 CE), in the context of the Thirty Year’s War. The excavations at the Hôtel-Dieu Saint-Jean-Baptiste of Amiens, France, were carried out in 2017 CE by the Service archéologie préventive d'Amiens Métropole, under the supervision of Richard Jonvel, who provided access to the city’s ancient plague cemetery (MNI = 259). The 34 individuals studied here come from multiple graves dating to 1630-1670 CE, when Amiens experienced notable military and epidemic turmoil. The site of the Maladrerie Saint-Lazare (Beauvais, France) was excavated between 2002 and 2014 CE by the Service archéologique municipal de Beauvais (MNI = 490). The discovery of a plague cemetery established since 1623 CE allowed the analysis of 14 individuals coming from the mass grave labelled 8559. The parish cemetery of La Major cathedral (La Major), excavated in 2008 CE was created at the beginning of September 1720 CE during the last great wave of plague in Marseille (1720-1722 CE). A total of 19 individuals were DNA tested, but showed poor macroscopic preservation, consistent with that observed across the many other individuals buried at the same site (MNI = 106). The site of Les Rayettes (Capucins de Ferrières trenches, Martigues, France, MNI = 208), was excavated in 2002 CE. It is located less than 40 kilometers away from La Major. The 12 individuals investigated here were found in multiple mass graves dated to 1720-1721 CE and characteristic of the epidemic ongoing locally at the time.

All molecular work, from dental pulp sampling to library construction and PCR mix preparation, was carried out in the ancient DNA facilities of the CAGT (CNRS UMR 5288/Université Paul Sabatier), following strict experimental standards to avoid and monitor contamination, including use of disposable personal protection equipment, positive air pressure, bleach/UV surface and instrument decontamination, and negative control blanks. PCR amplification, purification, quantitation as well as DNA capture and sequencing were carried out in post-PCR laboratories that are located in a building physically isolated from the ancient DNA facilities. Between 15 and 80 mg of dental pulp were collected for each of the 120 teeth analyzed (incisors, canines, premolars and molars), following the method developed by Neumann and colleagues, which is routinely used at CAGT. The extraction was performed following the protocol described by Seguin-Orlando and colleagues, which corresponds to a modified version of the Y2 protocol described by Gamba and colleagues. The standard protocol makes use of 963 μL of lysis buffer (0.45M EDTA, 0.25 mg/mL proteinase K and 0.5% sodium lauroyl sarcosinate) to carry out a first pre-digestion step of 1 hour at 37°C. The supernatant is then collected following centrifugation at 12,000 rotation per minute (rpm) for 30 seconds and a 963 μL volume of a fresh, pre-heated, lysis buffer is added to proceed to full digestion of the remaining pellets at 42°C overnight. After centrifugation at 12,000 rpm for 30 seconds, 200 μL of the lysis buffer is purified through QIAgen© Minelute columns and concentrated in a final elution buffer consisting of 23 μL EB supplemented with 0.05 % Tween 20 (later referred to as EB-Tween buffer). A total of 22.8 μL purified extract was treated with 7 μL of USER® enzyme (NEB®) for 3 hours at 37°C, following, to reduce spurious nucleotide mis-incorporation at those deaminated cytosine residues formed by post-mortem degradation. The pre-digestion supernatant was also purified and processed in parallel to the overnight digestion fraction, to contrast the DNA content following shotgun DNA sequencing and/or DNA capture. The standard extraction protocol was modified to assess the kinetics of DNA release in the lysis buffer. This was achieved through two main experiments, each focusing on two different teeth from four individuals that were previously tested positive for Y. pestis DNA (AMIENS268, tooth 2 (t2) and t4; LAR27 t2 and t4; SQ8723 t2 and t4, and; MAJ46 t9 and t11). The first experiment followed the standard extraction protocol, with the exception that the whole supernatant was collected every 10 min and replaced by 963 μL of a fresh lysis buffer. Supernatant collection was carried every 10 min until 1 hour of pre-digestion was completed. The second experiment followed the same procedure as the first, except that the supernatants were collected directly after suspending the dental pulp powder in the lysis buffer (0 min), and 0.5, 1, 2.5, 5, 7.5 and 10 min after. Each collected supernatant was further processed following the methodology described above to obtain 29.8 μL of USER-treated DNA extract.

Triple-indexed double-stranded DNA libraries compatible with Illumina sequencing were prepared following the protocol originally described by Rohland et al., 2015, as modified by Fages and colleagues. This protocol includes a unique combination of two barcodes, each located 3’ of the P5 and P7 library adapters. Due to their location immediately downstream of the sequencing primers, these “internal” barcodes form the first seven nucleotide positions within each sequencing read. Each library was amplified in eight parallel PCRs for eight cycles in the same conditions as described by Seguin-Orlando and colleagues: in a total reaction volume of 25 μL, using 0.4 μL of AccuPrime™ Pfx DNA polymerase (one unit), 2.5 μL 10X Accuprime™ Pfx reaction mix, 2 μL of DNA library, 1 μL of BSA (20 mg/mL), 17.9 μL of H2O, 0.2 μL of inPE1 primer (25 μM) and 1 μL indexing primer (5 μM) containing a unique 6-bp “external” barcode. The latter corresponds to a third barcode located within the PCR primer. The same external index was used in each of the eight parallel amplifications carried out on the same original library. Libraries were then amplified through a second additional PCR round to obtain sufficient DNA mass as input for capture enrichment with a custom myBaits® hybridization capture kit (Daicel Arbor Biosciences). To achieve this, two amplification products for a given library were merged and purified together using Agencourt Ampure XP beads (1:1.2 or 1:1.4 as DNA:beads ratio), and eluted in 11 μL of EB-Tween buffer. This provided four purified amplification products per library, which were further split and processed through two parallel PCR amplifications, using 1 μL of the eluate for each PCR reaction, for 10 cycles to obtain sufficient DNA material for DNA capture. These amplifications formed the second PCR round and were carried out in 25 μL reaction volumes, following the same conditions as for the first amplification round, except that IS5 and IS6 were used as PCR primers. The resulting eight PCR products were co-purified on a single MinElute column (QIAGEN©), and eluted in 12 μL of EB-Tween buffer. Library concentration and/or size distribution were measured on a Tapestation 4200 instrument using a High sensitivity D1000 ScreenTapeAssay (Agilent technologies) as well as on a Qubit HS dsDNA assay (Invitrogen).

A total of four to eight libraries from different extracts with unique internal and external barcodes were pooled and concentrated together on a MinElute column (QIAGEN©). This resulted in pools of 7 μL, containing between 2-12 μg of DNA, as recommended in the myBaits® High Sensitivity Protocol Version 5.00. DNA capture was performed according to the aforementioned protocol, consisting of two rounds of 16 to 24 hour hybridization with the probes described by Wagner and colleagues. Captured libraries were sequenced using the Paired-End mode on an Illumina MiniSeq instrument (2x81 cycles) at CAGT. A similar procedure was followed for shotgun sequencing, except that only one amplification round of 8-12 PCR cycles was performed on each library, before purification, quantification, pooling and DNA sequencing.

Read demultiplexing and collapsing, as well as adapter and poor-quality end trimming were performed using AdapterRemoval2 version 2.3.0, tolerating at most one mismatch in each internal barcode (--barcode-mm-r[12] 1 --minadapteroverlap 3 --mm 5). Read mapping against the human reference genome (hg19, GRCh37) and the revised Cambridge Reference Sequence (rCRS, Accession Number = NC_012920.1) were carried out using Bowtie2 version 2.3.4.1, and the recommended parameters from Poullet and Orlando, 2020. Read mapping against the Y. pseutotuberculosis reference genome (Accession Number NC_006155.1), the Y. pestis reference genome (strain CO92, Accession Number = NC_003143.1) and the three plasmids pCD1 (Accession Number NC_003131.1), pMT1 (Accession Number NC_003134.1) and pPCP1 (Accession Number AL109969.1) were performed using BWA version 0.7.17 “backtrack” mode, and the “stringent” parameters defined by Spyrou and colleagues. Read (re-)alignment, PCR duplicate removal, quality filtering and mapDamage profiling were carried out using the automated PALEOMIX pipeline v1.2.13.2. The same procedures were followed using previously published read sequence data (Table S2) to obtain the comparative panel of Y. pestis genomes and plasmids used for phylogenetic reconstructions and/or for assessing the presence/absence of specific virulence genes and lineage-specific mutations, which was summarized using snpToolkit v2.3 (see below). Basic alignment statistics for each sample, DNA library and experimental conditions are provided in Table S1. When comparing the respective performance of each experimental condition, Y. pestis and human DNA content of each library were estimated for normalized sequencing efforts. This was achieved by considering the fraction of collapsed reads only, and down-sampling the sequencing data 10 times for each library to the number collected for the least sequenced library. The investigation of pathogen and mitochondrial unique reads released through time was performed on the first set of extraction tests on four teeth, which were digested in lysis buffer during one hour, sampling and replacing the supernatant with fresh lysis buffer every 10 min. Results were normalized by down-sampling the number of reads for each digestion time by the lowest number of reads sequenced across the different digestion times for each sample. As the percentage of unique reads mapping against the pathogen genomes (CO92, pCD1, pMT1 and pPCP1) over digestion time approximatively fit an inverse exponential function, data have been log-linearized. Linear models of the logarithm of the percentage of unique reads mapping against pathogen and mitochondrial DNA (mtDNA) were performed using the R lm command. The remaining Y. pestis and human mtDNA to be released at a given time could be predicted using predict command in R, according to these models. The proportions of unique reads mapping against the circular chromosome and the different plasmids of the pathogen were calculated for each site to compare the plasmid load between 17th and 18th c. CE strains. These were contrasted to the proportion of probes targeting these different genetic materials. Calculations were based on libraries subjected to a one-hour digestion, as long as at least 500 high-quality alignments against the CO92 genome were found. Inter and intra-century differences in these proportions were compared for 17th and 18th c. CE samples using non-parametric Wilcoxon and Dunn tests in R. The read-to-reference edit distance distributions (Figure S2) were generated for each of the new genomes characterized in this study from the NM:I field obtained while running the samtools view command, considering either the Y. pestis (CO92, NC_003143.1), or the Y. pseudotuberculosis (IP 32953, NC_006155.1) as a reference. Depth-of-coverage variation was calculated using the PALEOMIX coverage command and non-overlapping sliding windows of 1,000 bp for the CO92 chromosome, or 100 bp for the pCD1, pMT1 and pPCP1 plasmids. GC content was estimated for the same windows using seqtk (https://github.com/lh3/seqtk), and visualized together with coverage estimates through circular plots produced using the R circlize library. The fraction of the positions that are covered at least once within a given gene was also calculated using the PALEOMIX depth command to assess the presence or absence of 235 plague virulence and pathogenicity loci, 190 of which located on the CO92 chromosome, and 45 on the three bacterial plasmids.,, Gene coverage plots for virulence genes were generated using the heatmap function from the R ggplot2 library (https://cran.r-project.org/web/packages/ggplot2/index.html). A 49-kb region of CO92 was found absent in most of the phylogenetically-close Y. pestis strains from New Churchyard (BED038, BED034, BED028, BED024, BED030), Azov (Azov38), CHE1, Rostov2033, Observance (OBS107, OBS116, OBS137, OBS110, OBS124), and La Major (MAJ75 and MAJ46). However, the coverage of some genes within this region was sporadically greater than zero in a subset of strains, suggesting the possible presence of spurious alignments, despite stringent alignment parameters. This was further confirmed using BlastN against the nucleotide collection (nt), showing best hits against environmental Yersinia strains than Y. pestis. Those cases are highlighted on the heatmap presented on Figure 7E by crossing the cell presenting the corresponding coverage estimate. We used snpToolkit v2.3 and the individual vcf files generated as described above to annotate the new Y. pestis genome and plasmid sequences characterized in this study for synonymous and non-synonymous sequence polymorphisms present in at least 90% of the reads. The minimal depth threshold was set to 3 (-d 3) for bases showing Phred scores above or equal to 30. The samples considered for this analysis were merged by archaeological site using samtools to provide sufficient coverage for the detection of strain-specific SNPs. Microbial taxonomic profiles were obtained using MetaPhlAn 4 on collapsed sequencing unique reads after filtering for those mapping against the human nuclear and mitochondrial genomes following the procedure presented above. To accommodate for variable sequencing efforts across DNA libraries, we random-sampled a number of 100,000 reads. This procedure was repeated 10 times to assess the robustness of the taxonomic assignments, given the limited of reads considered. Reads were mapped against the MetaPhlAn 4 biomarker database using Bowtie2 v.2.3.4.1 with default parameters, and further filtered for minimal mapping quality of 25, and PCR duplicates using the view and rmdup modules from samtools v1.11, respectively. Those species showing relative abundances <1% in a given dataset were defined as “low abundance”. Maximum Likelihood phylogenetic trees were constructed using IQ-TREE v1.6.12 and a panel of Y. pestis genomes, including 124 modern accessions and 61 second pandemic samples, 12 of which were sequenced here for the first time. Individual base calling was carried out using BCFtools (--ploidy 1), requiring minimum base, mapping and genotype Phred quality scores of 30, and filtering those variants located within 10 bp of indels. The best mutational model was selected based on AICc and BIC values with the Model Finder Plus module (MFP), as part of IQ-TREE. Node support was estimated from a total of 1,000 pseudo-replicates using the ultrafast bootstrap (UFBoot) approximation, optimized by nearest neighbor interchange (NNI), as well as the SH-like approximate likelihood ratio (aLRT) test. The final sequence alignments corresponded, on the first hand, to the entire CO92 chromosome, considering polymorphic sites only, or, on the other hand, were restricted to 29,609 genome-wide variant positions previously characterized to provide robust phylogenetic signal., Two sequence alignments were considered, requiring a conservative minimum 3-fold average depth-of-coverage per sample. Accessions for the sequence data underlying the comparative panel of 185 Yersinia pestis and 1 Yersinia pseudotuberculosis genomes are provided in Table S2.