OCTAL: Optimal Completion of gene trees in polynomial time

Species tree estimation from multi-gene datasets is now increasingly common. One challenge is that the evolutionary history for a single locus (called a “gene tree”) may differ from the species phylogeny due to a variety of different biological processes. Some of these processes, such as hybridization [1] and horizontal gene transfer [2], result in non-treelike evolution and so require phylogenetic networks for proper analysis [3–6]. However, other biological processes, such as gene duplication and loss, incomplete lineage sorting (ILS), and gene flow, produce heterogeneity across the genome but are still properly modeled by a single species tree [7, 8]. In the latter case, species tree estimation methods should be robust to heterogeneity across the genome.

Much of the recent focus in the mathematical and statistical phylogenetics literature has been on developing methods for species tree estimation in the presence of incomplete lineage sorting (ILS), which is modelled by the multi-species coalescent (MSC) model [9]. One popular approach for estimating species trees under the MSC model is to estimate trees on individual loci and then combine these gene trees into a species tree. Some of these “summary methods”, such as ASTRAL-II [10] and ASTRID [11], have been shown to scale well to datasets with many taxa (i.e., >1000 species) and provide accurate species tree estimates. (Summary methods share many features in common with supertree methods, but are based on mathematical properties of the MSC model and so can be proven statistically consistent under the MSC model; supertree methods, by contrast, assume conflict between source trees is due to estimation error rather than ILS, and so are generally not statistically consistent under the MSC model.)

A common challenge to species tree estimation methods is that sequence data may not be available for all genes and species of interest, creating conditions with missing data (see discussion in [12–14]). For example, gene trees can be missing species simply because some species do not contain a copy of a particular gene, and in some cases, no common gene will be shared by every species in the set of taxa [15]. Additionally, not all genomes may be fully sequenced and assembled, as this can be operationally difficult and expensive [13, 16].

Although summary methods are statistically consistent under the MSC model [17], the proofs of statistical consistency assume that all gene trees are complete, and so may not apply when the gene trees are missing taxa. Recent extensions to this theory have shown that some species tree estimation methods are statistically consistent under some models of missing data (e.g., when “every species is missing from each gene with the same probability \(p > 0\)”) [18]. However, missing data in biological datasets often violates such models (see discussion in [14]); for example, missing data may be biased towards genes with faster rates of evolution [19]. Furthermore, multi-gene datasets with missing data can be “phylogenetically indecisive”, meaning more than one tree topology can be optimal [20]. Because of concerns that missing data may reduce the accuracy of multi-locus species tree estimation methods, many phylogenomic studies have restricted their analyses to only include genes with most of the species (see discussion in [12, 13, 21]).

We approach the challenge of adding missing species into gene trees by formulating the Optimal Tree Completion problem, where we seek to add the missing species to a gene tree to minimize the distance (defined in some way) to another tree, called a “reference tree”. Since the Robinson–Foulds [22] distance is a common metric for comparing trees (where the Robinson–Foulds distance is the total number of unique bipartitions in the two trees), we specifically address the Robinson–Foulds (RF) Optimal Completion problem, which seeks a completion of the input gene tree with respect to a given reference tree that minimizes the RF distance between the two trees. We then present the Optimal Completion of Incomplete gene Tree Algorithm (OCTAL), a greedy polynomial time algorithm that we prove solves the RF Optimal Completion problem exactly. We also present results from an experimental study on simulated datasets comparing OCTAL to a heuristic for gene tree completion within ASTRAL-II. Finally, we conclude with a discussion of results and future research.

The algorithm begins with input tree t and adds leaves one at a time from the set \(S \setminus R\) until a tree on the full set of taxa S is obtained. To add the first leaf, we choose an arbitrary taxon x to add from the set \(S \setminus R\). We root the tree \(T|_{R \cup \{x\}}\) (i.e., T restricted to the leaf set of t plus the new leaf being added) at x, and then remove x and the incident edge; this produces a rooted binary tree we will refer to as \(T^{(x)}\) that has leaf set R.

We perform a depth-first traversal down \(T^{(x)}\) until a shared edge e (i.e., an edge where the clade below it appears in tree t) is found. Since every edge incident with a leaf in \(T^{(x)}\) is a shared edge, every path from the root of \(T^{(x)}\) to a leaf has a distinct first edge e that is a shared edge. Hence, the other edges on the path from the root to e are unique edges.

After we identify the shared edge e in \(T^{(x)}\), we identify the edge \(e'\) in t defining the same bipartition, and we add a new node \(v(e')\) into t so that we subdivide \(e'\). We then make x adjacent to \(v(e')\). Note that since t is binary, the modification \(t'\) of t that is produced by adding x is also binary and that \(t'|_R = t\). These steps are then repeated until all leaves from \(S \setminus R\) are added to t. This process is shown in Fig. 1 and given in pseudocode below.

Proof of correctness

In what follows, let T be an arbitrary binary tree on taxon set S and t be an arbitrary binary tree on taxon set R \(\subseteq\) S. Let \(T'\) denote the tree returned by OCTAL given T and t. We set \(r=RF(T|_R,t)\). As we have noted, OCTAL returns a binary tree \(T'\) that is an S-completion of t. Hence, to prove that OCTAL solves the RF Optimal Tree Completion problem exactly, we only need to establish that \(RF(T,T')\) is the smallest possible of all binary trees on leaf set S that are S-completions of t. While the algorithm works by adding a single leaf at a time, we use two types of subtrees, denoted as superleaves (see Fig. 2), to aid in the proof of correctness.

Definition 1

The backbone of T with respect to t is the set of edges in T that are on a path between two leaves in R.

Definition 2

A superleaf of T with respect t is a rooted group of leaves from \(S \setminus R\) that is attached to an edge in the backbone of T. In particular, each superleaf is rooted at the node that is incident to one of the edges in the backbone

Definition 3

There are exactly two types of superleaves, Type I and Type II:

1. 1

   A superleaf is a Type I superleaf if the edge e in the backbone to which the superleaf is attached is a shared edge in \(T|_R\) and t. It follows then that a superleaf X is a Type I superleaf if and only if there exists a bipartition A|B in \(C(t) \cap C(T|_R)\) where \(A|(B \cup X)\) and \((A \cup X)|B\) are both in \(C(T|_{R \cup X})\).

    
2. 2

   A superleaf is a Type II superleaf if the edge e in the backbone to which the superleaf is attached is a unique edge in \(T|_R\) and t. It follows that a superleaf X is a Type II superleaf if and only if for any bipartition A|B such that \(A|(B \cup X)\) and \((A \cup X)|B\) are both in \(C(T|_{R \cup X})\), \(A|B \not \in C(t)\).

    

Now we begin our proof by establishing a lower bound on the RF distance to T for all binary S-completions of t.

Lemma 4

Let Y be a Type II superleaf for the pair (T, t), and let \(x \in S \setminus R\). Let \(t^*\) be the result of adding x into t arbitrarily (i.e., we do not attempt to minimize the resulting RF distance). If \(x \not \in Y\), then Y is a Type II superleaf for the pair \((T,t^*)\). Furthermore, if \(x \in Y\), then \(RF(T|_{R \cup \{x\}}, t^*) \ge RF(T|_R,t) +2\).

Proof

It is easy to see that if \(x \not \in Y\), then Y remains a Type II superleaf after x is added to t. Now suppose \(x \in Y\). We will show that we cannot add x into t without increasing the RF distance by at least 2. Since Y is a Type II superleaf, it is attached to a unique edge in \(T|_{R \cup Y}\), and this is the same edge that x is attached to in \(T|_{R \cup \{x\}}\). So suppose that x is added to t by subdividing an arbitrary edge \(e'\) in t with bipartition C|D; note that we do not require that x is added to a shared edge in t. After adding x to t we obtain tree \(t^*\) whose bipartition set includes \(C|(D \cup \{x\})\) and \((C \cup \{x\})|D\). If C|D corresponds to a unique edge relative to t and \(T|_R\), then both of these bipartitions correspond to unique edges relative to \(t^*\) and \(T|_{R \cup \{x\}}\). If C|D corresponds to a shared edge, then at most one of the two new bipartitions can correspond to a shared edge, as otherwise we can derive that Y is a Type I superleaf. Hence, the number of unique edges in t must increase by at least one no matter how we add x to t, where x belongs to a Type II superleaf. Since t is binary, the tree that is created by adding x is binary, so that \(RF(T|_{R \cup \{x\}},t^*) \ge RF(T|_R,t) +2\). \(\square\)

Lemma 5

Let \(T^*\) be an unrooted binary tree that is a S-completion of t. Then \(RF(T^*,T) \ge r+2m\), where \(r=RF(T|_R,t)\) and m is the number of Type II superleaves for the pair (T, t).

Proof

We note that adding a leaf can never reduce the total RF distance. The proof follows from Lemma 4 by induction. \(\square\)

Now that we have established a lower bound on the best achievable RF distance (i.e., the optimality criterion for the RF Optimal Tree Completion problem), we show OCTAL outputs a tree \(T'\) that is guaranteed to achieve this lower bound. We begin by noting that when we add x to t by subdividing some edge \(e'\), creating a new tree \(t'\), all the edges other than \(e'\) in t continue to “exist” in \(t'\) although they define new bipartitions. In addition, \(e'\) is split into two edges, which can be considered new. Thus, we can consider whether edges that are shared between t and T remain shared after x is added to t.

Lemma 6

Let \(t'\) be the tree created by AddLeaf given input tree t on leaf set R and tree T on leaf set \(R \cup \{x\}\). If x is added to tree t by subdividing edge \(e'\) (thus creating tree \(t'\)), then all edges in t other than \(e'\) that are shared between t and T remain shared between \(t'\) and T.

Proof

Let \(T^{(x)}\) be the rooted tree obtained by rooting T at x and then deleting x. Let e be the edge in \(T^{(x)}\) corresponding to \(e'\), and let \(\pi _e=A|B\); without loss of generality assume A is a clade in \(T^{(x)}\). Note that C(T) contains bipartition \(A|(B \cup \{x\})\) (however, C(T) may not contain \((A \cup \{x\})|B\), unless e is incident with the root of \(T^{(x)}\)). Furthermore, for subclade \(A' \subseteq A\), \(A'|(R \setminus A') \in\) \(C(T|_R)\) and \(A'|(R \setminus A' \cup \{x\}) \in\) C(T). Now suppose \(e^*\) in t is a shared edge between t and \(T|_R\) that defines bipartition \(C|D \ne A|B\). Since A|B and C|D are both bipartitions of t, without loss of generality either \(C \subset A\) or \(A \subset C\). If \(C \subset A\), then C is a clade in \(T^{(x)}\), and so \(e^*\) defines bipartition \(C|(D \cup \{x\})\) within \(t'\). But since \(C \subset A\), the previous analysis shows that \(C|(D \cup \{x\})\) is also a bipartition of T, and so \(e^*\) is shared between T and \(t'\). Alternatively, suppose \(A \subset C\). Then within \(t'\), \(e^*\) defines bipartition \((C \cup \{x\})|D\), which also appears as a bipartition in T. Hence, \(e^*\) is also shared between T and \(t'\). Therefore, any edge \(e^*\) other than \(e'\) that is shared between t and T remains shared between \(t'\) and T, for all leaves x added by AddLeaf. \(\square\)

Lemma 7

OCTAL(T, t) preserves the topology of superleaves in T (i.e. for any superleaf with some subset of leaves \(Q \subseteq S\), OCTAL(T, t)\(|_Q\) equals \(T|_Q\)).

Proof

We will show this by induction on the number of leaves added. The lemma is trivially true for the base case when just one leaf is added to t. Let the inductive hypothesis be that the lemma holds for adding up to n leaves to t for some arbitrary \(n \in \mathbb {N}^+\). Now consider adding \(n+1\) leaves, and choose an arbitrary subset of n leaves to add to t, creating an intermediate tree \(t'\) on leaf set K using the algorithm OCTAL. Let x be the next additional leaf to be added by OCTAL.

If x is the first element of a new superleaf to be added, it is trivially true that the topology of its superleaf is preserved, but we need to show that x will not break the monophyly of an existing superleaf in \(t'\). By the inductive hypothesis, the topology of each superleaf already placed in \(t'\) has been preserved. Thus, each superleaf placed in \(t'\) has some shared edge in \(t'\) and \(T|_{K}\) incident to that superleaf. If x were placed onto an edge contained in some existing superleaf, that edge would change its status from being shared to being unique, which contradicts Lemma 6.

The last case is where x is part of a superleaf for the pair (T, t) that already has been added in part to t. AddLeaf roots \(T|_{K \cup \{x\}}\) at x and removes the edge incident to x, creating rooted tree \(T^{(x)}\). The edge incident to the root in \(T^{(x)}\) must be a shared edge by the inductive hypothesis. Thus, OCTAL will add x to this shared edge and preserve the topology of the superleaf. \(\square\)

Lemma 8

OCTAL(T, t) returns binary tree \(T'\) such that \(RF(T,T')=r+2m\), where m is the number of Type II superleaves for the pair (T, t) and \(r=RF(T|_R,t)\).

Proof

We will show this by induction on the number of leaves added.

Base Case Assume \(|S\setminus R|\) = 1. Let x be the leaf in S\(\setminus R\). AddLeaf adds x to a shared edge of t corresponding to some bipartition A|B, which also exists in \(T^{(x)}\).

1. 1.
   First we consider what happens to the RF distance on the edge x is attached to.

   • If x is a Type I superleaf, the edge incident to the root in \(T^{(x)}\) will be a shared edge by the definition of Type I superleaf, so AddLeaf adds x to the corresponding edge \(e'\) in t. The two new bipartitions that are created when subdividing \(e'\) will both exist in T by the definition of Type I superleaf so the RF distance does not change.

   • If x is a Type II superleaf, either \((A \cup \{x\})|\)B or \(A|(B\cup \{x\})\) must not exist in C(T). Since AddLeaf adds x to a shared edge, exactly one of those new bipartitions must exist in C(T).

    
2. 2.

   Now we consider what happens to the RF distance on the edges x is not attached to. Lemma 6 shows that AddLeaf (and therefore OCTAL) preserves existing shared edges between t and \(T|_R\), possibly excluding the edge where x is added.

    

Thus, the RF distance will only increase by 2 if x is a Type II superleaf, as claimed.

Inductive step Let the inductive hypothesis be that the lemma holds for up to n leaves for some arbitrary \(n \in \mathbb {N}^+\). Assume \(|S\setminus R|\) = \(n+1\). Now choose an arbitrary subset of leaves \(Q\subseteq S \setminus R\), where \(|Q|=n\), to add to t, creating an intermediate tree \(t'\) using the algorithm OCTAL. By the inductive hypothesis, assume \(t'\) is a binary tree with the RF distance between \(T|_{Q\cup R}\) and \(t'\) equal to \(r+2m\), where m is the number of Type II superleaves in Q. AddLeaf adds the remaining leaf x \(\in S\setminus R\) to a shared edge of \(t'\) and \(T|_{Q\cup R}\).

1. 1.

   Lemma 6 shows that AddLeaf (and therefore OCTAL) preserves existing shared edges between \(t'\) and \(T|_{Q\cup R}\), possibly excluding the edge where x is added.

    
2. 2.
   Now we consider what happens to the RF distance on the edge x is attached to. There are three cases: (i) x is not the first element of a superleaf (ii) x is the first element of a Type I superleaf or (iii) x is the first element of a Type II superleaf.

   • Case (i): If x is not the first element of a superleaf to be added to t, it directly follows from Lemma 7 that OCTAL will not change the RF distance when adding x.

   • Case (ii): If x is the first element of a Type I superleaf to be added, then x is attached to a shared edge in the backbone corresponding to some bipartition A|B existing in both C(t) and \(C(T|_R)\). Let \(e'\) be the edge in t s.t. \(\pi _{e'}=A|B\). Note there must exist an edge e in \(T|_{Q\cup R}\) producing A|B when restricted to just R. Hence, the bipartition \(\pi _e\) has the form M|N where \((M \cap R) = A\) and \((N \cap R) = B\). We need to show that \(M|N \in C(t')\).

   • By Lemma 6, any leaves from Q not attached to \(e'\) by OCTAL will preserve this shared edge in \(t'\).

   • Now consider when leaves from Q are added to \(e'\) by OCTAL. We decompose M and N into the subsets of leaves existing in either R or Q: let \(M = A \cup W\) and \(N = B \cup Z\). OCTAL will not cross a leaf from W with a leaf from Z along \(e'\) because this would require crossing the shared edge dividing these two groups: any leaf \(w \in W\) has the property that \((A\cup \{w\}) | B\) is a shared edge and any leaf \(z \in Z\) has the property that \(A | (B \cup \{z\})\) is a shared edge. Hence, any leaves added from Q that subdivide \(e'\) will always preserve an edge between leaves contained in W and Z on \(e'\).

   Thus, \(M|N \in C(t')\). Moreover, \((M\cup \{x\})| N\) and \(M | (N \cup \{x\})\) are bipartitions in C(T). AddLeaf roots T at x and removes the edge incident to x, creating rooted tree \(T^{(x)}\). We have shown that the edge incident to the root in \(T^{(x)}\) must be a shared edge, so adding x does not change the RF distance.

   • Case (iii): If x is the first element of a Type II superleaf to be added, we have shown in Lemma 4 that the RF distance must increase by at least two. Since AddLeaf always attaches x to some shared edge \(e'\), the RF distance increases by exactly 2 when subdividing \(e'\).

   Thus, OCTAL will only increase the RF distance by 2 if x is a new Type II superleaf.

    

\(\square\)

Combining the above results, we establish our main theorem:

Theorem 9

Given unrooted binary trees t and T with the leaf set of t a subset of the leaf set of T, OCTAL(T, t) returns an unrooted binary tree \(T'\) that is a completion of t and that has the smallest possible RF distance to T. Hence, OCTAL finds an optimal solution to the RF Optimal Tree Completion problem. Furthermore, OCTAL runs in \(O(n^2)\) time, where T has n leaves.

Proof

To prove that OCTAL solves the RF Optimal Tree Completion problem optimally, we need to establish that OCTAL returns an S-completion of the tree t, and that the RF distance between the output tree \(T'\) and the reference tree T is the minimum among all S-completions. Since OCTAL always returns a binary tree and only adds leaves into t, by design it produces a completion of t and so satisfies the first property. By Lemma 8, the tree \(T'\) output by OCTAL has an RF score that matches the lower bound established in Lemma 5. Hence, OCTAL returns a tree with the best possible score among all S-completions.

We now show that OCTAL can be implemented to run in \(O(n^2)\) time, as follows. The algorithm has two stages: a preprocessing stage that can be completed in \(O(n^2)\) time and a second stage that adds all the leaves from \(S \setminus R\) into t that also takes \(O(n^2)\) time.

In the preprocessing stage, we annotate the edges of T and t as either shared or unique, and we compute a set A of pairs of shared edges (one edge from each tree that define the same bipartition on R). We pick \(r \in R\), and we root both t and T at r. We begin by computing, for each of these rooted trees, the LCA (least common ancestor) matrix for all pairs of nodes (leaves and internal vertices) and the number \(n_u\) of leaves below each node u; both can be computed easily in \(O(n^2)\) time using dynamic programming. (For example, to calculate the LCA matrix, first calculate the set of leaves below each node using dynamic programing, and then calculate the LCA matrix in the second step using the set of leaves below each node.) The annotation of edges in t and T as shared or unique, and the calculation of the set A, can then be computed in \(O(n^2)\) time as follows. Given an edge \(e \in E(T)\), we note the bipartition defined by e as X|Y, where X is the set of leaves below e in the rooted version of T. We then let u denote the LCA of X in t, which we compute in O(n) time (using O(n) LCA queries of pairs of vertices, including internal nodes, each of which uses O(1) time since we already have the LCA matrix). Once we identify u, we note the edge \(e'\) above u in t. It is easy to see that e is a shared edge if and only if e and \(e'\) induce the same bipartition on R, and furthermore this holds if and only if \(n_u = |X|\). Hence, we can determine if e is a shared edge, and also its paired edge \(e'\) in t, in O(n) time. Each edge in T is processed in O(n) time, and hence the preprocessing stage can be completed in \(O(n^2)\) time.

After the preprocessing, the second stage inserts the leaves from \(S \setminus R\) into t using AddLeaf, and each time we add a leaf into t we have to update the set of edges of t (since it grows through the addition of the new leaf) and the set A. Recall that when we add \(s \in S \setminus R\) into t, we begin by rooting T at s, and then follow a path towards the leaves until we find a first shared edge; this first shared edge may be the edge incident with s in T or may be some other edge, and we let e denote the first shared edge we find. We then use the set A to identify the edge \(e' \in E(t)\) that is paired with e. We subdivide \(e'\) and make s adjacent to the newly created node. We then update A, the set of bipartitions for each tree, and the annotations of the edges of t and T as shared or unique. By Lemma 6, AddLeaf preserves all existing shared edges other than the edge the new leaf x is placed on, and these specific edges in E can each be updated in O(1) time. Furthermore, OCTAL places x on a shared edge, bifurcating it to create two new edges. Thus, just two edges need to be checked for being shared, which again can be done in O(n) as claimed. Thus, adding s to t and updating all the data structures can be completed in O(n) time. Since there are at most n leaves to add, the second stage can be completed in \(O(n^2)\) time. Hence, OCTAL runs in \(O(n^2)\) time, since both stages take \(O(n^2)\) time. \(\square\)

These results show that OCTAL was generally at least as accurate as ASTRAL-II at completing gene trees, and can be more accurate; this trend does not appear to be sensitive to the distance measure used to evaluate the accuracy of the completed gene trees. Within the scope of our study, the degree and frequency of improvement depended on the level of ILS, but not so much on the number of genes or on the reference tree, as long as the reference tree was estimated from the gene trees. Furthermore, using several techniques to produce the reference tree from the gene trees, including even a greedy consensus tree, produced reference trees that were as good as the true species tree in terms of the impact on the accuracy of the completed gene tree. However, a random tree produced very poor results. We also noted that OCTAL provided a clear advantage over ASTRAL-II under low to moderate ILS, but the improvement was smaller and less frequent under the high to very high ILS conditions. We offer the following as a hypothesis for the reason for these trends. Under low to moderate ILS, the true species tree is close to the true gene tree, and the estimated species trees (computed using ASTRID or the greedy consensus) are reasonably close to the true species tree; by the triangle inequality, the estimated species tree is close to the true gene trees. Therefore, when ILS is at most moderate, completing the estimated gene trees using the estimated species tree as a reference can be beneficial. However, under higher ILS, the true species tree is farther from the true gene trees, which makes the true species tree (or an estimate of that tree) less valuable as a reference tree. Despite this, we also saw that using estimated species trees as reference trees produced comparably accurate completions as using the true species tree as a reference, and that this held for both moderate and high ILS levels. Hence, OCTAL was robust to moderate levels of error in the estimated species tree. However, OCTAL is not completely agnostic to the choice of reference tree, since the random reference tree (which has close to 100% RF error) resulted in very poor performance.

OCTAL is a greedy polynomial time algorithm that adds species into an estimated gene tree so as to provably minimize the RF distance to a given reference tree. In our study, OCTAL frequently produced more accurate completed gene trees than ASTRAL-II under ILS conditions ranging from moderate to very high; however, the improvement under high ILS conditions was much lower and less frequent than under moderate ILS conditions.

There are many directions for future work. First, we compared OCTAL to ASTRAL-II, but ASTRAL-III [35] has recently been developed, and the comparison should be made to this new version of ASTRAL. OCTAL could also be compared to gene tree completion methods that are designed to handle gene tree heterogeneity resulting from gene duplication and loss [36], and these comparisons could be made on datasets that have evolved under multiple causes of gene tree discord (e.g., gene duplication and loss, horizontal gene transfer, and incomplete lineage sorting).

The current approach only adds missing species to the estimated gene tree, and so implicitly assumes that the gene tree is accurate; since estimated gene trees have some error, another approach would allow the low support branches in gene trees to be collapsed and then seek a complete gene tree that refines the collapsed gene tree that is close to the reference tree. This is similar to approaches used in [37–39], each of which aims to improve gene trees that use reference species trees, but are primarily (or exclusively) based on gene duplication and loss (GDL) distances. The optimal completion problem or the accuracy of the completed gene trees could also be based on other distances between trees besides the RF distance, including weighted versions [40] of the RF distance (where the weights reflect branch lengths or bootstrap support values), quartet tree distances, geodesic distances [41], or the matching distance. It is likely that some of these problems will be NP-hard, but approximation algorithms or heuristics may be useful in practice.

We did not evaluate the impact of using OCTAL on downstream analyses. Since missing data (i.e., incomplete gene trees) are known to impact species tree estimation methods using summary methods [21], this would be a natural next analysis. As an example, if the input includes some incomplete gene trees, a species tree could be estimated from the full set of gene trees and then OCTAL could use that estimated species tree as a reference tree to complete the gene trees. Then, the species tree could be re-estimated (using a good summary method) on the new set of gene trees, all of which are complete. This two-step process (completing gene trees using an estimated species tree, then re-estimating the species tree) could then iterate. It would be interesting to determine whether this improves the species tree, and if so under what conditions. It would also be helpful to evaluate the impact of completing incomplete gene trees when the genes are missing due to true biological loss rather than data collection issues, and hence also to see if OCTAL provides any helpful insight into gene evolution (such as better estimating the duplication/loss/transfer parameters).

Finally, there can be multiple optima to the RF Optimal Tree Completion problem for any given pair of trees, and exploring that set of optimal trees could be important. An interesting theoretical question is whether the set of optimal solutions admits a compact representation, even when it is large. From a practical perspective, the set of optimal completions could be used to provide support values for the locations of the missing taxa, and these support values could then be used in downstream analyses.