Identifiability and Reconstructibility Under a Modified Coalescent

# Identifiability and Reconstructibility of Species Phylogenies Under a Modified Coalescent

Colby Long and Laura Kubatko
###### Abstract.

Coalescent models of evolution account for incomplete lineage sorting by specifying a species tree parameter which determines a distribution on gene trees. It has been shown that the unrooted topology of the species tree parameter of the multispecies coalescent is generically identifiable. Moreover, a statistically consistent reconstruction method called SVDQuartets has been developed to recover this parameter. In this paper, we describe a modified multispecies coalescent model that allows for varying effective population size and violations of the molecular clock. We show that the unrooted topology of the species tree for these models is generically identifiable and that SVDQuartets is still a statistically consistent method for inferring this parameter.

## 1. The multispecies coalescent

The goal of phylogenetics is to reconstruct the evolutionary history of a group of species from biological data. Most often, the data available are the aligned DNA sequences of the species under consideration. The descent of these species from a common ancestor is represented by a phylogenetic tree which we call the species tree. However, it is well known that due to various biological phenomena, such as horizontal gene transfer and incomplete lineage sorting, the ancestry of individual genes will not necessarily match the tree of the species in which they reside [23, 27, 18]. There are various phylogenetic reconstruction methods that account for this discrepancy in different ways. One approach is to reconstruct individual gene trees by some method and them utilize this information to infer the original species tree [16, 15, 32, 20, 21].

The multispecies coalescent model incorporates incomplete lineage sorting directly. The tree parameter of the model is the species tree, an -leaf rooted equidistant tree with branch lengths. The species tree yields a distribution on possible gene trees along which evolution is modeled by a -state substitution model. For a fixed choice of parameters, the multispecies coalescent returns a probability distribution on the possible -tuples of states that may be observed. In order to infer the species tree from data, one searches for model parameters yielding a distribution close to that observed, using, for example, maximum likelihood.

In [5], the authors show that given a probability distribution from the multispecies coalescent model it is possible to infer the topology of the unrooted species tree parameter. Restricting the species tree to any 4-element subset of the leaves yields an unrooted 4-leaf binary phylogenetic tree called a quartet. For a given label set, there are only three possible quartets which each induce a flattening of the probability tensor. Given a probability distribution arising from the multispecies coalescent, the flattening matrix corresponding to the quartet compatible with the species tree will be rank or less while the other two will generically have rank strictly greater than this value. Since the topology of an unrooted tree is uniquely determined by quartets [24], these flattening matrices can be used to determine the unrooted topology of the species tree exactly. Of course, empirical and even simulated data produced by the multispecies coalescent will only approximate the distribution of the model. Therefore, the same authors also proposed a method called SVDQuartets [4], which uses singular value decomposition to infer each quartet topology by determining which of the flattening matrices is closest to the set of rank matrices.

The method of SVDQuartets offers several advantages over other existing phylogenetic reconstruction methods. It is statistically consistent, accounts for incomplete lineage sorting, and is computationally much less expensive than Bayesian methods achieving the same level of accuracy. It is often under-appreciated that this reconstruction method can be used to recover the species tree for several different underlying nucleotide substitution models without any modifications. It was shown in [5] that the method of SVDQuartets is statistically consistent if the underlying model for the evolution of sequence data along the gene trees is the 4-state general time-reversible (GTR) model or any of the commonly used submodels thereof (e.g., JC69, K2P, K3P, F81, HKY85, TN93). Thus, the method does not require any a priori assumptions about the underlying nucleotide substitution process other than time-reversibility.

In this paper, we show that the method of SVDQuartets has more theoretical robustness even than has already been shown. We will specifically focus on the case where the underlying nucleotide substitution model is one of the the 4-state models most widely used in phylogenetics. We describe several modifications to the classical multispecies coalescent model to allow for more realistic mechanisms of evolution. For example, we remove the assumption of a molecular clock by removing the restriction that the species tree be equidistant. We also allow the effective population size to vary on each branch and show how these two relaxations together can be interpreted as a scaling of certain rate matrices in the gene trees. Remarkably, we show that the unrooted species tree of these modified models is still identifiable and that SVDQuartets is a statistically consistent reconstruction method for these modified models. Thus, despite the introduction of several parameters, effective and efficient methods for reconstructing the species tree of these modified coalescent models are already available off the shelf and implemented in PAUP [25].

In Section 2, we review the classical multispecies coalescent model and discuss some of its limitations in modeling certain biological phenomena. We then describe several modifications to the classical model to remedy these weaknesses. In Section 3, we establish the theoretical properties of identifiability for these families of modified coalescent models. Finally, in Section 4, we describe why SVDQuartets is a strong candidate for reconstructing the species tree of these models and propose several other modifications that we might make to the multispecies coalescent.

## 2. The multispecies Coalescent

### 2.1. Coalescent Models of Evolution

In this section, we briefly review the multispecies coalescent model and explain how the model yields a probability distribution on nucleotide site patterns. As our main results will parallel those found in [5], we will import much of the notation from that paper and refer the reader there for a more thorough description of the model.

The Wright-Fisher model from population genetics models the convergence of multiple lineages backwards in time toward a common ancestor. Beginning with lineages from the current generation, the model assumes discrete generations with constant effective population size . In each generation, each lineage is assigned a parent uniformly from the previous generation. For diploid species, there are copies of each gene in each generation, and thus the probability of selecting any particular gene as a parent is . Two lineages are said to coalesce when they select the same parent in a particular generation.

As an example, if we begin with two lineages in the same species, the probability they have the same parent in the previous generation, and hence coalesce, is and the probability that they do not coalesce in this generation is . Therefore, the probability that two lineages coalesce in exactly the th previous generation is given by

 (12N)(1−12N)i−1.

For large , the time at which the two lineages coalesce, , follows an exponential distribution with rate , where time is measured in number of generations. Every generations is called a coalescent unit and time is typically measured in these units to simplify the formulas for time to coalescence. However, in this paper, we will introduce separate effective population size parameters for each branch of the species trees. So that our time scale is consistent across the tree we will work in generations rather than coalescent units. In these units, for lineages, the time to the next coalescent event has probability density,

 f(tj)=j(j−1)2(12N)exp(−j(j−1)2(12N)tj), tj>0.

This is typically referred to as Kingman’s coalescent [12, 13, 11, 28, 30, 29].

The multispecies coalescent is based on the same framework, but we assume that the species tree of the sampled taxa is known. We let denote the topology (without branch lengths) of the -leaf rooted binary phylogenetic species tree. We then assign a positive edge weight to each edge of . The weights are chosen to given an equidistant edge weighting, meaning that the sum of the weights of the branches on the path between any leaf and the root is the same. To specify an equidistant edge-weighting, it is enough to specify the vector of speciation times . The weights give the length of each branch and so we let denote the species tree with branch lengths. The condition that be equidistant implies that the expected number of mutations along the path from the root to each leaf is the same. This is what is known as as imposing a molecular clock.

Once this species tree is fixed, the multispecies coalescent gives a probability density on possible gene trees. All of the same assumptions as before apply, except that it is now impossible for two lineages to coalesce if they are not part of the same population. Hence, lineages may only coalesce if they are in the same branch of . We denote the topology of a rooted -leaf binary phylogenetic gene tree by and let denote the vector of coalescent times. As for the species tree, then represents a gene tree with weights corresponding to branch lengths. A coalescent history refers to a particular sequence of coalescent events as well as the branches along which they occur. Note that there are only finitely many possible coalescent histories compatible with and we call the set of all such histories . For a particular history, we can compute the probability density for gene trees with that history explicitly. We write the gene tree density of gene trees conforming to a particular coalescent history in as .

###### Example 2.1.

Let be the 4-leaf species tree depicted in Figure 1 and let refer to the coalescent history in which

1. Lineages and coalesce in .

2. Lineages and coalesce in .

3. Lineages and coalesce in .

The gene tree depicted inside of is compatible with .

The probability of observing a gene tree with history for the species tree under the multispecies coalescent model is

 ∫∞0∫t30∫τabcd−τab0(12N)3exp(−t12N)exp(−(τabcd−τcd)2N)exp(−3t22N)exp(−t32N) dt1dt2dt3.

The term involving is the probability that lineages and do not coalesce in the branch .

Once the gene tree for a locus is specified, we model the evolution along this gene tree as a continuous-time homogenous Markov process according to a nucleotide substitution model. The model gives a probability distribution on the set of all possible -tuples of observed states at the leaves of . We can write the probability of observing the state as . Precisely how this distribution is calculated is described in [5]. Here, we sketch the relevant details needed to introduce the modified multispecies coalescent model described in the next section.

For a -state substitution model, there is a instantaneous rate matrix where the entry encodes the rate of conversion from state to state . The distribution of states at the root is where is the stationary distribution of the rate matrix . To compute the probability of observing a particular state at the leaves, we associate to each vertex a random variable with state space equal to the set of possible states. To the root, we associate the random variable where . Letting be the length of edge , is the matrix of transition probabilities along that edge. That is, . Given an assignment of states to each vertex of the tree we can compute the probability of observing this state using the Markov property and the appropriate entries of the transition matrices. To determine the probability of observing a particular state at the leaves, we marginalize over all possible states of the internal nodes.

In this paper, we are primarily interested in 4-state models of DNA evolution where the four states correspond to the DNA bases. Different phylogenetic models place different restrictions on the entries of the rate matrices. The results that we prove in the next section will apply when the underlying nucleotide substitution model is any of the commonly used 4-state time-reversible models. We note here that because these models are time-reversible the location of the root in each gene tree is unidentifiable from the distribution for each gene tree [7]. In subsequent sections, we will introduce and describe similar results for the JC+I+ model, that allows for invariant sites and gamma-distributed rates across sites.

Now, given a species tree and a choice of nucleotide substitution model, let be the probability of observing the site pattern at the leaves of . To account for gene trees with topologies different from that of , we label each leaf of a gene tree by the corresponding label of the leaf branch of in which it is contained. To compute , we must consider the contribution of each gene tree compatible with . So that we may write the formulas explicitly, we first consider the contribution of gene trees matching a particular coalescent history,

 pi1…in|h,(S,τ)=∫tp∗i1…in|(G,t)fh((G,t)|(S,τ)) dt

To compute the total probability, we sum over all possible histories so that

 pi1…in|(S,τ) =∑h∈H∫tp∗i1…in|(G,t)fh((G,t)|(S,τ)) dt =∑h∈Hpi1…in|h,(S,τ)

Note that the bounds of integration will depend on the history being considered.

### 2.2. A Modified Coalescent

In this section, we introduce various ways that we might alter the multispecies coalescent to better reflect the evolutionary process. Recall that the length of the path from the root of the species tree to a leaf is the total number of generations that have occurred between the species at the root and that at the leaf. Since the length of a generation may vary for different species [19], it may be desirable to allow the lengths of the paths from the root to each leaf to differ. Therefore, we first consider expanding the allowable set of species tree parameters to include trees that are not equidistant. Such a tree cannot be specified by the speciation times alone, so we use the vector to denote the vector of branch lengths.

Fix a nucleotide substitution model. Let be the set of probability distributions obtained from the standard multispecies coalescent model on the -leaf topological tree . Let denote the set of all distributions obtained by allowing the species tree to vary over all equidistant -leaf trees. If the species tree of the model is not assumed to be equidistant, this removes the assumption of a molecular clock and we refer to this model as the clockless coalescent. The set of distributions obtained from a single tree topology in the clockless coalescent is and the set of distributions obtained by varying the species tree over all rooted -leaf trees is .

We can also account for the fact that the effective population size, , may vary for different species [3] by introducing a separate effective population size parameter for each internal branch of the species tree. We call this model the -coalescent and let be the set of probability distributions obtained from the -coalescent on . In analogy to our notation from above, we use to denote the set of all distributions obtained from the -coalescent and use and for the clockless -coalescent.

Notice that since we assume coalescent events do not occur within leaves, changing the effective population size on the leaf edges does not change the gene tree distribution or the site pattern probability distribution from the gene trees. In this section, we will show that remarkably, the unrooted species tree parameter of the clockless coalescent, -coalescent, and the clockless -coalescent are all generically identifiable. Conveniently, from the perspective of reconstruction, we also show that the method of SVDQuartets [4] can be used to reconstruct the unrooted species tree based on a sample from the probability distribution given by the model.

It is worth pausing for a moment to compare the mathematical and biological meaning of these changes. To illustrate, we will consider how the model on the species tree depicted in Figure 1 changes as we modify the parameters. First, consider the effect of increasing the length of the branch in . This is essentially modeling an increase in the generation rate of the common ancestor of species and . Notice that this changes the gene tree distribution and that this is as it should be. Indeed, thinking in terms of the Wright-Fisher process, an increase in the number of generations along this branch should correspond to an increase in the probability of two lineages coalescing in this branch.

Similarly, we may consider the -coalescent on the equidistant tree . There is now a parameter representing the effective population size of the immediate ancestor of species and . Altering also changes the gene tree distribution, since an increase or decrease of will have the opposite effect on the probability that branches and coalescence along . However, we emphasize that there is no reason to assume that we can obtain the same distribution obtained by altering the length of by adjusting the value of for . When we alter branch length, both the formulas for time to coalescence and for the transition probabilities change, whereas, when we alter the effective population size on a branch, only the formulas for time to coalescence are changed.

Finally, rather than altering branch lengths, one might think to model the increase in the number of generations along a branch by scaling the rate matrix by some factor along . To compute the transition probabilities for a specific branch of a gene tree, we would then use the rate matrices for the branches of the species tree in which it is contained. Thus, multiple transition matrices may have to be multiplied to compute the transition matrix on a single branch of a gene tree. For example, if we scale the rate matrix for , the transition matrix along the leaf edge labeled by in the gene tree depicted in Figure 1 will be . This idea also allows us to consider a multispecies coalescent model with an underlying heterogenous Markov process which we explore in a forthcoming paper.

Since the transition probabilities are given by entries of a matrix exponential, and , it would seem that scaling the rate matrix for a single branch or scaling the length of that branch give the same model. However, while this is true for a nucleotide substitution model, it is not necessarily true for the multispecies coalescent since scaling the rate matrix does not change the gene tree distribution, but changing the length of the branch does. Still, it may be desirable to scale a rate matrix to model an actual change in the per generation mutation rate. However, there is actually no need to introduce scaled rate matrices, since, as Example 2.2 demonstrates, the clockless -coalescent effectively models this situation as well. In fact, we can obtain the same model by fixing any one of the branch length (e.g., set ), effective population size, or rate matrix scaling factor in the species tree and introducing distinct parameters on each branch for the other two. We have chosen to introduce parameters for effective population size and branch length since these have the most natural biological interpretations.

###### Example 2.2.

Let and be two lineages entering a branch of a species tree as in Figure 3. Let be the length of this branch and be the effective population size parameter. The probability that and do not coalesce in is

 (1) 1−∫τ0(12N)exp(−t2N) dt=exp(−τ2N).

If and coalesce, then we can compute the probability of observing the state at and under a homogenous Markov model model where the rate matrix on the branch is scaled by a factor . We assume the distribution of states at the vertex is the vector . Thus, we have,

 (2) pxy=∑z1,z2∫τ0(12N)exp(−t2N)πz1exp(ρQ(τ−t))z1,z2exp(ρQ(t))z2,xexp(ρQ(t))z2,y dt.

Consider now if instead of scaling the rate matrix by , we scale the length of and the effective population size by . Then the probability that and do not coalesce remains unchanged since

 exp(−ρτ2ρN)=exp(−τ2N).

Likewise, the probability of observing state is given by the following formula, where we make the substitution .

 pxy =∑z1,z2∫ρτ0(12ρN)exp(−t2ρN)πz1exp(Q(ρτ−t))z1,z2exp(Q(t))z2,xexp(Q(t))z2,y dt =∑z1,z2∫τ0(12N)exp(−T2N)πz1exp(ρQ(τ−T))z1,z2exp(ρQ(T))z2,xexp(ρQ(T))z2,y dT.

Since the last line is equal to (2), the distribution at the leaves remains unchanged and the probability of coalescence also remains unchanged. Generalizing this example, we can obtain the distribution from a species tree with scaled rate matrices by adjusting population size and branch lengths across the tree.

## 3. Identifiability of the Modified Coalescent

One of the most fundamental concepts in model-based reconstruction is that of identifiability. A model parameter is identifiable if any probability distribution arising from the model uniquely determines the value of that parameter. For the purposes of phylogenetic reconstruction, it is particularly important that the tree parameter of the model be identifiable in order to make consistent inference.

In the following paragraphs, we will use the notation for the set of distributions obtained by varying the -leaf tree parameter in the standard multispecies coalescent model, though the discussion applies equally to , , and . To uniquely recover the unrooted tree parameter of the -leaf multispecies coalescent model we would require that for all -leaf rooted trees and that are topologically distinct when the root vertex of each is suppressed, . This notion of identifiability is unobtainable in most instances and much stronger than is required in practice. Instead we often wish to establish generic identifiability. A model parameter is generically identifiable if the set of parameters from which the original parameter cannot be recovered is a set of Lebesgue measure zero in the parameter space. In our case, although we cannot guarantee that , we will show that that if we select parameters for either model the resulting distribution will lie in with probability zero.

In [5], it was shown that the unrooted topology of the tree parameter for the coalescent model is generically identifiable using the machinery of analytic varieties. Let

 ψS:ΘS↦Δκ4−1

be the map from the continuous parameter space for to the probability simplex with . Label the states of the model by the natural numbers . Given any two species tree parameters and , the strategy is to find a polynomial

 f∈R[qi1…in:1≤i1,…,in≤κ]

such that for all , , but for which there exists such that . If we can show that the function

 f∘ψS2:ΘS2→R

is analytic, the theory of real analytic varieties and the existence of then imply that the zero set of is a set of Lebesgue measure zero [22]. Doing this for all pairs of -leaf trees establishes the generic identifiability of the unrooted tree parameter of .

We will show that the modified models introduced above are still identifiable using the same approach. In the discussion proceeding [5, Corollary 5.5], it was shown that for the standard multispecies coalescent, to establish identifiability of the tree parameter of the coalescent model for trees with any number of leaves it is enough to prove the identifiability of the tree parameter for the 4-leaf model. Essentially the same proof of that theorem applies to the clockless coalescent giving us the following proposition.

###### Proposition 3.1.

If the topology of the unrooted species tree parameter of is generically identifiable then the topology of the unrooted species tree parameter of is generically identifiable for all .

A similar proposition holds for the -coalescent but a slight modification is required. The subtlety is illustrated in Figure 4 where a species tree and its restriction to a 4-leaf subset of the leaves are shown. Notice that on the restricted tree, the effective population size may now vary within a single branch. Therefore, to show the identifiability of the unrooted species tree parameter of the -coalescent for -leaf trees, we must show the identifiability of the unrooted species tree parameter of a model on 4-leaf trees that allows for a finite number of bands on each branch with separate effective population sizes. We will revisit this point after the proof of Theorem 3.5, though it turns out to be rather inconsequential.

### 3.1. The analyticity of ψS

As mentioned above, to apply the results for analytic varieties, we need each of the parameterization functions to be a real analytic function in the continuous parameters of the model. That this is so may seem obvious to some and was stated without proof in [5]. However, the coordinates of are functions of the parameters defined by improper integrals. We note that it is not true in general that if is a real analytic function that is real analytic – even if is defined for all (e.g., if , then ). Further complicating matters is the fact that the entries of the transition matrices are defined by a matrix exponential. Thus, while they can be expressed as convergent power series in the rate matrix parameters, it is difficult to work with improper integral over these series. For the models JC69, K2P, K3P, F81, HKY85, and TN93 these issues become irrelevant, as we can diagonalize the rate matrices and obtain a closed form expression for the entries of the transition matrices. The entries are then seen to be exponential functions of branch length and we can solve the improper integrals from the multispecies coalescent and obtain exact formulas for each coordinate of that are clearly analytic. Thus, we have the following proposition.

###### Proposition 3.2.

Let be a rooted -leaf species tree. The parameterization map is analytic when the underlying nucleotide substitution model is any of JC69, K2P, K3P, F81, HKY85, or TN93.

The rate matrix for the 4-state general time-reversible model is similar to a real symmetric matrix and is thus also diagonalizable. However, actually writing down a closed form for the entries of the transition matrix is not possible due to the large number of computations involved. Instead, we will explain the steps necessary to diagonalize the GTR rate matrix. While we will not be able to write the entries explicitly, we will see that they must be expressions involving only elementary operations of the rate matrix parameters, roots of the rate matrix parameters, and exponential functions. This enables us to prove the following proposition about .

###### Proposition 3.3.

Let be a rooted -leaf species tree. Let be the parameterization map for the multispecies coalescent model when the underlying nucleotide substitution model is the 4-state GTR model. For a generic choice of continuous parameters , there exists a neighborhood around on which is a real analytic function.

###### Proof.

Let be a rate matrix for the 4-state GTR model. Then

 A=diag(π1/2)Qdiag(π−1/2)

is a real symmetric matrix and hence there exists a real orthogonal matrix and a diagonal matrix such that As similar matrices, and have the same eigenvalues and so is an eigenvalue of . Factoring the degree four characteristic equation of , we can use the cubic formula to write the other eigenvalues and explicitly in terms of the rate matrix parameters. All eigenvalues of will be real and less than or equal to zero [1, Section 2.2].

Suppose now that all eigenvalues of are distinct. By the Cayley-Hamilton theorem, we have

 A(A−λ2I)(A−λ3I)(A−λ4I)=0.

Since all eigenvalues of are distinct, the minimal and characteristic polynomials of are equal. Thus, there exists a non-zero column of that is an eigenvector of with eigenvalue . Generically, the first column of each of these matrices will be non-zero and so we can use the normalized first columns to construct the matrix so that .

Therefore,

 P(t)=eQt=(diag(π−1/2)U)diag(1,eλ2t,eλ3t,eλ4t)(UTdiag(π1/2)),

and each entry of the matrix exponential can be written as

 (3) Pij(t)=∑1≤k≤4f(ij)k(q)eλkt,

where the are rational functions of the rate matrix parameters and roots of the rate matrix parameters coming from the cubic formula. Thus, for a generic choice of parameters, the functions are all sums of products of these functions which are exponential in the branch length . The formulas coming from the coalsecent process, , are also exponential functions in . Because each is guaranteed to be less than or equal to zero, when we integrate with respect to branch length each of the integrals converges. Therefore, can be written in closed form as an expression involving rational functions of the model parameters, roots of the model parameters, and exponential functions of both of these.

The expression for given above only fails to be analytic at points where one of the vectors used to construct is actually the zero vector so that the formula for normalizing the columns involves dividing by zero. This includes points where the eigenvalues of the rate matrix are not distinct. Since each of the rate matrices is diagonalizable, if all of the eigenvalues are not distinct then at least one of the vectors we have constructed will certainly be zero [10, §6.4, Theorem 6]. Therefore, the formula given above, which we will denote , agrees with everywhere that it is defined. Our choice of rate matrix parameters from the K3P model in Theorem 3.5 demonstrates that the analytic functions for the entries of the columns of are not identically zero. Therefore, for a generic choice of parameters , agrees with on a neighborhood around and is a real analytic function in that neighborhood. ∎

We can extend the result in Proposition 3.2 to include the generalized -state Jukes-Cantor model since the entries of the transition matrices are given by

 Pij(t)={1κ+κ−1κe−μti=j,1κ−1κe−μti≠j.

For the -state GTR model, because zero is an eigenvalue of , finding all the eigenvalues of amounts to finding the roots of a degree polynomial with coefficients that are polynomials of the rate matrix parameters. Thus, the same arguments also extend to any -state time-reversible model with . Since there is no general algebraic solution for finding the roots of polynomials of degree five or higher (Abel-Ruffini Theorem), for , a different argument is required. We do not pursue these arguments here as the -state model is by far the most commonly used in phylogenetics.

### 3.2. Identifiability of the modified multispecies coalescent for 4-leaf trees

We may encode the probability distribution associated to a -state phylogenetic model on an -leaf tree as an -dimensional tensor where the entry is the probability of observing the state . In [5, Section 4], the authors explain how to construct tensor flattenings according to a split of the species tree. Our first result is the analogue of [5, Theorem 5.1] for the modified coalescent models. We use the notation to denote the probability tensor that results from choosing a species tree with vector of edge lengths and continuous parameters .

###### Theorem 3.4.

Let be a four-taxon symmetric or asymmetric species tree with a cherry . Let be one of the splits , , or . Consider the clockless coalescent when the underlying nucleotide substitution model is any of the following: JC69, K2P, K3P, F81, HKY85, TN93, or GTR. If is a valid split for , then for all ,

 rank(FlatL1|L2(P(S,τ,θ)))≤10.
###### Proof.

Let be the split that is valid for and consider the distribution . Increase the length of one of the leaf edges in the cherry of to create a new vector of edge lengths . Without loss of generality, suppose is extended so that as in Figure 5.

First, we claim that

 rank(FlatL1|L2(P(S,ξ,θ)))≤rank(FlatL1|L2(P(S,τ,θ)))

Notice that since coalescent events do not happen in the leaves of the species tree, the gene tree histories and the formulas for the gene tree distributions from and are identical. The only difference is that the leaf edge labeled by in the gene tree from is longer by than the same edge in the gene tree from . The probability of observing the state from these two gene trees will not be the same, so let us express this probability as when the species tree is and when the species tree is . Extending the branch of a gene tree is equivalent to grafting a new edge onto the leaf edge to create an internal vertex of degree two. To compute the probability of observing a particular state at the leaves of the extended gene tree, we sum over all possible states of this vertex. For clarity of notation, let us represent the matrix of transition probabilities along the grafted edge by . Thus,

 q∗i1i2i3i4|(G,t,θ)=∑1≤j≤4(Mji3)p∗i1i2ji4|(G,t,θ).

Therefore, the total probability for a particular history is given by

 pi1i2i3i4|h,(S,ξ,θ) =∫tq∗i1i2i3i4|(G,t,θ)fh((G,t)|(S,ξ)) dt =∫t(∑1≤j≤4(Mji3)p∗i1i2ji4|(G,t,θ))fh((G,t)|(S,τ)) dt =∑1≤j≤4(Mji3)(∫tp∗i1i2ji4|(G,t,θ)fh((G,t)|(S,τ)) dt) =∑1≤j≤4(Mji3)pi1i2ji4|h,(S,τ,θ)

Summing over all histories, we also obtain

 pi1i2i3i4|(S,ξ,θ)=∑1≤j≤4(Mji3)pi1i2ji4|(S,τ,θ).

Now consider the column of indexed by . The formula above shows that this column is a linear combination of the columns of indexed by and . Therefore,

 rank(FlatL1|L2(P(S,ξ,θ)))≤rank(FlatL1|L2(P(S,τ,θ))).

Given a 4-leaf species tree with a cherry, consider the tree , where all entries of are the same as for except that . The tree may not be equidistant, but it is still clear from the symmetry in the cherry that for any choice of continuous parameters we will have

 pi1i2i3i4|(S,τ0,θ)=pi1i2i4i3|(S,τ0,θ),

which implies that Now we can extend each of the cherry edges and in to construct . Since we have shown that extending these branches does not change the rank of the flattening, the result follows.

###### Theorem 3.5.

Let be a four-taxon symmetric or asymmetric species tree with a cherry . Let be one of the splits , , or . Consider the clockless coalescent when the underlying nucleotide substitution model is any of the following: JC69, K2P, K3P, F81, HKY85, TN93, or GTR. If is not a valid split for , then for generic distributions ,

 rank(FlatL1|L2(P(S,τ,θ)))=16.

We verify the statement computationally by finding a specific choice of parameters for the symmetric tree where and and doing the same for the asymmetric tree. In fact, we can address both cases with one tree by letting be the symmetric tree and setting . There is now only one internal edge and we choose , and . We also need only look at one of the flattening matrices, since as observed in [5, Theorem 5.1], the flattening matrices for the invalid splits will be the same.

We choose continuous parameters from the Jukes-Cantor model by setting the off-diagonal entries of the rate matrix to and letting the root distribution be uniform. With these choices, both flattening matrices for the invalid splits are rank 16. Since the Jukes-Cantor model is contained in JC69, K2P, K3P, F81, HKY85, and TN93, this choice of parameters establishes the result for each of these.

The Jukes-Cantor model is of course also contained in the 4-state GTR model. However, recall from Proposition 3.3 that we were not able to verify that each coordinate of is a real analytic function on the entire parameter space. Instead, in Proposition 3.3, we described the function that agrees with everywhere it is defined, which is to say, almost everywhere. Therefore, we will take a slightly different tack in proving the result for the GTR model. As before, we only need to look at either invalid split, so denote by the determinant of the flattening matrix where the formulas for are used in place of the variables . Since they agree almost everywhere, if is zero on a set of positive measure in the parameter space, the same is true of . Recall that when constructing the coordinate functions of , we normalize the vectors of , and that the resulting coordinate functions fail to be defined only when one of these four vectors is identically zero. Therefore, it is possible to multiply by these vector norms to clear denominators and obtain the function that is analytic on the entire parameter space. Importantly, is zero anywhere that is. Therefore, if is zero on a set of positive measure, so is , and hence so is . But since it is analytic, this would imply that is identically zero.

Therefore, we will establish the result by verifying that is not identically zero. Choosing parameters from the K3P model with , and (Figure 2), the rate matrix has eigenvalues and and the matrix of normalized eigenvectors from Proposition 3.3 is

 U=12⎛⎜ ⎜ ⎜⎝1−11111−1−1111−11−1−11⎞⎟ ⎟ ⎟⎠.

Specifying the other parameters of the model as above, we find that is not identically zero, and since none of the columns of are zero, neither is . ∎

A matrix has rank less than or equal to if and only if all of its minors vanish. Therefore, by Theorem 3.4, if is a 4-leaf tree displaying the split then all of the minors of the flattening matrix for that split vanish on the set . Theorem 3.5 shows that if is a -leaf tree that does not display the split , then one of these minors will not vanish on the set . Hence, as per the discussion at the beginning of Section 3, the set of parameters for mapping into is a set of measure zero. Thus, the unrooted species tree parameter of the clockless coalescent is generically identifiable.

Since coalescent events do not occur in leaves of the species tree, Theorem 3.4 applies equally to the -coalescent and clockless -coalescent models. Following Proposition 3.1, we observed that showing the identifiability of the unrooted species tree parameter of the -coalescent requires proving the same result for 4-leaf trees in a model that allows multiple effective population size parameters on a single edge. Specifically, a model where the maximum number of effective population size parameters on each edge is bounded by the total number of edges in an -leaf binary tree, , is sufficient to prove identifiability of the unrooted species tree parameter in . We have essentially already proven this for all .

Both distributions in the proof of Theorem 3.5 are still contained in the model where we allow multiple effective population size parameters on each edge since we can just choose the multiple parameters on to all be equal to one. We must still verify that the parameterization map for this model is analytic, but the argument from Section 3.1 remains unchanged when we allow multiple effective population size parameters on each edge. Thus, the same choices of parameters from the proof of Theorem 3.5 establish the result for the clockless -coalescent. We also intentionally chose a point corresponding to an equidistant tree so that it applies to the -coalescent. Thus, we have the following corollary.

###### Corollary 3.6.

The unrooted species tree parameters of the clockless coalescent, the -coalescent, and the clockless -coalescent models on an -leaf tree are generically identifiable for all when the underlying nucleotide substitution model is any of the following: JC69, K2P, K3P, F81, HKY85, TN93, or GTR.

### 3.3. Identifiability with invariant sites and gamma distributed rates

It is well known that the rate of evolution may vary across sites [33, 34]. One way to account for this is to let each site evolve according to the same model but where the rate matrix at each site is scaled by a random factor drawn from a specified gamma distribution. If the underlying nucleotide substitution model is assumed to be the GTR model, this is what is known as the GTR+ model.

In practice, the gamma distribution is approximated using rate categories, each with probability . Let the density of the gamma distribution be given by . The rate for each category for is determined by first setting and then computing the points for such that

 ∫ai+1aig(ρ,α,β) dρ=1m.

For , the rate is then given by the mean over each interval,

 ρi=m⋅∫ai+1aiρ⋅g(ρ,α,β) dρ=(mαβ)(I(bβ,α+1)−I(αβ,α+1))

where

 I(z,α)=(1Γ(α))∫z0exp(−x)xα−1 dx.

The last rate is then . We reproduce these formulas from [34] here only to show that the rates can be expressed as analytic functions in the parameters and and consequently that the distributions from the GTR+ model are given by real analytic functions of the parameters.

It is also common to account for invariant sites by using the GTR+I+ model, where is the proportion of sites that do not evolve at all. The multispecies coalescent with the -discrete -state GTR+I+ model was shown to exhibit the same flattening ranks as the multispecies coalescent with the -state GTR model in [5]. This is not terribly surprising as a probability distribution from the former is the sum of distributions each satisfying the same linear relations. Explicitly,

 pΓ+δi1i2i3i4|(S,τ,θ)=(1−δ)m(pρ1i1i2i3i4|(S,τ,θ)+…+pρmi1i2i3i4|(S,τ,θ))+δ(zi1i2i3i4|θ),

where is the probability of observing from the multispecies coalescent model with scaling factor and is the probability of observing this state an an invariant site. If has a cherry as above, then each summand is contained in the linear space defined by the linear relations of the form in the distribution space. The sum satisfies these relations as well, so we have

 rank(Flatab|cd(PΓ+δ(S,τ,θ)))≤(κ+12).

For a non-equidistant tree, the same result no longer applies. If we extend a leaf in the cherry of as we did in Theorem 3.4, we will still have

 rank(Flatab|cd(Pρj(S,τ,θ)))≤(κ+12),

but the particular linear relationships satisfied by the columns of each flattening matrix will depend on the entries of the transition matrix on the extended edge, which in turn depend on the . However, we can obtain an analogous result for JC+I+, where the JC refers to the -state Jukes-Cantor model. When , we prove the result for and , as four is the most common number of categories used in actual phylogenetic applications [14].

###### Theorem 3.7.

Let be a four-taxon symmetric or asymmetric species tree with a cherry . Let be one of the splits , , or . For , consider the -state -discrete JC+I+ model under the coalescent with species tree and .

1. If is a valid split for , then for all ,

 rank(FlatL1|L2(PΓ+I(S,τ,θ)))≤κ2−(κ−12)
2. If is not a valid split for , then for generic distributions ,

 rank(FlatL1|L2(PΓ+I(S,τ,θ)))>κ2−(κ−12)
###### Proof.

Let be the split that is valid for and consider the distribution from the Jukes-Cantor model. Without loss of generality, suppose . Construct the vector with all entries equal to those of but with . Again, by symmetry, we have

 p⋆⋆i3i4|(S,τ,θ)=p⋆⋆i4i3|(S,τ,θ).

As in Theorem 3.4, we will identify the tree as an extension of . For the JC model, there are only two distinct entries of . Let if and otherwise. Therefore, we have

 p⋆⋆i3i4|(S,ξ,θ)=ap⋆⋆i3i4|(S,τ,θ)+∑i3≠jbp⋆⋆ji4|(S,τ,θ),

and one can check that for distinct , the distribution , satisfies

 p⋆⋆k1k2−p⋆⋆k1k3−p⋆⋆k2k1+p⋆⋆k2k3+p⋆⋆k3k1−p⋆⋆k3k2=0.

We obtain such a relation for any 3-element subset of . Moreover, since this linear relation does not depend on or , it is satisfied by , by , and hence, by any distribution from the -discrete JC+I+ model. Consider the relations that come from choosing 3-element subsets of the form . For all , exactly one of these relations involves the variable . Therefore, these relations are linearly independent, and so (1) follows.

In [5, Theorem 5.1] the authors show that for all , when , if is not a valid split for , then

 rank(FlatL1|L2(PΓ+I(S,τ,θ))>κ2−κ.

When , we have

 κ2−κ ≥κ2−(κ−12),

which establishes our result. For , we must produce a choice of parameters to prove that the claim holds for and and for both the symmetric and asymmetric trees. Choosing and the same continuous JC69 parameters from Theorem 3.5 establishes the result. ∎

Since all of the parameterization functions involved are analytic, this is enough to prove the identifiability of the unrooted species tree parameter of the JC+I+ model. Thus, we have the following corollary.

###### Corollary 3.8.

The unrooted species tree parameters of the clockless coalescent, the -coalescent, and the clockless -coalescent models on an -leaf tree are generically identifiable for all when the underlying nucleotide substitution model is the -discrete -state JC+I+ model with and and with and or .

Moreover, the parameters that we used to demonstrate that the invalid flattenings are full rank come from an exponential distribution. Therefore, the same result holds for a model where the rates are constructed from an exponential distribution. In fact, this also applies to a more general variable rates model where the rates are free parameters.

## 4. Conclusions

In the previous section, we have proven that the unrooted species tree parameters of several more generalized versions of the multispecies coalescent model are generically identifiable. Moreover, the means by which we have proven identifiability give us the necessary framework for reconstructing the unrooted species tree from data. In each case, we showed that we can reconstruct the unrooted quartets of the species tree parameter if we know the distribution exactly by taking ranks of the flattening matrices. Specifically, for a 4-state model and generic choices of parameters, we showed that the rank of the flattening matrix for the quartet compatible with the species tree will be less than or equal to while the other two flattening matrices will both be rank 16.

This gives a natural method for inferring the unrooted species tree from actual biological data. Specifically, for each quartet, we infer the unrooted quartets of the species tree by determining which of the three flattening matrices is closest to being rank 10. The method of singular value decomposition from linear algebra already provides a means of determining how close a matrix is to being of a certain rank [8]. This is exactly the procedure used by the method SVDQuartets, which is already fully implemented in the PAUP software [26]. Hence, there is strong theoretical justification for applying SVDQuartets for phylogenetic reconstruction even when effective population sizes vary throughout the tree or when the molecular clock does not hold.

The model presented in Section 2, as well as that presented in Chifman and Kubatko (2015), describes the situation in which gene trees are randomly sampled under the multispecies coalescent model, and then sequence data for a single site evolve along each sampled gene tree according to one of the standard nucleotide substitution models. Data generated in this way have been termed “coalescent independent sites” [31] to distinguish them from SNP data. Although coalescent independent sites and SNP data refer to observations of single sites that are assumed to be conditionally independent samples from the model given the species tree, SNP data are generally biallelic, while coalescent independent sites may include three or four nucleotides at a site, or may be constant. Thus SNP data are a subset of coalescent independent sites data, and the results derived here apply to these data as well.

The other situation in which one might wish to apply these results is to multilocus data. Multilocus data are data in which individual genes are sampled from the species tree under the multispecies coalescent, but for each sampled gene tree, many individual sites are observed. Typical genes observed in phylogenomic studies range from 100 base pairs (bp) to 2,000bp in size, though most are less than 500bp. The site patterns observed within a gene are not independent observations under the model because they share the same gene tree, and thus it is not immediately obvious that the results presented here apply to this case. However, consider the case in which a large sample of genes, say , is obtained, and for each gene, sites are observed. Then, the flattening matrices of site pattern counts constructed from such data will be times the flattening matrix of site pattern counts that would have been observed if only a single site had been observed from each gene tree, which does not change the matrix rank. It is clear that as , the correct theoretical distribution will be well-approximated by the observed site pattern frequencies, and the results presented here will hold. In practice, the genes will vary in their lengths and a more careful argument is required. We have elsewhere carried out thorough simulation studies to show that the methods used in SVDQuartets hold for multilocus data as well as for SNP data and for coalescent independent sites for the original model [4], and we verify these results under the model presented here in comparison with other species tree estimation methods elsewhere [17].

We note two possible criticisms of this method. The first is that, while we showed that generically the flattening matrices for the invalid splits will be rank 16, we have no theoretical guarantees that they are not arbitrarily close to the set of rank 10 matrices. Therefore, we do not know a priori that this method will provide any insight with a finite amount of either simulated or actual biological data. Along the same lines, determining that a flattening matrix is close to the set of rank 10 matrices does not necessarily mean that it is close to the set of distributions arising from a coalescent model, as the latter is properly contained in the former. While both valid considerations, they appear to be academic, as SVDQuartets has already been shown to be an effective reconstruction method on several data sets, both real and simulated [4, 6]. As mentioned above, in a forthcoming paper, we will demonstrate that SVDQuartets also works well in practice by simulating data from these modified coalescent models and applying the method to real biological data sets known to violate the molecular clock.

In recent years, the amount of sequence data available for species tree inference has increased rapidly, presenting significant computational challenges for most model-based species tree inference methods that accommodate the coalescent process. The SVDQuartets method, however, is fully model-based but inference using this method is much more computationally efficient than other methods. This is because, for each quartet considered, all that is required is construction of the three flattening matrices, which involves the simple task of counting site patterns, and computation of singular values from these matrices. In addition, increases in sequence length benefit the performance of the method (because site pattern probabilities are estimated more accurately) with almost no increased computational cost, in direct contrast to other species tree inference methods, such as *BEAST [9] and SNAPP [2]. In the work presented here, we show that the theory underlying the SVDQuartets method holds in much more general settings than originally suggested. In particular, the method can be applied to data that violate the molecular clock and to the case in which each branch of the species has a distinct effective population size. Thus, this work is a significant advance that will contribute meaningfully to the collection of methods available to infer species-level phylogenies from phylogenomic data in very general settings.

## References

• [1] E. S. Allman, C. Ané, and J. A. Rhodes. Identifiability of a Markovian model of molecular evolution with gamma-distributed rates. Advances in Applied Probability, 40:228–249, 2008.
• [2] D. Bryant, R. Bouckaert, J. Felsenstein, N. Rosenberg, and A. Roy Choudhury. Inferring species trees directly from biallelic genetic markers: bypassing gene trees in a full coalescent analysis. Mol. Biol. Evol., 29(8):1917–1932, 2012.
• [3] B. Charlesworth. Effective population size and patterns of molecular evolution and variation. Nature Reviews Genetics, 10:195–205, 2009.
• [4] J. Chifman and L. Kubatko. Quartet inference from SNP data under the coalescent model. Bioinformatics, 30(23):3317–3324, 2014.
• [5] J. Chifman and L. Kubatko. Identifiability of the unrooted species tree topology under the coalescent model with time specific rate variation and invariable sites. Journal of Theoretical Biology, 374:35–47, 2015.
• [6] J. Chou, A. Gupta, S. Yaduvanshi, R. Davidson, M. Nute, S. Mirarab, and T. Warnow. A comparative study of SVDQuartets and other coalescent-based species tree estimation methods. BMC Genomics, 16(Suppl 10):S2, 2015.
• [7] J. Felsenstein. Evolutionary trees from DNA sequences: A maximum likelihood approach. J. Mol. Evol., 17(6):368–76, 1981.
• [8] G. H. Golub and C. F. V. Loan. Matrix Computation. Johns Hopkins University Press, 4th edition, 2013. Section 2.4.
• [9] J. Heled and A. J. Drummond. Bayesian inference of species trees from multilocus data. Mol. Biol. Evol., 27(3):570–580, 2010.
• [10] K. Hoffman and R. Kunze. Linear Algebra. Prentice Hall, second edition, 1971.
• [11] J. F. C. Kingman. Exchangeability and the evolution of large populations. Pp. 97-112 in G. Koch and F. Spizzichino, eds. Exchangeability in probability and statistics. North-Holland: Amsterdam., 1982.
• [12] J. F. C. Kingman. On the genealogy of large populations. J. Appl. Prob., 19A:27–43, 1982.
• [13] J. F. C. Kingman. The coalescent. Stoch. Proc. Appl., 13:235–248, 1982.
• [14] P. Lio and N. Goldman. Models of molecular evolution and phylogeny. Genome Research, 8:1233–1244, 1998.
• [15] L. Liu, L. Yu, and S. Edwards. A maximum pseudo-likelihood approach for estimating species trees under the coalescent model. BMC Evol. Biol., 10(302), 2010.
• [16] L. Liu, L. Yu, D. Pearl, and S. Edwards. Estimating species phylogenies using coalescence times among sequences. Syst. Biol., 58(5):468–477, 2009.
• [17] C. Long and L. Kubatko. Extension of SVDQuartets to general multispecies coalescent models, and comparison with other coalescent-based species tree estimation methods. in preparation, 2017.
• [18] W. P. Maddison. Gene trees in species trees. Syst. Biol., 46:523–536, 1997.
• [19] A. P. Martin and S. R. Palumbi. Body size, metabolic rate, generation time, and the molecular clock. Proceedings of the National Academy of Sciences USA, 90:4087–4091, 1993.
• [20] S. Mirarab, R. Reaz, M. D. Bayzid, T. Zimmermann, M. S. Swenson, and T. Warnow. Astral: Genome-scale coalescent-based species tree. Bioinformatics (ECCB special issue), 30(17):i541–i548, 2014.
• [21] S. Mirarab and T. Warnow. Astral-ii: Coalescent-based species tree estimation with many hundreds of taxa and thousands of genes. Bioinformatics (ISMB special issue), 31(12):i44–i52, 2015.
• [22] B. Mityagin. The zero set of a real analytic function. arXiv:1512.07276, 2015.
• [23] P. Pamilo and M. Nei. Relationships between gene trees and species trees. Mol. Biol. Evol., 5(5):568–583, 1988.
• [24] C. Semple and M. Steel. Phylogenetics. Oxford University Press, Oxford, 2003.
• [25] D. Swofford. PAUP. Phylogenetic Analysis Using Parsimony (and Other Methods). Version 4. Sinauer Associates, Sunderland, Massachusetts, 2002.
• [26] D. Swofford. PAUP. Phylogenetic Analysis Using Parsimony (and Other Methods). Version 4a150. 2016.
• [27] M. Syvanen. Horizontal gene transfer: evidence and possible consequences. Annu. Rev. Genet., 28:237–261, 1994.
• [28] F. Tajima. Evolutionary relationship of DNA sequences in finite populations. Genetics, 105:437–460, 1983.
• [29] N. Takahata and M. Nei. Gene genealogy and variance of interpopulational nucleotide differences. Genetics, 110:325–344, June 1985.
• [30] S. Tavaré. Line-of-descent and genealogical processes, and their applications in population genetics models. Theor. Popul. Biol., 26:119–164, 1984.