Author Response:
Reviewer #1:
Summary:
Moody et al. presented a comprehensive investigation into the choice of marker genes and its impact on the reconstruction of the early evolution of life, especially regarding the length of the branch that separates domains Bacteria and Archaea in the phylogenetic tree. Specifically, this work attempts to resolve a debate raised by a previous work: Zhu et al. Nat Commun. 2019, that the evolutionary distance between the two domains is short as estimated using an expanded set of marker genes, in contrast to conventional strategies which involve a small number of "core" genes and indicate a long branch.
Through a series of analyses on 1000 genomes, Moody et al. defended the use of core genes, and reinforced the conventional notion that the inter-domain branch (the AB branch) is long, as inferred by the core gene set. They proposed that with the 381 marker genes (the "expanded" set) used by Zhu et al., the observed short branch length is an artifact due to inter-domain gene transfer and hidden paralogy. Through topology tests, they ranked the markers by "verticality", and showed that it is positively correlated with the AB branch length. They also conducted divergence time estimation and showed that even the most vertical genes led to an implausible estimate of the origin of life.
In parallel, Moody et al. surveyed the best marker genes using a set of 700 genomes. They recovered 54 markers, and demonstrated that ribosomal markers do not indicate a longer AB branch than non-ribosomal markers do. With the better half (27) of these marker genes, they conducted further phylogenetic analyses, which shows that potential substitutional saturation and the use of site-homogeneous models could contribute to the underestimation of the AB branch. Using this taxon set and marker set, they reconstructed the prokaryotic tree of life, which revealed a long AB branch, a basal placement of DPANN in Archaea, and a derived placement of CPR in Bacteria.
Prokaryotic tree of life:
The scope(s) of the manuscript is somehow split. First, it is posed as a point-to-point rebuttal to the Zhu et al. paper, on the long vs. short AB branch question. Second, it introduces a new phylogeny of prokaryotes using 27 "good" marker genes, and demonstrates that DPANN is basal to Archaea, and CRP is derived within Bacteria.
Thanks for the summary. The two aspects of the manuscript identified by the reviewer are closely related, because the different issues boil down to the same underlying question: which genes should we use to infer the deep structure of the tree of life? The provocative work of Zhu et al. acted as an impetus to compare and evaluate the properties of several published marker gene sets, and then to identify (what our analyses suggest are) the subset best-suited for deep phylogeny, which we then use to infer an updated tree of life. We have clarified this logical structure in the revised manuscript, writing (at the end of the Introduction):
“Here, we investigate these issues in order to determine how different methodologies and marker sets affect estimates of the evolutionary distance between Archaea and Bacteria. First, we examine the evolutionary history of the 381 gene marker set (hereafter, the expanded marker gene set) and identify several features of these genes, including instances of inter-domain gene transfers and mixed paralogy, that may contribute to the inference of a shorter AB branch length in concatenation analyses. Then, we re-evaluate the marker gene sets used in a range of previous analyses to determine how these and other factors, including substitutional saturation and model fit, contribute to inter-domain branch length estimations and the shape of the universal tree. Finally, we identify a subset of marker genes least affected by these issues, and use these to estimate an updated tree of the primary domains of life and the length of the stem branch that separates Archaea and Bacteria.”
The second scope has inadequate novelty. A recent paper (Coleman et al. Science. 2021), which was from a partially overlapping group of authors, was dedicated to the topic of CPR placement, and indicated the same conclusion (CPR being derived and sister to Chloroflexi) as the current work does, albeit using more sophisticated approaches. The paper also addressed the debate of CPR placement (including citing the Zhu et al. paper). Additionally, the basal placement of DPANN has also been suggested by previous works (such as Castelle and Banfield. Cell. 2018). Therefore, re-addressing these two topics using a largely well-established and repeatedly adopted method on a relatively small taxon set does not constitute a significant extension of current knowledge.
We disagree. Resolving the deep structure of the tree of life is an important topic --- this is what we, Zhu et al. (2019), and of course many others have been trying to achieve, in different and sometimes conflicting ways. Most of the published work is based on limited or biased taxon sampling (see Figure 1 Figure Supplement 14,15,16) or else focused on just one of the two prokaryotic domains of life. Furthermore, deep phylogeny is uncertain, and new results become convincing only when they receive support from multiple datasets and approaches. For instance, Coleman et al. (2021) recently found support for the placement of CPR as a sister clade to Chloroflexota rather than as a basal branch within the Bacteria. Notably, this work focused only on Bacteria, and made use of a different rooting method (with its own strengths and limitations) and taxon sampling. Most previous analyses using Archaea as an outgroup to root the bacterial tree recovered CPR as a deeply branching lineage within Bacteria, a placement likely resulting from LBA. In turn, our present findings represent an important confirmation of the CPR+Chloroflexi clade. Similarly, the basal placement of DPANN within Archaea remains controversial despite a number of studies on the topic, and our study also contributes to that ongoing debate.
The debate:
The first scope appears to be the more important goal of this manuscript, as it extensively discusses the claims made by Zhu et al. and presents a point-to-point rebuttal, including counter evidence. This may narrow the interest of this work to a small audience of specialists. Nevertheless, to best evaluate the current work, it is necessary to review the Zhu et al. paper and compare individual analyses and conclusions of the two studies.
In doing so, I found that the two articles have distinct scopes that appear similar but not actually inline. To a large extent, the current work does not constitute actual rebuttal to the points made by Zhu et al. In contrast, some of the analyses presented in the current work support those by Zhu et al., despite being interpreted in a different way. For the claims that directly contest Zhu et al., I do not see sufficient evidence that they are supported by the analyses.
Below is a summary of the comparison, which I will explain point-by-point in later paragraphs.
- Moody et al. assessed AB branch length, while Zhu et al. assessed AB evolutionary distance (which is different).
- Moody et al. evaluated the phylogeny indicated by a small number of core markers, while Zhu et al. evaluated the genome average using hundreds of global markers.
- Zhu et al.'s results also showed that gene non-verticality, substitutional saturation, and site-homogeneous models shorten the AB distance, which is consistent with Moody et al.'s.
- However, Zhu et al. found that some core markers are outliers in the genome-wide context, and the long AB distance indicated by them cannot be compensated for by the aforementioned effects. Moody et al. hasn't addressed this.
Therefore, the novelty and potential impact of the current work is less compelling: It used a classical method (a few dozen core genes) and found a pattern that has been found many times by some of the same authors and others (including Zhu et al., who also analyzed core genes).
Thanks for this detailed comparison of the two studies --- the points raised here and elaborated on below have prompted us to perform additional analyses which provide further insight into the properties and behaviour of the various marker gene sets analyzed. We nonetheless disagree that “the current work does not constitute actual rebuttal to the points made by Zhu et al.”: our finding that ribosomal and other “core” proteins are among the best phylogenetic markers for resolving both within- and between-domain relationships, estimating the length of the AB stem, and performing divergence time estimation, challenges an important claim of Zhu et al.’s study, and will be of broad interest to the community of researchers working on early life/early evolution.
That said, we do also agree that one aspect of the disagreement between our study and that of Zhu et al. has to do with what is meant by evolutionary distance, and we have now discussed these issues in detail in the revised manuscript (as detailed below). In revising the manuscript, we have also sought to avoid a reductive focus on rebuttal, have revised the text to acknowledge important strengths and interesting features of the Zhu et al. analyses, and have made text revisions to ensure a consistent constructive tone: these are fundamental and challenging questions, and different perspectives and analyses are valuable in making
progress. We also note that there has been an ongoing debate about the suitability of ribosomal genes for deep phylogeny in the literature (e.g. Petitjean et al. 2014, discussed in more detail below). Our analyses, and those of Zhu et al. (2019) previously, contribute to that broader discussion.
Detailed responses to each of the above points follow below.
AB distance metric:
There is a subtle but critical difference between the scopes of the two papers: The Zhu et al. paper "reveals evolutionary proximity between domains Bacteria and Archaea". By stating "evolutionary proximity", it investigated two metrics:
The length of the branch separating Archaea from Bacteria in the phylogenetic tree, i.e., the "AB branch". This was the main focus of the current work.
The average tip-to-tip distance (sum of branch lengths) between pairs of Archaea and Bacteria taxa in the tree. A significant proportion of the Zhu et al. work was discussing this metric, and it led to several important conclusions (e.g., Figs. 4F, 5). The current work has not explored this metric.
Thanks for raising the point about relative AB distance. In our revised manuscript, we have expanded Figure 1 and the associated analyses to include this metric. These analyses demonstrate that relative AB distance behaves similarly to AB branch length: they are positively correlated with each other; both are reduced by inter-domain HGT, and both are negatively correlated with ΔLL and with split score, an additional metric of within- and between-domain marker gene verticality which we have included in the revised Figure 1. Taken together, these results suggest that high-verticality marker genes (as judged both by the recovery of reciprocal AB monophyly, and of established within-domain relationships) support a longer AB branch and show a higher relative AB distance.
These two metrics implicate distinct research strategies: For 1), HGTs and paralogy are usually considered problematic (as the current and many previous works argued). However, 2) is naturally compatible with the presence (and prevalence) of HGTs and paralogy.
Authors of the current work equate "genetic distance" to "branch length" (line 70), and only investigated the latter. This equation is misleading. If organism groups A and B diverged early, but then exchanged many genes post-divergence, then this is indisputable evidence that their "genetic distance" is close. This point needs to be clearly explained in the manuscript.
We agree with the reviewer that various definitions of evolutionary distance are possible, and some may be more useful than others for particular applications. The reviewer’s argument that “If organism groups A and B diverged early, but then exchanged many genes post-divergence, then this is indisputable evidence that their "genetic distance" is close” makes the case for a kind of phenetic distance: a distance based on overall similarity, regardless of how that similarity was brought about in terms of evolutionary process. We appreciate the democratic appeal of such a metric, and we have no desire to impose any particular philosophy of classification on the reader. However, the key point here is that methods that rely on concatenation for branch length or divergence time estimation (as used by Zhu et al., and in our current study) make the assumption that all of the sites in the concatenate evolved on the same underlying tree and if this assumption is not met, analyses can be misled. Thus, the shorter AB branch length and the more recent Archaea-Bacteria divergence times estimated from concatenations of incongruent marker genes result from unmodelled gene transfers which are misinterpreted as evidence for more recent common ancestry. Gene transfer is an important aspect of genome evolution, but none of the currently available methods, including those used by Zhu et. al., allow for genome-scale comparisons to be made in a way that accounts for our understanding of the underlying evolutionary processes.
The point about different possible definitions of evolutionary distance made by the reviewer is valid, and we have now revised the opening of our conclusion to discuss these issues in more detail, writing:
“We note that alternative conceptions of evolutionary distance are possible; for example, in a phenetic sense of overall genome similarity, extensive HGT will increase the evolutionary proximity (Zhu et al., 2019) of the domains so that Archaea and Bacteria may become intermixed at the single gene level. While such data can encode an important evolutionary signal, it is not amenable to concatenation analysis.”
Core vs genome:
This difference between "AB distance" and "AB branch length" is relevant to a more fundamental question: What defines the "evolutionary distance" between two groups of organisms? Both papers did not explicitly discuss this topic. It likely cannot be resolved in one article (as many scholars have continuously attempted on related topics in the past decades). But the discordance in understanding led to very different research strategies in the two papers, and rendering them incongruent in methodology.
Specifically, the current work (and multiple previous works) based phylogenetic inference on only genes that demonstrate a strong pattern of vertical evolution. HGTs were considered deleterious, and needed to be excluded from the analysis. This left a few dozen genes at most, and many are spatially syntenic and functionally related (e.g. ribosomal proteins). In this work, the final number is 27. Previous critiques of this methodology have suggested that this is not a tree of life, but a "tree of one percent" (Dagan and Martin, Genome Biol. 2006).
In contrast, Zhu et al. (and related previous works) attempted to evaluate the evolution of whole genomes by "maximizing the included number of loci.". They used a "global" set of 381 genes. They faced the challenge of "reconciling discordant evolutionary histories among different parts of the genome", because "HGT is widespread across the domains". To resolve this, they adopted the gene tree summary method ASTRAL.
Therefore, the "AB distance" estimated by Zhu et al. is a genome-level distance, calculated by merging conflicting gene evolutions (which itself can be disputed, see below). Whereas the "AB branch" evaluated in this work is strictly the branch length in the core gene evolution. Therefore, the results presented in the two papers do not necessarily conflict, because of the different scopes.
This point is closely related to the previous one, and the new section (final paragraphs of the Conclusion, quoted directly above) goes some way to addressing this comment. Regarding the issue of a focus on just a small proportion of vertically-evolving genes, the critical point is as above: current methods for branch length and divergence time estimation (including those used by Zhu et al.) require such vertically-evolving genes, because they make the assumption that all of the sites evolve on the same tree, i.e. trace back to the same origin via vertical evolution. We agree that most prokaryotic gene families do not evolve under these restrictive assumptions and therefore cannot be analysed using concatenation methods for branch length estimation. Indeed, one of the main points of our study is that most of the genes in the 381-gene set of Zhu et al. do not meet these assumptions and are thus unsuited for estimating evolutionary distance and divergence times.
There is much ongoing method development which will allow more of the genome to be used in deep-time comparative analyses; Astral-Pro, FastMulRFS and SpeciesRax, among others, are recent promising steps in this direction. However, our central critique of Zhu et al. is that inferences under concatenation-based methods can be misled by HGT and other sources of incongruence, and indeed our analyses show that these unmodelled signals underlie the difference between the conclusions of Zhu et al. and other studies (e.g. (Liu et al., 2021; Spang et al., 2015; Williams et al., 2020) that have instead supported a deep divergence between Archaea and Bacteria. In our revised manuscript, we have shown that the relative AB distance, like the AB branch length, is shortened by unmodelled gene transfers (Figure 1), and that estimates of the AB stem length from different studies are similar when the congruent subset of the data is analysed with the best available substitution models (Figure 6). We therefore disagree that the scopes are distinct: richer, broader measures of genomic diversity can be proposed and, with the development of new methods, estimated; but so far, the vertical signal is the only signal that can be harnessed to infer divergence times using concatenations.
The expanded marker set:
The authors made a valid critique (line 121-135) that many of the 381 genes in the "expanded marker set" adopted by Zhu et al., are under-represented in Archaea. According to the PhyloPhlAn paper (Segata et al. Nat Commun. 2013) which originally developed the 400 markers (a superset of the 381 markers), these genes were selected from ~3,000 bacterial and archaeal genomes available in IMG at that time time (note that it was 2013). Zhu et al. also admitted, in the discussion section, that this marker set falls short in addressing some questions (such as the placement of DPANN). What is important in the current context, is that they were not specifically selected to address the AB distance question.
We agree that the taxon sampling of archaea and the choice of marker genes in the Zhu et al. study were not ideal for estimating the evolutionary distance between the domains. However, we note that this distance (or proximity), and the hypothesis that traditional core genes over-estimate the Archaea-Bacteria divergence, was one of the main results of the paper (c.f. the title of that paper, “Phylogenomics of 10,575 genomes reveals evolutionary proximity between domains Bacteria and Archaea”).
However, note that Zhu et al.'s Fig. 5A, B presented the AB distance informed by 161 out of the 381 genes. These genes have at least 50% taxa represented in both domains - the same threshold discussed in the current work (line 132).
While the 50% sampling criterion indeed enriches for the genes of the expanded set that were present in LUCA and on the AB branch, we note that the 50% criterion represents a minimum of 4953 bacteria and 335 archaea; that is, it still reflects the unbalanced sampling of the dataset overall. For example, 30 of the genes had fewer than two archaeal homologues, and in 100 of the trees there were fewer than 50 archaea reflecting the large disparity in taxon sampling (Supplementary Information Table S1). The phylogenetic signal in these genes is discussed in more detail below. Looking at the subsampled versions of these 161 genes, we found the majority of these genes (123/161) to have no discernible AB branch length. The 38/161 genes which had an arguable AB branch length (but still with transfers/paralogs) possessed a range of AB lengths: 0.0814:5.26, with a mean AB length of 1.03 and a median of 0.635.
Even with those sufficiently represented genes, they still found that ribosomal proteins and a few other core genes are "outliers" in the far end of the AB distance spectrum.
The reviewer raises an interesting point about outliers with high relative AB distances, which gets to the heart of the debate about how best to estimate the evolutionary distance between Archaea and Bacteria. The new analyses of relative AB distance introduced in our revised manuscript (Figure 1) demonstrate that this metric is affected by HGT in a similar manner to AB branch length (that is, high-verticality marker genes have a greater relative AB distance (relative AB vs ΔLL: p = 0.0001051 & R = -0.2213292, relative AB vs between-domain split score: p = 2.572e-06 & R = -0.2667739). Thus, core genes can be viewed as “outliers” compared to other prokaryotic genes in the sense that they have experienced an unusually low amount of HGT. This high verticality makes them among the few prokaryotic gene families that can be analysed by concatenation methods, which make the assumption that all sites evolve on the same underlying tree topology.
Domain monophyly in gene trees:
The authors' efforts in manually checking the gene trees are appreciable (Table S1), considering the number and size of those trees. They found (line 147) "Archaea and Bacteria are recovered as reciprocally monophyletic groups in only 24 of the 381 published (Zhu et al., 2019) maximum likelihood (ML) gene trees of the expanded marker set."
The domain monophyly check was valid, however the result could be misleading because any sporadical A/B mixture was considered evidence of non-monophyly for the entire gene tree. As the taxon sampling grows, the opportunity of observing any A/B mixture also increases. For example, in Puigbò et al. J. Biology. 2009, 56% (a much higher ratio) of nearly universal genes trees had perfect domain monophyly based on merely 100 taxa. This is because even the "perfect" marker genes (such as ribosomal proteins) are not completely free from HGTs (e.g., Creevey et al. Plos One. 2011), let alone the fact that there are many artifacts in the published reference genomes (Orakov et al. Genome Biol. 2021).
Therefore, to have an objective assessment of this topic, it would be better to have a metric that allows some imperfection and reports an overall "degree" of separation (also see below).
We agree that complementing the monophyly check with a more nuanced metric is useful. In our revised manuscript, we now also evaluate the split score (Dombrowski et al. (2020) Nat Commun) of each marker, which reflects the degree to which a gene recovers the monophyly of established taxonomic ranks (a higher score reflects the splitting of monophyletic groups into a number of smaller clades in the gene tree, and so the metric permits a degree of “imperfection”, as suggested; in addition, the metric is averaged over bootstrap replicates, so that lack of resolution or poorly-supported disagreements with the reference taxonomy do not disproportionately affect the score). This expanded analysis (Figure 1) indicates that both within- and between-domain split score and ΔLL are significantly positively correlated (R = 0.836679, p < 2.2✕10-16), and that phylogenetic markers that more strongly reject domain monophyly (higher delta-LL) also perform worse at recovering between-domain (and within-domain) relationships (higher split score) and support a shorter AB branch length.
AB branch by gene: correlation and outliers
Figure 1 is the single most important result in this work, because it argues that the short AB branch observed in Zhu et al. is an artifact due to "inter-domain gene transfer and hidden paralogy" (line 202). This argument is based on the observation that the indicated AB branch length is negatively correlated with "verticality" (measured by ΔLL and split score) of the gene.
Our argument that the short AB branch results from inter-domain gene transfer and hidden paralogy is based on three main lines of evidence: (i) documentation of extensive transfers and intermixing of paralogues in the gene trees for the 381 gene set; (ii) the analyses in Figure 1, which demonstrate that verticality positively correlates with AB branch length and AB distance; (iii) the demonstration that the incremental addition of low-verticality markers to a concatenate results in a concomitant decrease in AB branch length.
However, Zhu et al. also investigated the impact of verticality on AB distance, and they also found that they are negatively correlated (Fig. 5E). Therefore, the current result does not appear to deliver new information (as do multiple other analyses, see below).
Zhu et al. indeed identified a weak positive relationship between gene verticality and AB distance. Our analyses go beyond that work by showing, using a variety of complementary metrics of verticality, that AB branch length and relative AB distance are strongly positively correlated with verticality (see Figure 1), and that the low verticality of the genes in the 381 gene set largely explains the difference in stem length inference between that dataset and earlier analyses (Figure 6). An additional factor not considered in the analyses of Zhu et al. was the question of whether a gene was present in LUCA, and so can provide information on the AB branch length. Our analyses (detailed below) suggest that the majority of genes (317) in the 381 gene set do not contain an unambiguous AB branch, and so do not contribute interpretable signal to estimates of the AB branch length.
An important finding in Zhu et al., which is largely not discussed in the current work, is that a handful of "core" genes are outliers in the spectrum of AB distance, as compared to the majority of the genome (Fig. 5A). The AB distance indicated by these core genes is so long compared with the genome average that it cannot be compensated for by the impact of non-verticality, substitutional saturation, site-homogeneous model, etc (see below).
Fig. 1A of the current work also clearly shows that many long-AB branch genes are outliers compared with the majority of the genome (the bottom of the blue bar).
Figs. 3 and 4 attempted to show that ribosomal proteins are not outliers, but that analysis was based on a very small set of core genes, and the figures clearly show that there are outliers even in this small set (to be further discussed below).
This comment re-iterates the reviewer’s earlier points about “core” genes as outliers compared to the majority of the genome. The key issue is that “most of the genome”, and a significant portion (317 genes) of the 381-gene set, contain features that make them unsuitable for estimation of AB branch length by concatenation, or indeed estimation of an interpretable relative AB distance. We have documented the cases of HGT and mixing of paralogues in the 381-gene dataset; this information is summarised in the main text and presented in more detail in Supplementary Information Table S1.
Focusing on the 161 genes with >50% representation in both Archaea and Bacteria, manual inspection of gene trees inferred on the 1000-species subsample under the LG+G+F model indicate that 123/161 do not have a clear AB branch (that is, a branch that separates most or all Archaea from Bacteria). While distinguishing such cases from early gene transfers is not straightforward, there is no compelling reason to think that these genes were present in LUCA. The simplest explanation for these gene phylogenies is instead an origin within Bacteria and subsequent transfer on one or multiple occasions into Archaea. As a result, estimates of AB branch length or relative AB distance inferred from these genes cannot be straightforwardly compared to those of the traditional “core” or other genes for which the evidence of a pre-LUCA origin is stronger. Considering only the 38/161 genes for which a LUCA origin appears, from the gene phylogeny, to be likely, the mean AB branch length is 1.03, greater than that estimated from the concatenation of the most vertical genes in the expanded set (0.56), and suggesting that phylogenetic incongruence, combined with (for some families) a more recent origin explains the shorter AB distances inferred from the 381 gene set. Thus, it is not the case that the AB branch lengths (or relative distances) estimated from the majority of genes form a null distribution” against which “core” genes can be seen as outliers; instead, our analyses suggest that “core” genes are among the limited number of genes that trace vertically to LUCA.
Regarding Figures 3 and 4, see the more detailed discussion below.
Verticality is not causative of short AB branch:
In spite of the outlier question, there is an important logic problem in these analyses: The authors observed that gene verticality (measured by negative ΔLL) is correlated with AB branch length (Fig. 1), and concluded that HGTs and paralogy shortened the AB branch (line 202). However, they did not directly assess the rate of evolution in this model. It is totally possible that the most vertical genes happen to be those that evolved faster at the AB split. In order to support the claim made in this work, it is important to separate the effect of the rate of evolution from the effect of HGT / paralogy.
The ideal solution would be to include ALL genes (not just "good" ones), build gene trees, identify parts of the gene trees that once experienced HGT or paralogy, and prune off these PARTS, instead of excluding the entire gene tree. The remaining data are thus free of HGT / paralogy, based on which one can quantify the "true" AB branch length, and further assess how much it is correlated with "verticality", and whether there are still "outliers". This solution is likely not trivial in implementation, though. However, without such assessment, the observed short AB branch still only applies to the "tree of one percent", not the "tree of life".
Thanks for this comment --- the reviewer raises a subtle and valid point. Our analyses indicate that vertically-evolving genes have longer AB branch lengths, but in the first version of our manuscript we did not test the alternative hypothesis that this relationship might simply result from a faster rate of evolution in vertically-evolving genes. To evaluate the relationship between evolutionary rate, verticality, AB branch length and relative AB distance on as broad a set of genes as possible, we took the 302 genes from the 381-gene expanded set, excluding 56 genes for which the 1000-species subsample included no archaea, and another 23 which included only 1 archaeon. To estimate per-gene evolutionary rate, we rooted each gene tree using MAD (Tria et al. 2017) and calculated the mean root- to-tip distance on the MAD-rooted gene tree, then evaluated the relationship between rate and verticality. This analysis indicated that vertically-evolving genes evolve more slowly (have shorter mean root-to-tip distances) than less vertical genes (using deltaLL and between-domain split score as proxies for marker verticality, with a Pearson’s product- moment correlation: MAD rooted mean root-to-tip distance against deltaLL: R = 0.1397803, p = 0.01506 or against split score: R = 0.1902056 p = 0.000893), despite having longer AB branches and relative AB distances (using a Pearson’s product moment correlation of MAD rooted mean root-to-tip distance against AB length: p = 0.2025, R= 0.1143076, or against relative AB distance p = 0.007435, R=0.1537479). Thus, the longer AB branches of vertically evolving genes do not appear to be the indirect result of faster evolution of those genes. These analyses are reported in the main text, where we write:
An alternative explanation for the positive relationship between marker gene verticality and AB branch length could be that vertically-evolving genes experience higher rates of sequence evolution. For a set of genes that originate at the same point on the species tree, the mean root-to-tip distance (measured in substitutions per site, for gene trees rooted using the MAD method (Tria et al., 2017)) provides a proxy of evolutionary rate. Mean root-to-tip distances were significantly positively correlated with ∆LL and between-domain split score (∆LL: R = 0.1397803, p = 0.01506, split score: R = 0.1705415 p = 0.002947; Figure 1 Figure Supplement 5,6, indicating that vertically-evolving genes evolve relatively slowly (note that large values of ∆LL and split score denote low verticality)). Thus, the longer AB branches of vertically-evolving genes do not appear to result from a faster evolutionary rate for these genes. Taken together, these results indicate that the inclusion of genes that do not support the reciprocal monophyly of Archaea and Bacteria, or their constituent taxonomic ranks, in the universal concatenate explain the reduced estimated AB branch length.
Differential metric for verticality:
In spite of the similarity between the current result and Zhu et al.'s (see above), the two works approached this goal using different metrics.
First, the authors attempted to quantify the AB branch length in individual gene trees, including those that do not have Archaea and Bacteria perfectly separated. To do so they performed a constrained ML search (line 210). I am wary of this treatment because it could force distinct sequences (due to HGT or paralogy) to be grouped together, and the resulting branch length estimates could be highly inaccurate.
We agree with the reviewer that estimating AB branch lengths in this way might lead to inaccuracy. We note that this is, in effect, what was done in the published analysis (Zhu et al. 2019): a topology in which Archaea and Bacteria were reciprocally monophyletic was inferred using ASTRAL (a reasonable analysis, given the robustness of ASTRAL to
some degree of HGT/gene tree incongruence), and then the AB branch length was estimated from the concatenation of these 381 genes, fixing on the ASTRAL topology. We
performed this experiment (inferring AB branch length on constrained trees) in order to evaluate how incongruence between the gene and species trees might affect AB branch
length inference.
In contrast, Zhu et al. quantifies the average taxon-to-taxon phylogenetic distance between the two domains, regardless of the overall domain monophyly. That method was free of this concern, although it computed a different metric.
Thanks for raising this point. As described above, in the revised manuscript we have also evaluated the relative AB distance metric used by Zhu et al., and show that it behaves similarly to the AB branch metric we evaluated in the first version of the manuscript (see revised Figure 1).
Second, the authors assessed "marker gene verticality" using two metrics: a) AU test result (rejected or not) (Fig. 1A), c) ΔLL, the difference in log likelihood between the constrained ML tree and ML gene tree (line 222, Fig. 1B, C). I am concerned that they are sensitive to taxon sampling and stochastic events, as I explained above regarding domain monophyly. It is possible that a single mislabeling event would cause the topology test to report a significant result. In addition, they evaluate how severely domain monophyly is violated, but they do not account for intra-domain HGTs and other artifacts, which are also part of "verticality", and they can potentially distort the AB branch as well.
In the revised manuscript, we also evaluate a complementary metric for marker gene verticality, the split score (see above), which measures the extent to which marker genes recover established relationships at a given taxonomic level (we computed both within-domain and between-domain split scores). The split score is a more granular measure than ΔLL and, by summing over bootstrap replicates, it also better accommodates phylogenetic uncertainty. The two metrics (ΔLL and split score) are positively correlated and the analyses come to the same conclusions regarding the impact of HGT and other sources of incongruence on estimates of the AB branch length and relative AB distance.
I did not find the ΔLL values of individual markers in any supplementary table. I also did not find any correlation statistics associated with Fig. 1B.
The ΔLL values for individual markers can be found in the data supplement in:
“Expanded_Bacterial_Core_Nonribosomal_analyses/Individual_gene_tree_analyses/Expanded//Expanded_AB_AU.csv”
We have now updated the readme.txt file for clarity and included all the new results from the analyses which we have undertaken as part of the review process in the latest version of the supplemental available on figshare (10.6084/m9.figshare.13395470) as well as updated the directory and file names for clarity. We also have added the statistics associated with the correlations in Figure 1 to the Figure Legend.
Statistical test:
Line 157: "For the remaining 302 genes, domain monophyly was rejected (p < 0.05) for 232 out of 302 (76.8%) genes." Did the authors perform multiple hypothesis correction? If not, they probably should.
Thanks for this suggestion. We have now used a Bonferonni correction to account for multiple testing. As a result, fewer marker genes are rejected at the 5% level (151/302), although the overall conclusions are unaffected.
Line 217: "This result suggests that inter-domain gene transfers reduce the AB branch length when included in a concatenation." and Fig. 1A. If I understand correctly, this analysis was performed on individual gene trees, rather than in a concatenated setting. Therefore, the result does not directly support this conclusion.
Thanks for pointing this out. The reviewer is correct that this inference depends not only on the single gene analyses, but also on the subsequent concatenation results presented in this section. We have therefore moved this sentence later in the section, after the concatenation analysis.
Line 224: "Furthermore, AB branch length decreased as increasing numbers of low-verticality markers were added to the concatenate (Figure 1(c))". While this statement is likely true, Zhu et al. also presented similar results (Fig. 5) despite using a different metric, and they concluded that the impact is moderate and cannot explain the status of some core genes as outliers.
Zhu et al. did identify some of these trends, as we acknowledge in our manuscript ("The original study investigated and acknowledged (Zhu et al., 2019) the varying levels of congruence between the marker phylogenies and the species tree, but did not investigate the underlying causes.“ --- line 178; “These results are consistent with (Zhu et al., 2019), who also noted that AB branch length increases as model fit improves for the expanded marker dataset.” --- line 337) and as discussed above. Our analysis (Figure 6,
Table 1) goes further in showing that the most vertical subset of the 381-gene set supports an inter-domain branch length closely similar (2.4 subs/site compared to e.g. 2.5 subs./site
for the 27-gene dataset) to analyses using the traditional marker gene set.
Concatenation and branch length:
The authors pointed out that "Concatenation is based on the assumption that all of the genes in the supermatrix evolve on the same underlying tree; genes with different gene tree topologies violate this assumption and should not be concatenated because the topological differences among sites are not modelled, and so the impact on inferred branch lengths is difficult to predict." (line 187).
This argument is valid. In my opinion, this is the one most important potential issue of Zhu et al.'s analysis. In that work, they inferred genome tree topology through ASTRAL, which resolves conflicting gene evolutions. However ASTRAL does not report branch lengths in the unit of number of mutations. Therefore, they plugged the concatenated alignment into this topology for branch length estimation, hoping that it will "average out" the result. That workaround was apparently not ideal.
Yes, we agree --- this is our main critique of the Zhu et al. analyses.
However, the practice of molecular phylogenetics is complicated. Theoretically, every gene, domain, codon position and site may have its unique evolutionary process, and there have been efforts to develop better partition and mixture models to address these possibilities. But there is a trade off; these technologies are computationally demanding and have the risk of overfitting. It is plausible that in some scenarios, the gain of concatenating many loci (despite conflicting phylogeny) may outweigh the cost of having unpredictable effects.
But this dilemma needs to be analyzed rather than just being discussed. The Zhu et al. paper did not assess the impact of such concatenation on branch length estimation. The best answer is to conduct an analysis to show that concatenating genes with conflicting phylogeny would result in an AB branch that is shorter than the mean of those genes, and the reduction of AB branch length is correlated with the amount of conflict involved. The current work has not done this.
Thanks for raising this point. We agree that phylogenetics is complicated and that
we lack methods that can account for all possible factors. With respect to the impact of gene transfers on the AB branch length, and as touched on above, there are two issues here.
The first is with the analysis actually performed by Zhu et al: of the 381 extended set genes, 79 have one or no archaea in the 1000-taxon subsample, and a further 176 have an AB branch length close to 0 (<0.00001) in the constrained analyses. To investigate further, we manually inspected ML gene trees for the 381 genes (1000 taxon subsample). Allowing for recent gene transfer, we nevertheless identified only 64 genes with an unambiguous branch separating most Archaea from most Bacteria that might correspond to the ancestral AB divergence (Supplementary File 1).
Taken together, these analyses suggest that there is no strong evidence that these genes were present in LUCA or evolved along the AB branch, and so they do not provide information on its length. Since the branch length in the concatenation is an average over the branch length per site, the inclusion of this set of genes in the analysis did reduce the AB branch length, as demonstrated by our analyses (Figure 1(H)).
The second issue is: for genes which were likely present in LUCA and evolved on the AB branch, does gene transfer cause a reduction in the AB branch length inferred from their concatenation? To test this, we initially tested iterative concatenations of increasing numbers of non-vertical markers (Figure 1H), as well as a comparison of the most vertical genes to the whole expanded marker set (Figure 3 Figure Supplement 2). This revealed that as more markers were added (with lower verticality), the inferred AB branch length from the concatenate was reduced. We also found an increased AB branch length when only the 20 most vertical markers were used as opposed to the whole (381 marker) dataset (0.56 vs 0.16 substitutions/site, Figure 6).
The reviewer proposes an additional test of the impact of marker gene incongruence on branch length inference from concatenations: to compare the AB branch length before and after pruning of HGTs from individual marker gene alignments. To do this, we took the 54 marker genes from our new dataset and concatenated them before and after pruning of unambiguous HGTs. The AB branch length inferred from the concatenation with HGTs removed was 1.946 substitutions/site, compared to 1.734 substitutions/site without pruning HGTs, demonstrating the impact of even a relatively small number of HGTs on branch length estimation from concatenates.
Divergence time estimation:
The manuscript dedicates one section (line 230-266) to argue that the divergence time estimation analysis performed by Zhu et al. was not good evidence for marker gene suitability. Zhu et al. showed congruence of the expanded marker set with geological records whereas ribosomal proteins were conflicting with the geologic record.To support their argument, the authors estimated divergence times using the top 20 most "vertical" genes measured by ΔLL.
It would be good to clarify which genes they are, and it would be important to check whether they include some of the most "AB-distant" ones found by Zhu et al. Their Fig. 5A shows that there are genes that divide the two domains several folds further than the ribosomal proteins (such as rpoC). If they are among the 20 genes, it will not be surprising that the estimated AB split is older than it should be.
We now include the annotations for these 20 genes in Supplementary File 5a. The 20 most vertical genes include two of the “AB-distant” outliers identified by Zhu et al., tuf and infB, and one ribosomal marker, rpsG.
Overall, I think this section is logically questionable. Zhu et al. suggested that "They show the limitation of using core genes alone to model the evolution of the entire genome, and highlight the value in using a more diverse marker gene set.". The current work showed that using another set of a few genes (I do not know if they include multiple "core" genes, as discussed above, but it is plausible) also did not work well. This does not refute Zhu et al.'s claim.
What's important in Zhu et al.'s analysis is this: they demonstrated that using a small set of genes in DTE may cause artifacts due to them significantly violating the molecular clock at certain stages of evolution. Instead, using a larger set of markers that represent a portion of the entire genome would help to "smooth out" these artifacts. This of course is not the ideal solution, likely because concatenating conflicting genes and modelling them uniformly is not the best idea (see above). But as an operational workaround, it was not challenged by the analysis in the current work.
Finally, I agree with the authors' statement that more and reliable calibrations are the best way to improve divergence time estimation.
The dating analyses presented in the first version of our manuscript demonstrated that the apparent agreement between molecular clock estimates using the 381-gene set and the fossil record was the result of artifactual shortening of the AB branch, as discussed in detail above. Once the subset of the data least affected by these issues (that is, the most vertical subset) was used, the limitations of current clock methods, particularly with few calibrations, for dating deep nodes became clear.
That said, we agree with the reviewer (and also R3) that the dating section in the first version of our manuscript was somewhat unsatisfactory: it identified an important limitation of the published analysis, but did not explore the underlying question of why molecular clock methods infer unrealistically old divergence times from vertically-evolving genes. In the revised manuscript we have reworked and improved this section extensively, including new analyses on the 27-gene dataset, with more fossil calibrations, that help to diagnose how and why clocks struggle to date the archaeal and bacterial stems from the available data. We now show that the old ages result from a combination of low rates of molecular evolution across the tree inferred from “shallow” calibrations, combined with a lack of age maxima for nodes other than the root of the tree; when the rate distribution is informed in this way, the long AB branch is interpreted as representing a long period of time and estimates of LUCA age are strongly influenced by prior assumptions about root maximum age. These analyses now suggest how the difficulties might be overcome in the future, for example using better calibrations (particularly maximum ages, and indeed any fossil calibrations within the Archaea), or alternatively other sources of time information such as from gene transfers. Reflecting the new, broader focus, we have moved this section to the end of the manuscript.
AB branch by ribosomal and non-ribosomal genes:
Two figures (Figs. 3 and 4) are two sections (line 270-303) dedicated to the argument that ribosomal markers do not indicate a longer AB branch than a non-ribosomal one. However, this is a small scale test (38 ribosomal markers vs. 16 non-ribosomal markers) compared with the similar analysis in Zhu et al. (30 ribosomal markers vs. 381 global markers). A closer look at Figs. 3 and 4 suggests that while the AB lengths indicated by the ribosomal markers are within a relatively narrow range, those by the non-ribosomal ones are very diverse, including ones that are several folds longer than the ribosomal average. This result is in accordance with that of Zhu et al.'s Fig. 5A, although the latter was describing a different metric. Do these genes also overlap the ones found by Zhu et al.?
Nevertheless, this analysis does not falsify Zhu et al.'s, because it compared a different, much smaller, and deliberately chosen group of genes.
As the reviewer indicates, the purpose of the analyses presented in Figures 3-4 is to evaluate the hypothesis of accelerated ribosomal protein evolution: that is, the idea that ribosomal proteins over-estimate the AB branch length due to accelerated evolution during the divergence of Archaea and Bacteria. Although this hypothesis was independently proposed in Zhu et al., to our knowledge it actually originates with Petitjean et al. (2014) GBE (https://academic.oup.com/gbe/article/7/1/191/601621; see their Figure 2), and has been at play in analyses of deep evolution and in particular the position of DPANN Archaea in the phylogeny since that time. Thus, this section of our manuscript (indeed, all but the first section) is not a critique of Zhu et al.’s work, but a contribution to the broader ongoing discussion about which marker genes are best to use in deep phylogeny. We compare only vertically-evolving genes in Figures 3-4 so as to distinguish the impact of gene function (ribosomal versus non-ribosomal) from confounding factors such as HGT, paralogy, and gene origination time.
To clarify this point, we have modified our main text discussion to make it clear that we are making a comparison between ribosomal genes and other vertically-evolving members of the traditional “core” gene set, rather than a broader genome-wide claim. We now write:
“If ribosomal proteins experienced accelerated evolution during the divergence of Archaea and Bacteria, this might lead to the inference of an artifactually long AB branch length (Petitjean et al., 2014; Zhu et al., 2019). To investigate this, we plotted the inter-domain branch lengths for the 38 and 16 ribosomal and non-ribosomal genes, respectively, comprising the 54 marker genes set. We found no evidence that there was a longer AB branch associated with ribosomal markers than for other vertically-evolving “core” genes (Figure 2(b); mean AB branch length for ribosomal proteins 1.35 substitutions/site, mean for non-ribosomal 2.25 substitutions/site).”
Substitutional saturation:
The comparative analysis of slow- and fast-evolving sites is interesting. The result (Fig. 5) is visually impactful. In my view, this analysis is valid, and the conclusion is supported. It would be better to explain the rationale with more detail to facilitate understanding by a general audience.
Thanks for this assessment. We have now expanded on the rationale of this analysis in the main text, writing:
“It is interesting to note that the proportion of inferred substitutions that occur along the AB branch differs between the slow-evolving and fast-evolving sites. As would be expected, the total tree length measured in substitutions per site is shorter from the slow-evolving sites, but the relative AB branch length is longer (1.2 substitutions/site, or ~2% of all inferred substitutions, compared to 2.6 substitutions/site, or ~0.04% of all inferred substitutions for the fastest-evolving sites). Since we would not expect the distribution of substitutions over the tree to differ between slow-evolving and fast-evolving sites, this result suggests that some ancient changes along the AB branch at fast-evolving sites have been overwritten by more recent events in evolution --- that is, that substitutional saturation leads to an underestimate of the AB branch length.”
Zhu et al. also tested the impact of substitution saturation on the AB branch, using a more traditional approach (Fig. S19). They also found that the inter-domain distance is more influenced by potential substitution saturation, but the difference is minor. They concluded that (AB distance) "is not substantially impacted by saturation."
Like other analyses, these two analyses involved very different locus sampling (27 most "vertical" genes vs. 381 expanded genes). They also differ by the metric being measured (AB branch length vs. average distance between AB taxa). Therefore, the analysis in the current work does not falsify the analysis by Zhu et al. In contrast, it is inline with (though not in direct support of) Zhu et al. and others' suggestion that there was "accelerated evolution of ribosomal proteins along the inter-domain branch" (line 25) in the 27 core genes (of which 15 are ribosomal proteins).
We disagree that our analysis is consistent with the hypothesis of accelerated ribosomal protein evolution. The analysis that directly addresses this point is Figure 3, where we show that the distributions of AB branch lengths in single gene trees are not significantly different between ribosomal and non-ribosomal datasets (Figure 3; mean AB branch length for ribosomal proteins 1.35 substitutions/site, mean for non-ribosomal 2.25 substitutions/site).
Evolutionary model fit:
The authors compared the AB branch length indicated by the standard, site-homogeneous model LG+G4+F vs. the site-heterogeneous model LG+C60+G4+F, and found that the latter recovered a longer AB branch (2.52 vs. 1.45). The author's reasoning for using a site-heterogeneous model is valid, and this analysis is sound.
However, Zhu et al. also analyzed their data using the site-heterogeneous model C60 -- the same as in this work, but through the PMSF (posterior mean site frequency) method. Zhu et al. also compared it with two site-homogeneous models (Gamma and FreeRate). The results were extensively presented and discussed (Figs. 3, 4E, F, S23, S24, Note S2). They also found that C60+PMSF elongated the AB branch compared with the site-homogeneous models (Fig. S24A). As for the average AB distance (another metric evaluated by Zhu et al., as discussed above), C60+PMSF increased this metric when using ribosomal proteins, but not much when using the expanded marker set (Fig. S25A). And overall, the elongation by C60+PMSF with the expanded markers cannot compensate for the long branch indicated by the ribosomal proteins.
Therefore, similar to the point I made above, this analysis is sound but it does not logically falsify the conclusion made by Zhu et al., as it only concerns a small set of markers, and it recovered a previously described pattern.
Thanks for this comment. As above, note that the second part of our manuscript presents a general analysis of the issues around marker gene and model selection using our meta-analysis and new dataset, and is not a direct response to Zhu et al’s work. On reflection, we agree that this was not sufficiently clear in the first version of the paper, and we have now modified the text to acknowledge the model fitting analyses of Zhu et al.
The manuscript also did not clarify what the phrase "poor model fit" refers to (line 34 and line 304). If this is addressing the Gamma model evaluated by the authors, then this claim is valid though not novel (but see my previous comment on the trade-off). If that is a general reference to Zhu et al.'s methodology, then the authors should at least include the C60+PMSF model in the analysis, and show that C60 indicates a significantly longer AB branch than C60+PMSF does (if that's the case, which is doubtful). Admittedly, C60+PMSF is cheaper than the native C60 in computation, but "In some empirical and simulation settings PMSF provided more accurate estimates of phylogenies than the mixture models from which they derive." (Wang et al. Syst Biol. 2018).
Thanks for this comment. We did not intend the phrase “poor model fit” to imply a critique of Zhu et al.’s work; as the reviewer notes, those authors carried out a range of analyses to investigate the impact of model choice on their inferences. Rather, the title of the section is intended to summarise its main conclusion, which is that substitutional saturation and poor model fit (on any dataset, and even with the best available models) can lead to under-estimation of the AB branch length. Note that the analyses in Table 1 illustrating the impact of model fit are from the new dataset that is assembled and analysed in the second part of the manuscript. As above, we agree that this was not sufficiently clear in the first version of the paper. We think the title of this section is accurate and so we have not changed it, but we have changed the final two paragraphs of the section (as quoted immediately above) so as to acknowledge the model fitting analyses of Zhu et al., and to clarify that the results are general (and based on our new dataset), rather than a critique of Zhu et al’s work.
Finally, Zhu et al. also performed an analysis using the native C60 model on a further reduced taxon set. That result was not presented in the published paper, but it can be found in the "Peer Review File" posted on the Nature Communications website. That tree also recovered a short AB distance, and placed CPR at the base of Bacteria, and showed that this placement was not impacted by the removal of Archaea.
Thanks for pointing us to this additional analysis. The unrooted, bacteria-only tree referred to by the reviewer (panel B) recovers a clan (that is, a cluster of branches on the unrooted tree) comprising CPR+Chloroflexi, in agreement with the analysis on the new marker dataset we present here (Figure 6). The disagreement between that analysis and the new tree presented here relates to the position of the archaeal outgroup, which in the Peer Review File panel A connects to the bacterial tree between CPR and Chloroflexi. If, as recently suggested, the bacterial root lies between Gracilicutes and Terrabacteria (Coleman et al. 2021), then CPR and Chloroflexi represent monophyletic sister lineages. We note that the CPR+Chloroflexi relationship recovered here and in Peer Review File Panel (B) has also been obtained in several other recent analyses (Taib et al. 2020, Coleman et al. 2021, Martinez-Gutierrez and Aylward 2021), as cited in the main text.
Taxon sampling:
My final comment is about taxon sampling. Zhu et al. developed an algorithm for less biased taxon sampling, and they argued that extensive taxon sampling is important in resolving the early evolution of life. They presented evidence showing that reduced taxon sampling changed overall topology and basal relationships (Figs. S13, S14, S23, Note S2). The analyses were performed in combination with the assessment of site sampling, locus sampling, substitution model and other factors. The importance of less biased and/or extensive taxon sampling was also noted by previous works, especially in a phylogenomic framework (e.g., Hedtke et al. Syst Biol. 2006; Wu and Eisen. Genome Biol. 2008; Beiko. Biol Direct. 2011). The current work is based on a smaller set of taxa, and it has not addressed the impact of taxon sampling. As I suggested above, some results may be sensitive to taxon sampling.
We agree that taxon sampling is important for phylogenetics. While the analyses of Zhu et al. (2019) included a very large number of genomes, sampling of genomes (and indeed marker genes) was biased, both towards Bacteria compared to Archaea, and also within Bacteria. In our revised manuscript, we now compare the taxon sampling between Zhu et al.’s work and our new analyses (see Figure 1 Figure Supplements 13,14,15 and Figure 4 Figure Supplements 1,2). Balanced sampling is important for phylogenetic inference (Heath et al., 2008; Hillis, 1998) and, by this criteria, the taxon sampling in the analyses of Zhu et al. was not ideal. Our new analyses made use of fewer genomes (700), but these sample the known diversity of Archaea and Bacteria in a more representative way (Figure 4 Figure Supplement 1,2).
Reviewer #3:
Moody and coworkers principally address a recent paper presented by Zhu et al. (Nature Communications, 2019). In their paper, Zhu and coworkers claim that (i) ribosomal protein genes, commonly used in resolving deep phylogenies, have experienced an increased rate of evolution right after LUCA, and (ii) that an expanded set of markers show that the branch separating archaea from bacteria (AB-branch) is 10-fold shorter than previously thought. Moody et and coworkers first demonstrate flaws in the Zhu et al. analysis: first, the expanded gene set is biased towards bacteria, with 25% of the single-gene trees having very few archaeal counterparts. Second, that over 75% of the single-gene trees from Zhu et al are not monophyletic at domain level, suggesting a large influence of horizontal gene transfers (HGT), inter-domain exchanges, and inclusion of paralogous sequences in the original datasets. Third, they show that genes with fewer HGT display longer AB-branches. Fourth, they show that the argument by Zhu et al. that the longer AB-branch yields absurd LUCA datation is not relevant. Fifth, and maybe most important, they show that the shorter AB-branches recovered by Zhu et al in their expanded dataset result from inadequate substitution models, which lead to underestimating rates and thus branch lengths.
Going further, they select a set of 54 manually curated markers (showing mostly monophyletic archaea and bacteria), both from ribosomal proteins (36) and non-ribosomal proteins (18) and retrieve these in a balanced set of 350 archaea and 350 bacteria. With this set, they show that ribosomal protein markers do not display longer AB-branches than non-ribosomal ones. They also show that diversity among Archaea and Bacteria, as measured as the total tree length within each domain, is very similar, when sampling equal number of genomes in both domains.
Strengths:
The paper is well-written and well structured. In general, the methodology chosen here is adapted to the question at hand and very rigorously followed. The balanced dataset (with equal amounts of bacteria and archaea) of 54 carefully selected genes is also appropriate to explore diversity differences between the two domains of life.
Although all arguments presented in Zhu et al are carefully re-evaluated, the part where Moody et al show that substitutional saturation and poor model fit is artifactually producing short AB-branches is quite compelling and elegantly presented.
Weaknesses:
One potential weakness, more in terms of significance than in terms of scientific soundness is that the paper is mostly "reactive", responding to a single other paper. The authors might have used the data and methodology presented here to give the paper a broader scope. An example would be to provide the audience with a solid protocol or general guidelines on how to avoid artifacts in making deep phylogenies. I believe that the authors have demonstrated that they have the authority to do that.
Thanks for this suggestion. We considered including guidelines of this type in the first version of the manuscript, but we were --- and remain --- wary of attempting to promote one particular way of doing deep phylogeny over others. These are difficult and slippery questions, and different approaches and perspectives (including ones we might disagree with) are, in a broader sense, useful in refining ideas and helping the field to make progress as a whole. That said, a recurring issue appears to be the question of the fit between model and data, both in terms of substitution model fit (as with the impact of site-heterogeneous models on branch length inferences) and the broader issue of using models that, for example, account for gene duplication or transfer. There are several recent reviews (including one by some of us) which treat these topics in detail and provide detailed advice. We have now raised and discussed these issues in our conclusion. We have also updated Figure 6 to illustrate the approach we used in assembling the new 27-gene dataset, which may be of use to others, and goes some way towards the suggestion of providing guidelines for future analyses. We now write:
“Our analysis of a range of published marker gene datasets (Petitjean et al., 2014; Spang et al., 2015; Williams et al., 2020; Zhu et al., 2019) indicates that the choice of markers and the fit of the substitution model are both important for inference of deep phylogeny from concatenations, in agreement with an existing body of literature (reviewed in (Kapli et al., 2021, 2020; Williams et al., 2021). We established a set of 27 highly vertically evolving marker gene families and found no evidence that ribosomal genes overestimate stem length; since they appear to be transferred less frequently than other genes, our analysis affirms that ribosomal proteins are useful markers for deep phylogeny. In general, high-verticality markers, regardless of functional category, supported a longer AB branch length. Furthermore, our phylogeny was consistent with recent work on early prokaryotic evolution, resolving the major clades within Archaea and nesting the CPR within Terrabacteria. Notably, our analyses suggested that both the true Archaea-Bacteria branch length (Figure 6A), and the phylogenetic diversity of Archaea, may be underestimated by even the best current models, a finding that is consistent with a root for the tree of life between the two prokaryotic domains.”
In the figure 6 legend, we also expand on guidelines for future analyses, writing:
“(B) Workflow for iterative manual curation of marker gene families for concatenation analysis. After inference and inspection of initial orthologue trees, several rounds of manual inspection and removal of HGTs and distant paralogues were carried out. These sequences were removed from the initial set of orthologues before alignment and trimming. For a detailed discussion of some of these issues, and practical guidelines on phylogenomic analysis of multi-gene datasets, see (Kapli et al., 2020) for a useful review.”
The authors use the difference in log-likelihood between the constrained and unconstrained gene trees as a proxy for verticality and thus marker gene quality (Figure 1b). However, they don't demonstrate that that metric is actually appropriate. Could the monophyly (or split score) be also involved here? The authors might want to comment on that.
Thanks for this suggestion, which has substantially improved our analysis of Archaea-Bacteria distance and marker gene verticality (see the revised Figure 1 and associated text). We have now evaluated the relationship between AB branch length and split score (both within- and between-domain level relationships) for the expanded marker set and have updated our results and discussion accordingly. We found that deltaLL and split score (both within- and between-domains) are positively correlated with each other, and negatively correlated with AB length (that is, high-verticality markers have longer AB branch lengths). These analyses also revealed that within-domain and between-domain split scores are strongly positively correlated, implying that genes that recover domain monophyly also do better at resolving within-domain relationships.
The argument about the age of LUCA an ad absurdum one, showing that using better suited genes one gets impossible time estimates. However, the argument presented by Zhu et al is also a "just so" argument (if we get a time estimate that doesn't make sense then the phylogeny must be wrong), which doesn't give it much weight. The authors themselves note well that this part of the paper is more revealing of the limitations of the strict clock method, or of the relaxed clock with one single calibration point, than of the quality or appropriateness of the dataset.
We agree that the dating section in the first version of our manuscript was somewhat unsatisfactory. We have now expanded it to include new analyses on our 27-gene dataset, using more fossil calibrations, in order to diagnose why current clock methods struggle to estimate evolutionary rate near the root of the tree, and how this impacts on the age of LUCA and other deep nodes. These analyses add substantial value to this section, which has been moved to the end of the manuscript to reflect its expanded focus.
Another small weakness (or loose end) is that manual curation of the 95 genes dataset is not consistently reducing the percentage of non-monopyhletic genes (e.g. 62 to 69% from the 95 to the 54 genes dataset for non-ribosomal genes; 21 to 33% from the 95 to the 27 genes dataset for ribosomal genes). The author don't discuss how this impacts the manual curation they perform on the datasets; however, they state that "manual curation of marker genes is important". The authors might want to discuss that aspect further.
Thanks for raising this point. We were not sufficiently clear in describing the logic of our approach in the first version of the manuscript, and have now revised the text to clarify. In this analysis, we used a strict binary definition of monophyly --- that is, even a single inter-domain transfer leads to non-monophyly (note that this is in contrast to the re- analysis of the expanded set, where we considered whether each marker statistically rejected domain monophyly). For some genes scored as non-monophyletic in this way, manual removal of a small number of unambiguous recent transfers is sufficient restore domain monophyly; for others, HGT is extensive and it is difficult to know how to filter the sequences so as to obtain a reliable marker gene alignment; it was these latter cases that we set aside. We have now revised this section to make the logic of the approach clear, writing:
“Prior to manual curation, non-ribosomal markers had a greater number of HGTs and cases of mixed paralogy. In particular, for the original set of 95 unique COG families (see ‘Phylogenetic analyses’ in Methods), we rejected 41 families based on the inferred ML trees, either due to a large degree of HGT, paralogous gene families or LBA. For the remaining 54 markers, the ML trees contained evidence of occasional recent HGT events. Strict monophyly was violated in 69% of the non-ribosomal and 29% of the ribosomal families. We manually removed the individual sequences which violated domain monophyly before re-alignment, trimming, and subsequent tree inference (see Methods). These results imply that manual curation of marker genes is important for deep phylogenetic analyses, particularly when using non-ribosomal markers. Comparison of within-domain split scores for these 54 markers indicated that markers that better resolved established relationships within each domain also supported a longer AB branch length (Figure 2A).”
In summary and despite the small weaknesses listed above, my opinion is that the authors reach their goal of showning that the AB-branch is indeed a long one, and that the results support the conclusion.
Impact:
The main point addressed by the authors here, the time of divergence between Archaea and Bacteria, is crucial to our understanding of early evolution. The long branch separating Bacteria and Archaea has long been thought to be a long one, and the paper by Zhu et al casted a doubt about the validity of this long-standing hypothesis. Here, Moody et al convincingly establish that the divergence between archaea and bacteria is a profound one. The paper also has profound implications on the validity of the commonly used core-gene phylogenies, particularly those based on ribosomal protein genes. Indeed, it shows that the these proteins are appropriate for deep phylogenies. They also show the impact of model violations on deep phylogenies, and how to avoid them.
We thank the reviewer for this positive assessment of impact.