We generated two sets of nucleotide sequence data by using computer simulations. In one, a 66-taxa tree representing the phylogenetic relationships among mammals was used (see in Rosenberg and Kumar 2001
). We simulated DNA evolution for 448 hypothetical genes along this tree, each with an independent set of evolutionary parameter values (base frequencies, sequence length, mean evolutionary rate, and transition–transversion rate ratio) estimated from the real sequence data (Rosenberg and Kumar 2003
). For each set of evolutionary parameters (448 different sets), the branch lengths of the model tree were estimated using the corresponding evolutionary rate. Sequence alignments were generated under the Hasegawa–Kishino–Yano (HKY) (Hasegawa et al. 1985
) model of nucleotide substitution at four different levels of rate variation among sites (α
= 0.25, 0.5, 1.0, and 2.0) that were implemented during computer simulations via a discretized gamma distribution with a very large number of categories. This resulted in a total of 1,792 alignment sets.
We also generated DNA sequence alignments containing 20–765 taxa, which were based on the corresponding sized trees derived from a master phylogeny of 765 taxa (see supplementary fig. S1
, Supplementary Material
online, in Battistuzzi et al. 
). This master phylogeny was obtained by pruning taxa and groups from the tree of 1,671 families in the Timetree of Life (Hedges and Kumar 2009
), such that the final tree was strictly bifurcating. The resultant tree of 765 taxa was scaled to time and spanned 4.2 billion years of evolution. This master topology was subsampled to produce model trees used to generate the sequence alignments containing varying number of taxa (20, 40, 60, …, 500), with one set containing all 765 taxa. Sequences were simulated using SeqGen (Rambaut and Grassly 1997
) under the HKY (Hasegawa et al. 1985
) model of nucleotide substitution with a G + C content of 48% and a transition–transversion rate ratio of 1.05, which were estimated from an alignment of small subunit rRNA sequences from 800 taxa of animals, fungi, plants, and archaebacteria. In order to make the evolutionary rate heterogeneous among tip and internal lineages, rates were varied randomly by drawing them from a uniform distribution with boundaries ±5% of the expected rate in each branch independently. We used substitution rates of 0.025, 0.050, and 0.100 per base pair per billion years to establish branch lengths. In total, 530 data sets were generated in this way and the results are presented in the main text. We also conducted 290–765 taxa simulations in which sequences evolved four times faster (0.4 substitutions per site per billion years) and found the differences between methods were very similar to those reported here (results not shown).
A benchmark comparison of ML phylogenetic inference between MEGA5, RAxML7, and PhyML3 was performed for all simulated datasets by collecting the computational and phylogenetic performance of these programs, including execution time (in seconds), the estimate of a gamma shape parameter, ML values, and topological accuracy. Because Windows is MEGA5’s native operating system, Windows executables were used for PhyML (version 3.0) and RAxML (version 7.04). All analyses were conducted on computers with identical hardware (Intel Q8400 2.66 GHz Quad Core processor and 6 GB RAM) and operating systems (64-bit Windows 7 Enterprise Edition). For direct comparison, each program was executed serially in a single thread of execution with one core utilized per dataset.
In order to generate comparable results on time and accuracy, we used identical substitution models and discrete gamma options across all programs. Because the fastest heuristic search in RaxML, MIX, assumes a general time reversible (GTR) model with four discrete gamma rate categories, we used these options in all cases, unless noted otherwise. For all three programs, analyses were conducted using the automatically generated initial trees and selecting the default options. And, heuristic searches starting with the initial trees were conducted with two different levels of branch rearrangements: quick searches (Nearest Neighbor Interchange [NNI] for MEGA5 and PhyML and MIX for RaxML) and slow searches (Close Neighbor Interchange [CNI] for MEGA5, Subtree–Pruning–Regrafting [SPR] for PhyML, and GTRGAMMA for RaxML). The accuracy of phylogenetic tree of n taxa was estimated from the topological distance (dT) between the inferred tree and the true topology was given by (n − 3 − ½dT)/(n − 3).