COMBINATORICS OF DISTANCE-BASED TREE INFERENCE

Several popular methods for phylogenetic inference (or hierarchical clustering) are based on a matrix of pairwise distances between taxa (or any kind of objects): the objective is to construct atreewithbranc hlengthssothatthedistancesbet w eentheleaves in that tree are as close as possible to the input distances. If we hold the structure (topology) oft h et r e efi x e d ,i ns o m er e l e v a n tc a s e s (e.g. ordinary least squares) the optimal values for the branch lengths can be expressed using simple combinatorial formulae. Here we define a general form for these formulae and show that they all have two desirable properties: first, the common tree reconstruction approaches (least squares, minimum evolution), when used in combination with these formulae, are guaranteed to infer the correct tree when given enough data (consistency ); second, the branch lengths of all the simple (nearest neighbor interchange) rearrangements of a tree can be calculated, optimally, in quadratic time in the size of the tree, thus allowing the efficient application of hill climbing heuristics. The study presented here is a continuation of that by Mihaescu and Pachter on branch length estimation (Mihaescu R, Pachter L( 2 0 0 8 )P r o cN a t lA c a dS c iU S A1 0 5 : 1 3 2 0 6 -1 3 2 1 1 ) . T h ef o c u sh ere is on the inference of the tree itself, and on providing a basis for novel algorithms to reconstruct trees from distances.

Introduction

At a s kw i t hs e v e r a lr e l e v a n ta p p l i c a t i o n si st h eu s eo fam a t r ix of distances to construct a tree whose leaves' relative positions somehow reflect the given distances. This is useful both in evolutionary biology, where the tree is intended to represent the evolution of a set of species, populations or genes, and in cluster analysis, where the tree shows the similarities in a collection of objects. In evolutionary biology, the distances are typically estimated from molecular sequences using probabilistic models of sequence evolution [1,2]. The resulting distances can be expected to be approximately additive;t h a t is, there exists a (phylogenetic)t r e ew i t hb r a n c hl e n g t h s ,s ot h a tt h el e n g t h so ft h ep a t h sb e tween its leaves (sequences) are approximately equal to the input distances. Finding this tree is the goal of several popular distance-based tree reconstruction methods.

For phylogenetic reconstruction -which this pap er concentrates on -the main advantage of distancebased methods is their speed of execution, which is orders of magnitude faster than that of other (potentially more accurate) approaches. As a consequence, distance methods are used whenever computational efficiency is of critical importance: for the reconstruction of very large trees, or -as in the case of bootstrapping -large collections of trees, or even toc o n s t r u c ti n i t i a lp h y l o g e n i e sf o rs e a r c h heuristics based on more sophisticated approaches. In fact,agen eralt ren di nb i oi n f orm at i csan dcomputational biology is the growing demand for methods that canc o p ew i t hm a s s i v ed a t a s e t so fD N A sequences. Distance-based methods are a possible answer to this demand, not only for phylogenetic inference, but also for related tasks such as sequence identification (e.g. in metagenomics) and gene orthology inference (e.g. in functional genomics). A proof of this demand is the continuing success of neighbor-joining (NJ) [3], which to-date remains the most cited algorithm in phylogenetics.

The advantage in speed of distance methods is counterbalanced by a lower accuracy than methods that take full sequence information into account [4], such asm a x i m u ml i k e l i h o o d( M L ) ,a l t h o u g hi t has recently been shown that under a certain measure of statistical efficiency, some distance methods are essentially as good as ML [5]. A limitation of distance-based methods lies in the fact that if the distances are estimated from pairwise sequence comparisonso n l y ,t h e ni tm a yb ei m p o s s i b l et oi n f e r some parameters common to the evolution of all the sequences [6]. However, it is still possible to estimate these parameters from limited sequence samples by ML, and then use distance methods for the whole sample.

If we consider the estimation of distances as a separate task,v i r t u a l l ya l ld i s t a n c em e t h o d sa r eb a s e d on two components, corresponding to the two main unknowns in ap h y l o g e n e t i ct r e e : b r a n c hl e n g t h s and topology. First, (i) we must define a method to assign lengths to the branches of any tree of fixed topology, so that the distances between leaves are as close asp o s s i b l et ot h ei n p u td i s t a n c e s . S e c o n d , (ii) we must choose a criterion to discriminate among the trees with different topologies obtained with the step above. Distance-based algorithms then look for the tree that optimizes this criterion, typically using heuristics such as successive agglomeration (e.g. NJ [3]) and hill climbing (e.g. FastME [7]).

For comp onent (i), a weighted least squares (WLS) approach is usually adopted: the lengths of the branches in a tree T are set to the values that minimize (1) i,j

w ij (δ ij -d T ij ) 2 ,
where the δ ij are the distance estimates, the d T ij are the distances between the leaves of T ,d et erm i n ed by the lengths assigned to its branches, and the weights w ij > 0 are intended to account for the variances of the δ ij :t h eh i g h e rt h ev a r i a n c e ,t h el o w e rt h ew e i g h ta n dt h ei m p o r t ance given to the corresponding residual δ ij -d T ij .I d e a l l y ,w ij should be proportional to Var[δ ij ] -1 (see Relationship with WLS and the M&P formulae below), but in practice setting the weights is a delicate art,becau se the variances are not known; for example, one trap to avoid is to assume zero variance (and therefore an infinite weight) for the distance between two identical sequences [8]. An even more ideal approach would be to also consider the covariances between distances for different pairs of taxa, which leads a generalized least squares (GLS) optimization criterion [START_REF] Bulmer | Use of the method of generalized least squares in reconstructing phylogenies from sequence data[END_REF]. The optimal branch lengths with respect to (1), and even GLS, can be expressed succinctly in matrix form.H o w e v e r ,d e s p i t es o m ep r o g r e s s[ 1 0 ] , the matrix calculations involved are computationally expensive and remain a limiting factor for the efficiency of the algorithms that use them (such as those implemented in PAUP* [START_REF] Swofford | PAUP* -phylogenetic analysis using parsimony (*and other methods)[END_REF]).

As for component (ii), distance methods fall into two broad categories: (pure) least squares (LS) methods [START_REF] Cavalli-Sforza | Phylogenetic analysis: models and estimation procedures[END_REF][START_REF] Fitch | Construction of phylogenetic trees[END_REF] use again a least squares criterion such as (1) to score trees; on the other hand, minimum evolution (ME) methods [START_REF] Kidd | Phylogenetic analysis: concepts and methods[END_REF][START_REF] Saitou | Relative efficiencies of the Fitch-Margoliash, maximum-parsimony, maximumlikelihood, minimum-evolution, and neighbor-joining methods of phylogenetic tree construction in obtaining the correct tree[END_REF] aim to find the tree with minimum total length (which can be defined in a number of different ways [START_REF] Kidd | Phylogenetic analysis: concepts and methods[END_REF][START_REF] Saitou | Relative efficiencies of the Fitch-Margoliash, maximum-parsimony, maximumlikelihood, minimum-evolution, and neighbor-joining methods of phylogenetic tree construction in obtaining the correct tree[END_REF]8], see Statistical consistency for details), among those whose branch lengths are fitted with component (i). The intuition underlying ME is the same as that of maximum parsimony for character-based tree reconstruction: simpler (i.e. shorter) explanations are preferable to more complicated ones.

An important realization has been that in some relevant casest h eb r a n c hl e n g t h st h a tm i n i m i z e( 1 ) can be expressed using simple "combinatorial" formulae, which allow to avoid slow matrix calculations. The best-known cases are that with constant weights w ij (ordinary least squares,O L S )[ 1 6 ,1 7 ]a n d that with weights proportional to 2 -tij ,w h e r et ij is the number of branches in the path between i and j in T (the balanced case) [START_REF] Desper | Theoretical foundation oft h eb a l a n c e dm i n i m u me v o l u t i o nm e t h odo fp h y l o g e n e t i c inference and its relationship to weighted least-squares tree fitting[END_REF]. These formulae allow to efficiently calculate branch lengths and to efficiently update the tree length while performing a local search for the optimal tree with respect to ME. For example, the balanced branch length formulae [START_REF] Pauplin | Direct calculation of a tree length using a distance matrix[END_REF] can be used to calculate in O(n 2 ) time, for any tree with n leaves, not only all its branch lengths, but also the total lengths of all of its NNI (nearest neighbor interchange) [ 7 ]a n dS P R( subtree pruning and regrafting)r e a r r a n g e m e n t s[ 2 0 , 2 1 ] .

Ak e yw o r ko ns u c hc o m b i n a t o r i a lf o r m u l a ef o rl e a s ts q u a r e sb ranch lengths has appeared recently [START_REF] Mihaescu | Combinatorics of least squares trees[END_REF]. The authors show that all the known formulae are particular cases of a more general framework: whenever the weights w ij (or, equivalently, the assumed variances w -1 ij )h a v eap a r t i c u l a r" m u l t i p l i c ative" form (see Relationship with WLS and the M&P formulae), then the optimal branch lengths with respect to (1) can be calculated using simple formulae -such as those for OLS or the balanced case -w h i c hh e r ew er e f e rt oa st h eM&P formulae (from the authors Mihaescu and Pachter or the word "multiplicative"). The multiplicative model is biologically and mathematically meaningful, as it can be shown that the variances of the distance estimates are approximately multiplicative for large distances [START_REF] Nei | Variances of the average numbers ofnucleotidesubstitutionswithinandbetweenpopulations[END_REF][START_REF] Bulmer | Use of the method of generalized least squares in reconstructing phylogenies from sequence data[END_REF][START_REF] Mihaescu | Combinatorics of least squares trees[END_REF].

In the following, we use the seminal work by Mihaescu and Pachter [START_REF] Mihaescu | Combinatorics of least squares trees[END_REF] as a starting point. Whereas these authors focused on the problem of branch length estimation, here we switch the focus to tree reconstruction itself -namely the statistical and algorithmic consequences of the use of combinatorial branch length formulae on tree reconstruction. Our results can be summarized as follows:

(1) We define a class of formulae for fitting branch lengths thatgeneralizestheM&Pform ulaeand consequently also all known combinatorial formulae. (2) We prove the statistical consistency of the main distance-based tree reconstruction principles (LS and ME), when combined with our formulae. In other words the optimal tree with respect to any of these principles converges to the correct tree as thei n p u td a t ab e c o m em o r ea n d more abundant and the estimated distances converge to their correct values. Particularly in the case of ME, where it is problematic, this issue has received much attention (e.g., [START_REF] Rzhetsky | Theoretical foundation of the minimum-evolution method of phylogenetic inference[END_REF][START_REF] Atteson | The performance of neighbor-joiningm e t h o d so fp h y l o g e n e t i cr e c o n s t r u c t i o n[END_REF][START_REF] Gascuel | Strengths and limitations of the minimum evolution principle[END_REF][START_REF] Denis | On the consistency of the minimum evolution principle of phylogenetic inference[END_REF][START_REF] Willson | Consistent formulas for estimating the total lengths of trees[END_REF][START_REF] Bordewich | Consistency of topological moves based on the balanced minimum evolution principle of phylogenetic inference[END_REF]). This addresses the question by Mihaescu and Pachter ([22], p.13211) of "what classes of semimultiplicative [a minor generalization of multiplicative] variance matrices result in consistent tree estimates", by showing that all multiplicative variance matrices have this property.

(3) We investigate the computational efficiency of local search heuristics in combination with our class of formulae. In particular, we describe an algorithm that calculates the branch lengths determined by the adopted formulae not only for a fixed tree T ,b u tal sof oral lt reesob t ai n ed by performing one NNI on T .T h ee n t i r e c a l c u l a t i o n o p t i m a l l y r e q u i r e s O(n 2 ) time. This algorithm can be used as the basic component for local searches, and can be combined with any classic tree reconstruction principle.

Preliminaries: Branch length formulae

We employ the standard terminology used in the phylogeneticsl i t e r a t u r e[ 4 , 2 8 ]( p h y l o g e n e t i ct r e e , topology, branch lengths, internal and external branches etc.). For simplicity, we identify the leaves of a phylogenetic tree with a set of taxa {1, 2,...,n} and we choose to consider only binary trees. We say that two subsets of taxa A and B in a tree are separated by a branch e if any path between an element of A and an element of B passes through e. A and B are k-separated when they are separated by exactly k distinct branches. A proper subset of taxa A {1, 2,...,n} is a clade if A and {1, 2,...,n}\A are separated by some branch e;i nf a c t ,e is unique and is called the root branch of A;t h ee n d po i n to fe to the side of A is called the root node of A.Ab r a n c hbelongs to clade A if it lies in the path between two elements of A.

We also adopt the following standard conventions for distance-based methods: δ denotes the n × n input distance matrix and δ ij its element expressing the distance between taxa i and j (in the following, indices i and j are always assumed to be elements of the set of taxa {1, 2,...,n}). The distances do not necessarily form a metric, as only δ ij = δ ji and δ ii =0are assumed. Given a tree T with branch lengths, d T denotes the distance matrix where d T ij coincides with the length of the path between i and j in T .W h e nδ = d T for some T ,w es a yt h a tδ is additive (with respect to T ) [ 2 9 ] .

In the rest of this section, we introduce a new class of formulae that express the branch lengths of a generic topology T over {1, 2,...,n} as a function of δ.T h i sc l a s si sp a r a m e t e r i z e db ys o m eq u a n t i t i e s that we present using a probabilistic interpretation (see SIA p p e n d i x1f o rm o r ed e t a i l s ) . L e tT be a binary tree topology. Assume that the rules for a random walk on T are defined in the following way: if we enter an internal node from a branch e,w ec a nt h e ne x i tt h i sn o d ef r o mi t st w oo t h e ra d j a c e n t branches, f and g,w i t hp r o b a b i l i t i e sγ ef and γ eg =1-γ ef ,r e s p e c t i v e l y .W er e q u i r e0 < γ ef , γ eg < 1 (note the strict inequalities). These parameters define a (non-zero) probability of reaching any branch of T from any other branch of T .

! " # $ " # # $ $ % & ' ( ) * "+ $ # # $ ( & % *
This also defines a probability distribution over the leaves of any clade A:f o ra n yi ∈ A,l e tp i|A = γ e0e1 • γ e1e2 • ...• γ e k-1 e k ,w h e r ee 0 is the root branch of A and e 1 ,e 2 ,...,e k are the branches on the path between the root of A and i.( S e eS IA p p e n d i x1f o r a fi g u r e i l l u s t r a t i n g t h i s . )C l e a r l y ,the probabilities {p i|A | i ∈ A} form a distribution over A.W ec a nt h e nd e fi n et h ea v e r a g ed i s t a n c eδ AB between any two clades A and B as the expected distance between two taxa chosen at random from A and B according to the distributions defined above:

δ AB = i∈A j∈B p i|A p j|B δ ij .
Note that the δ AB so defined depend on the γ ef parameters, as well as on the underlying topology T , but for simplicity we do not indicate this in the chosen formalism. Also note that δ AB = δ BA .F o r simplicity, we write δ iA (or δ Ai )i n s t e a do fδ {i}A .

In addition to the γ ef probabilities defined for each pair of adjacent branches (e, f ),w ei n t r o d u c ea parameter λ XY for each unordered pair {X, Y } of 3-separated clades in T (recall the definition of kseparated clades above). We constrain these parameters so that, for every internal branch separating clades

A = A 1 ∪ A 2 and B = B 1 ∪ B 2 (see Fig. 1(b)), λ A1B1 = λ A2B2 > 0, λ A1B2 = λ A2B1 > 0 and λ A1B1 + λ A1B2 =1
,m e a n i n gt h a to n l yo n ep a r a m e t e ra m o n gλ A1B1 , λ A2B2 , λ A1B2 and λ A2B1 determines all the others. A possible interpretation for λ XY is as the probability of drawing T so that X and Y are consecutive in a clockwise ordering of the taxa (which explains why λ A1B1 = λ A2B2 , λ A1B2 = λ A2B1 and λ A1B1 + λ A1B2 =1;s e eS IA p p e n d i x1f o rd e t a i l s ) .

In summary, we have three free parameters per internal node of T (γ ef , γ fg and γ ge determine γ eg , γ fe and γ gf )a n do n ef r e ep a r a m e t e rp e ri n t e r n a lb r a n c h( λ A1B1 determines λ A2B2 , λ A1B2 and λ A2B1 ). These parameters determine a set of formulae to estimate the length le of any branch in T :

(γ T , λ T )-formulae. Let the vectors γ T =( γ ef ) and λ T =( λ XY ) be defined for binary topology T , under the constraints described above. Then, for any branch e in T :

le (δ)= 1 2 (δ iA + δ iB -δ AB ) if e is external, 1 2 [λ A1B1 (δ A1B1 + δ A2B2 )+(1-λ A1B1 )(δ A1B2 + δ A2B1 ) -δ A1A2 -δ B1B2
] if e is internal, where, if e is external, we define A, B, i as in Fig. 1(a) and, if e is internal, we define A 1 , A 2 , B 1 , B 2 as in Fig. 1

(b).

Note that because λ A1B2 = λ A2B1 =1-λ A1B1 =1-λ A2B2 ,t h ef o r m u l aa b o v ef o ri n t e r n a lb r a n c h lengths is (as desired) independent of how we assign names A 1 , A 2 to the two subclades of A = A 1 ∪A 2 , and how we assign B 1 , B 2 to the two subclades of B = B 1 ∪ B 2 .A ni n t e r p r e t a t i o no ft h e(γ T , λ T )formulae as averages of simpler formulae is given in SI Appendix 1.

These formulae are a generalization of all the combinatorialf o r m u l a ep r o p o s e di nt h ep a s tt ofi tt h e branch lengths of a tree of fixed topology. In particular, the OLS branch lengths [START_REF] Vach | Conceptual and numerical analysis of data,e d .O p i t z[END_REF][START_REF] Rzhetsky | Theoretical foundation of the minimum-evolution method of phylogenetic inference[END_REF] can be obtained by setting λ A1B1 = |A1||B2|+|A2||B1| |A||B| (same clade naming as above), and by setting γ ef = |A1| |A1|+|A2| ,f o r every pair of adjacent branches e and f in the configuration of Fig. 1(b). Note that the γ ef parameters thus defined ensure that {p i|X | i ∈ X} is uniform for any clade X (in fact the word "unweighted" is often associated to OLS). Similarly, the balanced branch lengths [START_REF] Pauplin | Direct calculation of a tree length using a distance matrix[END_REF] at the basis of the balanced minimum evolution principle [7,[START_REF] Desper | Theoretical foundation oft h eb a l a n c e dm i n i m u me v o l u t i o nm e t h odo fp h y l o g e n e t i c inference and its relationship to weighted least-squares tree fitting[END_REF][START_REF] Gascuel | Neighbor-joining revealed[END_REF] are obtained by setting all parameters to 1 2 .T h en e x ts e c t i o n shows that the (γ T , λ T )-formulae also generalize the M&P formulae by Mihaescu and Pachter [START_REF] Mihaescu | Combinatorics of least squares trees[END_REF]. It is easy to see that the (γ T , λ T )-formulae still satisfy the independence of irrelevant pairs (IIP) property introduced by those authors [START_REF] Mihaescu | Combinatorics of least squares trees[END_REF] as a basic requirement for their formulae.

Finally, we show that the (γ T , λ T )-formulae above are correct,t h a ti s ,t h e yc a l c u l a t et h ec o r r e c t values of the branch lengths of any given tree whenever the distances are additive with respect to that tree (proof in SI Appendix 1). Naturally, because the input distances are only estimates of the real evolutionary distances, they are usually only approximately additive. However, this property is an important prerequisite of any branch length formula, as it ensures the statistical consistency of the branch lengths assigned to the correct topology (see Statistical consistency below).

Theorem 1. Let T be a binary topology. For any given branch e in T ,a s s i g nl e n g t hℓ e to e,a n d let δ be additive with respect to the resulting tree. Let le (δ) be the length that is assigned to e by a (γ T , λ T )-formula. Then, le (δ)=ℓ e .

Relationship with WLS and the M&P formulae

The choice of the weights in (1) is a key factor for the accuracyo fl e a s ts q u a r e st r e ee s t i m a t e s . T h e weights w ij should be proportional to Var[δ ij ] -1 ,ast h i si m p l i est h att h eb ran c hl en gt h st h atm i n i m i ze (1) have minimum variance among all linear unbiased estimators of the branch lengths (under the assumption that the distance estimates are unbiased and uncorrelated for different pairs of taxa) [START_REF] Aitken | On least squares and linear combinations of observations[END_REF]. In this section, we consider the case where the weights (and therefore the assumed variances) are "multiplicative": given a tree topology T and a collection of weights w =( w ij ) associated to pairs of taxa in T ,w es a yt h a tt h e s ew e i g h t sa r emultiplicative with respect to T ,i fw ec a na s s i g nt oe a c h branch e of T aw e i g h tw e > 0,s ot h a t ,f o re v e r yp a i ro ft a x ai and j, w ij = e∈Pij (T ) w e , where P ij (T ) denotes the set of branches in the path between i and j in T .T h i sc o n d i t i o ng e n e r a l i z e ss e v e r a l well-known cases: that of constant weights (coinciding withO L Sa n do b t a i n e db ys e t t i n gw e to 1 for internal branches and to a constant for external ones), that of taxon-specific weights [START_REF] Denis | On the consistency of the minimum evolution principle of phylogenetic inference[END_REF] (obtained like for OLS but with w e free to vary for external branches) and also that of weights exponentially related to the number of branches separating each pair of taxa( w h i c h ,w h e nt h eb a s eo ft h ee x po n e n t is b = 1 /2,c o i n c i d ew i t ht h eb a l a n c e dw e i g h t s[ 1 9 ]a n da r eo b t a i n e db ysetting w e = b for internal branches and to a constant for external ones). [START_REF] Mihaescu | Combinatorics of least squares trees[END_REF] have shown that if the assumed weights w are multiplicative with respect to T ,thentheoptimalbranc hlengthsofT with respect to the WLS criterion (1) are given by their M&P formulae. We refer to SI Appendix 2 for a description of these formulae. The following theorem shows that the class of the M&P formulae is contained in that of the (γ T , λ T )-formulae and, conversely, it characterizes the values of γ T and λ T corresponding to M&P formulae.

Mihaescu and Pachter

Theorem 2. Let T be a binary topology. (i) Given any w multiplicative w.r.t. T , the corresponding M&P formulae are also (γ T , λ T )-formulae for some choice of γ T and λ T satisfying the properties P1 and P2 below. (ii) Given any γ T and λ T satisfying the properties P1 and P2 below, the corresponding (γ T , λ T )-formulae are also M&P formulae for some choice of w, multiplicative w.r.t. T . P1. For every internal node of T ,i fe, f and g are the three branches incident to it, then

γ ef γ fg γ ge =(1-γ ef )(1 -γ fg )(1 -γ ge ).
P2. For every pair of clades A and B separated in T by three branches a, e and b (with a being the root branch of A,a n db being the root branch of B), λ AB = γ ea + γ eb -2γ ea γ eb .

Theorem 2, proved in SI Appendix 2, not only shows that the M&P formulae are particular types of (γ T , λ T )-formulae, but it also provides an alternative set of parameters to represent the M&P formulae: instead of the branch-associated weights w e ,o n ec a nu s eas e to fγ ef parameters satisfying P1. This condition implies that any of γ ef , γ fg and γ ge can be determined from the other two and P2 implies that all the λ XY parameters are determined by the γ ef parameters. This reduces the number of free parameters needed to describe the (γ T , λ T )-formulae that are also M&P formulae to 2 per internal node, that is 2n -4.T h i si se x a c t l yo n el e s st h a nt h e2n -3 branch-associated parameters w e describing multiplicative weightings, which corresponds to the fact that multiplying all the w e for external branches by any positive constant results in equivalent weightings with respect to (1).

Theorem 2 establishes that the (γ T , λ T )-formulae have enough "expressive power" to optimize the least squares criterion (1), when the weights are multiplicative. It is therefore important to discuss this assumption. First, multiplicative weights generalizet h eb a l a n c e dw e i g h t s ,w h i c hh a v eb e e ne xperimentally demonstrated to behave well in combination with ME [START_REF] Desper | Theoretical foundation oft h eb a l a n c e dm i n i m u me v o l u t i o nm e t h odo fp h y l o g e n e t i c inference and its relationship to weighted least-squares tree fitting[END_REF][START_REF] Vinh | Shortest triplet clustering: reconstructing large phylogenies using representative sets[END_REF]. Second, in the case of distances estimated from molecular sequences, we note that for many models of sequence evolution (for instance Jukes-Cantor [33]; see Chapter 13 in [4] or Appendix B in [2] for the general technique), the variance of δ ij can be approximated by a function of the correct evolutionarydistanced ij that, for small values of d ij ,b e h a v e sa sal i n e a rf u n c t i o no fd ij ,a n d ,f o rm o d e r a t e -t o -l a r g ed ij ,a sa ne x p o n e ntial of d ij .T h i sm e a n st h a t ,f o rp a i r so ft a x as e p a r a t e db ys m a l ld ij ,t h ev a r i a n c e so ft h e i rd i s t a n c e estimates will tend to be additive, whereas for pairs of taxa separated by moderate-to-large d ij ,t h e variances will tend to be multiplicative. The additive modelf o rt h ev a r i a n c e s[ 3 4 ] ,o ri t sv a r i a n tw i t h variances proportional to d 2 ij [START_REF] Fitch | Construction of phylogenetic trees[END_REF], are used in practice with δ ij in place of d ij ,ast h el at t eri su n kn o wn . As a result, these approaches need some precautions for very small distance estimates, so as to avoid an overconfidence in these estimates (for δ ij tending to 0, also the assumed variance tends to 0, and w ij tends to infinity): for example, one possibility is to add pseudocounts to the numbers of observed differences between sequences [8] (known as "Laplace smoothing"). In this context, the multiplicative model provides a simple and robust alternative for small distances (for δ ij → 0,t h ea s s u m e dv a r i a n c e tends to a constant), and is mathematically justified for moderate-to-large distances.

The other important assumption here, common to all WLS methods, is that the δ ij are uncorrelated for different pairs of taxa, which is clearly not true for distances estimated from molecular sequences [START_REF] Bulmer | Use of the method of generalized least squares in reconstructing phylogenies from sequence data[END_REF]. As mentioned above, covariances between different distance estimates can be accounted for by adopting aGLScriterion. Ho w ev er,settingtheco v ariancesandcalculating the resulting branch lengths [START_REF] Bryant | Rapid evaluation of least-squares and minimum-evolution criteria on phylogenetic trees[END_REF] are difficult problems, which explains the lack (to the best of our knowledge) of practical implementations of GLS for phylogenetic reconstruction.

Statistical consistency

Am e t h o df o rp h y l o g e n e t i ci n f e r e n c ei ss a i dt ob e(statistically) consistent if the probability that it reconstructs the correct tree (within any given accuracy)c o n v e r g e st o1a sm o r ea n dm o r ed a t a are analyzed. For distance-based methods, the consistency of tree inference usually depends in turn on the consistency of the distance estimates, that is, the assumption that δ converges to a matrix d T containing the distances in the correct phylogenetic tree for the taxa under consideration. Even though in reality the precise consistency of distance estimates cannot be expected to hold -because the models used to obtain these estimates are only approximations of reality -the ability to infer the correct tree in such a best-case scenario is an essential property of any phylogenetic inference method: it is a prerequisite for robust inference of the correct topology with real distance estimates, subject to sampling errors and not perfectly consistent [START_REF] Atteson | The performance of neighbor-joiningm e t h o d so fp h y l o g e n e t i cr e c o n s t r u c t i o n[END_REF][START_REF] Susko | On inconsistencyo ft h en e i g h b o r -j o i n i n g ,l e a s ts q u a r e s ,a n dm i n i m u m evolution estimation when substitution processes are incorrectly modeled[END_REF][START_REF] Pardi | Robustness of phylogenetic inference based on minimum evolution[END_REF].

In this section, we state our main results on the statistical consistency of the tree reconstruction methods using the (γ T , λ T )-formulae. We leave the proofs to SI Appendix 3. We assume that, for any binary topology T over the taxa of interest {1, 2,...,n},ac o l l e c t i o no fp a r a m e t e r sγ T =( γ ef ) and λ T =(λ XY ) is defined, thus defining in turn, for any such T ,asetof(γ T , λ T )-formulae for estimating the branch lengths of T .W ec a l lt h i sabranch length estimation scheme based on (γ, λ)-formulae. (Note the absence of superscript.) We stress that, for the consistency results here, no connection between (γ T , λ T ) and (γ T ′ , λ T ′ ) for different topologies T and T ′ needs to be assumed; in other words, completely unrelated formulae can be used for any pairo ft o p o l o g i e s . Now combine a branch length estimation scheme with an optimization principle, such as LS or ME, that allows us to choose among all the topologically-distinct fitted trees over {1, 2,...,n}.W eh a v e already described LS (but also see SI Appendix 3). As for ME, three variants of this principle have been proposed, essentially differing for how tree length is defined in the presence of negative branch lengths (which are allowed by many branch length estimation schemes, including those based on (γ, λ)-formulae). We call them ME -1 [START_REF] Kidd | Phylogenetic analysis: concepts and methods[END_REF], ME +1 [START_REF] Saitou | Relative efficiencies of the Fitch-Margoliash, maximum-parsimony, maximumlikelihood, minimum-evolution, and neighbor-joining methods of phylogenetic tree construction in obtaining the correct tree[END_REF][START_REF] Rzhetsky | A simple method for estimating and testing minimum-evolution trees[END_REF], and ME 0 [8]. Assuming that a tree has been assigned the branch lengths le , ME i defines its length as e:ℓe>0 le + e:ℓe<0 i • le .

The three versions of ME then differ in how they deal with negative branch lengths when calculating tree length: ME +1 adds together all branch lengths irrespective of their sign,w h e r e a sME 0 ignores negative branch lengths and ME -1 takes their absolute value. Gascuel et al. [START_REF] Gascuel | Strengths and limitations of the minimum evolution principle[END_REF] previously named ME +1 , ME 0 and ME -1 ,"all-BL","positive-BL"and"absolute-BL",respectively. The following theorem shows that for these three versions of ME, as well as for LS, tree inference is consistent when (γ T , λ T )formulae are used.

Theorem 3. Assume that the input distances δ are consistent estimates of the correct evolutionary distances d T *

,w h e r eT * is a binary tree with positive branch lengths. Adopt a branch length estimation scheme based on (γ, λ)-formulae. Then, the optimal trees with respect to LS, ME +1 , ME 0 and ME -1 are statistically consistent estimates of T * .

Whereas the consistency of LS is a simple consequence of the correctness of the (γ T , λ T )-formulae, and is included here for sake of completeness, the result for ME issomewhatsurprising,giv enthatMEhas been proven to be inconsistent when combined with WLS branch lengths (for some particular values of the weights w ij ) [START_REF] Gascuel | Strengths and limitations of the minimum evolution principle[END_REF]. F urt herm ore,Theorem3generalizesallpreviouslyknown cases of consistency for the ME principle [START_REF] Rzhetsky | Theoretical foundation of the minimum-evolution method of phylogenetic inference[END_REF][START_REF] Denis | On the consistency of the minimum evolution principle of phylogenetic inference[END_REF][START_REF] Desper | Theoretical foundation oft h eb a l a n c e dm i n i m u me v o l u t i o nm e t h odo fp h y l o g e n e t i c inference and its relationship to weighted least-squares tree fitting[END_REF]. In particular, it demonstrates the statistical consistency of tree reconstruction when using the formulae by Mihaescu and Pachter, thus answering their fundamental question mentioned in the Introduction.

Computational efficiency

While the statistical consistency results above provide a theoretical basis for the use of (γ T , λ T )formulae, we now consider a more practical advantage of thesef o r m u l a e : t h ef a c tt h a tt h e yc a nb e efficiently combined with hill climbing heuristics, a pervasive and successful tool for tree reconstruction. Hill climbing consists of repeatedly applying small changesthatimpro v ethescoreofacandidatetree, until no such change is possible anymore. The behavior of hillclimbingisessentiallydeterminedbythe changes allowed at each step, or in other words by a notion of neighborhood defined over tree space. Here, we consider the simplest such changes, known as nearest neighbor interchanges (NNIs), which consist of swapping the positions of two 3-separated subtrees in a topology: for example, the topology in Fig. 1(c) can be obtained from that in Fig. 1(b) by swapping clades A 1 and B 1 .W h e nt o p o l o g y T ′ can be obtained from topology T in this way, we say that T and T ′ are NNI-neighbors.A nN N I transforming T into T ′ is around e,i fe is the middle branch among the three branches separating the subtrees being swapped in T . While simple, NNIs can be used to efficiently implement more complex changes (such as SPRs) that can be obtained via a series of NNIs[ 2 8 ,2 1 ] .

Clearly, the computational efficiency of a hill climbing heuristic depends crucially on the ability to efficiently evaluate some/all neighbors of any candidate topology. For all distance-based optimization principles, the evaluation is essentially done on the basis of some function of the assigned branch lengths. It is then important to calculate efficiently the branch lengths of the neighbors that are considered at each iteration. Here, we show that if (γ T , λ T )-formulae are used for computing branch lengths, and a natural relation between the γ T parameters for NNI-neighbors is assumed, then the O(n 2 ) branch lengths of all the NNI-neighbors of a candidate topology can be calculated in O(n 2 ) time. This is optimal, because these O(n 2 ) branch lengths depend on all the O(n 2 ) input distances.

In order to express the required relation between the γ T parameters for NNI-neighbors, we assume that when performing an NNI around a branch e,a l lo t h e rb r a n c h e sk e e pt h e i rn a m e s . ( F o re x a m p l e , see branches f , g, h and l in Fig. 1(b) and (c).) Then, when T ′ is obtained from T with an NNI around branch e,w es a yt h a tp a r a m e t e rs e t sγ T =( γ e1e2 ) and γ T ′ =( γ ′ e1e2 ),d e fi n e df o rT and T ′ , respectively, are almost identical,i fγ e1e2 = γ ′ e1e2 ,f orev eryp ai rofad j acen tb ran c h es(e 1 ,e 2 ) in T such that their common endpoint is not also an endpoint of e (in which case e 1 and e 2 are also adjacent in T ′ ). The intuitive idea is that γ T and γ T ′ may only differ locally around the location of the NNI. This requirement is a prerequisite for the efficient evaluation of T ′ from that of T .N o t et h ed i fference here with the approach in the previous section, where we assumed no relationship between parameter sets for different topologies. Our result can now be stated as follows:

Theorem 4. Let T 0 be a binary topology over taxa {1, 2,...,n} and T 1 ,T 2 ,...,T 2(n-3) all its NNIneighbors. For all i ∈ {0, 1,...,2(n -3)},a s s u m et h a tt h eb r a n c hl e n g t h so fT i are defined by the (γ Ti , λ Ti )-formulae, with the constraint that γ Ti and γ T0 are almost identical. Then, (i) the branch lengths of T 0 can be calculated in O(n 2 ) time;

(ii) the branch lengths of all the NNI-neighbors of T 0 can be calculated in O(n 2 ) time.

We leave the pro of of this result to SI App endix 4. While p oint (i) merely generalizes further a property already known for all M&P formulae [START_REF] Mihaescu | Combinatorics of least squares trees[END_REF], the result in (ii) is novel. It is related to, and somehow explains the existence of a number of efficient hill climbing algorithms for distance-based tree reconstruction. In particular, it predicts the efficiency of hill climbing for balanced minimum evolution (BME), which assumes γ T parameters always equal to 1 /2 and therefore clearly having the property of being almost identical for NNI-neighbors. The existing hill climbing algorithm for BME [7] directly updates the total tree length, rather than the lengths of each branch, but the worst-case time complexity for each iteration is still O(n 2 ) and results in one of the most accurate and fast distance-based methods [START_REF] Desper | Theoretical foundation oft h eb a l a n c e dm i n i m u me v o l u t i o nm e t h odo fp h y l o g e n e t i c inference and its relationship to weighted least-squares tree fitting[END_REF][START_REF] Vinh | Shortest triplet clustering: reconstructing large phylogenies using representative sets[END_REF]. Theorem 4 also predicts thee fficiency of hill climbing for OLS: the γ ef parameters for OLS depend in fact on the sizes of the three clades to the sides of e, f and g (where the latter is the branch adjacent to both e and f ), and these do not change when performing an NNI around a branch other than e, f and g,w h i c hi m p l i e st h ea l m o s ti d e n t i t yo ft h eγ ef parameters for NNI-neighbors.

Note that Theorem 4 has very wide applicability, not only because of the generality of the formulae it assumes, but also because it makes no assumption on the optimization criterion used to score trees (apart from its dependence on the branch lengths). This is unlike the hill-climbing algorithms we mentioned above, which were only applicable to the classic version of ME (the one we call ME +1 ), where all branch lengths are added together, irrespective oft h e i rs i g n .

Discussion

We presented here a framework unifying some of the most succe ssful approaches for distance-based tree reconstruction: for example, ordinary least squares methods for clustering [START_REF] De Soete | A least squares algorithm for fitting additive trees to proximity data[END_REF] and balanced minimum evolution (BME, the optimization principle behind Neighbor-Joining [START_REF] Gascuel | Neighbor-joining revealed[END_REF]) for phylogenetic inference. We have shown that all the metho ds that fit into this general framework have highly desirable statistical properties (the consistency of the tree estimates) and algorithmic properties (efficiency of hill climbing heuristics).

Our study opens the way for improvements of existing methods and the development of new ones. Novel combinations of branch length formulae and tree optimization principles can be envisaged. For example, our results enable the efficient implementation of hill climbing for the versions of ME discouraging negative branch lengths (or at least not favoring them: see ME -1 and ME 0 above), in combination with any of the classic branch length estimationschemes(e.g.OLSorthatusedinBME). Alternatively, our framework enables to explore novel, biologically-motivated ways of estimating branch lengths, for example assuming multiplicative variance models based on the current tree estimate.

We conclude by noting that although the class of branch lengthf o r m u l a ew ec o n s i d e rh e r ei si n s p i r e d by previous work on multiplicative variance models [START_REF] Mihaescu | Combinatorics of least squares trees[END_REF], nothing excludes that it may be applicable to least squares criteria other than WLS with multiplicativew e i g h t s . I nf a c t ,i ti se a s yt oc o n s t r u c t covariance models with nonzero covariances that result in GLS branch length estimators coinciding with (γ T , λ T )-formulae. Future research should aim to elucidate the full potential of our class of formulae.

SI Appendix 1

Here, we present an interpretation of the (γ T , λ T )-formulae as averages over a large number of simple branch length formulae, which allows us to prove Theorem 1: because these simple formulae are correct, that is, they correctly provide the branch lengths of a tree whenever the input distances are additive with respect to that tree, it follows that also the (γ T , λ T )-formulae are correct. The detailed arguments follow below.

Let le (δ) denote the length that is assigned to branch e by an adopted length estimation method. Let T * be a tree where e has length ℓ e .T h e a d o p t e d m e t h o d i s correct if, for any such tree T * , le (d T * )=ℓ e . Suppose e is an external branch, and define A, B, i as in Fig. 1(a) in the main text. Choose a taxon a from A and a taxon b from B.T h e nc a l c u l a t et h el e n g t ho fe with:

lab e (δ)= 1 2 (δ ia + δ ib -δ ab ).
If instead e is an internal branch, let A 1 , A 2 , B 1 and B 2 be the four clades surrounding it, as in Fig.

1(b) in the main text. Choose taxa a 1 , a 2 , b 1 , b 2 from A 1 , A 2 , B 1 , B 2 ,
r e s p e c t i v e l y . I ti sc l e a rt h a t any drawing of the tree on the plane either places A 1 to the side of B 1 ,a n dt h e r e f o r eA 2 to the side of B 2 (as in Fig. S1(ii), top) or alternatively A 1 to the side of B 2 and A 2 to the side of B 1 (Fig. S1(ii), bottom). We associate the former drawing with the following formula for the length of e:

la1b1a2b2 e (δ)= 1 2 (δ a1b1 + δ a2b2 -δ a1a2 -δ b1b2 ).
The alternative drawing is associated to the formula la1b2a2b1

p i|A = γ e0e1 • γ e1e2 • ...• γ e k-1
e k is the probability of ending up in i when entering clade A from its root and following the random walk rules described in the main text; (ii) λ A1B1 = λ A2B2 can be seen as the probability of drawing the tree in the top configuration, whilst λ A1B2 = λ A2B1 =1λ A1B1 can be seen as the probability of drawing the tree in the bottom configuration.

The γ T and λ T parameters introduced in the main text can be interpreted as controlling a probability distribution over all possible such formulae for calculating the length of a given branch e in T .F i r s t , as illustrated in Fig. S1(i), the γ T parameters determine the probability of choosing a given taxon out of any given clade (a from A, b from B in the case of external branches, a 1 from A 1 , a 2 from A 2 , b 1 from B 1 , b 2 from B 2 in the case of internal branches). Second, as illustrated in Fig. S1(ii), the λ T parameters determine the probability of choosing either of the two possible drawings for the clades around an internal branch e,a n dt h e r e f o r ee i t h e ro f la1b1a2b2 e (δ) or la1b2a2b1 e (δ) for the length of e: λ A1B1 = λ A2B2 is the probability of drawing A 1 to the side of B 1 and A 2 to the side of B 2 ,w h e r e a s its complement λ A1B2 = λ A2B1 =1-λ A1B1 is the probability of drawing A 1 to the side of B 2 and A 2 to the side of B 1 .

Given this probability distribution, let us take the resulting expected value of the length assigned to e. In the (harder) case of an internal branch, this is given by

a 1 ∈A 1 , a 2 ∈A 2 , b 1 ∈B 1 b 2 ∈B 2 p a1|A1 p a2|A2 p b1|B1 p b2|B2 λ A ! B1 la1b1a2b2 e (δ)+λ A ! B2 la1b2a2b1 e (δ) = 1 2 a 1 ∈A 1 , a 2 ∈A 2 , b 1 ∈B 1 b 2 ∈B 2 p a1|A1 p a2|A2 p b1|B1 p b2|B2 [λ A ! B1 (δ a1b1 + δ a2b2 )+(1-λ A1B ! )(δ a1b2 + δ a2b1 ) -δ a1a2 -δ b1b2 ] = 1 2 [λ A1B1 (δ A1B1 + δ A2B2 )+(1-λ A1B1 )(δ A1B2 + δ A2B1 ) -δ A1A2 -δ B1B2 ] .
Thus what we obtain (also in the easier case of an external branch; not shown) are precisely the (γ T , λ T )-formulae. In other words, these formulae can be seen as providing the expected length of a branch when this is assigned following the random procedure described above. Given this observation, the correctness of the (γ T , λ T )-formulae follows trivially from the correctness of the basef o r m u l a e lab e (δ) and la1b1a2b2 We note that the approach of expressing a length estimator as the combination of several simple formulae has already been considered by Willson [1]. His basef o r m u l a e ,h o w e v e r ,e x p r e s st h el e n g t h of a path in the tree (as a function of the distances between three taxa) rather than a single branch (which we express as a function of the distances between threeo rf o u rt a x a ,f o re x t e r i o ra n di n t e r i o r branches, respectively). Moreover, the combination of his base formulae provide an estimate of the total length of the tree (in the ME +1 sense).

Here, we prove the relationship between the M&P formulae [2] and our (γ T , λ T )-formulae, as stated in Theorem 2. We start by formally defining the M&P formulae (A2.1); then we introduce a few additional formalisms and a useful observation (A2.2) and then prove separately the two parts of Theorem 2 ((i) in A2.3 and (ii) in A2.4). This requires to show how to derive the parameters of each class of formulae from the parameters of the other class (i.e. (γ T , λ T ) from w,a n dv i c ev e r s a ) .

A2.1 The M&P formulae. We assume that the weights w =(w ij ) are multiplicative w.r.t. a binary topology T .T h e n ,f o ra n yt w oc l a d e sA and B of T ,d e fi n e

Z AB = i∈A j∈B w ij and δ w AB = 1 Z AB i∈A j∈B w ij δ ij .
Mihaescu and Pachter [2] have shown that the optimal branch lengths of T with respect to the WLS criterion (1) are then given by the following formula, applicable to any branch e in T :

le (δ)=            1 2 (δ w iA + δ w iB -δ w AB ) if e is external, 1 2 ZA 1 B 2 +ZA 2 B 1 ZAB (δ w A1B1 + δ w A2B2 ) + ZA 1 B 1 +ZA 2 B 2 ZAB (δ w A1B2 + δ w A2B1 ) -δ w A1A2 -δ w B1B2 if e is internal,
where, if e is external, we define A, B, i as in Fig. 1(a) in the main text and, if e is internal, we define where we recall that the w e are the branch-associated weights that compose the pairwisew e i g h t sw ij , and P xy (T ) is the set of branches on the path between x and y in T .T h e nd e fi n et h emultiplicative weight of X (a clade with root x)a s :

A 1 , A 2 , A = A 1 ∪ A 2 , B 1 , B 2 , B = B 1 ∪ B 2 as in
(2)

Z X = i∈X w ix .
We assume w xx =1for any node x,w h i c hi m p l i e sZ {i} =1 ,f o ra n yo n e -t a x o nc l a d e{i}.I ti st h e n easy to check that, if X and Y are any two disjoint clades in T ,w i t hr oo t sx and y respectively, then

Z XY = w xy Z X Z Y . (3)
A2. 3 The M&P formulae are also (γ T , λ T )-formulae.

Proof. (i) Equation (3) allows us to express δ w

XY in a very similar form to that of δ XY :

δ w XY = i∈X j∈Y w ij Z XY δ ij = i∈X j∈Y w ix Z X w jy Z Y δ ij ,
where x and y are the root nodes of X and Y ,r e s p e c t i v e l y . I no r d e rt oh a v eδ XY = δ w XY ,i ti st h e n sufficient to prove that, for any clade X with root x,a n da n yt a x o ni ∈ X,

p i|X = w ix Z X .
Let the path from x to i traverse branches e 1 ,e 2 ,...,e k (in this order), and let X j be the subclade of X having e j as root branch (0 ≤ j ≤ k,w i t hX 0 = X and e 0 being the root branch of X); then, X k = {i} and

p i|X = γ e0e1 • γ e1e2 • ...• γ e k-1 e k = w e1 Z X1 Z X • w e2 Z X2 Z X1 • ...• w e k Z X k Z X k-1 = w e1 • w e2 • ...• w e k Z X Z {i} = w ix Z X .
(ii) Using again (3),

Z A1B2 + Z A2B1 Z AB = w f w l Z A1 Z B2 + w g w h Z A2 Z B1 Z A Z B = γ ef γ el + γ eg γ eh = γ ef + γ eh -2γ ef γ eh .
Points (i) and (ii) are thus both verified.

Setting γ ef as in (4) has a simple intuitive meaning: if we call g the root branch of A 2 (as in Fig. 1

(b)), then Z A = w f Z A1 + w g Z A2 .
T h eγ ef above can then be seen as the relative multiplicative weight of the subtree corresponding to clade A 1 in the subtree corresponding to clade A.T h er a n d o mw a l kd e fi n e d by these parameters is then "attracted" by the heavier subtrees, in a way that is directly proportional to the weights of the subtrees. It is easy to check that this implies λ A1B1 = λ A2B2 > 0, λ A1B2 = λ A2B1 > 0 and λ A1B1 + λ A1B2 =1 , and therefore the definition of the λ T parameters is also admissible.

It is now easy to verify that the resulting (γ T , λ T )-formulae coincide with the M&P formulae corresponding to w:f o rt h ee x t e r n a lb r a n c h e st h i si sa ni m m e d i a t ec o n s e q u e n c eof δ XY = δ w XY ,w h i l e for the internal branches we also use the above definition of λ A1B1 and the fact that 1 -λ A1B1 = (Z A1B1 + Z A2B2 )/Z AB .F i n a l l y ,t h eγ T and λ T defined above satisfy properties P1 and P2: the first can be verified by using (4) in P1 and the second is a direct consequence of Lemma 1, part (ii).

A2.4 Characterization of the (γ, λ)-formulae that are also M&P formulae.

Proof of Theorem 2, part (ii).

The proof has the following structure: as an intermediate step, we introduce -for every three clades A 1 , A 2 and B whose respective root branches f , g and e are incident to the same internal node (as in Fig. 1(b)) -three values ϕ A1 , ϕ A2 and ϕ B such that

ϕ A1 + ϕ A2 + ϕ B =1, 0 < ϕ A ! , ϕ A2 , ϕ B < 1 and (5) γ ef = ϕ A1 ϕ A1 + ϕ A2 , γ fg = ϕ A2 ϕ A2 + ϕ B , γ ge = ϕ B ϕ B + ϕ A1 .
The existence of such values is guaranteed by property P1. In intuitive terms, we are requiring that each clade has somewhat a "preference value", such that the probabilities γ ef and γ eg are proportional to the preference values of the clades that f and g lead to. On the basis of these values, we then define asetofbranc h-associatedw eigh tsw e specifying a multiplicative model such that the preference values can be obtained as ( 6)

ϕ B = w e Z B w e Z B + w f Z A1 + w g Z A2 .
(With A 1 ,A 2 ,B,f,g,e as above.) It is easy to see (shown below) that this implies that condition (4) of Lemma 1, and therefore its conclusions, hold. This, together with the fact that the λ T parameters satisfy P2, implies that the M&P formulae for w coincide with the given (γ T , λ T )-formulae.

Let us now look in detail at each step. First of all, property P1i m p l i e st h a tt h ed e t e r m i n a n to ft h e coefficient matrix for the system of linear equations in ϕ A1 , ϕ A2 and ϕ B corresponding to ( 5) is equal to 0; this system is then solved by the following subspace of solutions:

(ϕ B , ϕ A1 , ϕ A2 ) ∈ x, 1 -γ ge γ ge x, 1 -γ ef γ ef 1 -γ ge γ ge x | x ∈ R .
If furthermore we impose ϕ A1 + ϕ A2 + ϕ B =1,i ti se a s yt os e et h a tau n i q u es o l u t i o ni sd e t e r m i n e d , such that ϕ A1 , ϕ A2 , ϕ B > 0 (and therefore also < 1).

Given the ϕ X parameters, we now show how to define the branch-associated weights w e .F o r a n y internal branch e separating clades A and A,l e t ( 7)

w e = ϕ A 1 -ϕ A ϕ A 1 -ϕ A .
As for the external branches, which for simplicity we call with the same names as the taxa they are incident to (e.g. branch i being the one incident to taxon i), we assign their weight in the following way: choose arbitrarily w 1 > 0 and then, for any other external branch i ∈ {2, 3,...,n} define ( 8)

w i = w 1 γ * 1i γ * i1 1 -ϕ 1 ϕ 1 ϕ i 1 -ϕ i ,
where γ * ef = γ ee1 • γ e1e2 • ...• γ e k f ,f o ra n yp a i ro fb r a n c h e se, f linked by a path composed of the ordered sequence of branches (e 1 ,e 2 ,...,e k ),an df orsi m p l i ci t yw ewri t eϕ 1 and ϕ i instead of ϕ {1} and ϕ {i} . The weights thus defined determine a multiplicative weighting w =( w ij ) that satisfies (6). In order to show this, we first show that, for every clade A (whose root branch we call e)suc ht hat|A| <n-1 (i.e. the endpoint of e on the other side of A is not a leaf), ( 9)

w e Z A = w 1 γ * 1e γ * e1 1 -ϕ 1 ϕ 1 ϕ A 1 -ϕ A .
We prove this by induction on the size of A.I fA consists of a single taxon, then either this is taxon 1, in which case both sides of the equation reduce to w 1 ,o rt h i si sa n o t h e rt a x o ni,i nw h i c hc a s e( 9 ) coincides with (8). In both cases (9) trivially holds. If |A| > 1,l e tA 1 , A 2 , B, f and g be as in Fig. 1(b). Then, by inductive hypothesis, (9) holds for w f Z A1 and w g Z A2 and we have:

Z A = w f Z A1 + w g Z A2 = w 1 1 -ϕ 1 ϕ 1 γ * 1f γ * f 1 ϕ A1 1 -ϕ A1 + γ * 1g γ * g1 ϕ A2 1 -ϕ A2 Now note that (10) γ * 1f γ * f 1 = γ * 1e γ * e1 ϕ A1 (1 -ϕ A1 ) ϕ B (1 -ϕ B ) and γ * 1g γ * g1 = γ * 1e γ * e1 ϕ A2 (1 -ϕ A2 ) ϕ B (1 -ϕ B ) ,
which can be verified by noting that, depending on the positionoftaxon1(inA )(γ fg γ ef )/(γ fe γ gf ).T h ee q u a t i o n si n( 1 0 )c a nt h e nb eo b t a i n e df r o m these expressions by making the substitutions γ xy = ϕ Y /(1-ϕ X ) (equivalent to (5)) for x, y ∈ {e, f, g}, where X, Y are the clades separated by and having x, y as root branches, respectively. If now we use [START_REF] Bryant | Rapid evaluation of least-squares and minimum-evolution criteria on phylogenetic trees[END_REF] in the expression above for Z A ,w eo b t a i na f t e ro b v i o u ss i m p l i fi c a t i o n s (11)

Z A = w 1 γ * 1e γ * e1 1 -ϕ 1 ϕ 1 1 -ϕ B ϕ B .
If now we use (7) and [START_REF] Swofford | PAUP* -phylogenetic analysis using parsimony (*and other methods)[END_REF] to express w e and Z A in w e Z e ,w h a tw eo b t a i ni sp r e c i s e l y( 9 ) ,w h i c h therefore is proven.

We are now ready to prove (6). Note that (11) holds for any 'composite' clade A (i.e. one that that can be decomposed into two other clades A 1 and A 2 ). Then,

w e Z B w e Z B + w f Z A1 + w g Z A2 = w e Z B w e Z B + Z A = ϕB 1-ϕB ϕB 1-ϕB + 1-ϕB ϕB = ϕ B ,
where for the second equality we have used both [START_REF] Swofford | PAUP* -phylogenetic analysis using parsimony (*and other methods)[END_REF] and [START_REF] Bulmer | Use of the method of generalized least squares in reconstructing phylogenies from sequence data[END_REF].

But this implies that, for every composite clade

A = A 1 ∪ A 2 in the configuration of Fig. 1(b), (12)
w f Z A1 Z A = w f Z A1 w f Z A1 + w g Z A2 = ϕ A1 ϕ A1 + ϕ A2 = γ ef .
Condition (4) of Lemma 1 is therefore verified. But this ensures that δ XY = δ w XY ,f o ra n yt w od i s j o i n t clades X and Y (Lemma 1, part (i)), while the fact that the λ T parameters satisfy P2 implies (Lemma 1, part (ii)) that λ A1B1 =( Z A1B2 + Z A2B1 )/Z AB for every pair of 3-separated clades A 1 , B 1 in the configuration of Fig. 1(b). That is, the M&P formulae for w coincide with the given (γ T , λ T )-formulae, which is what we set out to prove.

But, because δ

p → d T *
,t h ep r o b a b i l i t yt h a tδ belongs to this neighborhood, and consequently M(δ) has the same topology as T * ,c o n v e r g e st o1 . F i n a l l y ,w h e nM(δ) has the same topology as T * ,t h e continuity of the branch length estimation scheme implies that the branch lengths in M(δ) converge in probability to those in M(d T * )=T * .W e c a n t h e n c o n c l u d e t h a t b o t h t h e t o p o l o g y a n d b r a n c h lengths of M(δ) consistently converge to T * .

A3.2 Consistency of Least Squares.

We have briefly defined LS methods in the Introduction. Here we assume the most general form for LS and prove its consistency under very general conditions (Proposition 4 below). We define LS methods as those that use a tree score function with the following form:

(13) Q(T,δ)=(δ -d T ) t W T,δ (δ -d T ),
where T is the tree fitted using the assumed branch length estimation scheme, δ is a column vector with the n 2 input distances and W T,δ is a n 2 × n 2 matrix which may depend on the topology of T and, continuously, on δ and the branch lengths of T .A d d i t i o n a l l y ,w ea s s u m et h a t ,f o ra n yT and δ,t h em a t r i xW T,δ is positive-definite. (W T,δ should be interpreted as the inverse of the assumed covariance matrix for δ.)

Note that whereas the dependence of W T,δ on δ is common (e.g. the version of WLS by Fitch and Margoliash [4] uses a diagonal matrix with W ij,ij = δ -2 ij ), the dependence on T is non-standard, and we have included it here for completeness. (But for example the balanced version of WLS [5] at the basis of the balanced branch lengths [6] does assume a variance model that depends on tree topology.) Also recall that the criterion Q(T,δ) above is used to score trees with already-fitted branch lengths, so the dependence on the branch lengths does not cause any computational problem. Proposition 4. Adopt any correct and continuous branch length estimation scheme. Then, LS is consistent.

Proof. We prove that for LS, the hyp otheses of Prop osition 3 are satisfied and therefore LS is consistent. First, the branch length estimation scheme is continuous (byhypothesis)andthescorefunctionin( 13) is a continuous function of both δ and of the branch lengths assigned to T (note that d T is linear, and thus continuous, in the branch lengths). It remains then to show that, for δ = d T * ,L Su n i q u e l y identifies T * as optimal. Because W T,δ is positive definite, Q(T,d T * )=0if and only if d T * d T = 0, that is, if and only if T = T * ,w h e r e a sf o ra l lo t h e rt r e e sW = T * , Q(W, d T * ) > 0.M o r e o v e r ,b e c a u s e the branch length estimation scheme is correct, T * is precisely what is obtained when fitting its branch lengths. Therefore T * uniquely minimizes the score function Q and is returned by LS.

Corollary 5. Adopt a branch length estimation scheme based on (γ, λ)-formulae. Then, LS is consistent.

A3.3 Dependency between the consistencies of different variants of Minimum Evolution.

Recall that ME i reconstructs the fitted tree that minimizes the following tree score function, where l =( le ) denotes the branch lengths in the fitted tree: In the main text, we assume i ∈ {-1, 0, +1},butherew econsider,moregenerally ,ME x with x being any real number. We do this not only for the sake of mathematical completeness, but also to include variants of ME that may be considered in the future (e. Proof. We adapt a line of reasoning app eared elsewhere [7,3]. Because of the correctness of the branch length estimation scheme, when the branches of the topology of T * are fitted using d T * ,t h e i rl e n g t h s are set to their correct values. Because these are all positive, the scores assigned by ME x and ME y to T * equal precisely the sum L * of all branch lengths in T * .N o wl e t l denote the branch lengths assigned to an incorrect topology using d T *

; x<yimplies that L x ( l) ≥ L y ( l),an dbecau seME y (d T * ) is unique and equal to T * ,w ea l s oh a v eL y ( l) >L * .B u tt h e nL x ( l) >L * for any incorrect topology. Because L * coincides with the score assigned to T * by ME x ,w ec a nt h e nc o n c l u d et h a tME x (d T *

) is unique and coincides with T * .

A3.4 Consistency of the classic version of Minimum Evolution.

We now concentrate on the consistency of ME +1 ,w h i c hbe c a u s eo fL e m m a6 ,i m p l i e st h a to fME x for any x<+1.W eu s eas t a n d a r df r a m e w o r kf o ri n v e s t i g a t i n gt h ec o n s i s t e ncy of ME +1 (e.g. [1]), which consists in verifying a property ("Willson's condition") ofb r a n c hl e n g t he s t i m a t i o ni nt h ep r e s e n c eo f as p e c i a lk i n do fb i n a r yd i s t a n c em a t r i x . I nt h ef o l l o w i n g ,w ei n t r o d u c ea n ds t a t eW i l l s o n ' sc o n d i t i o n (A3.4.1), then prove some properties of the (γ T , λ T )-formulae that are useful to verify it (A3.4.2), and finally prove the consistency of ME +1 via Willson's condition (A3.4.3). For simplicity, we write LT (δ) as a shorthand for L +1 ( lT (δ)),t h a ti s ,t h et r e el e n g t h( s e n s uME +1 )r e s u l t i n gf r o mfi t t i n gt h eb r a n c h lengths of T using δ.

A3.4.1. Willson's condition [1].

In order to state it (Lemma 7 below), we denote by d S|S (where S ⊆ {1, 2,...,n} and S = {1, 2,...,n}\ S)t h ef o l l o w i n gc o l l e c t i o no f n 2 distances, indexed by i = j:

d S|S ij = 1 if |S ∩ {i, j}| =1, 0 otherwise.
Lemma 7. Adopt any linear branch length estimation scheme such that, for any binary topology T and any proper and nonempty subset S of taxa from T ,

LT (d S|S ) =1 if S is a clade in T, > 1 otherwise.
Then ME +1 (d T * )=T * .

Informally, this holds because d T * is a weighed sum of all d S|S corresponding to the clades in T * ,t h e coefficients of this weighted sum being the branch lengths of T * .T h el i n e a r i t yo f LT then implies that, in turn, LT (d T *

) is a weighted sum of the branch lengths of T * ,w h e r et h ec o e fficients are now either 1o rs t r i c t l yg r e a t e rt h a n1 ,d e pe n d i n go nw h e t h e ro rn o tt h ec orresponding S is a clade in T .C l e a r l y this weighted sum is minimized when all the clades in T * are also clades in T ,t h a ti sf o rT = T * .

A3.4.2 Tools to verify Willson's condition.

Given any subset of taxa S ⊆ {1, 2,...,n} and any clade X in a binary topology T for which a set of (γ T , λ T )-formulae is defined, define p S|X as the probability of picking an element of S from the random distribution over X defined by the γ T parameters. That is,

p S|X = i∈X∩S p i|X .
Moreover, a clade X is monochromatic (w.r.t. S)i fe i t h e rX ⊆ S or X ⊆ S = {1, 2,...,n}\S. In this case it is clear that p S|X ∈ {0, 1}. (ii) Branch e is either internal or external. We consider here only the casew h e r ei ti si n t e r n a l ,a st h e external case is analogous (and simpler). We assume clades A 1 , A 2 , B 1 and B 2 are defined as in Fig. 1(b). Because e belongs to a monochromatic clade, at least three clades out of A 1 , A 2 , B 1 and B 2 are all subsets of S or all subsets of S.W i t h o u tl o s so fg e n e r a l i t y ,w ea s s u m et h a tA 1 , A 2 , B 1 are all subsets of S.A p p l y i n gp a r t( i ) ,i ti se a s yt os e et h a tt h i si m p l i e sd

S|S A1A2 = d S|S A1B1 = d S|S A2B1 =0and d S|S A1B2 = d S|S A1B2 = d S|S A1B2 = p S|B2 .B u tt h e n , le = 1 2 λ A1B1 (d S|S A1B1 + d S|S A2B2 )+(1-λ A1B1 )(d S|S A1B2 + d S|S A2B1 ) -d S|S A1A2 -d S|S B1B2 = 1 2 λ A1B1 p S|B2 +(1-λ A1B1 )p S|B2 -p S|B2 =0.
(iii) Let B = {1, 2,...,n}\(A 1 ∪ A 2 ),a si nF i g .1 ( b ) . I naw a ya n a l o g o u st op a r t( i i ) ,i ti se a s y to check that, independently of f being internal or external, lf = Proof. We show that any branch length estimation scheme based on (γ, λ)-formulae satisfies Willson's condition (i.e. the hypotheses of Lemma 7), and therefore we must have

ME +1 (d T * )=T * .
First, it is trivial to see that any such branch length estimation scheme is linear. Second, it is correct (Theorem 1), which implies LT (d S|S )=1whenever S is a clade of T :i nt h i sc a s e ,i nf a c t ,d S|S is additive with respect to a tree with topology T and with all branches of length 0 except the root branch of S,whic hhaslength1;becauseoftheircorrectness,the(γ T , λ T )-formulae result in assigning T precisely these branch lengths and therefore a total length LT (d S|S )=1.

It remains to prove that LT (d S|S ) > 1 whenever S is not a clade of T -f o ra n yb r a n c hl e n g t h estimation scheme based on (γ, λ)-formulae. We do this by induction on the size of T .

For n =3taxa, this is trivially true, as all proper, nonempty sets of {1, 2, 3} are clades of T .

For n>3,i fS is not a clade of T ,i ti sa l w a y sp o s s i b l et ofi n dap a i ro f2 -s e p a r a t e dc l a d e sA 1 ,A 2 such that A 1 S and A 2 S,i . e .s u c ht h a tA 1 and A 2 are monochromatic, but A 1 ∪ A 2 and B = {1, 2,...,n}\(A 1 ∪ A 2 ) are not monochromatic. To see this, consider the tree T (S) that is obtained by substituting every monochromatic clade in T with a taxon; because S is not a clade of T , then T (S) must have at least two cherries (i.e., pairs of 2-separated taxa); any of these corresponds to ap a i ro fc l a d e sA 1 , A 2 in the original tree T with the required properties. Let e, f and g be the root branches of B, A 1 and A 2 ,r e s p e c t i v e l y ,a n da their common endpoint, as in Fig. 1(b). Because A 1 and A 2 are monochromatic, it is clear that all branches belonging tot h e s ec l a d e sa r ea s s i g n e dl e n g t h 0w h e nt h ei n p u td i s t a n c e sa r ed = d S|S (Lemma 8, part (ii)). Now let T ′ be the topology that is obtained by deleting from T all branches belonging to A = A 1 ∪ A 2 , so that a is a leaf of T ′ .I ti sc l e a rt h a tt h e r ei so n e -t o -o n ec o r r e s p o n d e n c eb e t w e e nthebranches/clades of T ′ and a subset of the branches/clades of T .W h e nc a l c u l a t i n gt h eb r a n c hl e n g t h si nT ′ ,w eassu m e that the γ ef and λ XY parameters are the same as those for the corresponding branches/clades in T . Now define the following distances over {1, 2,...,n} ∪ {a}\A,i . e .t h et a x ai nT ′ : (2) .

d (1) = d (S∪{a}\A1)|(S\A2) , d (2) = d (S\A1)|(S∪{a}\A2) , d ′ = γ ef d (1) + γ eg d
Note that d ′ coincides with d except for the distances that involve a.F o rt h e s e ,w eh a v e

d ′ aj = γ ef d (1
d ′ XY = d XY if a/ ∈ X ∪ Y, d X ′ Y if a ∈ X and where X ′ = X ∪ A \{a},
that is, the average distances between disjoint clades remain the same when going from (T,d) to (T ′ , d ′ ).T os e et h i s ,n o t et h a tt h efi r s tc a s ei sas i m p l ec o n s e q u e n c eo ft h ef a c tt h a td ′ coincides with d for the distances that do not involve a.A sf o rt h es e c o n dc a s e ,i ti sac o n s e q u e n c eo fc o m b i n i n g d ′ aj = d Aj (shown above) with the first case (and it can be easily proved byi n d u c t i o no nt h es i z eo f X).

The important consequence of ( 15) is that the lengths of branches belonging to B -w h i c ho n l y depend on average distances between disjoint clades in T ′ -r e m a i nc o n s t a n tw h e ng o i n gf r o m(T,d) to (T ′ , d ′ ).T h e r e f o r e ,t h eo n l yd i fference between the two tree lengths will come from the lengths of branches e, f and g: 2) ).

LT ′ (d ′ )=γ ef LT ′ (d (1) )+γ eg LT ′ (d (
Because d (1) and d (2) are equal to some d S ′ |S ′ where S ′ is a proper and nonempty subset of the n ′ (<n) taxa in T ′ ,w ec a na p p l yt h ei n d u c t i o nh y p o t h e s i sa n di n f e rt h a t LT ′ (d (1) ) ≥ 1 and LT ′ (d (2) ) ≥ 1 (where equality is achieved when S ′ is a clade of T ′ ). Therefore LT ′ (d ′ ) ≥ 1.

In order to prove that the tree length in ( 16) is strictly greater than 1, we then just need to prove that lf + lg + le -l′ e > 0. First, because of Lemma 8, part (iii), lf + lg =1.I no r d e rt oc a l c u l a t e le and l′ e ,i ti su s e f u lt on o t e that d A1A2 =1and, for any clade

B ′ ⊆ B, d A1B ′ = p S|B ′ =1-p S|B ′ and d A2B ′ = p S|B ′ (Lemma 8, part (i)). Then, le = 1 2 λ A1B1 (1 -p S|B1 + p S|B2 )+(1-λ A1B1 )(1 -p S|B2 + p S|B1 ) -1 -d B1B2 = 1 2 (1 -2λ A1B1 )(p S|B1 -p S|B2 ) -d B1B2 .
Similarly,

l′ e = 1 2 d ′ aB1 + d ′ aB2 -d ′ B1B2 = 1 2 [d AB1 + d AB2 -d B1B2 ] = 1 2 [γ ef (d A1B1 + d A1B2 )+(1-γ ef )(d A2B1 + d A2B2 ) -d B1B2 ] = 1 2 γ ef (2 -p S|B1 -p S|B2 )+(1-γ ef )(p S|B ! + p S|B2 ) -d B1B2 = γ ef + 1 2 (1 -2γ ef )(p S|B ! + p S|B2 ) -d B1B2 . Then, lf + lg + le -l′ e =1-γ ef + p S|B1 (γ ef -λ A1B1 )+p S|B2 (γ ef -(1 -λ A1B1 ))
. But this is a linear function of (p S|B1 ,p S|B2 ) in the square [0, 1] 2 ,a n di st h u sm i n i m i z e di no n eo fi t s four vertices. In (0, 0), (0, 1), (1, 0) and (1, 1),t h ef u n c t i o nh a sv a l u e s1 -γ ef , λ A1B1 , 1 -λ A1B1 , γ ef , respectively. Since these are all strictly greater than 0 by hypothesis, then so is lf + lg + le -l′ e and therefore LT (d) > 1.

A3.5 Wrapping it all together.

By applying Proposition 9, Lemma 6 and Proposition 3, we then conclude: Corollary 10. Adopt a branch length estimation scheme based on (γ, λ)-formulae. Then, for any x ≤ +1, ME x is consistent.

Which, together with Corollary 5 and Assumption 2, completeso u rp r oo fo fT h e o r e m3 .

SI Appendix 4

Here, we prove the efficiency of calculating branch lengths with our formulae in hill climbing heuristics, as stated in Theorem 4. We start by showing that efficient branchlengthcalculationsessentiallydepend on the availability of the average distances between (some) clades in the current tree, and that these can be calculated in quadratic time, which allows us to prove Theorem 4, part (i) (A4.1). When performing an NNI, calculating the new branch lengths can be done efficiently by recalculating only some of these average distances, which leads us to prove Theorem 4, part (ii) (A4.2). Finally, we show that updating the average distances following an NNI can alsobedoneefficiently (A4.3), which is not a claim of Theorem 4, but is nevertheless potentially useful. The results and proofs here are inspired by those of Desper and Gascuel [8]. However, their results were specific to the balanced and OLS branch lengths in combination with the ME +1 optimization principle. A key property of these estimators is that, when performing an NNI, the sum of the branch lengths in each of the four corner subtrees around the location of the NNI remains constant. Thanks to this property, the difference between the ME +1 lengths of any two NNI neighbors T and T ′ can be efficiently calculated using simple formulae. This property does not hold in general for (γ T , λ T )-estimators, and so we have to recalculate all branch lengths every time we perform an NNI. The good news is that (1) the complexity of each iteration in ah i l lc l i m b i n gh e u r i s t i c sf o rM E( c o m p u t i n gt h el e n g t ho fa l lN N In e i g h b o r so fag i v e nt o p o l o g ya n d updating the data structures for the new best topology), which for BME was quadratic in the worst case, remains quadratic in the size of the tree, and that (2) recalculating all branch lengths makes it possible to use optimization principles such as ME 0 , ME -1 and ME -∞ .

In the following, T , T ′ and T i always denote binary topologies, and γ T , γ T ′ and γ Ti ,c o l l e c t i o n so f γ ef parameters defined for them, in the way described in the main text.

A4.1 Computing the branch lengths of a fixed topology. Lemma 11. Adopt a set of (γ T , λ T )-formulae for the branch lengths of T .G i v e nδ and δ XY for every pair of 3-separated clades X, Y in T ,t h el e n g t ho fa n yb r a n c hi nT can be calculated in O(1) time.

Proof. The (γ T , λ T )-formulae are simple linear combinations of average distances δ XY between 2-and 3-separated clades X and Y and can be computed in O(1) once these average distances are available. Since we assume that the average distances between 3-separated clades are given, it remains to show that δ XY between any pair X, Y of 2-separated clades can be obtained in O(1).B u tt h i si st r i v i a l : either both X and Y consist of one taxon only, i.e. X = {i} and Y = {j},i nwh i c hcaseδ XY = δ ij ,o r at least one of the two clades, say X,i ss u c ht h a tX = X 1 ∪ X 2 ,w h e r eb o t hX 1 and X 2 are clades, in which case δ XY = γ ef δ X1Y + γ eg δ X2Y ,w h e r ee, f and g are the root branches of X, X 1 and X 2 , respectively, and both δ X1Y and δ X2Y are known, as X i and Y are 3-separated (for i ∈ 1, 2}).

Although the one above is a straightforward observation, it determines the minimum amount of information necessary to determine any branch length in T in constant time. Motivated by it, we define ∆ T (γ T ) as a data structure holding all the average distances δ XY between pairs of disjoint clades in T ,andmakeexplicititsdependenceonγ T .N o t et h a t∆ T (γ T ) specifies the average distances between 3-separated clades as a particular case. Lemma 12. Given δ, T and γ T ,t h ec a l c u l a t i o no f∆ T (γ T ) requires O(n 2 ) time.

Proof. Consider any total ordering A 1 ,A 2 ,...,A 2(2n-3) of the clades in T ,s u c ht h a ti fA k = A i ∪ A j then i<kand j<k .F i n d i n go n es u c ho r d e r i n gi st r i v i a la n dc a nb ed o n ei nan u m b er of ways, for example by sorting the clades in ascending order of depth [2],orb yroot i n gt h et reei non eofi t sl ea v es and then performing a postorder traversal, listing the clades oriented away from the root, followed by a preorder traversal, listing the clades oriented towardst h er o o t . T h ef o l l o w i n gp r o c e d u r et h e n calculates δ XY for all pairs of clades (including non-disjoint ones):

For i =2,...,2(2n -3), in T i ,r e q u i r e sO(n) time (Lemma 13). Combining this to Lemma 11 yields that all O(n) branch lengths in T i can be calculated in O(n) time. Since there are O(n) neighbors of T and each is treated in O(n) time, the whole calculation requires O(n 2 ) time.

A4.3 Updating the accessory information when performing an NNI. The proof of Lemma 13 above suggests a related result which may also be useful forh i l lc l i m b i n g ,w h e nt h eb e s tN N Ineighbor T ′ of T has been identified and we need to calculate ∆ T ′ (γ T ′ ) in order to explore efficiently the NNI-neighborhood of T ′ .D e fi n ediam(T ),thediameterofT ,asthemaxim umn um berofbranc hes separating any two leaves of T . Proposition 14. Let T and T ′ be NNI-neighbors and let γ T and γ T ′ be almost identical. Given ∆ T (γ T ),i t su p d a t ei n t o∆ T ′ (γ T ′ ) requires O(n • diam(T )) time.

Proof. Let T be as in Fig. 1(b) and T ′ as in Fig. S2. Let δ XY and δ ′ XY denote average clade distances from ∆ T (γ T ) and ∆ T ′ (γ T ′ ), respectively. In order to obtain ∆ T ′ (γ T ′ ) from ∆ T (γ T ),o n en e e d st o calculate the entries of ∆ T ′ (γ T ′ ) that have no corresponding entry in ∆ T (γ T ) or those that have changed. These are the δ ′ XY for all pairs of clades X, Y in T ′ such that some of the branches f, g, h, l belong to X or Y .T h eo n l yc a s ew h e r eb o t hX and Y have at least one of f, g, h, l belonging to them is that where X = A 1 ∪ B 2 and Y = A 2 ∪ B 1 .I nt h i sc a s e ,δ ′ XY can be obtained from ∆ T (γ T ) with

δ ′ XY = γ ′ e ′ f γ ′ e ′ g δ A1A2 + γ ′ e ′ f γ ′ e ′ h δ A1B1 + γ ′ e ′ g γ ′ e ′ l δ A2B2 + γ ′ e ′ h

Figure 1 .

 1 Figure 1. Standard naming of clades and branches when (a) e is external, and (b) e is internal. (c) An NNI-neighbor (around e)o ft h et r e ei n( b ) .

e

  (δ).N o t et h a t labe (δ), la1b1a2b2 e (δ) and la1b2a2b1 e (δ) are all trivially correct.

Figure S1 .

 S1 Figure S1. (i) p i|A = γ e0e1 • γ e1e2 • ...• γ e k-1e k is the probability of ending up in i when entering clade A from its root and following the random walk rules described in the main text; (ii) λ A1B1 = λ A2B2 can be seen as the probability of drawing the tree in the top configuration, whilst λ A1B2 = λ A2B1 =1λ A1B1 can be seen as the probability of drawing the tree in the bottom configuration.

e

  (δ).T h e o r e m1i st h e r e f o r ep r o v e d .

Fig. 1 (

 1 b) in the main text. A2.2 Decomposition of Z XY . Extend the w ij notation to any pair of nodes x and y (possibly internal) in T :

Lemma 1 .

 1 Given weights w =(w ij ) multiplicative w.r.t. a binary topology T ,d e fi n e ,f o re a c hp a i ro f adjacent branches e and f : (4) γ ef = w f Z A1 Z A , where A 1 and A = A 1 ∪ A 2 are the clades having f and e as root branches, respectively, as in Fig. 1(b). Then, (i) The resulting average distances between clades are such that δ XY = δ w XY ,f o ra n yt w od i s j o i n tc l a d e s X and Y in T . (ii) For any internal branch e,l e tA = A 1 ∪ A 2 and B = B 1 ∪ B 2 be the clades in the configuration of Fig. 1(b); then (Z A1B2 + Z A2B1 )/Z AB = γ ef + γ eh -2γ ef γ eh .

Proof of Theorem 2 ,

 2 part (i). Given w multiplicative w.r.t. T ,d e fi n et h eγ T parameters as in(4). Note that because w e > 0 for any branch e,t h e nZ X > 0 for every clade X.M o r e o v e r ,f o ra n yt h r e e adjacent branches e, f and g in the configuration of Fig.1(b), we have w f Z A1 + w g Z A2 = Z A and w f Z A1 ,w g Z A2 > 0,w h i c hi m p l yγ ef + γ eg =1and 0 < γ ef , γ eg < 1.T h e r e f o r et h ed e fi n i t i o no ft h e γ T parameters is admissible and implies that δ XY = δ w XY (Lemma 1, part (i)). As for the λ T parameters, set λ A1B1 =( Z A1B2 + Z A2B1 )/Z AB ,f o re v e r yp a i ro fc l a d e sA 1 and B 1 separated by 3 branches, and being in the configuration of Fig.1(b)withA = A 1 ∪A 2 and B = B 1 ∪B 2 .

Lemma 8 .

 8 Assume that the branch lengths of T are assigned with (γ T , λ T )-formulae using the input distances in d S|S ,f o rs o m eS ⊆ {1, 2,...,n}.T h e n , (i) For any two clades X and Y in T ,t h ea v e r a g ed i s t a n c eb e t w e e nt h e md S|S XY = p S|X p S|Y + p S|X p S|Y . (ii) If branch e belongs to a monochromatic clade, it is assigned length le =0. (iii) If adjacent branches f and g separate (and are the only ones to separate) two monochromatic clades A 1 and A 2 ,w i t hA 1 ⊆ S and A 2 ⊆ S,t h e n lf + lg =1. Proof. (i) d S|S XY = i∈X j∈Y p i|X p j|Y d S|S ij = i∈X∩S j∈Y ∩S p i|X p j|Y + i∈X∩S j∈Y ∩S p i|X p j|Y = p S|X p S|Y + p S|X p S|Y .

1 2 ( 1 +Proposition 9 .

 219 p S|B -p S|B ) and, similarly, lg = 1 2 (1 + p S|B -p S|B ).T h e r e f o r e , lf + lg =1. Points (i), (ii) and (iii) are thus all verified. A3.4.3 Consistency of ME +1 with perfect data. Adopt a branch length estimation scheme based on (γ, λ)-formulae. Then, ME +1 (d T * )= T * .

( 16 )

 16 LT (d)= LT ′ (d ′ )+ lf + lg + le -l′ e , where lb and l′ b represent the lengths assigned to branch b for (T,d) and (T ′ , d ′ ),r e s p e c t i v e l y . Because for ME +1 the tree length is a linear function of the branch lengths and the branch lengths themselves are linear functions of the input distances, LT ′ (d ′ ) is then linear in d ′ and we can write

  γ ′ e ′ l δ B1B2 . All the other cases correspond to a pair of clades X, Y such that one of them, say X,i si n c l u d e di n one of the four corner clades A 1 ,A 2 ,B 1 ,B 2 and the other, Y ,i n c l u d e st w oo rt h r e eo ft h eo t h e rc l a d e s (see e.g. Fig.S2, whereX ⊂ A 1 and Y ⊇ B 1 ∪ B 2 ∪ A 2 ).It is clear that for any such X,t h en u m b e r of possible choices for Y equals the number of branches in the path starting with e ′ and ending in the root of X.I no t h e rw o r d s ,t h e r ea r eO(n) possible choices for X,e a c ho fw h i c hc o r r e s p o n d st o diam(T ) choices for Y .T h e r e f o r ew en e e dt oc o n s i d e rO(n diam(T )) pairs of clades. For each of these pairs, we now prove that δ ′ XY can be calculated in O(1) time from ∆ T (γ T ),o n c es t e p s( 1 ,1 a ,1 b ,1 c ) from the proof of Proposition 13 have been executed (inO(n) time): if Y = A 2 ∪ B 1 or Y = A 1 ∪ B 2 , then it is straightforward to obtain δ ′ XY as γ ′ e ′ g δ XA2 + γ ′ e ′ h δ XB1 or as γ ′ e ′ f δ XA1 + γ ′ e ′ l δ XB2,respectively; otherwise, if Y includes three of A 1 ,A 2 ,B 1 ,B 2 ,i ti se a s yt os e et h a t( 2 0 ) ,a n ds i m i l a re q u a t i o n sf o r X ⊂ A 2 ,B 1 ,B 2 ,s t i l lh o l d( w i t h o u tt h ea s s u m p t i o n ,m a d ei nt h ep r o o fo fL e m ma 13, that X and Y are 3-separated). It is then possible to calculate each δ ′ XY in constant time.

  Lemma 6. Adopt any correct branch length estimation scheme and let x<y . If ME y (d T *)=T * , then also ME x (d T

g. ME -∞ ,w h i c hc o r r e s p o n d st oa v o i d i n ga t all costs trees which are assigned negative branch lengths).T h ef o l l o w i n gp r o p o s i t i o ns h o w st h a ti fw e can prove the consistency of ME y using Proposition 3, then the same can be done for any ME x with x<y. * )=T * .

  A2j ,w h i c hi m p l yd ′ aj = d Aj . This in turn implies that, for any disjoint clades X and Y in T ',

	) aj + γ eg d (2) aj ,
	for every j ∈ {1, 2,...,n}\A.N o t ea l s ot h a td (1) aj = d A1j and d (2)
	(15)

aj = d

SI Appendix 3

Here, we prove that the main criteria to score trees, LS and ME (in all their common variants), are statistically consistent when used in combination with our branch length formulae, as stated in Theorem 3. We start by showing that the consistency of any distance-based principle is essentially determined by its behavior on perfect data (A3.1). Next, we move on to proving the consistency of LS (A3.2) and then that of ME: for the latter, first we show a useful dependency property between different variants of ME (A3.3), and then we prove the consistency on perfect data of the classic version of ME (A3.4), which is the key nontrivial result of this appendix and allows us to conclude the proof of Theorem 3 (A. 3.5).

We recall that a branch length estimation scheme is a method that, for any binary topology over the set of taxa under consideration {1, 2,...,n} (n ≥ 3), determines how to fit the length of its branches on the basis of an n×n distance matrix δ.W es a yt h a tab r a n c hl e n g t he s t i m a t i o ns c h e m ei scontinuous [linear ] if, for any branch e in any binary topology, the function le (δ) giving its fitted length is continuous [linear] in δ.Ab r a n c h l e n g t he s t i m a t i o n s c h e m ei scorrect if, for any branch e in any binary tree with branch lengths, le (δ) returns the length of e whenever δ is additive with respect to that tree (as in Theorem 1). The branch length estimation schemes that we consider here are those based on (γ, λ)-formulae,w h e r e b yac o l l e c t i o no fγ T and λ T parameters is chosen for each binary topology T , thus determining a set of (γ T , λ T )-formulae for T 's branch lengths (with no assumed relation between the values of these parameters across different topologies). It is clear that the resulting branch length estimation schemes are linear (thus continuous) and correct( T h e o r e m1 ) .

Any branch length estimation scheme can be combined to a number of principles identifying an optimal tree among all the topologically-distinct fitted trees. The optimization principles we consider here are defined by a tree score function,w h i c hc a nd e p e n do nt h et o p o l o g yo ft h et r e e ,t h ea s s i g n e db r anch lengths and (in the case of LS, but not ME) the input distances.A no p t i m i z a t i o np r i n c i p l eM then consists of seeking the fitted tree(s) that minimize this function, and we denote this tree (or set of trees) with M(δ).W es a yt h a tM is statistically consistent if M(δ) converges (in probability) to the correct tree.

The following assumption (the consistency of the distance estimates and the positive additivity of the correct evolutionary distances) applies to all the propositions that follow, and we state it here so that we do not have to repeat it in every statement. Assumption 2. Let the correct phylogenetic tree for the taxa under consideration, T * , be a binary tree with positive branch lengths. Assume that the input distances δ converge (in probability) to d T * .

A3.1 Consistency for perfect data implies statistical consistency.

The following is a well-known sufficient condition for consistency, which has been proven for ME with the same continuity arguments (e.g. [3]). It can be applied tom o s tt r e eo p t i m i z a t i o np r i n c i p l e s( L S , ME, possibly combinations of the two or even totally differentc r i t e r i a ) .

Proposition 3. Adopt a continuous branch length estimation scheme. Let M be an optimization principle based on a tree score function that is continuous in all its continuous parameters (i.e. all but the topology). If M(d T *

) is unique and coincides with T * ,t h e nM is statistically consistent.

Proof. To any binary top ology, M assigns a score by first assigning branch lengths to it, and then applying the adopted tree score function. Note that because both the branch length estimation scheme and the tree score function are continuous, then also the score associated to any particular topology is continuous in δ.B e c a u s eM(d T * )=T * is unique, when δ = d T * the score of the topology of T * must be strictly smaller than the score of all other binary topologies. But then, because the scores of topologies are continuous in δ and finite in number, this must still hold for every δ in a neighborhood of d T * ;t h a ti s ,f o re v e r yδ in this neighborhood, M(δ) is unique and has the same topology as T * .

for j =1,...,i-1,

In the second case of [START_REF] Rzhetsky | Theoretical foundation of the minimum-evolution method of phylogenetic inference[END_REF] we assume that e, f, g are the root branches of A i ,A i1 ,A i2 ,r e s p e c t i v e l y , whereas in the third case they are the root branches of A j ,A j1 ,A j2 ,r e s p e c t i v e l y . N o t et h a tt h e s et w o cases are not mutually exclusive, and the result is the same independently of which case is applied. Moreover, because of the way the ordering is defined, we must have i 1 ,i 2 <i ,i nt h es e c o n dc a s e ,o r j 1 ,j 2 <j,i nt h et h i r dc a s e ,w h i c hm e a n st h a tδ Ai 1 Aj and δ Ai 2 Aj (second case), or δ AiAj 1 and δ AiAj 2 (third case) have already been calculated and are available when we calculate δ AiAj .B e c a u s ee a c h δ AiAj can be calculated in constant time, the whole calculation requires O(n 2 ) time.

The complexity we obtain in Lemma 12 is optimal. Even if we restrict the calculation to 3-separated clades, we still cannot do better than O(n 2 ),ast h ea v eraged i st an cesbet w eensu c hp ai rsofcl ad esst i l l depend on O(n 2 ) input distances.

Proof of Theorem 4, part (i).C o m b i n i n gL e m m a1 2w i t hL e m m a1 1y i e l d st h a tt h eb r a n c hl e n g ths determined by a set of (γ T , λ T )-formulae for a binary topology T can be calculated in O(n 2 ) time.

A4.2 Computing the branch lengths of the NNI neighbors of a given topology. Lemma 13. Let T and T ′ be NNI-neighbors and let γ T and γ T ′ be almost identical. Let δ XY and δ ′ XY denote the average clade distances in ∆ T (γ T ) and ∆ T ′ (γ T ′ ),r e s p e c t i v e l y . T h e n ,g i v e n∆ T (γ T ), the calculation of δ ′ XY for every pair of 3-separated clades X, Y in T ′ ,r e q u i r e sO(n) time.

Proof. We assume that T is as in Fig. 1(b) and T ′ as in Fig. S2. Let the elements of γ T and γ T ′ be denoted by γ e1e2 and γ ′ e1e2 ,r e s p e c t i v e l y . First, we show that δ ′ XY is straightforward to obtain in the case of pairs of 3-separated clades in T ′ such that none or 1 of the 3 branches separating X and Y belongs to one of the corner clades

w h e r ea l lt h es e t si n v o l v e da r e also clades of T ′ (and therefore T ). It is trivial to verify that (18)

, as γ T and γ T ′ are the same within all the clades in the subscripts above. Nowo b s e r v et h a tA 1 ∪ B 2 and A 2 ∪ B 1 are clades in T ′ but not in T .T h e i ra v e r a g ed i s t a n c e sw i t ho t h e r3 -s e p a r a t e dc l a d e sm u s t then be obtained with expressions such as

We have therefore proved that these δ ′ XY can be obtained from one or two corresponding entries in ∆ T (γ T ) in O(1) time.

We still have to show how to derive δ ′ XY when 2 or all 3 of the branches separating X and Y belong to a corner clade. Without loss of generality, we assume this clade to be A

.L e tc l a d e sY 1 ,Y 2 ,...,Y k be defined as in Fig. S2 (andnote

as the probability that the random walk defined by the γ T parameters reaches Y ′ ,a s s u m i n gt h a ti te n t e r sY from its root branch:

, where the last equality uses the almost identity of γ T and γ T ′ .S i m i l a r l y ,

Therefore, (20)

It is easy to derive similar equations for the cases where (a)

(2) Use (20), or similar equation, to derive δ ′ XY for all 3-separated clades X, Y in T ′ such that 2 or all 3o ft h eb r a n c h e ss e p a r a t i n gX and Y belong to a corner clade A 1 ,A 2 ,B 1 ,B 2 .

(3) Use the simple equations in ( 18) and ( 19) to calculate δ ′ XY for the remaining 3-separated clades.

Step (1) can be done in O(n) time, by starting with the smallest clades including B 1 ∪ B 2 ∪ A 2 and using the derived values to calculate those for the larger clades. The same holds for steps (1a, 1b, 1c). Then, each δ ′ XY can be calculated in O(1) time. Since there are O(n) 3-separated pairs of clades, the entire algorithm runs in O(n) time and thus Lemma 13 is proved.

We are now ready to complete the pro of of Theorem 4. Theorem 4,part (ii).R e c a l l t h a t a l l t h e b r a n c h l e n g t h s i n T 0 and its NNI-neighbors T 1 ,T 2 ,...,T 2(n-3) are defined by (γ Ti , λ Ti )-formulae, with the constraint that γ Ti and γ T0 are almost identical. We wish to prove that the branch lengths of T 1 ,T 2 ,...,T 2(n-3) can be calculated in O(n 2 ) time. Let δ (i) XY denote the average clade distances in ∆ Ti (γ Ti ).B e c a u s eo fL e m m a1 2 ,∆ T0 (γ T0 ) can be calculated in O(n 2 ) time. From this, the calculation of δ

Proof of