{'fraction_non_alphanumeric': 0.03658805752219902, 'fraction_numerical': 0.011069797449871814, 'mean_word_length': 4.402182885444373, 'pattern_counts': {'":': 0, '<': 0, '<?xml version=': 0, '>': 4, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 5, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Impacts of terraces on phylogenetic inference

Terraces are potentially large sets of trees with precisely the same likelihood or parsimony score, which can be induced by missing sequences in partitioned multi-locus phylogenetic data matrices. The set of trees on a terrace can be characterized by enumeration algorithms or consensus methods that exploit the pattern of partial taxon coverage in the data, independent of the sequence data themselves. Terraces add ambiguity and complexity to phylogenetic inference particularly in settings where inference is already challenging: data sets with many taxa and relatively few loci. In this paper we present five new findings about terraces and their impacts on phylogenetic inference. First we clarify assumptions about model parameters that are necessary for the existence of terraces. Second, we explore the dependence of terrace size on partitioning scheme and indicate how to find the partitioning scheme associated with the largest terrace containing a given tree. Third, we highlight the impact of terraces on bootstrap estimates of confidence limits in clades, and characterize the surprising result that the bootstrap proportion for a clade can be entirely determined by the frequency of bipartitions on a terrace, with some bipartitions receiving high support even when incorrect. Fourth, we dissect some effects of prior distributions of edge lengths on the computed posterior probabilities of clades on terraces, to understand an example in which long edges"attract"each other in Bayesian inference. Fifth, we show that even if data are not partitioned, patterns of missing data studied in the terrace problem can lead to instances of apparent statistical inconsistency when even a small element of heterotachy is introduced to the model generating the sequence data. Finally, we discuss strategies for remediation of some of these problems.

Abstract.-Terraces are potentially large sets of trees with precisely the same likelihood or parsimony score, which can be induced by missing sequences in partitioned multi-locus phylogenetic data matrices. The set of trees on a terrace can be characterized by enumeration algorithms or consensus methods that exploit the pattern of partial taxon coverage in the data, independent of the sequence data themselves. Terraces add ambiguity and complexity to phylogenetic inference particularly in settings where inference is already challenging: data sets with many taxa and relatively few loci. In this paper we present five new findings about terraces and their impacts on phylogenetic inference. First we clarify assumptions about model parameters that are necessary for the existence of terraces.
Second, we explore the dependence of terrace size on partitioning scheme and indicate how to find the partitioning scheme associated with the largest terrace containing a given tree.
Third, we highlight the impact of terraces on bootstrap estimates of confidence limits in clades, and characterize the surprising result that the bootstrap proportion for a clade can be entirely determined by the frequency of bipartitions on a terrace, with some bipartitions receiving high support even when incorrect. Fourth, we dissect some effects of prior distributions of edge lengths on the computed posterior probabilities of clades on terraces, to understand an example in which long edges "attract" each other in Bayesian inference.
Fifth, we show that even if data are not partitioned, patterns of missing data studied in the terrace problem can lead to instances of apparent statistical inconsistency when even a small element of heterotachy is introduced to the model generating the sequence data.
Finally, we discuss strategies for remediation of some of these problems. Among the most promising is the strategic deletion of a minimal number of taxa from the data in order to reduce terrace sizes.
(Keywords: phylogenetics, terrace, partitioned model, bootstrap, posterior probability) Phylogenetic trees with thousands to tens of thousands of species are becoming increasingly commonplace (Rabosky et al. 2013;Zanne et al. 2014). They can serve at least two purposes: quantifying and conveying the true scale and breadth of biodiversity, and providing statistical power to distinguish between alternative models of evolution (Wiens 2011;Boettiger et al. 2012;Chamberlain et al. 2012;Goldberg and Igic 2012;Marazzi et al. 2012;Christin et al. 2013;Davis et al. 2013). Reconstruction of large trees entails many challenges (Sanderson 2007;Izquierdo-Carrasco et al. 2011;Liu et al. 2012), including a recently discovered one: "terraces" (Sanderson et al. 2011). A terrace is a region in tree space in which all trees have precisely the same likelihood and parsimony score, which adds ambiguity to the "landscape" of trees ( Fig. 1) and complexity to tree inference. Terraces may have been overlooked among the inevitable small numerical differences that arise in computing the likelihood score on different trees, especially in large data sets. Indeed, the issue of whether two trees have exactly the same likelihood rarely arises in modern phylogenetic inference for this reason.
In trying to improve the efficiency of likelihood calculations in RAxML, Stamatakis and Alachiotis (2010) pointed out one important context in which likelihoods of different trees could be precisely identical: when the subtrees for different loci in a multi-locus partitioned analysis are the same. We suggested the term "terrace" for the resulting trees and exploited a variety of results concerning subtrees and supertrees to characterize these terraces, which in some data sets can be quite large (Sanderson et al. 2011). In general, trees with equal or nearly equal optimality scores can arise for many reasons, including lack of variable sites, homoplasy, and missing data. Terraces can arise in parsimony analysis when there are missing data or in likelihood-based analyses when missing data are distributed differently in different elements or blocks of a specified data partition.
Data may be missing for many reasons, ranging from sampling biases inherent in studies that mine GenBank (Driskell et al. 2004;Sanderson 2008) to more biological causes, as in the loss of plastid genes transferred to the nuclear genomes of some plants (Sabir et al. 2014), or the differential expression of genes found in EST libraries or transcriptomes (Letsch et al. 2012). Missing data affect phylogenetic tree reconstruction in many ways (Kearney 2002;Burleigh et al. 2009;Lemmon et al. 2009;Cho et al. 2011;Simmons and Freudenstein 2011;Wiens and Morrill 2011;Crawley and Hilu 2012;Simmons 2012a,b;Simmons and Goloboff 2013;Hinchliff and Roalson 2013;Roure et al. 2013;Simmons 2014;Simmons and Goloboff 2014), but many properties of terraces depend only on the pattern of "taxon coverage," the set of taxa for which data are present, in the different elements of the data partition Sanderson et al. 2010Sanderson et al. , 2011.
In practice, terraces are generally small or absent in phylogenomic studies in which the number of loci greatly exceeds the number of taxa (Hejnol et al. 2009: 94 taxa × 1487 loci; Salichos and Rokas 2013: 23 taxa × 1070 loci; Zwickl et al. 2014: 11 taxa × 473 loci), but when the number of taxa greatly exceeds the number of loci, the number of trees on each terrace can be extremely large and the resulting increase in ambiguity challenging to overcome (e.g., Pyron andWiens 2011: 2871 taxa × 12 loci;Smith et al. 2009: 55,473 taxa × 6 loci;Fabre et al. 2012Fabre et al. : 1265Rabosky et al. 2013: 7822 ray-finned fish taxa × 13 loci). Building on several basic mathematical properties of terraces we have characterized quantitatively (Sanderson et al. 2011), we extend our understanding of terraces in several directions, examining several new properties, the problems they induce, and strategies for overcoming them. We show that terraces can arise under more general conditions than we thought, and that they can be larger than believed before. We construct for parsimony a method to recover the maximal terrace for a given tree. Most importantly, we explore how terraces can affect confidence assessments and conventional views on how methods for characterizing ambiguity, such as consensus methods or bipartition support values, can be misled by them. Finally, we begin to extend results to the case in which the model partitioning scheme violates the sufficient conditions for terraces but still generates data sets with patterns of ambiguity related to our other results on terraces.

Definitions
Data.-Let D be a data matrix, such as a multiple sequence alignment, of n taxa ("rows") and l characters ("sites", "columns"), and let P be a partition of the columns into m elements or blocks, referred to colloquially throughout this paper as "loci," although blocks might be something like different sets of codon positions, etc. Let the taxon coverage matrix, C = C(D), be an n × m matrix where the ij-th element is 1 if any sequence data are present for taxon i and locus j, and 0 if the data are entirely missing (Fig. 2). This often happens simply because the locus was not sampled for that taxon. Note that there is no guarantee that taxa scored as "1" for a locus have phylogenetically informative sequence data-they merely have sequence data of some kind.
The taxon label set for D, L, is the set of all taxon names in D, and the label set, L j , for block j of D is the set of taxon names for which data are present for block j in partition P (i.e., the set of taxon names for which C ij = 1).
Trees and subtrees.-For any tree, T 0 , on the complete label set L (e.g., the best maximum likelihood tree found in a heuristic tree search based on D), the label set for any locus, j, L j (which may have taxa missing) induces a subtree of T 0 , which we write as T 0 |L j . Let Q(P, C, T 0 ) = T 0 |L 1 , T 0 |L 2 , ..., T 0 |L m , be the set of subtrees for this tree that are induced by the partitioning scheme and taxon coverage matrix-that is, T 0 |L j is the tree obtained from T 0 by removing any taxa that have data that are entirely missing for block L j , for each block L j of the partition.
A tree T is a resolution of T if T can be obtained from T by collapsing one or more edges of T , which transforms binary nodes to polytomies. A tree T on label set L(T ) displays another tree T on L(T ) (with L(T ) ⊆ L(T )) if T is equal to or a resolution of T |L(T ). Intuitively, this allows a larger tree to display a smaller tree even if that smaller tree is less resolved than the larger one. This extra generality is appropriate given the usual notion of polytomies as reflecting uncertainty rather than multiple speciation (Semple and Steel 2003).
Stands.-Here we define a new term not used in our previous work on terraces. A stand (of Taxon  coverage  matrix   A tree  Induced subtrees Stand/ Terrace Figure 2: Stands and terraces. The two-locus taxon coverage matrix at left (? indicates missing data) is partially filled and not decisive. For at least one tree (e.g., T ), and the two subtrees induced by the coverage matrix, there is more than one tree (the three at right) consistent with the subtrees. This collection of trees is a stand. In addition, because the parsimony score and likelihood score (with an unlinked model), will be identical on these three trees, the stand is also a terrace.
trees), S = S(P, C, T 0 ), is the set of all binary phylogenetic trees on leaf set L that display every subtree in Q(P, C, T 0 ). In the example in Figure 2, there are three trees that display the subtrees in Q, including the original T 0 .
Decisiveness.-The taxon coverage matrix is said to be decisive for tree T 0 and partitioning scheme P if T 0 is the only tree that displays all subtrees in Q (Steel and Sanderson 2010; Sanderson et al. 2010) -in other words, if the stand has only one tree (i.e. |S| = 1). In the example in Figure 2, the taxon coverage matrix is not decisive for T 0 .
Terraces.-Stands and decisiveness depend only on the coverage pattern, partition scheme and tree, {P, C, T 0 }; not on the actual sequence data, D, per se. However, under certain conditions, all trees in a stand, S, have precisely the same optimality score with respect to D, in which case we call the stand, S, a terrace, denoted T (Fig. 2). In particular, when parsimony is used as the score, all trees in S have the same score for any partitioning scheme, P, and S is always a terrace. The set S is also a terrace when likelihood is the optimality criterion if the likelihood is determined by a model for each locus in the partition, P, that is "unlinked" (Sanderson et al. 2011). An unlinked (also known as "partitioned", or "fully partitioned" in the literature) model has different sets of parameters, including edge lengths, for each locus (Hess and Goldman 2011;Hedin et al. 2012;Xi et al. 2013).

Properties of terraces
Because all terraces are stands, a number of results derived for stands, based only on trees, subtrees, and coverage patterns, are helpful for characterizing terraces. Some of these results hold only for rooted trees, but many problems can be effectively "rooted" as long as there is one taxon in C that is sampled for all blocks of the partition, in which case that taxon can serve as an "operational" root for the purposes of an algorithm. The number of trees on a stand (or terrace) can increase exponentially with the size of the tree (Semple 2003). Despite this, several properties of terraces make them more tractable than they might otherwise seem, mainly because several summary statistics can be obtained directly without any computation involving the data matrix, D (Sanderson et al. 2011). For example, all trees on a terrace can be enumerated without recalculation of optimality scores. This takes advantage of an algorithm due to Constantinescu (1995), with a running time that scales linearly with the size of the terrace (rather than, say, exponentially with the size of the tree).
The trees on a terrace can also be summarized by a strict consensus tree (Gordon 1986) or Adams consensus tree, either of which can be constructed in polynomial time.
This last claim (that we can sidestep the enumeration of trees on the terrace again, or any further search using the data) is not obvious, but it holds in general (for rooted trees). In the case of strict consensus this was shown in (Steel 1992) (using results from Aho et al. (1981)), while for Adams consensus, it relies on a particularly elegant result due to Bryant (1997) (Theorem 6.2) which states that the Adams consensus of a terrace is equal to the so-called BUILD supertree of the set of subtrees for each locus (i.e., the set Q(P, C, T 0 )), and this supertree can be computed quickly by the algorithm of Aho et al. (1981). Note that the Adams consensus tree displays each of those subtrees; that is, it is identical to (or resolves) them if extraneous taxa and edges are removed. Finally, testing if two trees, T 1 and T 2 , are on the same terrace can be done efficiently merely by checking for equality of the pair of induced subtrees, T 1 |L j = T 2 |L j ; for every locus j (cf. Day 1985).
Terraces are reminiscent of tree islands but are contained within them, at least for rooted trees. An island is a region of tree space with optimality score better than some threshold, separated from other such regions by regions of lower score (Maddison 1991;Salter 2001). Here a "region" is a set of trees that can be enumerated by a series of topological rearrangements that do not leave the region. Terraces are always wholly contained within a tree island, because all trees on a terrace of rooted trees can be reached by a series of nearest neighbor interchanges (NNIs) between trees of the same optimality score (Bordewich 2003;Sanderson et al. 2011). For unrooted trees this property does not necessarily hold. For rooted trees, then, tree space can be thought of as a collection of terraces (on islands), the size and topological position of which are determined entirely by the partition P and taxon coverage matrix C. Only the heights of the terraces depend on the data. To the extent that a tree "island" is a useful metaphor, it may be best to envision it as a rough landscape covered with terraces of different sizes and heights (Fig. 1).
New Results on Terraces 1. Terraces occur in likelihood inference under a less restrictive set of assumptions In our previous work, we showed that terraces can occur in maximum likelihood inference whenever the substitution models for separate loci in partition P have separate parameter sets (Sanderson et al. 2011). Such models are known as "partitioned" or "unlinked" models. To be precise, suppose the substitution model at locus i consists of two sets of free where M i is a set of free parameters associated with the substitution rate matrix alone, and E i is a set of free edge rate parameters. An edge-unlinked (EUL) model has free edge rate parameters E i for each locus. On the other hand if any constraints are imposed to reduce the number of free edge parameters (such as those of the form E i ∝ E j , for at least some i, j), we will say the model is edge-linked (EL); likewise for substitution matrix rate-linked (RL) and rate-unlinked (RUL) models. Later we will refer to two simple EL models: the "homotachy" model, in which E i = E j for all i, j, and the "proportional" model in which E i = αE j and α varies across loci to reflect locus specific overall rates.
In terms of these definitions, we showed (Sanderson et al. 2011) that terraces can arise when models are simultaneously EUL and RUL. This was a sufficient but not necessary condition. In fact, as we show now, terraces can arise even when models are EUL and the substitution rate matrices are constrained to be the same across loci.
Proposition 1. Let S(C, P, T 0 ) be a stand of trees containing tree T 0 for partitioning scheme, P and taxon coverage matrix, C. For any model M for computing the likelihood score that is edge-unlinked (EUL) and where the M i parameters are constrained to be identical (i.e. M i = M for all i), all trees in S have the same maximum likelihood score, ϕ(M, T 0 |D), as T 0 , and hence S is a terrace.
The proof of this result is presented in the Appendix. This finding implies that it is the lack of commonality between edge length parameters across loci that is necessary for a stand of trees to be a terrace under likelihood. Rate matrices between loci need not have different parameter sets. Henceforth we will refer to this broader class of maximum likelihood inference with just an EUL assumption as "ML-EUL."

Maximum size of stands and terraces
Given a coverage matrix, C and partition scheme, P, we can calculate the size of the stand, S containing some tree, T 0 , without reference to the data, D. If the optimality criterion is MP or ML-EUL, then S is also a terrace and the size of the terrace will be the size of the stand. Generally, the presumption is that P is chosen to reflect meaningful biological aspects of the data, such as different loci, or disjoint sets of codon positions.
Nonetheless, the choice is arbitrary and it is possible that this terrace for P is actually imbedded in a larger terrace (perhaps with a different likelihood score) under a different partitioning scheme. First, we define a special kind of partitioning scheme and then show that this will correspond to the largest terrace among all possible partitioning schemes.
Given a taxon coverage matrix, C, consider any partition P = {B 1 , . . . , B l } of the columns of the data matrix, D (i.e. the sets B i are disjoint subsets of columns, which cover every column). For a given block B i of P, let L i = L(B i ) be the taxa that are present for at least one column in B i . Now suppose we have another partition P = {B 1 , . . . , B l }. We say that P refines P if each block of P is a subset of some block of P (equivalently, each block of P is either equal to a block of P or is the union of two or more blocks of P ).
However, a refinement that includes breaking B j into say, K, smaller sets, should be disallowed if any of the labels sets, L(B i k ), are duplicated, as otherwise this would allow a trivial refinement of any partition into one in which each block consists of just a single column in the data matrix, D. With this restriction in place, there will be a unique (and usually nontrivial) maximal partition for any taxon coverage matrix, C, that cannot be refined further, we denote this as P max . Without this restriction the maximal partition would simply consist of l blocks, one for each column in the data matrix.
This maximal permitted partition P max can be described formally as follows. Let Proposition 2. If P refines P then S(C, P, T 0 ) ⊆ S(C, P , T 0 ); in particular, the former set is never larger than the latter.
Proof. Let L j = L(B j ), the taxa present in block B j for j = 1, . . . , l . Suppose that T is a tree in the stand of T 0 relative to P. This means that T |L i = T 0 |L i for i = 1, . . . , l.
Consider a block B j of P . Since P refines P, B j is a subset of some block of P, say block B i , and so L j ⊆ L i . In that case: Thus, T |L j = T 0 |L j for each block L j of P , and so T is in the stand of T 0 relative to P .
Corollary 3. For any permitted partition P we have: Proof. This follows from Proposition 2, since P max is a refinement of every permitted partition of the loci.
Thus, the largest stand containing a given tree can be obtained from the maximal partitioning scheme allowed for a given taxon coverage pattern, C. Under MP or ML-EUL inference, this will also be the largest terrace containing this tree.
For example, suppose a multiple sequence alignment is first partitioned into two loci, P I = {I 1 , I 2 }, where the loci are labeled by the partition name with a subscript for the locus number, as shown in Figure 3. This partition induces a stand of 13 trees containing the tree shown. However, perhaps locus I 2 actually consists of three biologically meaningful blocks: two exons and an intron, and suppose further that the label sets for the exons are the same but differ from that for the intron. Our second partition is then described by P II = {II 1 , II 2 , II 3 } (and has the indicated distribution of columns: Fig. 3). Formally, P II is a refinement of P I . The stand induced by this new three block partition, P II , has 23 trees ( Figure 3). This is the maximal partition for these data, and 23 is the size of the largest stand containing the given tree. Also, from Proposition 2, the smaller stand of trees is a subset of the larger.
Under MP and ML-EUL inference, these stands are also terraces, so the largest terrace containing the indicated tree has 23 trees on it. The sets of trees in these terraces arising from different partitioning schemes are the same for MP and ML-EUL inference, but there is an interesting distinction with respect to their optimality scores. For MP, the situation is particularly simple. For two partitioning schemes in which P II is a refinement of P I , the two terraces have the same parsimony score, so it is reasonable to visualize the smaller terrace is being literally imbedded in the larger one at the same "elevation" in the landscape of tree space ( Figure 4)  Partition elements/blocks/"loci" are labeled with 1,2, 2(1), etc. Label sets for each block are indicated in curly brackets. Partition II is a "maximal partition". There are 13 trees on the terrace for Partition I and 23 trees for Partition II (for parsimony), each containing T . partition, and if not, also check P max . The latter will be a more accurate representation of the actual extent of equally optimal trees around the given tree.
For ML-EUL, the situation is more complex. The likelihood scores of the trees in the terraces for partitioning schemes P I and P II could well be different, but all trees within either terrace will still have the same score once a partitioning scheme is fixed. This means changing a partitioning scheme to the maximal partitioning scheme, for example, not only expands the set of trees on the terrace, but potentially changes the likelihood score as well, unlike in parsimony. c b d a Figure 4: Refined view of the landscape of terraces with respect to different partitioning schemes. Regions 1-4 are different regions of tree space. Region 4 is a terrace with the highest likelihood. Region 2 is a subset of region 3 and both are terraces identified from two different partitioning schemes. Region 3 corresponds to a maximal partitioning scheme.

Impact on confidence assessment: bootstrap proportions
Terraces add ambiguity to phylogenetic tree inference directly by increasing the number of equally optimal solutions. In the next two sections we discuss how they also alter estimates of the quality of trees. The bootstrap (Felsenstein 1985) case is especially clear, as we demonstrate with detailed discussion of a simple example (Fig. 5). In keeping with many results related to terraces, much of this can be predicted by reference to the stand, S, induced by the taxon coverage matrix and partition scheme for a given tree.
Henceforth, we assume inference is with MP or ML-EUL and that results for stands carry over to terraces exactly. We simplify the problem further by letting sequence lengths within loci be so long that the only variation in bootstrap replicates has to be due to the existence of terraces . Specifically, let T * be a tree with four subtrees, X,Y,A,B, in a topology such that there is a bipartition given by XA|BY (Fig. 5). Assume that X and Y are leaf taxa, but A and B are subtrees with n A and n B leaves respectively. Further assume that sites are all evolving on T * according to a model with "reasonable" parameter values such that, together with the assumption that the number of sites, l, is large enough, if complete taxon coverage were available, both ML-EUL or MP would infer the correct tree with ∼ 100% accuracy. Assume there is a partition, P, of data D into two loci, 1 and 2.
The taxon coverage matrix, C is such that taxa in two of the subtrees, X and Y, are sampled for both loci, but taxa in subtree A are only sampled for locus 1, and subtree B only for locus 2 (i.e., the same as in Fig. 2, except no outgroup is present).
With complete sampling and "ideal" simulation conditions ("F84" substitution model with equal base frequencies in Seq-Gen [Rambaut and Grassly 1997]; l = 10,000 nt; all edge lengths = 0.1 substitutions/site; constant rates across sites), the bootstrap proportion (BP) for XA|BY are indeed 100%, whether inferred with MP or ML-EUL .
However, with the partial taxon coverage, C, outlined above, the true tree, T * , occurs in a stand, S, with other trees. The number of trees in |S| grows with n A and n B (Fig. 6), and many conflict with T * specifically by interleaving leaves from subtree A with leaves from subtree B, such that AB|XY is seen instead of XA|BY. In fact, perhaps surprisingly, these incorrect bipartitions rise in frequency in the set of bootstrap bipartitions with increasing n A and n B (Fig. 6). For example, if n A = 3 and n B = 3, then the stand has 107 binary trees, but only 2 of these are consistent with the XA|BY bipartition in the true tree, T * .
Figure 5: Impact of terraces on the bootstrap. True tree and taxon coverage as in Fig. 2 except outgroup removed and taxa A and B are now subtrees with n A taxa (present only for locus 1) and n B taxa (present for locus 2) respectively. As n A and n B get larger the bootstrap proportion for the incorrect bipartition at right goes to 100%.
Since topological error due to the sequence data has been factored out by long sequence length, bootstrapping will randomly sample the 107 trees on the terrace, and assuming one tree per replicate is kept (see below), the BP for the incorrect bipartition AB|XY will To explore this small example more generally, we need to know how many rooted binary trees correctly display the subtrees T * |L A and T * |L B , where we have changed notation a bit and L A refers to the subtree having just the leaf taxa in clade A, and likewise for B. Notice that the taxon sets L A and L B are disjoint. Let n = n A + n B , and R(k) = 1 × 3 × 5 × · · · × (2k − 3) be the number of rooted binary trees for k leaves. If n A or n B ≤ 3, then the number of trees on the terrace is a result that depends only on the two clade sizes. The "2" term corresponds to the two trees on the terrace with bipartitions of XA|BY and XB|AY. The term on the right corresponds to the set of trees all of which have the bipartition AB|XY, but there are several ways to interleave the A and B clades and still display the original two subtrees. If both A and B are larger than three leaves, the expression depends on the topology of the subtrees (i.e., the entire tree), not just their sizes, a result due to Constantinescu and Sankoff (1986), so a more complex calculation for Equation (1) would be necessary. On average, however, across a uniform random sample of clade topologies, the mean number of trees converges to the expression in Equation (1). Thus, for certain patterns of partial taxon coverage, the bootstrap proportion for a given clade may tend to 0 or 1 depending not on the sequence data but on the number of taxa (and possibly tree topology) and taxon coverage pattern in specific clades in the tree.
These results implicitly assume that bootstrap protocols sample trees from a terrace uniformly at random. They may not. Heuristic search strategies may not be constructed to guarantee this. Moreover, the very notion of "equality" of the likelihood score in numerically intensive likelihood calculations is problematic, and maximum likelihood reconstruction in RAxML (Stamatakis 2014) and GARLI (Zwickl 2006) operates on the implicit assumption that equally optimal trees either do not exist or cannot be unambiguously identified during tree search. For maximum parsimony inference, the situation is a little clearer. With integer value parsimony scores, maximum parsimony inference can certainly identify sets of trees, and there are then at least two choices for treating them. PAUP (Swofford 1999) uses a "frequency within replicates" (Davis et al. 2004) approach in which the the frequency of a bipartition in the equally optimal trees found in a given replicate is summed across replicates. This procedure would effectively treat the trees on a terrace in the same way as our analytical results above. However, an alternative is a more conservative "strict consensus" approach in which a set of equally optimal trees must all have the bipartition for it to be regarded as present in that replicate (Davis et al. 2004;Simmons and Freudenstein 2011). Were we to adopt such a rule here, no bipartitions involving A, B, X, and Y would be supported, a very different result.

Impact on confidence assessment: Bayesian posterior probabilities
Given the issues raised by terraces for MP and ML-EUL, it is reasonable to expect some effects on calculation of Bayesian posterior probabilities, but in this case the impact is determined by an interaction between the taxon coverage pattern, tree topology, and priors on edge lengths. Let us simplify the bootstrap example to a tree with four leaves (i.e., This phenomenon can be traced ultimately to how Bayesian inference handles missing data. We can begin to understand this by considering the the simplest context of partial taxon coverage: with a trivial "partition" consisting of only a single block and having one leaf taxon, x, with missing data throughout that block. In other words, this is a data set in which one leaf, x, has no data at all. Consider a tree, T , with branch lengths, containing taxon x as a leaf, and the tree, T −x (with its inherited branch lengths) obtained by deleting leaf x (and its incident branch) from T . If the sequence of x consists entirely of "?"s, the parsimony and likelihood scores of any tree, T , obtained by attaching x to some edge of T −x is the same. In other words, a method such as maximum likelihood, which relies only on these scores is unable to decide the position of x.
However, in Bayesian inference the position of x may also be influenced by prior probabilities. In fact, we will now show that for the usual exponentially distributed prior on edge lengths used in phylogenetic inference, the posterior probability for different placements of x in the tree is determined entirely by the length of the edge to which x attaches; if these differ, then some placements of x in this tree will have higher posterior probabilities. Ultimately this will affect posterior probabilities of trees on terraces in more complex partitioning schemes with partial taxon coverage.
Consider data generated by a reversible Markov process on a binary phylogenetic tree T with branch lengths assigned, and analyzed under a Bayesian approach in which (i) all rooted binary phylogenetic trees have the same prior probability, and (ii) an exponential prior applies independently across the branch lengths. Suppose we have a binary phylogenetic tree T that satisfies T −x = T −x (i.e. T has the same or different placements of leaf x in T −x ). The following result describes the posterior probability of T (in the limit of long sequences) as the proportion of the total tree length of T that comes from the edge of T −x to which x would attach in order to produce T . T −x , and e x is the edge of T −x that x must attach to in order to produce the tree topology T .
The proof of this result is presented in the Appendix. Intuitively, the theorem says that the relative chances of a leaf taxon attaching to two alternative edges of a tree are (for long sequences) given by their relative edge lengths: the longer edge will "attract" the taxon with missing data. This result may seem like something of a curiosity for a trivial partitioning scheme consisting of just a single locus, and a leaf taxon with no data, but it begins to explain why the posterior probability of trees on a terrace may be different even if their likelihoods are the same. In particular it pinpoints the important role that edge lengths can play in influencing subsets of the trees on a terrace to have higher posterior probabilities, as seen in the example in Figure 7.

Extensions to linked models
Our final result extends results on terraces to the setting in which there is some linkage among edge rate parameters across the loci in a data partition. We describe an example in which conditions are present such that terraces would occur for ML-EUL inference, but when an EL model is used instead, the accuracy of phylogenetic inference is negatively affected. Specifically, the relative likelihood scores can be re-ordered so that an incorrect tree is favored, much as in the Bayesian case described above.
The example we describe involves sequences that exhibit heterotachy; that is, the pattern of edge length variation among edges is different between loci. We describe a simple case of this in a 4-taxon data set for two loci. Assume all edges for both loci have evolved with the same rate, u = 0.1 substitutions/site, except for a single terminal edge for one of the loci, which has rate v, ranging from 0.1 to 1.0 substitution/site (Fig. 8). Each locus in the partition has 5000 sites and sequences are evolved in Seq-Gen (Rambaut and all three trees are the same. We examined two EL models. The first ("homotachy") model assumes the substitution rate matrix is the same for both loci and the edge rate parameters are the same for both loci; the second ("proportional") model assumes that the substitution rate matrix is the same for both loci but that edge rates for all edges at the second locus are strictly proportional to the corresponding edge rates for the first locus.
This allows rate variation between loci, but is not "heterotachy" per se. Likelihood calculations were carried out using PAUP* 4.0 (the "proportional" model implemented with the "siterates" command: Swofford 1999). Under these conditions the likelihood scores for the three possible binary trees are different for both models (Fig. 9). In particular, if v is more than about 4u, then an incorrect tree has the highest likelihood ( Fig. 9). In other words, use of an EL model for inference does not just avoid the ambiguity of the terrace phenomenon that would have been seen using EUL inference, it is positively misleading. Under these conditions, using an EUL model for inference and dealing with the resulting ambiguity might be a preferable (if more conservative) procedure.
This case clearly involves model misspecification, since the EL model described above is not the model actually generating the sequences. Interestingly, if we "fix" the data set by filling in the missing entries with data generated under our heterotachous model, the correct ranking is restored under EL inference, even with this model misspecification. By the same token, eliminating heterotachy by setting u = v in the generating model, but retaining missing data according to the specified partial coverage pattern, also restores the correct ranking of trees when using the EL model. Thus, the re-ranking of trees we observe in this example is not a product of either partial coverage alone or heterotachy alone. Both are required, and the correct ranking of the three trees' likelihood scores can be rescued by eliminating one or the other.
Although these results imply strongly that ML inference can be statistically inconsistent (i.e., converging to the wrong tree) under the combination of missing data and heterotachy, we have not been able to formally prove this in any sort of generality.
However, some aspects of the converse are provable. For data evolving according to our heterotachy model, both parsimony and likelihood inference using a homotachy model are statistically consistent. The proof that parsimony will be consistent for the heterotachy model with all taxa present, and the branch lengths indicated (for any u, v) follows from a straightforward application of Theorem 8.7.1 in Semple and Steel (2003), which shows that parsimony is consistent for tree 1 and for tree 2 with those branch lengths (for any u and also any v). Since parsimony is a linear scoring scheme, if it is consistent on each block of a partition it is consistent on all the data.
To prove the same for likelihood, consider first a single site -if in any tree just one branch length is changed, say a pendant branch length from u to v, and if z is the state at this pendant leaf, x the state at the other end (interior node) of this edge, and Y is the collective set of states at the rest of the tree, then by the Markov property, if proportion α of sites evolve under tree 1 (and 1 − α under tree 2) then for the heterotachy model we have where P r s (Y, z|x) denotes the conditional probability of leaf states Y and z, given that interior node is x, and s is the length of the pendant edge in question.
However, if we set the length of this pendant edge to θ, where (for the 2-state symmetric model) for all x, y, Z, and so P r(Y, z) = P r θ (Y, z) and in particular, P r(Y, z) (=site pattern distribution under the heterotachy model) matches exactly the probability distribution of a homotachy model in which the length of the pendant edge is a particular value θ chosen by (3) and lying somewhere between u and v. In other words, this heterotachy model produces site pattern probabilities the same as a homotachy model on the same topology and with an intermediate branch length. These results are based on a single site, but it follows by Theorem 5.1 of Chang (1996), that if we now use ML (on sequences from the model) under the homotachy model this will be a statistically consistent estimator of tree topology.

Discussion
The new results on terraces presented here are relevant to an increasingly common setting for phylogenetic analysis in which large multi-locus data sets with significant numbers of missing sequences are combined with models that partition the data in various ways during likelihood or Bayesian inference. They are also relevant to analyses using maximum parsimony irrespective of partitioning schemes, which remains a computationally attractive method for analyses of large data sets (e.g., Yoder et al. 2013). Previously we had shown that the terraces of equally optimal trees surrounding any specific tree could be quite large, although it was possible to exploit certain mathematical results on subtrees and supertrees to help characterize these sets of trees. In the present paper we found it useful to define "stands" of trees and distinguish between stands and terraces. Stands are collections of trees that all display the subtrees associated with the separate blocks within a partitioned data set. Under some conditions the trees in a stand can all have the same optimality score, in which case they are called terraces, to reflect their constant "elevation" along a vertical axis corresponding to that score.
If the score is the tree's parsimony score, a stand is always a terrace. The first new result we described specifies conditions under which this will also be true for the likelihood score. We showed that a sufficient condition was that the parameters describing edge lengths in different partition blocks must be independent ("edge unlinked"). Previously we had also assumed parameters of the rate matrix to be unlinked (Sanderson et al. 2011).
This is sufficient but not necessary. Because practice in current phylogenetics of multi-locus super matrices ranges from using completely linked models to highly partitioned ones that unlink both rate matrices and edge parameters, impacts of these divergent strategies should be assessed. Our results indicate it is possible that edge unlinked models can induce potentially large terraces around any given tree found in a reconstruction. This raises the deeper question of whether this ambiguity accurately reflects phylogenetic uncertainty or is best ameliorated by adjustments to the model (setting aside the possibility of simply filling in the missing data).
The second new result indicates that the stand associated with one partitioning scheme can actually be part of a larger stand under a different partitioning scheme. For parsimony, this leads to strong conclusions: it is possible to easily identify the extent of the largest terrace associated with any given partitioning scheme simply by constructing the so-called maximal partitioning scheme for the data set (a simple product of the distribution of missing data in the matrix), and then characterizing the terrace for that. This places a useful upper bound on the extent of the ambiguity associated with partial taxon coverage, and therefore should serve as the default terrace reported for any specific tree when parsimony is the optimality score.
These results are fundamentally based on the properties of stands. However, in considering the optimality scores of the trees on a stand, it is worth remembering that even a "maximal" terrace may be imbedded in a yet larger collection of trees with the same optimality score owing to homoplasy, constant sites, or possibly other factors in the data not related to partial taxon coverage. The maximal number of trees on a terrace provides a lower bound on the actual number of equally parsimonious trees for T , and it can be calculated without reference to the sequence data proper, D. This may be relevant for tree search heuristics, which are clearly challenged by tree landscapes such as these (Goloboff 2014). For example, rather than doing tree rearrangements to find equally parsimonious trees, it should be possible to enumerate them directly. These can form "seeds" from which to continue searching using expensive tree score computations. This applies to likelihood inference as well, conditional on the partitioning scheme.
The remaining new results address various aspects of the impact of terraces on assessing tree accuracy. The third and fourth results refer to situations in which bootstrap proportions and Bayesian posterior probabilities, respectively, are determined in large part not by the information con tent of the sequence data but by the pattern of partial taxon coverage. These happen for different reasons. In bootstrapping, the trees sampled from a terrace end up reflecting the frequencies of trees on that terrace, which in turn reflects the pattern of taxon coverage and tree shape. One consequence is that a clade can be highly supported in such an analysis because most trees on the terrace (all with equal likelihood scores) have that clade. It is possible to view this as either a "feature" or a "bug", but the bottom line is that one might not be expecting such sensitivity to arise from factors other than the sequence data. Simmons and colleagues (Simmons and Freudenstein 2011;Simmons 2012b;Simmons and Goloboff 2014) have noted a number of impacts of missing data generally on bootstrap estimates and have been critical of its application in this context, suggesting several strategies to check well supported clades for spurious support.
We suspect some but not all of their observations in real data sets (e.g., Simmons and Goloboff 2014) are due to issues related to terraces, but other factors involving homoplasy and phylogenetic signal proper are no doubt also involved.
Similarly, Bayesian inference can be heavily influenced by posterior probabilities not being the same across trees on a terrace (again, despite the likelihood being the same). In this case, the result is largely explained by the influence of priors. With standard exponentially distributed prior probabilities on edge lengths, the probability that a single taxon with all its data missing will attach to a particular edge of the tree is determined by that edge length in relation to the length of the entire tree. Absent any other information, Bayesian inference will place a taxon on the longest edge of the tree (though not necessarily with high probability). Again, this may be viewed as desirable or not, but it has the downstream consequence that trees in a stand can have very different posterior probabilities, which are determined by a fairly opaque convolution of the priors, edge lengths, partial coverage patterns, and combinatorics of terraces and tree shapes.
A reasonable rejoinder to some of these concerns is that one should simply avoid edge-unlinked partitioning schemes in likelihood or Bayesian inference (and maximum parsimony) and enforce some more homogeneous model during statistical inference. Of course, this kind of underparameterization may cause its own problems, but we pursued briefly whether it at least can ameliorate some of the issues we have described with terraces. The results presented in the last section suggest reasons to proceed carefully.
When sequence data are generated on a tree with even a small amount of heterotachy between loci, an edge-linked model is used for inference, and there is only partial taxon coverage, then the likelihood score for the correct tree can actually be worse than for an incorrect tree. Simulations with long sequence lengths imply this is an instance of statistical inconsistency. As stated this might seem like simply another case in which model misspecification causes problems, but in this instance, it is clearly a negative interaction between model misspecification and missing data, because the problem can be avoided by fixing either aspect.
Much further work is needed to narrow down the precise effects of terraces on confidence estimation and accuracy of inference. In the meantime, however, there are options for remediation. Setting aside the strategy of acquiring the missing sequences to eliminate partial taxon coverage entirely, which may be expensive or impossible depending on availability of DNA samples, there are computationally promising approaches.
Previously we posed the "maximum defining label set" (MDLS) problem (Sanderson et al. 2011), in which we seek the smallest number of leaf taxa to delete from the coverage matrix such that the stand for a given tree is reduced to a single tree. This has an exact and efficient solution for two loci, and experiments with data sets indicate there are interesting instances in which elimination of relatively few taxa can solve the problem, even while leaving a significant amount of missing data. However, there is no known exact solution to the case of three or more loci. The good news is that simple heuristics in the case of three or more loci can find solutions that eliminate terraces (Sanderson et al. 2011); they just may not be optimal (it may have been possible to do the same thing and keep more taxa in the matrix).
To highlight the impacts clearly, most of our results were in the context of sequence data sets so large that the only error was due to partial taxon coverage. In real data sets there is also error from the finite sample taken from the substitution process. This translates into a broadening of the bootstrap or posterior distribution of trees. In addition, there may be distinct terraces associated with each sample tree taken from these distributions; and there may be a distinct MDLS solution for deleting some set of taxa for each of these trees. How do we integrate across this information to make headway in reducing the overall impact of terraces? A simple but conservative fix might be to replace any sampled tree in a bootstrap replicate or an MCMC run with the strict consensus of the maximal stand in which that tree is imbedded. Then any clade on that tree is present in all trees on the terrace. Although this would tend to reduce the false positive clades uncovered in an analysis, it might lack sensitivity. Clades present at an overwhelmingly high frequency on the terrace but not quite 100% will be missed by a strict consensus.
Another approach would be to rely on other sorts of summaries about trees than consensus. For example, a terrace could be characterized by the average dissimilarity among its trees, based on a measure of distance between trees, such as the Robinson-Foulds (RF: Robinson and Foulds 1981) distance. Then, perhaps a more synthetic assessment of a confidence set of trees sampled from bootstrap replicates or posterior distributions could be to explore the "distances" between terraces, using a measure of distance between sets, such as the Hausdorff distance (Yu et al. 2014). The Hausdorff distance is small when each tree on one terrace is close (measured here by RF distance) to some tree on the other terrace.
This might be true even if the average RF distance between trees within a terrace is much larger. Whether this approach ultimately proves promising or not, some means to characterize more fully the relationships between entire sets of trees seems to be a necessity when terraces are commonplace.
Appendix: Mathematical proof of Proposition 1 and Theorem 4 Proof of Proposition 1 Consider an analysis where for each locus i we are free to select E i -parameters, but the M i parameters are constrained to be identical (i.e. M i = M for all i). Let ϕ(T ) denote the log-likelihood of tree T (having leaf set L) for the data D. Then assuming the loci evolve independently (conditional on the parameter choices) we have: where k is the number of loci, 'sup' refers to supremum (i.e. maximum if it is attained, else its limiting value) as we search over M and the E i parameter spaces, and where (E i ) is short for (E 1 , . . . , E k ). Now, (notice that T and (E i ) on the left has been replaced by (T |L i ) and E i on the right).
Combining these last two equations gives: and so lies in the same terrace at T ). Then from Eqn. (4) we have: so T is an ML tree also. In other words, all trees on the same terrace as T are ML trees.
This completes the proof. 2

Proof of Theorem 4
We first begin by defining more formally some of the notions mentioned earlier.
• Given any binary phylogenetic X tree T , and any taxon x from X, consider the tree T −x that is obtained from T by deleting leaf taxon x and its incident edge e(x). Note that each edge of T −x corresponds to an edge of T , except for the edge e −x of T −x which corresponds to the two edges (e 1 , e 2 ) of T , that are incident with e(x) in T . In this way, if T comes equipped with an branch length assignment l (so l(e) is the length of edge e), then the induced branch length l −x function for T −x is thus given by: l(e 1 ) + l (e 2 ), if e = e −x .
• For any data set D that consists of a sequence of k aligned site patterns on X, and any taxon x ∈ X, let D −x denote sequence of k aligned site patterns on X − {x} obtained by deleting the sequence for x.
• Now suppose that the sequence sites in D have evolved i.i.d. on some fixed binary phylogenetic X−tree T with branch length assignment λ, under a reversible Markovian process. Thus the sites in D −x evolve i.i.d. on T −x with branch length assignment λ −x .
We wish to apply a Bayesian approach to compare different placements of the taxon x into T −x given the censored data D −x . We assume a prior probability distribution on the set of binary phylogenetic X-trees with branch lengths for which: (i) each binary tree has the same probability (i.e. the 'PDA distribution'); (ii) edge lengths are independent exponential random variables.
Without loss of generality (by rescaling) we may assume that the mean of the exponential distribution in (ii) is 1.
Suppose we have two binary phylogenetic X-trees T and T that satisfy . One (or neither) of these trees might be T . We are interested in the ratio of posterior probabilities P(T |D −x ) P(T |D −x ) . The following result states that for long sequences this ratio converges towards the ratio of the lengths of the two edges of T −x to which the missing taxon (x) is attached. To establish Theorem 4, we prove the following result.
Theorem 5. For data generated by a reversible Markov process on a phylogenetic X-tree T with branch length assignment λ, consider, for any x ∈ X, any two phylogenetic X-trees T and T obtained by attaching x to edges of T −x of length l and l respectively. Then the ratio P(T |D −x ) P(T |D −x ) converges in probability to l /l , as the sequence length k grows.
Notice that Theorem 5 implies Theorem 4 since the statistical consistency of Bayesian phylogenetics under identifiable models implies that any tree which is different from T −x when x is deleted has a posterior probability that converges to zero as the sequences length k grows (thus P(T |D −x ) sums to 1 as we sum over all binary trees T that agree with T up to the placement of x). Thus, the the remainder of our argument is tailored towards proving Theorem 5 under the same conditions stated for Theorem 4 (in particular, conditions (i) and (ii) in the preamble to that theorem).
From Bayes' identity we have: Now, P(T ) = P(T ) (by assumption (i)), and so, from Eqn. (5) and the analogous identity for P(T |D −x ): Moreover, (7) where l is the set of branch length assignment on T −x , and where f (l) refers to the density of the branch lengths on T −x that is induced by independent exponential prior branch lengths on T (Γ is the set of possible branch lengths of T −x ). Similarly, where l is the branch length assignment on T −x (= T −x ), and where f (l) refers to the density of the branch lengths on T −x induced by the independent exponential priors on T .
Letŝ denote the empirical frequency distribution of site patterns on X − x, and let p(l) denote the vector of site pattern probabilities generated by T −x with branch lengths l.
The proof of this last equation is given in a separate subsection below. We apply it as follows. Notice that for g = f or g = f , we have: Moreover, since g is continuous, and the nested sequence of sets B k convergences on the vector λ x as k → ∞ we have: B k exp(−kd(ŝ||p(l))g(l)dl B k exp(−kd(ŝ||p(l))dl attend to below.

Proof of Eqn. (12)
A classic result (e.g. Wilk's theorem) ensures that the following convergence in distribution holds: where χ 2 N −1 is a chi-square distribution with N − 1 degrees of freedom (here N is the number of possible site patterns). By the continuous mapping theorem, it now follows that: where W = exp(−χ 2 N −1 ) is a continuous and non-negative random variable.
Moreover, if a sequences of branch length vectors l k lies within (l ∞ ) distance 1 k of λ −x then we also have: exp(−kd(ŝ||p(l k ))) D − → W.
Next, Pinsker's inequality (see Cover and Thomas (2006)) gives for any l ∈ Γ: d(ŝ||p(l)) ≥ 1 2 ||ŝ − p(l)|| 2 1 , where · || 1 refers to the l 1 metric. The triangle inequality for this metric then gives: Now, for any l ∈ Γ − B k , we have l − λ −x ∞ ≥ k −1/4 , and so, by Theorem 2.1(2) of length between these leaves under branch lengths l and λ −x is at least 1 2 k −1/4 . Since the site substitution model is reversible, the probability two leaves are in the same state is a monotone decreasing function of the path length between them (a positive mixture of exponential functions). This in turn implies that the event that leaves i and j are in the same state differs in probability under the branch lengths l and λ −x by an amount that is at least k −1/4 times some constant (dependent on the model, and λ −x ). In particular, p(λ −x ) − p(l) 1 ≥ ck −1/4 , for some constant c > 0. Also, with probability converging to 1 as k grows. Consequently, by combining Eqns. (20), (21) and (22), the following inequality holds for all l ∈ Γ − B k with probability converging to 1 as k grows: exp(−kd(ŝ||p(l))) exp(−dk 1/2 ) ≤ 1, for some constant d > 0. Thus, with probability converging to 1 as k grows: Γ−B k exp(−kd(ŝ||p(l))g(l)dl exp(−dk 1/2 ) ≤ lim k→∞ Γ−B k 1 · g(l)dl = 1.
On the other hand, if we let B * k ⊂ B k be the set of branch length vectors that lie within (l ∞ ) distance at most 1/k from λ −x then B * k exp(−kd(ŝ||p(l))g(l)dl Now, for a value C > 0 that is independent to k (but dependent on λ −x ) we have B * k g(l)dl ≥ Ck −E , where E is the number of edges of T . Also, from Eqn. (19) sup l∈B * k {exp(−kd(ŝ||p(l))} converges in distribution to the random variable W = exp(−χ 2 N −1 ) as k grows. Thus, from (24), the following inequality holds with probability converging to 1 as k grows: B * k exp(−kd(ŝ||p(l))g(l)dl/Ck −E ≥ W, and thus, B k exp(−kd(ŝ||p(l))g(l)dl since B * k ⊂ B k and the integrand is non-negative. Combining Eqns. (23) and (25) the following inequality holds with probability converging to 1 with increasing k: Notice that the second term on the right ( 1 W ) is a continuous random variable, but since P(W > 0) = 1 and since the first term on the right converges to zero (absolutely) as k tends to infinity, this suffices to establish Eqn. (12).