The deep(er) roots of Eukaryotes and Akaryotes

It is challenging to review a reply without considering the original paper carefully. In this case it is doubly challenging because the most relevant portion of Williams et al. (2020)¹ is itself a reply to Harish and Kurland (2017)². This prompted me reread Harish and Kurland (2017) and Williams et al. (2020) in detail and evaluate the Harish and Morrison manuscript (which I will call HM-F1000 hereafter) in light of those publications. Thus, I have decided to provide a review of HM-F1000 along with a limited post-publication review of Harish and Kurland (2017) and Williams et al. (2020). I will divide this review into two sections: 1) a discussion of the larger philosophical questions and 2) a description of the changes to HM-F1000 that I believe to be necessary as well as a few minor issues with HM-F1000.
 
Before I provide that combined review, I want to answer two questions: 1) does HM-F1000 warrant publication as a peer-reviewed publication? and 2) did HM-F1000 convince me that the Eukaryote/Akaryote* two-domain tree of life (2D-ToL) represents an accurate placement of the root of the tree of life (ToL)? My answer to the first question is “yes” and my answer to the second question is “not at this time.”
 
* NOTE: Throughout this review I will use “akaryote” in because it is the terminology used in HM-F1000; however, I am agnostic regarding benefits of that term relative to “prokaryote.”
 
I have answered the first question in the affirmative despite expressing substantially greater caution when I answer the second because identifying the position for the root of the ToL is arguably one of the most difficult problems in evolutionary biology. The field needs more ideas regarding the best way to estimate a robust topology for the ToL and place the root, not fewer. Excluding the HM-F1000 from the peer-reviewed literature would exclude their concise defense of the idea that the KVR (Klopfstein et al. 2015)³ model of evolution can be used with protein fold presence/absence data to place the root of the ToL.
 
*****
Section 1:
 
Williams et al. (2020) argued that using the KVR model with protein domain composition data was inappropriate based upon their simulations; they chose instead to focus on analyses using models of protein sequence evolution. It is certainly true that the KVR model is a simplistic model of evolution. However, the fact that the KVR model is imperfect does not invalidate the use of the model; as Box (1979) famously stated “all models are wrong but some are useful.”
 
Simply showing that that a model has imperfect fit to empirical data (as Williams et al. 2020 did) does not mean the model is useless for inference. Indeed, it is very likely that the models of protein evolution, like the CAT, C60, and NDCH2 models, which Williams et al. (2020) used in their other analyses, also have an imperfect fit to the true underlying process of evolution. The central issue is whether analyses of protein fold data using the KVR model are more likely to recover true historical signal than analyses of aligned proteins using various models of protein sequence evolution.
 
HM-F1000 highlights a corollary of this fundamental issue in its abstract when they state that “(i)t is well known that different character types present different perspectives on evolutionary history that relate to different phylogenetic depths.” With that said, is clear that the models of protein sequence evolution that Williams et al. (2020) used are much more sophisticated than the KVR model. So why should we embrace the results of the KVR model over the results of analyses using those sophisticated models of protein evolution? Obviously, the purpose of HM-F1000 is to convince readers to accept the results of the KVR model (applied to protein fold presence/absence data) as more likely to be correct (in this context I will use “correct” to mean “closer to the truth”) than the results of the models of protein evolution used by Williams et al. (2020). Why might that be the case? I can think of two reasons that I will discuss below:
 
---
I. Historical signal might decay more rapidly in aligned protein sequences than in protein fold content data.
 
The simplest explanation for preferring the results of analyses using the protein fold data to those obtained using aligned proteins is the possibility that historical signal might have decayed in the latter. Mossel (2003)⁴ proved that “…it is impossible to reconstruct the topology of ‘deep’ trees with high mutation rates…” More accurately, Mossel (2003) identified a bound on the number of characters necessary for tree reconstruction, but that bound implies that impossible to reconstruct some past events (also see Sober and Steel 2002)⁵. Perhaps protein sequence alignments cannot provide accurate information about the deepest branches in the tree of life and we have to look to other data types, like protein fold content, to estimate the topology of the deepest branches in the tree of life.
 
In my opinion, two arguments are necessary to establish that data type is more important that model fit for reconstructing the deepest branches in the tree of life. HM-F1000 states that “(p)rotein structural domains, unlike amino acids, are biochemically non-redundant (see below) and have proven to be excellent ‘genomic characters’…” as a defense of the idea that protein fold data might be superior to aligned amino acids. I was expecting an explicit statement that it might be appropriate to view changes in the protein fold repertoire as rare genomic changes (RGCs; Rokas and Holland 2000; Bleidorn 2017)^6,7. I think that linking protein folds to RGCs is important because it provides an explicit link between protein fold data and the body of theory surrounding RGCs. Specifically, the fact that analyses using the maximum parsimony criterion are expected yield the correct tree when applied to RGC data (Steel and Penny 2004; 2005)^8,9. I believe this has implications for the idea that the relatively simple KVR model might be useful for rooting the tree of life.
 
Whether the maximum parsimony criterion should be viewed as a simple model (or any sort of model) has been a topic of philosophical debate in phylogenetics (Goloboff 2003’ Huelsenbeck et al. 2008)^10,11; I will accept the idea that maximum parsimony is “simple” for the sake of this argument (also see Yang 1995)¹². However, if we accept that a “perfect RGC” model (which I define as a process that results in some binary character that can only undergo a single transition on one edge in the gene tree associated with that genomic character) it allows us to pose a question about the KVR model: is the KVR model consistent for characters generated by a hypothetical “asymmetric perfect RGC” model? The asymmetric perfect RGC model modifies the perfect model so the ensemble state frequencies at the root differ from the tip frequencies. I recognize that, in addition to the treatment of the root state frequencies, the KVR model differs from parsimony in an important way (specifically, the treatment of branch lengths). However, this conjecture regarding the behavior of the KVR model might point the way toward a falsifiable hypothesis because it lends itself to testing by simulation.
 
The question of whether the KVR model is consistent given the asymmetric perfect RGC model is interesting from theoretical standpoint but there is a second (and more important question) that should be answered: is whether the true underlying model of fold content is sufficiently close to the asymmetric perfect RGC model for that model to be useful? The true underlying model of fold evolution includes fold origination (which is almost certainly a very rare event) and horizontal transfer (likely to be much more common). Williams et al. (2020) discuss this in their supplementary materials, where they state “…a change from 0 to 1 might indicate de novo origin of an existing fold by convergent evolution (which is likely to be rare), or the gain of an existing fold by [horizontal gene transfer]; if the latter, then the pattern of presences and absences for that fold cannot be reliably used to infer the underlying tree.” I agree with the first part of that sentence (which is a reason why I have invoked the idea of RGCs) but I disagree with the second; even when there is horizontal gene transfer novel fold acquisition might be sufficiently rare for that type of event to be considered an RGC.
 
Answering those questions will be challenging and outside the scope of a short note like HM-F1000. In that context, I think it would be good for HM-F1000 to express a little more caution. Statements like “[t]he KVR model is an optimal explanation of the evolution of clade-specific composition of homologous features” (first full paragraph on page 5 of HM-F1000). The point of my arguments above is that it might be reasonable to view the KVR model as an excellent approximating model for protein fold evolution. The first author has written multiple papers dealing with patterns of protein fold evolution over deep evolutionary time and I do not want to disrespect those efforts, but I do not think this is a settled issue at this time.
 
II. The best interpretation of the unrooted trees in Williams et al. (2020) is unclear.
 
The Harish and Kurland (2017) trees are the only intrinsically rooted trees under discussion in Williams et al. (2020). Since the position of the root of the 2D-ToL in Williams et al. (2020) ultimately represents the imposition of a root on an otherwise unrooted tree. Strictly from a logical standpoint there are four interpretations of the results of Harish and Kurland (2017) and Williams et al. (2020).
 

1. Both trees are inaccurate estimates of the ToL.
    
2. The Harish and Kurland (2017) topology is correct.
    
3. The unrooted Williams et al. (2020) topology is correct and the Harish and Kurland (2017) root is also correct (i.e., the root lies between eukaryotes and akaryotes).
    
4. The unrooted Williams et al. (2020) topology is correct and the root is not between eukaryotes and akaryotes.

 
Possibility 3 is important, and it would seem to be implicit in Figure 1b of HM-1000. The authors should make this point more explicitly; they use Figure 1b to make another point, using it to stress that “Williams et al. have overlooked important aspects of assessing phylogenetic signal in empirical data, and that it may be premature to reject a well-supported phylogeny based on simulated data.” I think it would be valuable to make the point that, at least in principle, the Williams et al. (2020) results would be consistent with a tree rooted between eukaryotes and akaryotes that also places the root of akaryotes within the Asgard archaea.
 
Obviously, embracing the unrooted Williams et al. (2020) topology would require rejecting the akaryote topology of the Harish and Kurland (2017). Specifically, the fact that the Harish and Kurland (2017) tree is consistent with archaeal monophyly (as long as one assumes the root of the tree of life is not within archaea) makes it fundamentally inconsistent with the Williams et al. (2020) topology. However, the possibility of archaeal non-monophyly would seem to be consistent with Figures 1 and 2 in Harish (2018)¹³, both of which show substantial uncertainty at the base of Archaea. Figure 2 in Harish (2018) is based on distances calculated using protein fold data, raising some questions regarding the strength of support for monophyly of Archaea.
 
One aspect of the Harish and Kurland (2017) tree that HM-F1000 should acknowledge is the fact plants are not monophyletic. Specifically, the root of the eukaryotic sub-tree of Figure 3 in that paper was placed between rice and all other eukaryotes. Obviously, this is troubling given that the Harish and Kurland (2017) tree includes other angiosperms. In fact, the Harish and Kurland (2017) dataset includes two other grasses; non-monophyly of both angiosperms and grasses is unreasonable. Even placing the eukaryotic root between the green plants and other eukaryotes seems unlikely given the best available information about the eukaryotic tree (reviewed by Burki et al. 2020)¹⁴.
 
One might wonder whether the root position for the ToL should be viewed as accurate given the unexpected position of the eukaryotic root. However, it is reasonable to postulate that the rice data were problematic in some way. Alternatively, it could reflect the observation that different data perform differently at different levels in the tree (Chen et al. 2015) (HM-F1000 already alluded to this). If I had reviewed Harish and Kurland (2017) I would have asked the authors to conduct a second set of analyses after excluding rice to see if that changed the root of the eukaryotic sub-tree. I do not think it would be reasonable to ask HM-F1000 to add a reanalysis of the Harish and Kurland (2017) after excluding rice, but it would be nice for HM-F1000 to acknowledge this issue.
---
 
Looking back, I realize that I have written this review as an advocate for the position articulated by HM-F1000. Given that tone it would be fair to ask why I not convince that their placement of the root between eukaryotes is accurate? I would answer that question I am not convinced that the tools exist to place of the root of the ToL is accurate exist at this point. As I discussed above, I have made an argument can be made that the KVR model might be able to yield an accurate estimate of the ToL. However, my reasons for that assertion rest on assumptions regarding the nature of the protein fold data. I believe that the Williams et al. (2020) analyses do show that the KVR model is imperfect; I simply disagree with their conclusion that this imperfect fit means that we should dismiss the Harish and Kurland (2017) result. Fundamentally, I view the use of the KVR model with protein fold data as similar to the Jukes-Cantor model. Felsenstein (2001)¹⁵ shared an anecdote regarding that model, stating that: “Tom Jukes once told me that the reason the Jukes–Cantor model was buried in the midst of a large empirical paper was that this was the only way to get it published. He felt that if he had attempted to publish it on its own, it would have been rejected by editors as idle and oversimplified speculation.” However, without the pioneering work of Jukes and Cantor (1969)¹⁶ and Neyman (1971)¹⁷ (or Felsenstein 1981)¹⁸ it is difficult to envision the development of the more sophisticated models of sequence evolution that developed over the subsequent five decades. Dismissing the use of protein fold data at this point will slow the development of those models. Will further model development support the Eukaryote/Akaryote 2D-ToL? I am uncertain whether it will, but I am interested to find out.
 
I would like to add a final discussion regarding model fit. Although there is a long history of model development for protein sequences and the models are now quite sophisticated, there is still much that we do not know about protein evolution. Williams et al. (2020) used the LG+C60 model. Presumably this is the C60 profile mixture from Le et al. (2008) (although I was unable to find Le et al. 2008 cited in Williams et al. 2020) combined with the LG (Le and Gascuel 2008)¹⁹ matrix. However, Pandey and Braun (2020)²⁰ reported that the Le et al. (2008)²¹ profile mixtures can exhibit surprising (and disturbing) behavior for at least one phylogenetic problem. The Le et al. (2008) mixture models are very similar to the CAT model (Lartillot and Philippe 2004); it is unclear whether these issues Pandey and Braun (2020) noted for the Le et al. (2008) models are general features of CAT-type models (or even more widespread than the particular case studied by Pandey and Braun 2020) but it is important to be careful regarding the use of any models of evolution so deep in the tree.
 
It is tempting to look at the sophistication of existing models of protein sequence evolution and conclude that the results obtained using those models trump other sources of information. Although I do not want to be overly dismissive of the Williams et al. (2020) analyses, which are state of the art, I do want to emphasize that I believe data type matters and that we should be looking at other sources of information. In my opinion, that is the message that HM-F1000 should convey; that is why I think some statements in HM-F1000, like the statement that the KVR model provides an optimal explanation for protein fold evolution, actually undercut the case. In my opinion, obtaining a strongly corroborated estimate of the deepest branches in the ToL, if it is possible, will require us to examine multiple sources of information and to be very careful regarding the models we use for analyses. That applies to analyses of aligned protein sequences and to analyses of protein fold content. 
 
 
*****
Section 2: Minor issues and description of necessary revisions:
 

1. I have written a fairly long review, but I feel the changes to HM-F1000 that are necessary are actually fairly minimal. I think HM-F1000 needs to walk back the claims that the KVR model is an optimal explanation for protein fold distribution and simply point out that it is likely to be a reasonable approximating model. I think HM-F1000 also needs to acknowledge that “both trees could be telling us part of the truth” (i.e., that a tree rooted between eukaryotes and akaryotes with a paraphyletic archaea might be a way to reconcile Harish and Kurland 2017 with Williams et al. 2020). HM-F1000 should also acknowledge the unexpected (and incorrect) rooting of the eukaryotic sub-tree. Finally, I hope the minor comments that follow are carefully considered.
    
2. I was surprised that the work of Poole et al. (1998; 1999)^22,23 was not cited. It provides another line of evidence supporting the placement of the root on the eukaryotic branch (i.e., it supports the Eukaryote/Akaryote 2D-ToL).
    
3. The first full line of the second column of the fourth page of HM-F1000 states “Support from fossils or other sources are not reliable, despite claims to the contrary [Williams et al. 2020].” I could not find an explicit statement in Williams et al. (2020) that makes this assertion. I agree with the basic point that the fossil record to establish the deep topology for the ToL provides, at most, limited information. However, HM-F1000 should be careful regarding the attribution of statements like this. I hope that I did not miss any such statement in Williams et al. (2020); if I have missed it, Harish and Morrison should point to the statement.
    
4. The legend of Figure 2 states “[t]he majority of proteins are multi-domain proteins formed by duplication and recombination of domain units.” However, Ekman et al. (2005)²⁴ reports that fewer than 50% of prokaryotic (akaryotic) proteins are multidomain. I was unable to find an explicit survey of proteins showing that the numbers of multidomain proteins is generally >50% in the literature. This statement should have an associated citation and, if the number is <50% in some lineage be a bit more cautious. Perhaps something like “a large proportion” would be a better statement.
    
5. The legend of Figure 2 also states that “[a]lthough it is common to suspect that the rooting between akaryotes and eukaryotes could be biased due to a larger domain cohort in eukaryotes [Williams et al. 2020], it is not the case.” Since the statement that the large domain cohort of eukaryotes is a source of bias only cites Williams et al. (2020) I don’t think it is valid to state that “it is common to suspect” unless there are additional citations. The statement “…it is not the case” cites three papers with Harish as an author, but the evidence that a large domain cohort cannot be a source of bias was not clear to me. The explanation should be expanded a bit and moved to the main text.