Plagued by a cryptic clock: insight and issues from the global phylogeny of Yersinia pestis

Plague has an enigmatic history as a zoonotic pathogen. This infectious disease will unexpectedly appear in human populations and disappear just as suddenly. As a result, a long-standing line of inquiry has been to estimate when and where plague appeared in the past. However, there have been significant disparities between phylogenetic studies of the causative bacterium, Yersinia pestis, regarding the timing and geographic origins of its reemergence. Here, we curate and contextualize an updated phylogeny of Y. pestis using 601 genome sequences sampled globally. Through a detailed Bayesian evaluation of temporal signal in subsets of these data we demonstrate that a Y. pestis-wide molecular clock is unstable. To resolve this, we developed a new approach in which each Y. pestis population was assessed independently, enabling us to recover substantial temporal signal in five populations, including the ancient pandemic lineages which we now estimate may have emerged decades, or even centuries, before a pandemic was historically documented from European sources. Despite this methodological advancement, we only obtain robust divergence dates from populations sampled over a period of at least 90 years, indicating that genetic evidence alone is insufficient for accurately reconstructing the timing and spread of short-term plague epidemics.

• Is there a specific reason why the oldest documented Y. pestis lineages were not included in this study? Currently, they are represented by the genomes Gok2 or RV2039. I suggest their addition as they might help improve the phylogenetic signal across all Y. pestis.
• One of the determining priors in Bayesian molecular dating analysis is the demographic model used. Yet, no testing of demographic models is performed here and the authors choose to proceed with a coalescent constant population size model that seems inappropriate for describing the evolutionary history of epidemic pathogens. The suitability of other models, such as the coalescent skyline and exponential growth models should be fully evaluated for the entire dataset, as well as all the subsets.
• The authors use different substitution models for the maximum likelihood and Bayesian phylogenetic analyses. If no specific reasoning exists, the model chosen by ModelFinder should be used for all analyses.
• The first section of this paper "Population Structure" includes a critical discussion of genomic nomenclature and population distinctions in Y. pestis. Yet the authors choose to use those arbitrary distinctions between populations for their subsequent analyses. How are the 12 "populations" defined in this study? Given their criticism regarding the distinction between the second and third pandemic, could all post-Big Bang lineages be considered as one population? Is a molecular dating analysis of combined Branches 1-4 supported by BETS?
• The authors mention that "BETS was inconclusive when attempting to fit a single clock to the updated global diversity of Y. pestis" but the accompanying Bayes Factor comparisons are not presented in the paper. To aid transparency, those should be presented as a supplementary display item.
• To what extent could false SNPs affect estimates of rate variation across different phylogenetic lineages? For example, I notice very long terminal branch lengths in some 2.MED and 0.PE genomes, which are probably responsibly for skewing the clock rate on these lineages. How are possible artefactual variants (SNPs) evaluated and accounted for in this study?
• The authors mention that "The inability to infer divergence dates due to sampling bias also has several historical implications." Could informative dates be deduced from 0.ANT lineages when co-analysed with their ancient representatives? The possibility of such an analysis should be evaluated with BETS. Moreover, the authors mention, "Perhaps the most significant concerns the emergence of plague in Africa which makes up 90% of all modern plague cases, yet for which there remains not a single ancient sequence." Given the inability to estimate realistic dates for 1.ANT alone, could these be better estimated when considering all Branch 1 genomes as a single group? Is such a setup supported through model-testing (BETS approach)?
• Within the section "Informative Rates and Dates" the resulting maximum clade credibility trees from all analyses should be presented. Moreover, the point of "origin" in each case should be indicated on the trees with its corresponding 95% HPD interval.
• Regarding the estimated 1214-1315 date for the origin of 1.PRE, what is the genetic distance of the most basal genome on this lineage from the "point of origin"? In Figure S6, it seems that the MRCA of 1.PRE is placed at a zero genetic distance from the most basal genome. What is the isolation location and archaeological date of this genome? Based on the cited literature, the earliest Second Plague Pandemic genomes are dated to the 14th century. If the observation of zero genetic distance to the MRCA is correct, what does an origin date of 1214 mean and is it realistic? Given the important historical implications discussed by the authors, I emphasize that feasibility should not be confused with statistical uncertainty and that the resulting dates should be discussed in a more balanced and critical manner in the Results and Conclusions sections.
• Similarly to the comment above, what is the genetic distance of the most basal genome from the MRCA of the 0.ANT4 lineage? Given the estimate of one mutation every ~9.5 years for Y. pestis, how realistic is the result of a 272 CE origin date? In terms of interpretation, does the origin of a pandemic always coincide with the origin of its associated lineage? Again, the authors should discuss these aspects in a more balanced and critical manner within the Results and Conclusions sections.
• Other studies using different approaches have reached comparable results when estimating the date intervals for the emergence of the First and Second Pandemic lineages. Such studies should be credited within the Results and Conclusions sections.
• The authors mention "…we independently fit three discrete migration models to the maximum likelihood phylogeny using the sampling locations by: (1) continent, (2) country, and (3) province" Should time not be also considered as a factor for these migration models?
• Migration analysis: Why were the First Pandemic and the Bronze Age lineages not considered for this analysis? The latter might be a particularly good candidate as it shows an almost perfect temporal signal and is sampled from diverse locations.
• The authors mention "We estimated that the Second Pandemic (1.PRE) diverged from an ancestral population located in China (probability: 0.93) as part of the "Big Bang" polytomy." This statement is rather confusing as to what exactly it means. Figure 7 shows that it was the 0.ANT lineage that diverged in China, not the 1.PRE lineage. Also, how is the 1.PRE divergence separated from the beginning of the Second Pandemic? Are those not referring to the same divergence? Given the low support of these models, a clear description is important here. Moreover, what do the yellow arrows in Figure 7 indicate and why are they not described within the corresponding results section? If these are inferred migration events of higher support then I would suggest they are appropriately discussed. Alternatively, the authors should give a reasoning why these results cannot be trusted. If the results cannot be trusted with regard to migration events, but are in other ways informative for the study, I would suggest this entire analysis is moved to the supplement.
Minor points: • "The most intensively researched events have been: (1) the first appearance of Y. pestis in human populations" This statement includes an incorrect citation and should be updated with more recent studies.
• "(3) the inter-pandemic or "quiescent" periods where Y. pestis recedes into wild rodent reservoirs and disappears from the historical record" Ancient Y. pestis evidence from the 3rd century AD could be included here, as it is a good example of a historical genome that does not fall within the defined "pandemic" periods.
• "As a result, considerable debate has emerged over whether Y. pestis has no temporal signal" Here only one citation is included. Please include more to highlight different views on the debate.
• "This uncertainty has resulted in radically different rate and date estimates between studies, with node dates shifting by several millennia." Here the addition of basal lineages contributing to a tightening of the estimated emergence dates should be included as a factor.
• "This contention concerns competing hypotheses about the relative importance of localized persistence versus long-distance reintroduction." This sentence is presented in the context of the pandemic's origins, yet it seems more relevant to the pandemic's progression. This precise association is not clear here and should be further clarified.
• "Among both sides of this issue, there is an expectation that genomic evidence will play a significant role, if not resolve the debate" The indicated citations within this sentence seem inappropriate as they do not present new genomic evidence or help resolved the debate regarding the pandemic's origin.
• It is unclear what is meant by "with 3,844 sites shared by at least two strains". What does this information contribute?
• Regarding the observation "The Pestoides (0.PE) and Medievalis (2.MED) biovars are informative examples, as these populations have co-existed in the Caucasus mountains since at least the 20th century", are these different biovar perhaps identified in distinct hosts? In this case a table of isolation hosts for all the presented strains would be appreciated.
• "Phylogenetic analysis reveals genetic continuity between these two events, as the Third Pandemic (1.ORI) is a direct descendant of the Second Pandemic (1.PRE)." While this statement is true, one of the major lineages responsible for the Second Plague Pandemic, is one that currently includes the majority of data produced from this period and is today extinct. In that regard, the consideration of Second and Third pandemic strains as being genetically distinct, at least to an extinct, is not completely arbitrary and should be mentioned here.
• Figure 1: A higher resolution phylogeny would be useful as a supporting display item, where the phylogenetic positioning of each individual strain, with its respective name is shown.
• Figure 1: A Bronze Age genome associated with 0.PE is missing from the "Time Period" part of panel A.
• Figure 1: Why is the polytomy shown as a separate event from the beginning of the pandemic? An explanation should be provided.
• Figure 4 is a bit confusing as to what the boxes with dates indicate. A temporal indication of all depicted strains would make the figure clearer.
• Figure 5: An accompanying table with isolation hosts from all genomes used for this paper should be included.
• Figure 6+7: The migration analysis maps shown here should be presented with their corresponding phylogenetic trees, to aid visualisation and interpretation.
• Figure 6+7: I suggest a removal of all non-significant migration arrows from the main figures as they can be confusing and lead to false interpretations.
• There is a typo in the legend of Figure 7, 0.PRE should be 1.PRE.
• There is a typo in the Conclusions section "Plague of Justinian (531CE)" should be 541.

Comments
Through published genomic data, the authors re-analyzed the population specific evolutionary rates of Y. pestis, and evaluated the justification of phylogeographical analysis in Y. pestis using solely genomic evidences. The improved estimation on divergent time of key nodes on genealogy would help to correlate the historical events with evolution of Y. pestis in future studies.

Major
1. ✅There is strong temporal signal presented in 0.PRE population, of which sampling time period lasted for ~1500 years, however, in the similar time scaffold, the temporal signal disappeared in the whole species level. Could the authors discuss the possible reason on this?
The mechanisms that govern rate variation within Y. pestis are currently not well understood. The two major hypotheses are (1) adaptations to new environments and hosts, and (2) changing bacterial replication rates between epidemic and endemic cycles (Cui et al., 2013). To date, little support for either hypothesis has been found using genomic data. Our paper focuses on an earlier stage of this analysis, which is obtaining more reliable estimates of the rate variation, and critiquing their interpretation and application. We hope this manuscript will serve as a foundation for future work to test which variables explain (or do not explain) the populationspecific patterns observed here.
2. ✅There are lot of strains revealed ultra long branch length (grey branch in the phylogeny). Are there any possible that these long branch strains were caused by poor assembly quality, and how can these outliers influence the estimation of rate variations?
We have added a new  (Cui et al., 2013).

✅Table 1, it would help a lot if the authors clarify the rule on definition of different categories.
The caption to Table 1 has been extended to include definitions for each category. In addition, categories were re-named to match the in-text paragraph header where they are discussed in detail.

✅Fig. 5, I like this plot-very informative, but could the author add a panel that use all genomes (including ones from human), as control?
The data underlying Figure 5 includes all genomes used in this study (N=601, including  We thank the reviewer for catching this and have adjusted the legend of Figure 7.
6. ✅Method, Sequence Alignment section, as author claimed they only used assembled genome in the analysis, it is quite weird to emphasize parameters like "depth of 10X, base quality of 20, mapping quality of 30, major allele frequency of 0.9".
To align assemblies to the reference, we use the snippy pipeline which first shreds contigs into 250 bp single-end reads before using a short-read aligner (BWA) for mapping. We have added this explanation to the section Materials and Methods: Sequence Alignment.

Comments
In this study, entitled "Plagued by a cryptic clock: Insight and issues from the global phylogeny of Yersinia pestis", Eaton and colleagues explore published Y. pestis genomic diversity to assess substitution rate variations across different phylogenetic lineages and their impact on the global clock rate. Using a published dataset of ~600 genomes, the authors identify an absence of temporal signal across the entire Y. pestis phylogeny and rather propose a lineagespecific approach for molecular dating. Moreover, the authors identify sufficient phylogenetic signal in five "populations" that encompass wider temporal sampling ranges. Individual molecular dating analyses were then used to re-estimate the dates of origin for the respective "populations" and to discuss possible historical implications of their inferred divergence dates. Finally, the authors explore possible dissemination patterns in relation to historical epidemics/pandemics, and identify significant migration events during the Third Plague Pandemic, which coincide with documented historical accounts.
While I do appreciate the authors effort towards a comprehensive characterisation of all Y. pestis clades using a methodology that has not been presented before, I have a number of major and minor comments regarding the analytical approach and interpretation of results. Moreover, I find that relevant supporting display items are missing for fully assisting the authors arguments. My specific comments can be found below.

✅The authors should provide a table of all 601 genomes used in this study together with their accession codes, isolation dates, geographic coordinates and isolation hosts in order to ensure transparency in the used dataset and aid reproducibility of the analysis.
Genomic metadata is provided in Dataset 1: Supplementary Table S8. 2. ✅Is there a specific reason why the oldest documented Y. pestis lineages were not included in this study? Currently, they are represented by the genomes Gok2 or RV2039. I suggest their addition as they might help improve the phylogenetic signal across all Y. Pestis.
Gok2 (Rascovan et al., 2019) was excluded as it did not meet our coverage filter of 70% of the genome covered at a minimum depth of 3X. RV2039 (Susat et al., 2021) was not used in this study as it was published on 2021 June 29, which post-dated our data collection period (2020 January 01). We have updated the section Methods: Data Collection to clarify our data collection cut-off.
3. ✅One of the determining priors in Bayesian molecular dating analysis is the demographic model used. Yet, no testing of demographic models is performed here and the authors choose to proceed with a coalescent constant population size model that seems inappropriate for describing the evolutionary history of epidemic pathogens. The suitability of other models, such as the coalescent skyline and exponential growth models should be fully evaluated for the entire dataset, as well as all the subsets.
This is an important point and we have explained our reasoning for our choice of tree prior in section Methods: estimating rates of evolutionary change.
We do agree with the reviewer that the constant coalescent tree prior is overly simplistic for these data. However, the only models appropriate here are those that account for population structure (e.g. structured coalescent or multitype birth-death), which are computationally prohibitive for our data. For this reason we chose a model that was as simple as possible and expected to have little impact in estimates of node times when the data are informative (see Möller et al., 2018;Ritchie & Ho, 2019

Artefactual variants are one possible explanation for the inflated evolutionary rate estimates in 2.MED and 0.PE. However, the relationship between the proportion of long branches and the rate variation in a population is also unclear. The populations with the greatest number of long branches (pestoides 0.PE and medievalis 2.MED) both had detectable temporal signal according to BETS. But of the two, informative rates and divergence dates could only be obtained from medievalis 2.MED. This suggests that numerous long-branches (possibly from artifactual SNPs)
does not necessarily preclude recovery of temporal sigal.

✅The authors mention that "The inability to infer divergence dates due to sampling bias also has several historical implications." Could informative dates be deduced from 0.ANT lineages when co-analysed with their ancient representatives? The possibility of such an analysis should be evaluated with BETS. Moreover, the authors mention, "Perhaps the most significant concerns the emergence of plague in Africa which makes up 90% of all modern plague cases, yet for which there remains not a single ancient sequence." Given the inability to estimate realistic dates for 1.ANT alone, could these be better estimated when considering all Branch 1 genomes as a single group? Is such a setup supported through model-testing (BETS approach)?
We agree that experimenting with population definitions (combining/splitting) is the intuitive next step in improving molecular clock analyses in Y. pestis. Similar to our response to Reviewer 2's Major point #5, however the complexity of the experimental design and scale of that analysis (i.e. testing all possible population combinations) is outside the scope of this work and would be better served as a standalone publication.
9. ✅Within the section "Informative Rates and Dates" the resulting maximum clade credibility trees from all analyses should be presented. Moreover, the point of "origin" in each case should be indicated on the trees with its corresponding 95% HPD interval.  We have added additional tMRCA estimates to Table 1 where they have been reported (Rasmussen et al., 2015;Spyrou et al., 2019). However, two studies (Namouchi et al., 2018), (Spyrou et al., 2018) do not provide dates. Namouchi et al. don't report any time estimates due to a lack of temporal signal in their dataset and Spyrou et al. 2019, only report the estimate of the origin of all plague, which we didn't not attempt and as such it is not directly comparable.We have left these as NA in the Table 1. 13. ✅The authors mention "…we independently fit three discrete migration models to the maximum likelihood phylogeny using the sampling locations by: (1) continent, (2) country, and (3) province" Should time not be also considered as a factor for these migration models?

We have created a new supplementary dataset (Dataset 2) which contains tree files (newick and nexus) as well as a formatted metadata file (tsv) for visualization. In addition, we have created a https://nextstrain.org/ community page to host interactive versions of all of our phylogenies
We agree that time is an important factor and is the preferred approach (i.e. jointly estimating the phylogeny, molecular clock, and geographic location). However, as we demonstrate here, determining a molecular clock model for Plague is difficult due to dramatic variation in the evolutionary rate and population dynamics. By avoiding a molecular clock model our approach effectively results in an 'overdispersed clock', where each branch is allowed an independent rate. We argue that this is the most conservative approach for phylogeography of this organism, until a more realistic molecular clock model is defined. We explained this point in the corresponding section: "We chose this approach, over one where time is coestimated (i.e. a molecular clock) because the ages of some samples have uncertainties associated with them and due to the large rate variation across the entire tree, which would be a large source of error in migration rates and events over time".
14. ✅Migration analysis: Why were the First Pandemic and the Bronze Age lineages not considered for this analysis? The latter might be a particularly good candidate as it shows an almost perfect temporal signal and is sampled from diverse locations.  Figure s1 illustrates the proximity of these two nodes and the clear evidence for ancestral lineages being predominantly based in China. The 1.PRE divergence is estimated back to a MRCA and the confidence around that estimate is given as the 95% CI.
16. ✅Are those not referring to the same divergence? Given the low support of these models, a clear description is important here. Moreover, what do the yellow arrows in Figure 7 indicate and why are they not described within the corresponding results section? If these are inferred migration events of higher support then I would suggest they are appropriately discussed. Alternatively, the authors should give a reasoning why these results cannot be trusted. If the results cannot be trusted with regard to migration events, but are in other ways informative for the study, I would suggest this entire analysis is moved to the supplement.
We now focus our discussion of migration events to those with high statistical support and only arrows with such high support are shown in Figures 6 and 7. We also added some text to explain why we only focus on those migration events with high support: "Importantly, we focus our interpretation of migration events to those that have high statistical support, with high bootstrap (topology inference) and discrete trait reconstruction (phylogeographic inference). In practice, low statistical support in these events means that the data are not sufficiently informative about migration pathways, and thus their interpretation can be misleading." Minor 1. ✅"The most intensively researched events have been: (1) the first appearance of Y. pestis in human populations" This statement includes an incorrect citation and should be updated with more recent studies.
We have updated the citation of "(1) the first appearance of Y. pestis in human populations" to (Susat et al., 2021).

✅"
(3) the inter-pandemic or "quiescent" periods where Y. pestis recedes into wild rodent reservoirs and disappears from the historical record" Ancient Y. pestis evidence from the 3rd century AD could be included here, as it is a good example of a historical genome that does not fall within the defined "pandemic" periods.
We have added a citation to the statement "(3) the inter-pandemic or "quiescent" periods where Y. pestis recedes into wild rodent reservoirs and disappears from the historical record" (Damgaard et al., 2018) to reflect the 3rd century strain DA101.
3. ✅"As a result, considerable debate has emerged over whether Y. pestis has no temporal signal" Here only one citation is included. Please include more to highlight different views on the debate.
We have updated the citations for this sentence in the Introduction, to include four sources that characterize the viewpoints that (1) Y. pestis has species-wide temporal signal, (2) Y. pestis has population-specific temporal signal, and (3) Y. pestis has no temporal signal.

4.
✅"This uncertainty has resulted in radically different rate and date estimates between studies, with node dates shifting by several millennia." Here the addition of basal lineages contributing to a tightening of the estimated emergence dates should be included as a factor.
We have edited that sentence in the Introduction to indicate that divergence dates are shifting and also narrowing, largely due to additional ancient Y. pestis genomes.

5.
✅"This contention concerns competing hypotheses about the relative importance of localized persistence versus long-distance reintroduction." This sentence is presented in the context of the pandemic's origins, yet it seems more relevant to the pandemic's progression. This precise association is not clear here and should be further clarified.
We have edited the top sentence of that paragraph to "The geographic origins and progression of past pandemics are similarly contentious…" to expand the context beyond just the pandemic's origins.
6. ✅"Among both sides of this issue, there is an expectation that genomic evidence will play a significant role, if not resolve the debate" The indicated citations within this sentence seem inappropriate as they do not present new genomic evidence or help resolved the debate regarding the pandemic's origin.
We have adjusted the positioning of the citations within that Introduction sentence. We hope to convey that it is not that the cited authors have solved the debate, but rather that they have stated in their articles that they hope/expect additional genomic data to solve it in the future. 9. ✅"Phylogenetic analysis reveals genetic continuity between these two events, as the Third Pandemic (1.ORI) is a direct descendant of the Second Pandemic (1.PRE)." While this statement is true, one of the major lineages responsible for the Second Plague Pandemic, is one that currently includes the majority of data produced from this period and is today extinct. In that regard, the consideration of Second and Third pandemic strains as being genetically distinct, at least to an extinct, is not completely arbitrary and should be mentioned here. 10. ✅Figure 1: A higher resolution phylogeny would be useful as a supporting display item, where the phylogenetic positioning of each individual strain, with its respective name is shown.
A high-resolution phylogeny with tips labelled as "Country (Date) Strain" has been included as Supplementary Figure S9. An interactive version is also available at: https://nextstrain.org/community/ktmeaton/yersinia-pestis/maximum-likelihood/all 11. ✅Figure 1: A Bronze Age genome associated with 0.PE is missing from the "Time Period" part of panel A.
Thank you for catching this error. We have updated Figure 1: Panel A to properly label the time period of the Bronze Age 0.PE genome.
12. ✅Figure 1: Why is the polytomy shown as a separate event from the beginning of the pandemic? An explanation should be provided.
The maximum-likelihood phylogeny displayed in Figure 1: Panel A, is a strictly bifurcating tree with no true polytomies. Because of this, the internal node representing the MRCA of Branches 1,2,3 and 4 is distinct from the MRCA of Second Pandemic strains (yellow star). We have updated the caption of Figure 1 to clarify the tree is strictly bifurcating.
13. ✅Figure 4 is a bit confusing as to what the boxes with dates indicate. A temporal indication of all depicted strains would make the figure clearer.
We have removed the date boxes from Figure 4 as the focus of this 14. ✅Figure 5: An accompanying table with isolation hosts from all genomes used for this paper should be included. 15. ✅Figure 6+7: The migration analysis maps shown here should be presented with their corresponding phylogenetic trees, to aid visualisation and interpretation.

Information about isolation hosts can be found in Supplementary
We pointed the reader back to Figure 1, which includes the complete tree, to avoid cluttering these figures (already very compressed). Moreover, the current version of Figures 6 and 7 only include migration events with high support, and should therefore be more informative than the trees.
16. ✅Figure 6+7: I suggest a removal of all non-significant migration arrows from the main figures as they can be confusing and lead to false interpretations.
Thank you for this visualization suggestion, we have removed all non-significant migration arrows from Figure 6 and Figure 7. 17. ✅There is a typo in the legend of Figure 7, 0.PRE should be 1.PRE.
Thank you for catching this error. We have updated the legend of Figure 7 to indicate 1.PRE instead of 0.PRE.
18. ✅There is a typo in the Conclusions section "Plague of Justinian (531CE)" should be 541.
Thank you for catching this error. We have updated the Conclusions to have the correct date (541 CE).