Species tree estimation using ASTRAL: how many genes are
enough?
Shubhanshu Shekhar1 , Sebastien Roch2 , and Siavash Mirarab∗1

arXiv:1704.06831v2 [q-bio.PE] 16 Sep 2017

1

Department of Electrical and Computer Engineering, University of California,
San Diego
2
Department of Mathematics, University of Wisconsin-Madison, Wisconsin
September 19, 2017

Abstract
Species tree reconstruction from genomic data is increasingly performed using methods
that account for sources of gene tree discordance such as incomplete lineage sorting. One
popular method for reconstructing species trees from unrooted gene tree topologies is ASTRAL. In this paper, we derive theoretical sample complexity results for the number of genes
required by ASTRAL to guarantee reconstruction of the correct species tree with high probability. We also validate those theoretical bounds in a simulation study. Our results indicate
that ASTRAL requires O(f −2 log n) gene trees to reconstruct the species tree correctly with
high probability where n is the number of species and f is the length of the shortest branch
in the species tree. Our simulations, some under the anomaly zone, show trends consistent
with the theoretical bounds and also provide some practical insights on the conditions where
ASTRAL works well.

Index terms— Phylogenetics, Species tree estimation, Incomplete lineage sorting, Sample
complexity, ASTRAL.

1

Introduction

Evolutionary relationships between organisms are typically represented in the form of a phylogeny
known as a species tree. Each branching point in a species tree represents a historical speciation
event. On the other hand, phylogenies of individual parts of genomes of extant species do not
always follow the same patterns as the speciation events. This discordance can be due to various
biological processes that include incomplete lineage sorting (ILS), duplication and loss, horizontal
gene transfer, and hybridization [34, 16, 8]. These processes can impact not only functional genes
but any region of the genome; however, following the convention, we use the term gene tree to refer
to a tree specific to a single region of the genome. The potential discordance between gene trees
and species tree has motivated a growing set of methods for inferring species trees while accounting
for these differences [20, 31, 35, 36, 30, 23, 12, 32, 37, 13, 2, 1]. Most existing techniques focus on
one of the biological causes of discordance (but there are exceptions [50, 9, 41, 54]). In particular
ILS, a ubiquitous cause of discordance found in many datasets [3, 39, 48], has been the subject of
much interest.
ILS is best understood by considering the interplay between phylogenetics and the population
genetic coalescent process [25]. The most widely used model to study ILS is the multispecies
∗ smirarab@ucsd.edu

Accepted for oral presentation at RECOMB 2017
Under review at IEEE TCBB (submitted on 12/26/2016)

1

Table 1: Versions of ASTRAL, their constraint set (X), and time complexity with n (the number
of species), m (the number of genes), and |X|.
Version

Constrain set (X)

Time

ASTRAL*

Unconstrained (i.e., exact)

O(nm3n )

ASTRAL-I

Gene tree bipartitions

O(n4 m3 )

ASTRAL-II

ASTRAL-I+heuristic additions

O(nm|X|2 )

coalescent model (MSC) [38, 40]. Under the MSC model, to generate a gene tree, the Kingman’s
coalescent process [25] is followed on each species tree branch, which represents a population and is
often assumed to have a fixed population size; all lineages that reach the top of a branch are copied
to the ancestral population, which follows its own coalescent process. When two lineages from two
different species fail to coalesce in their lowest common ancestral population, they both go back
to an earlier branch where alleles from other species are also present; there, the two lineages may
first coalesce with lineages from those other species before coalescing with each other. When this
happens, the resulting lineage tree differs from the species tree, and is said to experience ILS.
When only concerned with gene tree topologies, species tree branch lengths can be specified
in coalescent units (the number of generations divided by the effective population size) [16]. A
species tree with branch lengths in coalescent units uniquely defines a distribution on gene tree
topologies [17] under the MSC model.
Species tree estimation despite ILS is possible using an array of methods that use various
strategies, including co-estimation of gene trees and species tree [29, 23, 27], species tree estimation without gene trees [12, 10], and two-step approaches, which first estimate all gene trees
independently and then combine them using a summary method (e.g., [37, 13, 32, 31]). Since gene
tree inference typically produces unrooted trees, summary methods that operate on unrooted trees
are more useful. An early unrooted method was NJst [30], which can be considered [5] the unrooted version of STAR, and is implemented efficiently in a tool called ASTRID [52]. A newer
suite of methods called ASTRAL [35, 36] has also been designed and has been widely used on
biological datasets (e.g., [53, 28, 22, 44, 6, 3, 11, 45, 19, 55, 24]).
The ASTRAL suite of tools seek to solve the Maximum Quartet Support Species Tree (MQSST)
problem. Each quartet of leaves can have one of three unrooted tree topologies and, for each
quartet selected from the full set of leaves, the species tree will induce one of the three topologies.
The quartet score of a species tree with respect to a set of gene trees is the sum of the number of
induced quartet topologies shared between the species tree and each gene tree (see Appendix A.1
for a formal definition). The MQSST problem is to find the species tree with the maximum quartet
score with respect to a set of input gene trees. The problem is NP-hard [26], and this has led
to the development of several versions of ASTRAL (Table 1). The exact version of ASTRAL,
which we call ASTRAL*, uses dynamic programming to solve the problem exactly in exponential
time. ASTRAL* has limited scalability, motivating the definition of a constrained version of
the MQSST problem that restricts the species tree to draw its branches from a predefined set
of bipartitions. In ASTRAL-I, the constraint set is the set of bipartitions in the input gene
trees [35]. ASTRAL-II further expands the constraint set heuristically using several complicated
and empirically motivated strategies [36] that evade straightforward mathematical formulation.
ASTRAL*, ASTRAL-I, and ASTRAL-II are all statistically consistent estimators of the species
tree given true gene trees [35]. The proof of consistency follows from a result by Allman et al. that
shows that, for any quartet of leaves, the species tree topology has no less than a 1/3 probability
of matching each gene tree [4]; thus, the most likely quartet gene tree matches the species tree.
However, this result does not extend to more than four species. In fact, for five species or more,
there always exist species trees where the probability of an unrooted gene tree matching the
unrooted species tree is less than some other tree topology [14]. A species tree that does not
match the most likely gene tree is said to be in the “anomaly zone” [15, 14, 43].
Despite widespread use of the ASTRAL suite and its high accuracy in simulations [35, 36, 52,

2

46], little is understood about its sample complexity. In this context, the sample complexity is
simply the number of genes (m) required to guarantee correct species tree reconstruction with
high probability. Sample complexity is typically established asymptotically and with respect to
the number of species (n) or the length of the shortest branch in the species tree (f ).
Under the MSC model, shorter branches in coalescent units result in increased gene tree discordance and consequently increase the difficulty of reconstructing the species tree. Sample complexity results have been established for some methods that are not in wide use. For instance, for
sufficiently small f , the GLASS and STEAC algorithms, which use both gene tree branch length
and topology, need the number of true gene trees to scale linearly with f −1 and f −2 , respectively [42]. For summary methods, no better bound has yet been found. More specifically, for
algorithms that only use gene tree topology, the best demonstrated data requirement is f −2 . It
is worth noting that these results are all assuming that true gene trees are known. We revisit the
issue of gene tree error and review known results for methods that do not rely on input gene trees
in the discussion section.
The only previous sample complexity results specifically relevant to ASTRAL-II is the bipartition cover results established by Uricchio et al. [51]. They establish loose upper bounds on the
required number of genes such that each bipartition in the species tree is observed in at least one
of the gene trees with high probability. While this result does not help in analyzing the data
requirement of ASTRAL*, it can be used to limit the complexity of ASTRAL-I and ASTRAL-II
once the complexity of ASTRAL* is established.
In this paper, we establish the first theoretical bounds for the sample complexity of ASTRAL*.
We show that the sufficient number of true gene trees required by ASTRAL* to reconstruct the
correct species tree with high probability grows quadratically with the inverse of the shortest
branch length in the species tree and logarithmically with the number of species. We further show
that the necessary number of genes for ASTRAL* to reconstruct the species tree has a similar
asymptotic behavior. We then test the data requirements of ASTRAL-II in simulations and show
that its empirical data requirements match the theoretical results for ASTRAL*. Our simulation
study focuses on datasets with extremely short branches, including some conditions that fall under
the anomaly zone, and the results indicate that ASTRAL-II can be made arbitrarily accurate given
enough true gene trees even under very difficult conditions. We also test the data requirement of
NJst [30], as implemented in the ASTRID [52] software package, in simulations and show that the
relative performance of ASTRAL-II and NJst depends on the shape of the species tree.

2

Theoretical bounds

We first present theoretical bounds on the number of gene trees required by the exact ASTRAL*
algorithm to reconstruct the true species tree with high probability. We present the main results
here and leave the proofs for the appendix.
Let n be the number of leaves in the species tree, and let m denote the number of gene trees,
all generated from the species tree according to the MSC model. We use f to denote the length of
the shortest species tree branch (as an unrooted tree). For each quartet of leaves, there are three
possible unrooted tree topologies that can be induced by each gene tree. Under the MSC model,
we know that for each quartet qi , the probability that a gene tree will be congruent with the species
tree is pi = 1 − 23 e−di , and the other two configurations have equal probability (ri = (1 − pi )/2)
of occurring; here, di gives the length of the path between the two endpoints of the middle edge
of the quartet topology in the species tree, measured in coalescent units [4]. ASTRAL* returns
a fully binary species tree which shares the maximum number of induced quartets with the input
gene trees [35] (Appendix A.1). We say ASTRAL* has an error when there exists a bipartition in
the reconstructed tree that does not appear in the true species tree.

3

2.1

Upper bounds

We first ask how many gene trees are sufficient for ASTRAL* to reconstruct the true species tree
with high probability.
√
Theorem 2.1. Consider a model species tree with minimum branch length f < log( 2). Then,
for any  > 0, ASTRAL* returns the true species tree with probability at least 1 −  if the number
of input error-free gene trees satisfies
 n 
4 4
9
1
.
(1)
m > log
2

(1 − e−f )2
−f

Remark. In the limit of small f , we can use the fact that limf ↓0 1−ef ↑ 1, and thus for f small
√
−f
enough we have 1−ef ≥ 0.9, and so a sufficient condition on the number of gene trees required
is
 
n 1
.
(2)
m > 20 log
6 f 2
That is, for small values of f , m = Ω(f −2 log n) gene trees are sufficient for ASTRAL* to return
the correct species tree with an arbitrarily high probability.
Our proofs, shown in Appendix A.2, are based on the observation that quartet scores follow
a multinomial distribution and therefore, for a large number of genes, we expect the frequencies
to concentrate tightly around their means. We use Hoeffding’s inequality to get the required
concentration.

2.2

Lower bounds

We now show that there exist species trees for which Ω(f −2 log n) gene trees are also necessary.
2.2.1

Quartet species trees (n = 4)

It is simpler to first analyze the case of n = 4, i.e., reconstructing a single quartet. Let Y =
{q1 , q2 , q3 } be the set of three possible topologies, let {n11 , n12 , n13 } be their frequencies out of m
gene trees and let Q = {p, r, r} be the corresponding probability distribution on Y according to the
MSC model. Let E denote the event that ASTRAL* returns the wrong species tree. Observe that,
for n = 4, the error event of ASTRAL* is equivalent to n11 < max{n12 , n13 }. Let F denote the
event {n11 < m/3}, which requires max(n12 , n13 ) > m/3. ASTRAL* returns the tree associated
with the highest frequency; thus, P (E) ≥ P (F ).
Before proving the lower bounds, we first present a simple numerical experiment. Suppose we
want to find the m below which the error probability is at least  for any  ∈ [0, 0.5). We know
that P (F ) is a lower bound on P (E), and that it has the following closed form expression
bm/3c 

P (F ) =

X
i=0


m i
p (1 − p)m−i ,
i

(3)

which is simply the Cumulative Distribution Function (CDF) of a binomial calculated at bm/3c,
where p = 1 − 32 e−f . For given f and , we can compute the value of P (F ) for various choices of
m. In Figure 1, we show the largest m for which the computed P (F ) is larger than . As f −2
increases, the required number of binomial trials grows almost linearly with f −2 , with a slope
that increases with decreasing , as expected. This experiment suggests that for m = O(f −2 ),
the reconstruction error for one quartet can be made arbitrarily close to 1/2. Our next result
formalizes this intuition.
Theorem 2.2. For a species tree with n = 4 leaves, let f be the length of the single internal branch
in the unrooted species tree. Then, for any  ∈ [0, 0.5), there exists an f0 > 0 and a constant d ,
such that for all m ≤ d /f 2 with f ≤ f0 , the probability that ASTRAL* outputs the wrong quartet
is at least .
4

ε

140000

●

●

120000

●
●
●

100000

# genes

●

●

0.05
0.1
0.2
0.33
0.49

●
●
●
●
●
●
●

●

80000

●
●

●
●

●

●
●

●

60000

●
●

●
●

●

●
●

●

40000

●
●

●
●
●
●
●

●
●
●
●

●
●

●
●
●
●
●

●
●
●

●
●
●

●
●
●

●
●
●

●
●
●

●

●
●

●

●

●

●

●

●

●

●

●

●

●

●

●

●

●

●

●

●

●

●

●

●
●

●

●

●

●

●

●

●

●

●

●

50

15

10

10

0

(0

.
00 100
(0 00
00 .01 )
4
0
(0 14)
00 .01
00
0
0)
20 (0.
00 00
81
0
6)
25 (0.
00 00
7
0
07
30 (0.
)
00 00
63
0
2)
35 (0.
00 00
57
0
7)
40 (0.
00 00
53
0
5)
45 (0.
00 00
50
0
0)
50 (0.
00 00
47
0
1)
55 (0.
00 00
44
0
7)
60 (0.
00 00
42
0
6)
65 (0.
00 00
40
0
8)
70 (0.
00 00
39
0
2)
75 (0.
00 00
37
0
8)
80 (0.
00 00
36
0
5)
85 (0.
00 00
35
0
4)
90 (0.
00 00
34
0
3)
95 (0.
00 00
33
0
3
1e (0.
+0 00 )
32
5
(0
4)
.0
03
16
)

●

●
●

●
●

20000

1 f2 (f)

Figure 1: The number of genes below which a certain level of error () is unavoidable according
to the Binomial distribution. Values on the x axis show 1/f 2 where f is the branch length in
coalescent units (f is shown parenthetically). The y axis shows the largest value of m where the
CDF of a binomial with parameter 1 − 2/3e−f evaluated at m/3 is above the error level . A line
is fitted to points for each  (shown for each line).
This theorem shows that with only four species and m = O(f −2 ) genes, the probability of
error for ASTRAL* can be made arbitrarily close to 1/2. Thus, the lower bound and the upper
bound asymptotically match as f decreases.
The proof, shown in Appendix A.3, essentially uses the Berry-Esséen theorem to first approximate the binomial with a Gaussian distribution and then shows that for sufficiently small f , when
the number of observations is below O(f −2 ), the Gaussian has a sufficiently large probability of
deviating from its mean enough to become below 1/3.
2.2.2

Species trees of arbitrary size (n > 4)

If all quartets in a gene tree were generated independently, extending Theorem 2.2 to more species
would be trivial. However quartets have a complicated dependence structure. To establish a
complement to our upper bound, it suffices to show that there exists a tree where m = O(f −2 log n)
gene trees would result in a high probability of error. We start with a simpler claim that establishes
the lower bound with respect to f —but not n.
Claim 2.1. For any n and  ∈ [0, 0.5), there exists a species tree with n leaves and shortest branch
length f such that when ASTRAL* is used with m ≤ d0 /f 2 gene trees, for some constant d0 , the
probability that ASTRAL* reconstructs the wrong tree is at least .
The proof, detailed in Appendix A.4, simply uses trees of the form (((u, v), w), X ), where X
is an arbitrary tree with (n − 3) leaves (represented by the set X = [n − 3]) that is connected to
the root with a very long branch. By making this branch sufficiently long, we ensure that with
high probability all quartets with leaves {x, u, v, w} for x ∈ X have the same frequencies in the
gene trees; this observation enables us to reduce the problem to the case of n = 4 and make use
of Theorem 2.2.
By construction, our result in Claim 2.1 does not depend on n. To match the lower bound, we
strengthen this result to an upper bound that depends on both f and n.
Theorem 2.3. For any ρ ∈ (0, 1) and a ∈ (0, 1), there exist constants f0 and n0 such that the
following holds. For all n ≥ n0 and f ≤ f0 , there exists a species tree with n leaves and shortest
n
branch length f such that when ASTRAL* is used with m ≤ a5 log
f 2 gene trees, the event E that
ASTRAL* reconstructs the wrong tree has probability
P (E) ≥ 1 − ρ.

5

(4)

Remark. The above result tells us that, for all n large enough, there exists a species tree such
log n
gene trees, the probability of error is bounded away
that when ASTRAL* is used with m ≤ a 5f
2
from zero. This lower bound can be made arbitrarily close to 1 (as n goes to infinity).
The proof (Appendix A.5) uses ideas similar to Claim 2.1; instead of using one construct of
the form ((u, v), w) matched with x ∈ X (as we did for Claim 2.1), we use n/3 such constructs
(each called a triplet). We show that each triplet, when matched with one more leaf to build a
quartet, has a probability of error of at least Ω( n1a ). We then argue that when branches above all
triplets are long enough, error events for individual triplets can be made independent conditioned
on each triplet coalescing in its root branch. Finally, the conditional independence of Ω(n) error
events and the Ω( n1a ) lower bounds are used to derive the bound.

2.3

Summary of results

We showed that the upper bound on the required number of genes for ASTRAL* to reconstruct the
correct species tree with high probability grows as f −2 log n. We further proved that there exist
parts of the species tree space where this upper bound is tight in the sense that the probability
of error can be bounded strictly away from zero when the number of input gene trees is of order
f −2 log n. Thus, ASTRAL* requires Ω(f −2 log n) gene trees to universally guarantee correct
species tree reconstruction with high probability.

3

Simulations

We now study the performance of ASTRAL-II and the NJst approach [30] in simulations. We
seek to find the number of genes required by each method to recover the true unrooted species
tree with high probability. We test ASTRAL-II, version 4.10.10 [36] and NJst as implemented in
ASTRID [52].

3.1

Simulation procedures

Testing sample complexity in simulations requires care. We seek the minimum number of genes
with which methods stay under a specific level of error for specific model trees.
3.1.1

Model species trees

We study three different tree forms, shown in Figure 2(a), all with exactly n = 8 species and
varying branch lengths (in coalescent units). The three forms we test are a fully unbalanced
tree where all internal branches have length f , a fully balanced tree where all internal unrooted
branches have length f , and a balanced tree where the two branches incident on the root have
length 30 (i.e., very long) and the remaining branches have length f . We refer to these trees
as caterpillar, balanced, and double-quartet, respectively. We use nine values for f in the range
[0.005, 0.1], spaced so that 1/f 2 is divided into equal chunks.
The unrooted anomaly zone requires at least two consecutive short branches (often defined as
below 0.1) [14]. Based on calculations carried by the hybrid-coal software [56],the caterpillar tree
is in the anomaly zone, but the balanced tree and the double quartet trees are not (Table 2). More
importantly, most of the branch lengths we have used are extremely short and go beyond what one
might expect in most real biological datasets. With f = 0.1, assuming a diploid population size
of 100,000, the branch would represent 20,000 generations (or roughly 100,000 years assuming a
generation time of 5 years). In our simulations with f = 0.1, the percent of quartets that agree with
each branch (computed by ASTRAL-II [47]) is never more than 46%. Thus, f = 0.1 represents
short branches, but is within ranges observed in biological datasets. On the other end, using the
same assumptions, f = 0.005 corresponds to a branch of only 1,000 generations, which is perhaps
unrealistically short, and close to the range where biologists would consider a hard polytomy. In
our simulations, all branches except the one long internal branch in the double-quartet tree (note
6

Table 2: Anomaly zone Conditions. For each each species tree topology and the smallest and the
largest value of f , we show the rank of the species tree topology among the 10395 unrooted gene
trees, its probability, and the probability of the most likely gene tree other than the species tree.
Tree

f

ST rank

p(S) × 104

p(g t ) × 104
13.671

Caterpillar

0.1

50

8.708

Caterpillar

0.005

6028

0.698

2.420

Balanced

0.1

1

28.096

23.170

Balanced

0.005

1

2.5694

2.5690

Double-quartet

0.1

1

247.843

143.215

Double-quartet

0.005

1

128.457

124.679

that in the unrooted version, the double-quartet tree has only one long branch) result in high
levels of discordance. For f < 0.1, the quartet support for species tree branches is very close to
1/3 (i.e., random) and is never more than 36%.
Thus, we emphasize that our choice of branch lengths is motivated by a desire to explore the
ability to reconstruct the species tree under extreme conditions and establishing trends in data
requirement; we refer the reader to earlier publications for simulations seeking to emulate realistic
conditions [35, 36].
For each species tree, gene trees are generated according to the MSC model using the Dendropy
package [49].
3.1.2

Number of replicates and 

We seek to find the smallest number of gene trees, m, so that the probability of incorrect tree
reconstruction by each method is no more than  = 0.1. Computing the exact  for a given m in
simulations requires an infinite number of replicates. In our simulations, for each species tree, we
simulate 401 replicate datasets, with varying number of genes. To approximately achieve  = 0.1,
we search for the smallest m such that in no more than 40 replicates each method outputs an
incorrect species tree. Using normal approximation for binomial, observing 40 error events out of
401 replicates gives us a 90% confidence interval of (0.075, 0.124) for . Getting 90% confidence
intervals of 0.1 ± 0.01 would require close to 2,500 replicates, which is computationally infeasible
given the number of genes required per replicate.
3.1.3

Binary search for m

Finding m requires running ASTRAL-II and NJst with increasing numbers of genes until reaching
error levels below the specified threshold. Since our choices of f require tens of thousands of gene
trees, this approach would become computationally impractical. Instead, we use a binary search
to find a tight window around the true value of m. The binary search starts with an initial range
for m, picked as a guess by dividing results presented in Theorem 2.1 by 20. We set m to the
midpoint of the range and estimate  by counting the number of replicates where each method
outputs an incorrect tree. Based on the estimated , we then decide to narrow down the range
to the portion above or below the midpoint, and we repeat the process. We stop when the range
becomes narrow enough (1/20 of our initial guess). We depict final ranges and fit lines to their
midpoints. Additionally, we check that upper and lower halves have been each chosen at least in
one iteration of the binary search. This is to ensure that the true m is within our initial range,
and in occasions where this was not true, we expanded the initial range.

3.2

Results

Results of our simulations are shown in Figure 2. As expected, decreasing f results in increased
data requirements for both methods. The average number of genes required by ASTRAL-II ranges

7

f
f/2

f

f/2

f
f

f

f

f

f

30

f

H

F

G

E

D

Caterpillar

C

30

f

f

B

A

H

G

F

E

D

Balanced

B

C

A

H

G

f

F

f

E

D

C

Double-quartet

f

B

A

(a)
Astral

Astrid

1e+05

Balanced
Caterpillar
Double−quartet

9e+04
8e+04

# genes (m)

7e+04
6e+04
5e+04
4e+04
3e+04
2e+04
1e+04
0e+00

100

0)
6)
7)
2)
7)
5)
0)
0)
6)
7)
2)
7)
5)
0)
4)
4)
0)
0)
100
081
070
063
057
053
050
100
081
070
063
057
053
050
141
141
000
000
(0.1 00 (0.0 00 (0.0 00 (0.0 00 (0.0 00 (0.0 00 (0.0 00 (0.0 00 (0.0 00 (0.1 00 (0.0 00 (0.0 00 (0.0 00 (0.0 00 (0.0 00 (0.0 00 (0.0 00 (0.0
1
50
50
100
150
200
250
300
350
400
100
150
200
250
300
350
400

1 f2 (f)

ASTRAL/ASTRID

(b)
1.2
1.1
1.0
0.9
0.8

100

0)

000

(0.1

.014

0 (0

500

00)

14)
00 (0

100

.010

00 (0

150

.008

16)

07)

32)

.007

00 (0

200

.006

00 (0

250

1 f2 (f)

300

77)
35)
00)
.005
.005
.005
00 (0
00 (0
00 (0
350
400

(c)
Figure 2: Data requirement of ASTRAL-II and NJst (implemented in ASTRID) in simulations
with  = 0.1. For each of the three different species tree shapes (panel a) and values of f , the
length of the shortest branch, 401 replicate datasets are simulated using the MSC model, each with
up to 105 gene trees. (b) A binary search is used to find an approximate range for the smallest
number of genes with which ASTRAL-II (left) or NJst (right) recover the correct tree in at least
90% of the 401 replicates. Boxes show these ranges and the horizontal lines show the midpoint
of each range. A line is fitted to midpoints for each tree shape. (c) The ratio of the number of
genes required by ASTRAL-II to the number required by NJst; below 1.0 indicates ASTRAL-II
is better.
from 206, 255, or 297 genes, respectively for the balanced, caterpillar, or double-quartet trees with
f = 0.1, to 48948, 59528, and 93750 genes with f = 0.005. Matching our theoretical results, the
number of genes required by ASTRAL-II increased proportionally to 1/f 2 . Similarly, NJst seems
to have data requirements that scale with 1/f 2 under these model conditions. Interestingly, the
tree shape had a substantial and somewhat surprising impact on the exact number of required
genes.
With both methods, the caterpillar tree required slightly more genes than the balanced tree,
a result that we don’t find surprising as we will discuss. However, counter-intuitively, when
the first two branches in the balanced tree were made extremely large (i.e., the double-quartet
tree), the data requirement for both methods went up, and the increase was more dramatic for
ASTRAL-II. For example, for f = 0.00816, the balanced tree requires 21,096 gene trees, while
the double-quartet requires 36,115, an increase of about 70%. Comparing the data requirement of
ASTRAL-II and NJst reveals that ASTRAL-II typically required between 10 to 20% fewer genes
for the caterpillar and balanced trees, while NJst required fewer genes for the double quartet tree.
The results show that the data requirement of ASTRAL-II and its relative performance depend on the tree shape as well as on the exact pattern of branch lengths. We provide possible
explanations for the surprising results in the Discussion section. The results also suggest that the
sample complexity of NJst might be also quadratic with respect to 1/f .

8

4

Discussion

The results we presented heavily rely on properties of the MSC model, which assumes that genes
evolve in an i.i.d. fashion. This assumption, which is perhaps biologically unrealistic, enables us
to treat true gene trees as draws from a multinomial distribution over the discrete space of all tree
topologies, and hence provide a nice mathematical framework for studying sample complexity. The
main remaining difficulty in analyzing the MSC process is the fact the tree space is combinatorial
and hard to enumerate. Luckily, using quartets provides a way to address this challenge because
there are only three possibilities for an unrooted tree topology, and moreover, the anomaly zone
does not exist. The analysis of ASTRAL* sample complexity for quartets mostly reduces to analyzing the concentration of a multinomial distribution with three probabilities (1− 32 e−f , 13 e−f , 13 e−f );
thus, we can use the well-established Hoeffding concentration inequality and the standard union
bound technique to deal with dependencies for upper bounds. In the case of lower-bounds, even
though the quartet case reduces to analyzing the CDF of binomials (Fig. 1), the binomial CDF
function has a complex dependence on the parameter m, which motivates us to use the BerryEsséen theorem (which essentially is a stronger version of the Central Limit Theorem) and also
more strong bounds specific to binomials. Extending lower bound results from quartets to arbitrary n is not quite as easy as the upper bounds because simple methods like union bounds are
not available. Instead, we exploit the properties of the MSC model to design species trees with
sufficiently long branches, and this makes the quartets independent conditioned on events of an
arbitrarily high probability. To summarize, the key to obtaining our results is to reduce large n
cases to n = 4, and to then turn the problem into the analysis of binomial distributions defined
by MSC model.

4.1

Gene tree error

Throughout the paper, we assumed a random sample of true gene trees generated by the MSC
model is given. However, in practice, the gene trees are reconstructed from sequence data, and
hence might contain errors. Moreover, processes other than ILS may contribute to gene tree
discordance. These violations of assumptions have the potential to invalidate our results. One
may suspect that when violations are not severe, many of our conclusions may remain intact. We
can address this conjecture using a simplistic error model.
Assume that each gene tree quartet can be incorrectly estimated with probability not exceeding
some constant e0 , and also assume these events are independent across different quartets. In the
event of an error, each of the two alternative topologies is assumed to be equally likely (i.e., has 1/2
conditional probability). This model is not a realistic model of gene tree reconstruction error for a
variety of reasons. Quartets and their error events are not independent. Moreover, the assumption
of having unbiased error is very strong, and does not necessarily follow from properties of tree
reconstruction under any of the realistic models of sequence evolution. Nevertheless, the following
result provides insights about robustness of our results.
Let us denote by Qi the probability distribution according to the MSC model for the quartet i,
i.e., Qi = {pi , ri , ri } is the probability distribution on the set of the three possible topologies of the
ith quartet. Then under the above error assumption, the observed quartets will have a probability
0
0
0
0
0
distribution Q0i = {pi , ri , ri }, with pi ≥ pi (1 − e0 ) + ri e0 , and ri ≤ ri (1 − e0 /2) + pi e0 /2. As long
0
0
as pi > ri , the upper bound for m will have the same form as in Eq.5, but with larger constants.
A sufficient condition for this to happen is pi (1 − e0 ) + ri e0 > 1/3 or e0 < 2/3. Thus, with this
particular error model, up to 66% of gene trees can be incorrectly estimated and our sufficient
conditions for correct reconstruction using ASTRAL* will not asymptotically change, indicating
some level of robustness.
A rigorous analysis of the data requirement in the presence of gene tree error requires considering gene tree inference requirements. Since a sufficient condition for accurate gene tree inference
is having genes with length that scales as f −2 [21], one may expect that m and the length of
each gene both growing like f −2 will be sufficient to bound error. This will result in total data
requirement that scales with f −4 . The METAL [13] method that estimates species tree directly
9

from sequence data has f −2 overall data requirement (using short constant-length genes), but it
has remained theoretical and its performance on data remains untested.

4.2

ASTRAL-I and ASTRAL-II

The constrained versions, ASTRAL-I and ASTRAL-II (see Table 1), proceed similarly to ASTRAL*, with the added constraint that they restrict the species tree to draw its bipartitions
from a fixed set (X). ASTRAL-I sets X to the input set of gene trees, and ASTRAL-II further
expands the set of allowable bipartitions using various heuristics. Clearly, any sufficient condition
for ASTRAL-I to recover the true species tree will also be sufficient for ASTRAL-II, but necessary
conditions for ASTRAL-II need not be as restrictive as ASTRAL-I.
Analyzing the data requirement of ASTRAL-I requires answering the following question: how
many genes are needed before every bipartition of the species tree has an arbitrarily high probability of appearing in at least one gene tree. To our knowledge, a recent result by Uricchio et
al. is the only work on this “bipartition cover” problem [51]. Uricchio et al. give upper bounds
on the required number of genes and, as our Appendix B shows, their upper bound correspond
to an m that increases as Θ(f −(n−3) ). If that many genes were to be proved also necessary, we
would conclude that ASTRAL-I and ASTRAL* have dramatically different data requirements.
However, these analyses are very conservative and likely give only very loose upper bounds; the
data requirement of the bipartition cover problem may prove to be much less. Tighter bounds for
the bipartition cover problem are challenging to derive and need to be addressed in future work.

4.3

Discussion of simulation results

In simulations, we found the species tree to be consequential for the required number of genes
for both methods. Understanding the exact impact of tree shape would require a more extensive
study. Here we provide some conjectures.
The fact that caterpillar trees required more genes compared to balanced trees is not due to
increased discordance. For example, both trees had a quartet score of 44% with f = 0.1 or 34%
with f = 0.005 (computed based on 2000 simulated gene trees), indicating identical levels of gene
tree discordance. However, differences in data requirements of ASTRAL-II may be explained by
the dependency structure of quartets. We say a quartet is defined around an internal branch of
a tree if the middle path of the quartet corresponds exactly to that branch. Take the “centroid”
branch in each tree that divides the leaves into two sets of size four. In the caterpillar tree, there
are nine quartets around the centroid branch, and all these nine quartets share at least two leaves
(E and D). In the balanced tree, there are sixteen quartets around the centroid branch, and eight
pairs of quartets share no leaves in common. Thus, in the caterpillar case, we have fewer quartets
and quartet trees are highly dependent, providing little additional information with respect to each
other. On the other hand, in the balanced tree, we have more quartets and many of them share
few or no taxa and so, intuitively, may be less dependent. Note that having less dependence in the
quartets around a branch can potentially help ASTRAL-II because even if the binomial associated
with one of the quartets happens to give a score below m/3 for the species tree resolution, the
other quartets can potentially make up for it. In other words, having fully independent quartets
around a branch is like having more genes as input; balanced trees have more quartets and less
dependence among quartets compared to caterpillar trees.
We also observed that the double-quartet tree requires many more genes than the balanced
tree. This is counter-intuitive because the longer branch in the double-quartet tree reduces gene
tree discordance. In our simulations, the model balanced tree had quartet scores of 44% and 34%
for f = 0.1 and f = 0.005 while the double-quartet tree had scores of 71% and 68% (computed
based on 2000 simulated gene trees). Why does the double-quartet tree require more genes despite
having dramatically less discordance?
The pattern may be related to quartet dependencies. The long branch at the root of the
double-quartet tree results in a high likelihood that all lineages coalesce before reaching the root.
For example, take the parent of (H, G); quartets around that branch include H, G, E or F , and
10

one of {A, B, C, D} (Fig. 2); all quartets that include F will all have the same exact frequency
(same for E), hence providing no additional information with respect to each other. Thus, the
long branch at the top is causing the eight quartets around each of the four lower branches to have
very redundant information; this redundancy is unhelpful, as we discussed earlier. Thus, counterintuitively, reduced discordance at the basal branch may also be resulting in more uncertainty
about resolutions of the branches closer to the tips. The double-quartet tree that looks easier in
terms of gene tree discordance turns out to be more difficult for ASTRAL-II to resolve correctly.
Our results do not imply that the long branch itself is hard to recover. It simply indicates that
it may complicate the recovery of short branches. This pattern, if observed more broadly, would
incidentally resemble the long branch attraction problem in sequence-based analyses. We note
that the conditions where ASTRAL-II had surprisingly high data requirements resemble the types
of trees that we invoked to prove the lower bounds for arbitrary n.
Finally, our simulations start to shed light on differences between conditions where ASTRAL-II
and NJst can each be expected to work better. In model conditions where all branches where short,
ASTRAL-II outperformed NJst, especially with caterpillar trees. Both methods were sensitive to
the presence of long branches, but ASTRAL-II more so. More broadly, our research indicates that
a thorough examination of tree shapes for which ASTRAL and NJst perform well is required.
Finally, the simulation results lead to the conjecture that NJst also requires O(f −2 ) gene trees;
we hope future work will prove or disprove this conjecture.

4.4

Practical implications

Asymptotic results are not meant to provide practitioners with means to predict how many gene
trees are needed. Two issues hamper such a goal. On the one hand, some of the parameters of the
species tree (e.g., f ) are not known in advance. On the other hand, the lower and upper bound
values match only asymptotically and in practice, the effect of the unknown constants is important.
Practitioners trying to decide the number of required genes in advance may be only marginally
helped by our results. For example, given a large number of loci sampled from a few species (e.g.,
those with full genomes sequenced), they can empirically estimate the number of required genes
(e.g., judging accuracy by examining local posterior probabilities generated by ASTRAL [47]).
Then, to predict the number of required genes if they were to sample more species, they could
extrapolate using our results; this would still need guessing the reduction in f due to increased
sampling and such estimates will have to be taken as ballpark guesses.
Sample complexity results are useful to judge the effectiveness of methods. The best possible
sample complexity for the problem of inferring species trees from gene tree topologies is not known.
If future work proves O(f −2 ) to be the best possible, our results will give practitioners reasons
to have confidence in the efficiency of ASTRAL. On the other hand, if some other tool is ever
proved to have a better sample complexity, such asymptotic results should give practitioners some
indication that alternative algorithms may be preferable. In either case, asymptotic results are not
a replacement for careful empirical and simulation studies and should be used as a complementary
source of information when judging methods.

Acknowledgements
The authors thank Erfan Sayyari and anonymous reviewers for helpful comments. The work was
supported by National Science Foundation (NSF) grant IIS-1565862 to SM and SS and NSF grant
DMS-1149312 (CAREER) and NSF grant DMS-1614242 to SR. Computations were performed on
the San Diego Supercomputer Center (SDSC) and through the Extreme Science and Engineering
Discovery Environment (XSEDE), supported by NSF grant number ACI-1053575.

11

References
[1] DupTree: a program for large-scale phylogenetic analyses using gene tree parsimony. Bioinformatics, 24(13):1540–1541, 2008.
[2] Inferring species trees from incongruent multi-copy gene trees using the Robinson-Foulds
distance. Algorithms for Molecular Biology, 8:28, 2013.
[3] Phylotranscriptomic analysis of the origin and early diversification of land plants. Proceedings
of the National Academy of Sciences, 111(45):E4859–4868, 2014.
[4] E S Allman, James H. Degnan, and J A Rhodes. Identifying the rooted species tree from the
distribution of unrooted gene trees under the coalescent. J. Math. Biol., 62:833–862, 2011.
[5] Elizabeth Allman, James H. Degnan, and John Rhodes. Species tree inference from gene
splits by Unrooted STAR methods. IEEE/ACM Transactions on Computational Biology and
Bioinformatics, PP(99):1–7, 2016.
[6] Sónia C S Andrade, Marta Novo, Gisele Y Kawauchi, Katrine Worsaae, Fredrik Pleijel, Gonzalo Giribet, and Greg W Rouse. Articulating “archiannelids”: Phylogenomics and annelid
relationships, with emphasis on meiofaunal taxa. Molecular Biology and Evolution, 2015.
[7] Robert B. Ash. Information theory. Dover Publications, Inc., New York, 1990. Corrected
reprint of the 1965 original.
[8] Eric Bapteste, Leo van Iersel, Axel Janke, Scot Kelchner, Steven Kelk, James O. McInerney,
David A. Morrison, Luay Nakhleh, Mike Steel, Leen Stougie, and James Whitfield. Networks:
Expanding evolutionary thinking, 2013.
[9] Bastien Boussau, GJ J Szöllősi, and Laurent Duret. Genome-scale coestimation of species
and gene trees. Genome Research, 23(2):323–330, 2013.
[10] David Bryant, Remco Bouckaert, Joseph Felsenstein, Noah A. Rosenberg, and Arindam Roychoudhury. Inferring species trees directly from biallelic genetic markers: Bypassing gene
trees in a full coalescent analysis. Molecular Biology and Evolution, 29(8):1917–1932, 2012.
[11] Johanna Taylor Cannon, Bruno Cossermelli Vellutini, Julian Smith, Fredrik Ronquist, Ulf
Jondelius, and Andreas Hejnol. Xenacoelomorpha is the sister group to Nephrozoa. Nature,
530(7588):89–93, 2016.
[12] Julia Chifman and Laura S Kubatko. Quartet Inference from SNP Data Under the Coalescent
Model. Bioinformatics, 30(23):3317–3324, 2014.
[13] Gautam Dasarathy, Robert Nowak, and Sebastien Roch. Data requirement for phylogenetic
inference from multiple loci: a new distance method. IEEE/ACM Transactions on Computational Biology and Bioinformatics (TCBB), 12(2):422–432, 2015.
[14] James H. Degnan. Anomalous unrooted gene trees. Systematic Biology, 62:574–590, 2013.
[15] James H. Degnan and Noah A. Rosenberg. Discordance of Species Trees with Their Most
Likely Gene Trees. PLoS Genet, 2(5), 2006.
[16] James H. Degnan and Noah A. Rosenberg. Gene tree discordance, phylogenetic inference and
the multispecies coalescent. Trends in ecology & evolution, 24(6):332–340, 2009.
[17] James H. Degnan and Laura A Salter. Gene tree distributions under the coalescent process.
Evolution: International Journal of Organic Evolution, 59(1):24–37, 2005.
[18] Rick Durrett. Probability: theory and examples. Cambridge Series in Statistical and Probabilistic Mathematics. Cambridge University Press, Cambridge, fourth edition, 2010.
12

[19] Julien Y Dutheil, Gertrud Mannhaupt, Gabriel Schweizer, Christian M K Sieber, Martin
Münsterkötter, Ulrich Güldener, Jan Schirawski, and Regine Kahmann. A tale of genome
compartmentalization: the evolution of virulence clusters in smut fungi. Genome biology and
evolution, 8(3):681–704, 2016.
[20] Scott V Edwards, Liang Liu, and Dennis K Pearl. High-resolution species trees without
concatenation. Proceedings of the National Academy of Sciences of the United States of
America, 104(14):5936–5941, 2007.
[21] Peter Erdos, Mike Steel, L Szekely, and T Warnow. A few logs suffice to build (almost) all
trees: Part II. Theoretical Computer Science, 221(1-2):77–118, 1999.
[22] Thomas C Giarla and Jacob A Esselstyn. The Challenges of Resolving a Rapid, Recent
Radiation: Empirical and Simulated Phylogenomics of Philippine Shrews. Systematic biology,
page syv029, 2015.
[23] Joseph Heled and Alexei J Drummond. Bayesian inference of species trees from multilocus
data. Molecular Biology and Evolution, 27(3):570–580, 2010.
[24] Chien-Hsun Huang, Renran Sun, Yi Hu, Liping Zeng, Ning Zhang, Liming Cai, Qiang Zhang,
Marcus A Koch, Ihsan Al-Shehbaz, and Patrick P Edger. Resolution of Brassicaceae phylogeny using nuclear genes uncovers nested radiations and supports convergent morphological
evolution. Molecular biology and evolution, 33(2):394–412, 2016.
[25] J F C Kingman. On the genealogy of large populations. Journal of Applied Probability,
19(1982):27–43, 1982.
[26] Manuel Lafond and Céline Scornavacca. On the Weighted Quartet Consensus problem. 2016.
[27] Bret R Larget, Satish K Kotha, Colin N Dewey, and Cécile Ané. BUCKy: gene tree/species
tree reconciliation with Bayesian concordance analysis. Bioinformatics, 26(22):2910–2911,
2010.
[28] Christopher E Laumer, Andreas Hejnol, and Gonzalo Giribet. Nuclear genomic signals of the
’microturbellarian’ roots of platyhelminth evolutionary innovation. eLife, 4, 2015.
[29] Liang Liu. BEST: Bayesian estimation of species trees under the coalescent model. Bioinformatics (Oxford, England), 24(21):2542–2543, 2008.
[30] Liang Liu and Lili Yu. Estimating species trees from unrooted gene trees. Systematic Biology,
60:661–667, 2011.
[31] Liang Liu, Lili Yu, and Scott V Edwards. A maximum pseudo-likelihood approach for estimating species trees under the coalescent model. BMC Evolutionary Biology, 10(1):302,
2010.
[32] Liang Liu, Lili Yu, Dennis K Pearl, and Scott V Edwards. Estimating species phylogenies
using coalescence times among sequences. Systematic Biology, 58(5):468–477, 2009.
[33] Gábor Lugosi. Concentration-of-measure inequalities. 2004.
[34] Wayne P. Maddison. Gene Trees in Species Trees. Systematic Biology, 46(3):523–536, 1997.
[35] Siavash Mirarab, Rezwana Reaz, Md. Shamsuzzoha Bayzid, Théo Zimmermann, M Shel Swenson, and Tandy Warnow. ASTRAL: genome-scale coalescent-based species tree estimation.
Bioinformatics, 30(17):i541–i548, 2014.
[36] Siavash Mirarab and Tandy Warnow. ASTRAL-II: coalescent-based species tree estimation
with many hundreds of taxa and thousands of genes. Bioinformatics, 31(12):i44–i52, 2015.

13

[37] Elchanan Mossel and Sebastien Roch. Incomplete lineage sorting: consistent phylogeny estimation from multiple loci. IEEE/ACM Transactions on Computational Biology and Bioinformatics, 7(1):166–171, 2010.
[38] P Pamilo and M Nei. Relationships between gene trees and species trees. Mol Biol Evol,
5(5):568–583, 1988.
[39] Daniel A. Pollard, Venky N. Iyer, Alan M. Moses, and Michael B. Eisen. Widespread discordance of gene trees with species tree in drosophila: Evidence for incomplete lineage sorting.
PLoS Genetics, 2(10):1634–1647, 2006.
[40] Bruce Rannala and Ziheng Yang. Bayes estimation of species divergence times and ancestral
population sizes using DNA sequences from multiple loci. Genetics, 164(4):1645–1656, 2003.
[41] MD Rasmussen and Manolis Kellis. Unified modeling of gene duplication, loss, and coalescence
using a locus tree. Genome research, 22(4):755–765, 2012.
[42] Sebastien Roch. An analytical comparison of multilocus methods under the multispecies
coalescent: the three-taxon case. In Proc. Pacific Symposium on Biocomputing, volume 18,
pages 297–306, 2013.
[43] Noah A. Rosenberg. Discordance of species trees with their most likely gene trees: A unifying
principle. Molecular Biology and Evolution, 30(12):2709–2713, 2013.
[44] C. J. Rothfels, F.-W. Li, E. M. Sigel, L. Huiet, a. Larsson, D. O. Burge, M. Ruhsam, M. Deyholos, D. E. Soltis, C. N. Stewart, S. W. Shaw, L. Pokorny, T. Chen, C. DePamphilis, L. DeGironimo, L. Chen, X. Wei, X. Sun, P. Korall, D. W. Stevenson, S. W. Graham, G. K.-S.
Wong, and K. M. Pryer. The evolutionary history of ferns inferred from 25 low-copy nuclear
genes. American Journal of Botany, 102:ajb.1500089–, 2015.
[45] Greg W Rouse, Nerida G Wilson, Jose I Carvajal, and Robert C Vrijenhoek. New deep-sea
species of Xenoturbella and the position of Xenacoelomorpha. Nature, 530(7588):94–97, 2016.
[46] Erfan Sayyari and Siavash Mirarab. Anchoring quartet-based phylogenetic distances and
applications to species tree reconstruction. BMC Genomics, 17(S10):101–113, 2016.
[47] Erfan Sayyari and Siavash Mirarab. Fast Coalescent-Based Computation of Local Branch
Support from Quartet Frequencies. Molecular Biology and Evolution, 33(7):1654–1668, 2016.
[48] Sen Song, Liang Liu, Scott V Edwards, and Shaoyuan Wu. Resolving conflict in eutherian
mammal phylogeny using phylogenomics and the multispecies coalescent model. Proceedings
of the National Academy of Sciences of the United States of America, 109(37):14942–7, 2012.
[49] Jeet Sukumaran and Mark T Holder. DendroPy: a Python library for phylogenetic computing.
Bioinformatics, 26(12):1569–1571, 2010.
[50] G J Szöllõsi, E Tannier, V Daubin, and B Boussau. The inference of gene trees with species
trees. Systematic Biology, 64(1):e42–e62, 2014.
[51] Lawrence Uricchio, Tandy Warnow, and Noah Rosenberg. An analytical upper bound on the
number of loci required for all splits of a species tree to appear in a set of gene trees. 2016.
[52] Pranjal Vachaspati and Tandy Warnow. ASTRID: Accurate Species TRees from Internode
Distances. BMC Genomics, 16(Suppl 10):S3, 2015.
[53] Ya Yang, Michael J Moore, Samuel F Brockington, Douglas E Soltis, Gane Ka-Shu Wong,
Eric J Carpenter, Yong Zhang, Li Chen, Zhixiang Yan, Yinlong Xie, Rowan F Sage, Sarah
Covshoff, Julian M Hibberd, Matthew N Nelson, and Stephen A Smith. Dissecting molecular
evolution in the highly diverse plant clade Caryophyllales using transcriptome sequencing.
Molecular biology and evolution, pages msv081–, 2015.
14

[54] Yun Yu, J. Dong, Kevin Liu, and Luay Nakhleh. Maximum likelihood inference of reticulate
evolutionary histories. Proceedings of the National Academy of Sciences, 111(46):16448–16453,
2014.
[55] Zhi-Yong Yuan, Wei-Wei Zhou, Xin Chen, Nikolay A Poyarkov, Hong-Man Chen, NianHong Jang-Liaw, Wen-Hao Chou, Nicholas J Matzke, Koji Iizuka, Mi-Sook Min, Sergius L
Kuzmin, Ya-Ping Zhang, David C Cannatella, David M Hillis, and Jing Che. Spatiotemporal
Diversification of the True Frogs (Genus Rana): A Historical Framework for a Widely Studied
Group of Model Organisms. Systematic Biology, 2016.
[56] Sha Zhu and James H Degnan. Displayed Trees Do Not Determine Distinguishability Under
the Network Multispecies Coalescent. Systematic biology, 66(2):283–298, 2016.

A
A.1

Proofs
Notation and definitions


We denote the number of species by n and the number of gene trees by m. Let k = n4 denote
the total number of quadruples of leaves, also known as quartets, in the species tree. We use [k]
to represent the set {1, 2 . . . , k}.
For each quartet i ∈ [k], let qi represent the unrooted tree topology induced by the species tree
for quartet i. There are three possible unrooted tree topologies that can be induced by gene trees
and we represent these three configurations by qi1 , qi2 , and qi3 . We choose to use qi1 to denote the
configuration of the species tree, i.e., qi1 = qi . Also let nij , for j = 1, 2, 3, represent the number
P3
of gene trees that induce the qij topology. By definition, j=1 nij = m for all i ∈ [k].
Under the MSC model, for each qi , the probability that a gene tree has configuration qi1
is pi = 1 − 32 e−di where di is the length of the middle path in the quartet topology of the
species tree measured in coalescent units [4]. The other two configurations have equal probability
ri = (1 − pi )/2. Let us denote by δi the quantity pi − 1/3 and let δ be the minimum over all δi s.
All quartets defined around the shortest branch in the species tree have δi = δ and we use f to
denote the length of this shortest species tree branch.
Pk
ASTRAL* returns the species tree that maximizes the score s(T ) = i=1 ni,ti (T ) where ti (T )
gives the index of the
Pktopology for quartet i found in tree T . Hence observe that the score of the
true species tree is i=1 ni1 .

A.2

Proof of Theorem 2.1

√
Theorem 2.1. Consider a model species tree with minimum branch length f < log( 2). Then,
for any  > 0, ASTRAL* returns the true species tree with probability at least 1 −  if the number
of input error-free gene trees satisfies
 n 
4 4
9
1
m > log
.
(5)
2

(1 − e−f )2
To prove this result, we require an upper bound on the probability of error in terms of m, n
and f . For that, we observe that a sufficient condition for ASTRAL* to return the correct species
tree is that for all the quartets i ∈ [k], the true topology is observed in a majority of the m gene
trees, i.e., ni1 > max{ni2 , ni3 }. Thus, for ASTRAL* to make an error, at least one quartet must
have an alternate topology in a majority of the gene trees. We upper bound the probability of
this event for one quartet by using Hoeffding’s inequality, and then take a union bound over all
quartets to bound the overall probability of error.
Let w < δ/2. We use Ai to refer to the event ni1 /m ∈ (pi − w, pi + w), and use Bi and Ci to
refer to events ni2 /m ∈ (ri − w, ri + w) and ni3 /m ∈ (ri − w, ri + w), respectively.

15

Lemma A.1. For any w < δ/2, ASTRAL* is guaranteed to return the true species tree if the
event Di = Ai ∩ (Bi ∪ Ci ) occurs for all i ∈ [k].
Proof. Recall the definition of the score s(T ) and observe that ASTRAL* returns the true species
tree when the following condition holds:
ni1 > max{ni2 , ni3 },

∀i ∈ [k].

We claim that, for all values of w smaller than δ/2, we can ensure that ni1 is the largest among
nij , when the event Di occurs. Indeed, let Di hold and assume w.l.o.g. that Bi occurs. In the
worst case, ni1 = (pi − w)m and ni2 = (ri − w)m, which results in ni3 = (ri + 2w)m because the
sum of all three nij s is m. In this case, we need ri + 2w = (1 − pi )/2 + 2w < pi − w, or after
multiplying by 2 and rearranging 3pi − 1 > 6w. This is true for all w < 21 δi = 12 (1/3 − pi ).
Lemma A.2. Suppose the input to ASTRAL* is m true gene trees generated under the MSC
model. For every  > 0, ASTRAL* reconstructs the correct species tree with probability at least
1 −  if
 n 
4 4
1
.
(6)
m > 2 log

δ2
Proof. If E denotes the event that ASTRAL* returns the wrong species tree, then by Lemma A.1
we have by a union bound

P (E) ≤ P ∪ki=1 Dic
(7)
= P (∪ki=1 Aci ∪ (Bic ∩ Cic ))

≤ P ∪ki=1 (Aci ∪ Bic )
≤

k
X

[P (Aci ) + P (Bic )]

(8)
(9)
(10)

i=1

Thus we need a bound on P (Aci ) and P (Bic ). Let us define a random variable Xig , for i ∈ [k] and
g ∈ [m], which takes the value 1 if the configuration of the ith quartet in the g th gene tree is the
m
same as qi ; and is zero otherwise. Then
ig )g=1 is a sequence of independent Bernoulli random
P(X
m
variables with parameter pi and ni1 = g=1 Xig . The event Aci = {ni1 ∈
/ (m(pi − w), m(pi + w))}
δ
with w < 2 becomes increasingly unlikely as m grows and its probability can be bounded using
Hoeffding’s inequality (see e.g.[33]) to get
2

P (Aci ) ≤ 2e−2mw .
The same bound holds for P (Bic ). Plugging this into (7), we get
 
2
n
P (E) ≤
4e−2mw .
4
To make the Right Hand Side (RHS) smaller than , we take
 n 
4 4
1
m≥
log
.
2w2


(11)

(12)

The result follows from the fact that we can choose any w < δ/2.
Proof of Theorem 2.1. Lemma A.2 gives us the gene tree requirement in terms of the parameter
−f
)
δ. The result in Theorem 2.1 immediately follows by replacing δ with 2(1−e
. In the limit of
3
−f 2

small f , we note further that (1−ef 2 ) increases to 1. This gives the O(f −2 ) dependence of the
bound on m as mentioned in the Remark after the statement of the theorem.
16

A.3

Proof of Theorem 2.2

Theorem 2.2. For a species tree with n = 4 leaves, let f be the length of the single internal branch
in the unrooted species tree. Then, for any  ∈ [0, 0.5), there exists an f0 > 0 and a constant d ,
such that for all m ≤ d /f 2 with f ≤ f0 , the probability that ASTRAL* outputs the wrong quartet
is at least .
Proof. We divide the proof into two cases. First we show that the result holds for m ≤ c /f , for
some constant c . That follows from the fact that, in that case, a non-trivial fraction of genes
coalesce only in the ancestral root population. We then show separately that the result also holds
for c /f < m < d /f 2 , for some choice of d . There, we use the Berry-Esséen theorem to control
the probability that a majority of the genes return the wrong topology.
Case 1 (m ≤ c /f ) Let G be the event that, in all gene trees, all the lineages reach the root
population without any prior coalescence. Depending on the rooting, the quartet species tree can
either have the symmetric topology or the asymmetric topology. In both cases we have
P (G) = (e−f )m = e−mf .

(13)

(Recall that f represents the internal branch length in the unrooted tree. Thus it equals the
length of the single internal branch in the asymmetric case and it equals the sum of the two
internal branch lengths in the symmetric case.)
Under G, for every gene tree, all three topologies are equally likely. So n11 , n12 and n13 have
the same distribution. (There is only one quartet in this case.) Since ASTRAL* selects the
output tree randomly in case of a tie, by symmetry we can conclude that P (E|G) = 2/3 and
hence, by (13),
2
P (E) ≥ P (E|G)P (G) ≥ e−c .
3
2
) ensures finally that P (E) ≥ .
For all  ∈ [0, 0.5), the choice c = log( 3

Case 2 (c /f < m < d /f 2 )
We write

Recall that, when n = 4, we let F be the event that {n11 ≤ m/3}.
P (F ) = P (n11 ≤ m/3) = P (Y ≤ y)

where Y = √n11 −pm

p(1−p)m

√
and y = −m(p−1/3)
.
p(1−p)m

Then, by the Berry-Esséen theorem (as stated

e.g. in [18]), we have
cp(1 − p)[p2 + (1 − p)2 ]
[p(1 − p)]3/2 m1/2
cp
≤
,
[p(1 − p)]3/2 m1/2

P (Y ≤ y) − Φ(y) ≤

where Φ(·) denotes the normal CDF and c is a universal constant (which can be taken to be 3).
Thus
P (F ) = P (Y ≤ y) ≥ Φ(y) −

cp
.
[p(1 − p)]3/2 m1/2

(14)

It remains to show that the RHS above is bounded away from 0 when f is small. We estimate
the two terms separately. We know that as f → 0, p → 1/3 and p(1 − p) → 2/9, and p(1 − p)
is monotonic for 0 < p < 1/2. Thus, there exists an f1 such that p ≤ 10/27 and p(1 − p) ≥ 1/9
whenever f ≤ f1 . So, for all f ≤ f1 ,
y=

m1/2 (p − 1/3)
≤ 3m1/2 (p − 1/3),
(p(1 − p))1/2

17

and
10c
cp
≤ 1/2 .
(p(1 − p))3/2 m1/2
m
Thus we get the following relation

10c
P (F ) ≥ Φ − 3m1/2 (p − 1/3) − 1/2 .
m

(15)

Let us first consider the second term on the RHS. For m > c /f ,
10c
10cf 1/2
<
1/2
1/2
m
c
√
10c f¯
Let f¯2 be such that √c 2 = γ2 , for γ = 1/2 − . For all f ≤ f2 = min{f¯2 , f1 } and m > c /f , we
√
have 10c/ m < γ/2.
As for the first term on the RHS of (15), it is decreasing in m for a fixed f . Thus for any
f ≤ f2 , if we find a value of m as a function of f , say m0 (f ), which is larger than c /f and
for which this first term is no smaller than 1/2 − γ/2, then the required result also holds for all
m ∈ (c /f, m0 (f )].
Let us set m0 (f ) = d /f 2 , where d will be determined below. Then we have

 p  
γ
2 1 − e−f
−
(16)
P (F ) ≥ Φ −3 d
3
f
2
 p  γ
= Φ −2 d − ,
(17)
2
where we used that for f > 0,

1−e−f
f

≤ 1. Now, if we set


d =

Φ−1 (1/2 − γ/2)
2

2
(18)

we get that, for all f ≤ f2 and m ∈ (c /f, d /f 2 ], P (E) ≥ 1/2 − γ = . Finally, we define
f0 = min{f2 , d /c } to ensure that the above interval is non-empty.

A.4

Proof of Claim 2.1

Claim 2.1. For any n and  ∈ [0, 0.5), there exists a species tree with n leaves and shortest branch
length f such that when ASTRAL* is used with m ≤ d0 /f 2 gene trees, for some constant d0 , the
probability that ASTRAL* reconstructs the wrong tree is at least .
Proof. The idea is to reduce the proof to the case of Theorem 2.2 with an appropriate choice of
species tree. We consider a species tree with the topology (((a, b), c), X ), where X is an arbitrary
tree with (n − 3) leaves, which we denote by X = [n − 3] for convenience. We denote the edge
above (a, b) by e and set its length to f . We also set the length of the the branch above X to g, a
large value which we will specify later.
For x ∈ X, let nx1 , nx2 , nx3 be the frequencies of the quartet trees congruent with the unrooted
topology of (((a, b), c), x), (((a, c), b), x), and (((a, x), b), c) respectively. When nx1 < m/3 for all
x
P∈ X, ASTRAL* cannot output the
P true tree.P To see this, note that, under that condition,
n
<
m(n
−
3)/3;
thus,
max{
x1
x∈X
x∈X nx2 ,
x∈X nx3 } > m(n − 3)/3. W.l.o.g. assume that
P
x∈X nx2 > m(n − 3)/3. Then a tree that includes the branch e cannot have the maximal score
because the unrooted topology of (((a, c), b), X ) will have a higher score. Indeed, observe that all
quartets not of the form {a, b, c, x} for some x ∈ X contribute equally to the unrooted topologies
of (((a, b), c), X ) and (((a, c), b), X ).

18

Let Hi denote the event that, in gene i, all lineages from X coalesce before reaching the root
branch (i.e., within or below the long branch with length g) and let H = ∩m
i=1 Hi . Under the
event H, nxj = nyj for all x, y ∈ X and all j ∈ {1, 2, 3}; hence P (∩x∈X nxj < m/3|H) = P (n1j <
m/3|H). Thus, {n1j < m/3} is sufficient to produce the wrong topology. We now bound the
probability of that event.
Let us define


1 1
1 1
− >
1 = + , =  +
4 2
2 2
and
η = 1 − /1 .
From Theorem 2.2, we know that there exists an f0 such that, for all f ≤ f0 , if m ≤ d1 /f 2 , then
P (n1j < m/3|H) ≥ 1 . Let
m0 = m0 (f ) = d1 /f 2 ,
as defined in the proof of Theorem 2.2. If E denotes the event that ASTRAL* produces the wrong
output, then we have
P (E) ≥ P (n1j < m/3|H)P (H) ≥ P (H)1 .
(19)
Proving a lower bound of 1 − η for P (H) then implies the result, by our choice of η. For this, we
choose appropriately the value of g. Specifically, we pick g large enough to ensure that P (Hi ) ≥
1 − η1 for η1 = 1 − (1 − η)1/m0 . Due to the independence of the genes, we get that for all m ≤ m0
P (H) ≥ (1 − η1 )m ≥ (1 − η1 )m0 = 1 − η.

(20)

Combining this with (19), we get that P (E) ≥ . Thus, we have shown the existence of a tree
with n leaves for which ASTRAL* is wrong with probability at least  ∈ [0, 0.5) if the number of
input gene trees is less than or equal to d0 /f 2 , where d0 = d+1/2(1/2−) .

A.5

Proof of Theorem 2.3

Theorem 2.3. For any ρ ∈ (0, 1) and a ∈ (0, 1), there exist constants f0 and n0 such that the
following holds. For all n ≥ n0 and f ≤ f0 , there exists a species tree with n leaves and shortest
n
branch length f such that when ASTRAL* is used with m ≤ a5 log
f 2 gene trees, the event E that
ASTRAL* reconstructs the wrong tree has probability
P (E) ≥ 1 − ρ.

(21)

Proof. To simplify the presentation, we assume that n is divisible by three. (It is straightforward
to generalize the argument.) Consider a species tree constructed as follows. Make n/3 rooted
triplet trees all of them with an internal branch of length f and arbitrarily index the triplet trees
by i ∈ [n/3]. Then, pick any arbitrary tree with n/3 leaves labeled by [n/3] and all branch
lengths equal to exactly g. Finally, connect each leaf of this arbitrary tree with the root of the
corresponding triplet tree. Let Ji denote the event that, in all gene trees, all three lineages of the
n/3
triplet i coalesce in its root branch. Let J = ∩i=1 Ji . We proceed in two steps:
1. First we introduce a result analogous to Theorem 2.2 with an extra log n factor in the number
of genes. This gives us a smaller probability of error of the form n−a for some 0 < a < 1 on
a single quartet.
2. Then we consider the species tree above in which there are Ω(n) chances for a quartet-based
error in ASTRAL* to occur and show that at least one of these happens with sufficiently
large probability.

19

Step 1 Let us consider one of the triplets above i ∈ [n/3], with a topology ((u, v), w), as well
as a leaf x outside of the triplet. Let ni,1 be the number of genes with the correct species tree
topology and let Fi denote the event that {ni,1 ≤ m/3}. For reasons that will be made clear in
Step 2, we are interested in P (Fi |Ji ). Conditioned on Ji , let p0g and rg0 = (1/2)(1 − p0g ) be the
probabilities of observing ((u, v), (w, x)) and ((u, w), (v, x)) or ((v, w), (u, x)) respectively in the
unrooted gene trees. As g → +∞, p0g ↓ 1 − (2/3)e−f . Let g0 be large enough (as a function of f )
so that, for all g ≥ g0 , p0g − 1/3 ≤ (2.1/3)(1 − e−f ). Assume from now on that g ≥ g0 .
To deal with the extra log n factor in the number of genes, we revisit (14). This time, instead
of Berry-Esséen, we use a tailored estimate on the binomial. We will need the following bounds.
The Kullback-Liebler divergence between two Bernoulli random variables satisfies
1−a
a
+ (1 − a) log
b
1−b



a
1−a
≤a
− 1 + (1 − a)
−1
b
1−b
(a − b)2
,
=
b(1 − b)

D(akb) = a log

where we used the fact that log x ≤ x − 1. By [7, Lemma 4.7.1], we also have
 
m k
exp (−mD (k/mkp))
p (1 − p)m−k ≥ p
.
k
8m(k/m)(1 − k/m)
Using these two bounds we get
m/3 

P (Fi |Ji ) =

X m
(p0g )k (1 − p0g )m−k
k

k=0

m/3

≥

X
√
k=m/3− 2m
m/3

≥

X
√
k=m/3− 2m

 
m
(p0g )k (1 − p0g )m−k
k

exp −mD k/m p0g
p
8m(k/m)(1 − k/m)

m/3

≥

X
√
k=m/3− 2m

√


1
exp −mD k/m p0g
2m

m/3

(p0g − k/m)2
1
√
exp −m 0
≥
pg (1 − p0g )
√
2m
k=m/3− 2m
!
p
(p0g − 1/3 + 2/m)2
≥ exp −m
,
p0g (1 − p0g )

!

X

where we used that p0g > 1/3 on the last line. Above we assumed, to simplify the notation, that
√
2m is an integer. Using p0g −1/3 ≤ (2.1/3)(1−e−f ), 1−e−f ≤ f and m ≤ m1 (f ) = (a/5)f −2 log n,
we get
P (Fi |Ji )
!
p
(a/5) log n {2.1f /3 + f / (a/10) log n}2
≥ exp −
·
.
f2
p0g (1 − p0g )
Choose constants f0 and n0 depending on a ensuring that f ≤ f0 and n ≥ n0 imply that p0g (1 −
p
p0g ) ≥ 1/9 and that the expression in curly brackets in the previous display is less than f 5/9.
20

Then
P (Fi |Ji ) ≥ n−a .

(22)

Step 2 For each i ∈ [n/3], let Ei denote the event that all quartets containing the triplet i have
an alternative topology (i.e., different from the species tree) in the majority. We note that, when
Ji occurs, then Fi implies Ei . If E denotes the event that ASTRAL* outputs a wrong tree, then
a sufficient condition for an error to occur is that at least one of the Ei events occurs. Under the
MSC, the events {Fi ∩ Ji : i ∈ [n/3]} are independent. Thus,
Y
[1 − P (Fi ∩ Ji )]
P (E) ≥ 1 −
i∈[n/3]

Y

≥1−

[1 − n−a P (Ji )],

(23)

i∈[n/3]

by (22). It remains to bound P (Ji ).
Let 1 − ρ2 be the probability that the full coalescence for triplet i occurs in its root branch in
one gene tree. Then P (Ji ) = (1 − ρ2 )m ≥ (1 − ρ2 )m1 (f ) for m ≤ m1 (f ). Thus, for any ρ1 ∈ (0, 1),
we can choose g1 large enough to ensure that g ≥ g1 implies ρ2 ≤ 1 − (1 − ρ1 )1/m1 (f ) and therefore
P (Ji ) ≥ 1 − ρ1 . Assume from now on that g ≥ max{g0 , g1 }.
Going back to (23), using the fact that (1 − x1 )x ≤ e−1 , we have
P (E) ≥ 1 − [1 − n−a (1 − ρ1 )]n/3


(1 − ρ1 )n1−a
≥ 1 − exp −
3
≥ 1 − ρ,
for n ≥ n1 , for an appropriate n1 depending on a.

B

Bipartition cover bounds

A set of gene trees is called a bipartition cover of a species tree S if the set of bipartitions included
in the gene trees includes all the bipartitions of S. In [51], the authors show that input gene trees
form a bipartition cover with probability at least q if
m ≥ Mn,q =

log((1 − q)(n − 3))
Cn,q
=
,
log(1 − gn−2,1 (f ))
− log(1 − gn−2,1 (f ))

where gn−2,1 (f ) is the probability that n − 2 lineages coalesce to 1 within time f and Cn,q is
implicitly defined above. When n is fixed, we claim that Mn,q is Θ(f −(n−3) ) as f → 0.
Pn−3
Note that gn−2,1 (·) is the CDF of a random variable Tn−2 =  1 Wi where the Wi s are
independent exponential waiting times with parameters λi = n−1−i
. Thus,
2
P (Tn−2 ≤ f ) ≤

n−3
Y

n−3
Y

P (Wi ≤ f ) ≤ (

i=1

Since we have limf ↓0
So we can write

1−e−λi f
λi f

λi )f n−3 .

= 1, for f small enough we have

P (Tn−2

n−3
Y

(24)

i=1


P Wi ≤

1−e−λi f
λi f


f
≤ f) ≥
n−3
i=1
 n−3
Y  f n−3
≥ 0.9
λi
.
n−3
i=1
21

≥ 0.91/(n−3) for all i.

Thus for small f and fixed n, the probability gn−2,1 (f ) scales as f n−3 .
We now note that
Mn,q = Cn,q

(a)
1
Cn,q
≤
,
− log(1 − gn−2,1 (f ))
gn−2,1 (f )

(25)

where (a) follows from log(x) ≥ 1 − 1/x for positive x. Also, by using log(x) ≤ x − 1 and assuming
that f is small enough to ensure that the upper bound on gn−2,1 in (24) is smaller than 0.9, we
have
0.1
1 − gn−2,1 (f )
≥ Cn,q
.
(26)
Mn,q ≥ Cn,q
gn−2,1 (f )
gn−2,1 (f )
The bounds in (25) and (26) tell us the term Mn,q is Θ(f −(n−3) ), which suggests that the number of
trees required by heuristic ASTRAL-I is asymptotically significantly larger than that of ASTRAL*.
Thus our theoretical analysis shows that, consistently with simulation results presented in [36],
for large n there should be conditions where ASTRAL-I does not work well, but ASTRAL* has
high accuracy. Note that ASTRAL-II has sought to close this gap heuristically and seems to have
succeeded to a large extent based on simulation results [36].

22

