Phylogenomic evidence for multiple losses of flight in ratite birds

All phylogenetic analyses of our 20-locus dataset revealed strong support for the ostrich as the sister group of all other paleognaths, placing the volant tinamous within ratites and making ratites polyphyletic.^o When all loci were combined, the branch uniting tinamous with all ratites except the ostrich received 100% bootstrap support in maximum-likelihood (ML) analyses (Fig. 1), whether analyses were partitioned by locus (each locus assigned its own best-fit evolutionary model) or unpartitioned (a single model for all loci). Maximum-parsimony (MP) analyses produced similar results (Fig. S1A). The critical branch uniting non-ostrich paleognaths had a posterior probability of 1.00 in both partitioned and unpartitioned Bayesian analyses (Fig. 1) and highly significant support (P < 0.001) in the Shimodaira–Hasegawa (SH) test (see SI Methods and Fig. S2C).

Support for ratite polyphyly is robust to the assumptions of specific analyses. The critical branch was strongly supported in analyses using purine/pyrimidine (RY) coding (Fig. S1A), only protein-coding exons (Fig. S3) and different data partitioning schemes (see SI Methods and Fig. S1B). Furthermore, our conclusions do not reflect the specific set of outgroups used; analyses including crocodilians (Fig. 2) and/or up to 150 additional neognaths, representing all major living avian lineages, also support ratite polyphyly (25, 26, 30).

Separate analyses of individual loci show that 19 of 20 support paleognath monophyly in all analytical approaches (data not shown). The one that does not, BDNF, has a severe base compositional bias (see below). The great majority of single-locus trees (17 for ML and 15 for MP) support the ostrich as the sister group of all other paleognaths (Table 1). The probability of 15 or more of 20 independent gene trees agreeing by chance is extremely low (P = 2 × 10⁻⁹; binomial test using equiprobable trees null model). A number of loci that fail to support ratite polyphyly in individual analyses show hidden support (31, 32) in combined analyses (Table S3). Thus, the phylogenetic signal is widespread in the nuclear genome and any attempt to explain ratite polyphyly as an artifact must invoke a genome-wide systematic bias.

Beyond the strong signal present in nucleotide substitutions, three insertion/deletion events (indels) provide additional information regarding ratite phylogeny. Rare genomic changes like indels are thought to be valuable phylogenetic markers that may be free from a number of caveats that apply to nucleotide substitutions (33). The ostrich shares the ancestral character state with neognaths for two indels, an 8-bp deletion in ALDOB (Fig. 3) and a 9-bp insertion in MYC (Fig. S4). Tinamous share the derived character state with all other ratites, so both indels support the optimal tree found here. A single 1-bp indel in CLTC (Fig. S5) is found in ratites (including ostrich) but not tinamous, thus conflicting with the other indels and analyses of nucleotide substitutions. Vertebrate indels consistently exhibit less homoplasy than nucleotide substitutions (26, 34), and an examination of avian FGB intron 7 indels (34) revealed that 1-bp indels are more than twice as likely to exhibit homoplasy as longer (≥5-bp) indels. Thus, the likelihood of homoplasy in the 1-bp indel is much higher than the combined likelihood of homoplasy in the two longer indels that support ratite polyphyly.

Several distinct systematic biases that can mislead phylogenetic analysis have been characterized (35–38), and we determined whether any of these artifacts could explain ratite polyphyly.

Inconsistent phylogenetic estimation due to long-branch attraction (35) is a concern in phylogenomics. Superficially, it might appear problematic for paleognath phylogeny because tinamous have a high rate of molecular evolution (Table S4). Although the branch leading to tinamous is long (e.g., Fig. 1), it is not united with the other long branch on the tree, which leads to the outgroup. In fact, long-branch attraction would be expected to favor the conventional hypothesis of ratite monophyly rather than the surprising alternative of ratite polyphyly that we observed. We confirmed this expectation in two ways. First, a parametric bootstrap (Swofford–Olson–Waddell–Hillis, SOWH) test (39) rejected a null hypothesis of ratite monophyly (P < 0.001). Second, a series of 1,000 MP analyses in which the outgroup was replaced with a random sequence of similar base composition, simulating the longest possible branch (40), all rooted the tree within tinamous or along the branch leading to them. This confirmed our expectations and rejected long-branch attraction as an explanation for ratite polyphyly (P < 0.001).

Convergence in base composition can also create artifactual relationships (23). If the observed topology reflects compositional convergence, we would expect a tree derived from base-compositional information alone to have tinamous nested within ratites. When we clustered taxa using base-compositional distances, however, tinamous fell within ratites for few loci (Table S5). Instead, the most common grouping (10 of 20 loci) clustered ratites, a signal expected to reinforce the conventional topology. Furthermore, ratite polyphyly was supported both by analyses conducted after RY coding (Fig. S1A), which increases historical signal relative to compositional bias (37, 41), and ML analyses using a model allowing base compositional change (Table S6 and Fig. S1B). Base compositional convergence therefore cannot explain our results.

The one locus that failed to support paleognath monophyly was BDNF, which places tinamous within neognaths in ML analyses. However, BDNF has the greatest base-compositional variation and longest tinamou branch of all loci studied (Table S4). ML analyses of BDNF after RY coding do not place tinamous within neognaths (data not shown), suggesting that the anomalous results for that locus are artifactual.

Another potential problem occurs when individual gene trees differ from the species tree. Under some conditions, stochastic lineage sorting can produce a bias toward gene trees more symmetrical than the species tree (38), and the gene tree signal can predominate when many loci are analyzed (42). However, this is unlikely to occur unless the relevant branches are short relative to the coalescence time for the genes examined. The expected coalescence time for nuclear genes is twice the effective population size in generations. Because the relevant branch length and the ancestral paleognath population size are unknown, we assessed the potential impact of lineage sorting empirically. The gene-tree bias is weak (42), so it is unlikely to produce the consistency of single-gene trees we observed (Table 1). The aggregate probability of all trees with non-ostrich paleognath monophyly is 1/10, given the pure birth model for gene trees (42), so the probability of finding this result in 15 or more of 20 single-gene analyses is very small (P = 9 × 10⁻¹²). This calculation is conservative; it reflects the minimum number of loci supporting the relevant branch in single-gene MP analyses. More loci support non-ostrich paleognath monophyly using ML (Table 1 and Table S6), and all loci that fail to support this group in single-gene ML analyses have positive partitioned hidden branch support in combined analyses (Table S3).

Phylogenetic analyses can be influenced by alignment (43), although our alignment protocol was conservative, and we excluded ambiguously aligned sites. We further tested the possibility of an alignment artifact by performing two analyses using only length-conservative regions: first, those regions that could be reliably aligned with distantly related crocodilian sequences (Fig. 2) and second, protein-coding exon sequences (Fig. S3). Both analyses strongly supported the non-ostrich paleognath clade, like the full data. In sum, our analyses indicate that the signal evident in the 20-locus dataset does not reflect any known phylogenetic artifact.

Other relationships strongly supported by our results include monophyly of rheas, tinamous, and an emu–cassowary clade (e.g., Fig. 1), in agreement with previous studies. Monophyly of extant Australasian ratites (kiwis, emus, and cassowaries) is also strongly supported (Fig. 1, Fig. S1), in agreement with previous molecular analyses (9–11, 14) but contrary to some morphological analyses (9, 12). Kiwis may be difficult to place in morphological phylogenies because of the striking differences between kiwis and other ratites in life history, behavior, and size (44). The molecular support for nesting kiwis within ratites is especially interesting given suggestions that kiwis were derived directly from volant paleognaths (45).

Although placement of tinamous within the ratites is strongly supported, the sister group of tinamous is unclear. ML and Bayesian analyses place tinamous sister to Australasian ratites (Fig. 1, Fig. S1B), whereas MP and RY-coded ML analyses place tinamous sister to rheas (Fig. S1A). The alternative trees do not differ significantly when the SH test is used (Fig. S2), suggesting that the topology is sensitive to the model of sequence evolution applied. The second alternative is most parsimonious from a biogeographic standpoint, because both tinamous and rheas are exclusively Neotropical.

Our large nuclear dataset allowed us to address a controversy regarding the root of the avian tree and examine the impact of the root on our conclusions. Although the traditional view is that the root lies between paleognath and neognath clades (12, 27, 46), early analyses of mtDNA placed the root either between passerines and all other birds (47) or within passerines (48, 49), contradicting neognath monophyly. More sophisticated analyses (e.g., RY-coding) of mtDNA data strongly support the traditional rooting (22, 23, 50, 51), unlike the early analyses. Some morphological studies also suggest nonmonophyly of paleognaths (8, 45). Given these questions, we analyzed our dataset using several different methods and confirmed that the position of the root lies between paleognath and neognath clades (Fig. 2).

The strong support for ratite polyphyly in our study raises an important question: Why have most modern analyses supported ratite monophyly? We cannot fully answer this question, but we can offer plausible hypotheses. A study based on a single, short nuclear intron (52) showed limited support for ratite monophyly, probably reflecting limited power. DNA–DNA hybridization (14) required extrapolation of distances far beyond the useful range of the method. Early studies using mtDNA strongly supported ratite monophyly (9–11), but recent analyses using more sophisticated methods revealed that support for a ratite clade is weak (22, 23). Avian mtDNA evolves rapidly and exhibits high among-site rate variation (10, 11), making analyses sensitive to the evolutionary model and taxon sample used (23). Furthermore, tinamous have a much higher mtDNA evolutionary rate than other paleognaths (11), so long-branch attraction uniting tinamous with the long branch leading to the outgroup is likely (as described above). Thus, we do not consider mtDNA to be in strong conflict with ratite polyphyly; indeed, the ambivalent signal evident in recent analyses (22, 23) suggests more sophisticated analyses of mtDNA may ultimately support polyphyly.

Several morphological studies strongly support ratite monophyly (9, 12) in conflict with our nuclear genetic data. Simply collecting more data or combining molecular and morphological data are unlikely to resolve the impasse generated by these incongruent signals (53). Both signals cannot represent evolutionary relationships; at least one must be nonphylogenetic. Cranial morphology offers a potential solution, because two studies using cranial characters (13, 19) agree with our data. There is no reason to expect convergence in the cranial characters of tinamous and any particular subset of ratites. However, convergence is common in the postcranial skeletons of other flightless birds (54), and morphological convergence certainly can mislead phylogenetic analysis (55). Flightlessness is expected to result in reduction or loss of the sternal keel; reduction in size, complexity, and number of wing bones; increase in size of leg bones; and nonaerodynamic changes in plumage structure. Because the volant ancestors of each ratite lineage may have been morphologically similar, parallel evolution could have produced some ratite traits identical in state but not by descent.

Any topology that nests the volant tinamous within the flightless ratites requires either multiple losses of flight or a loss of flight in the ancestral paleognath and a regain in tinamous. Although loss and regain is more parsimonious if both transitions are equally probable, multiple losses of flight are more likely. Flight has been lost in members of 18 extant bird families, many more times in extinct groups, and hundreds of times in the family Rallidae alone (21, 54, 56). Thus, the loss of flight is much more probable than gain. Given the position of tinamous in either optimal tree based on the complete dataset (Fig. 1, Fig. S1A), flight must have been lost independently at least three times, in ostriches, rheas, and Australasian ratites. A scenario in which tinamous regained flight would be even more interesting, but there are no examples of avian lineages that have lost and regained flight.

There have been many proposals relating divergences among living ratites to the breakup of Gondwana (10, 11, 14, 16), but all assume ratite monophyly. Our phylogeny actually fits a strictly vicariant hypothesis better than previous phylogenies. Africa was the first piece of Gondwana to separate (16), and the African ostriches are the first clade to diverge (e.g., Fig. 1). Nevertheless, no proposed phylogeny, ours included, can be explained entirely by the order of separation of Gondwanan fragments.

Multiple losses of flight, with the implication of greater dispersal capability for ancestral paleognaths, make a strictly vicariant model less compelling. The existence of volant paleognaths in the Paleogene of Europe and North America also suggests that dispersal must be considered (45, 57). Dispersal of ratites is further suggested by phylogenies in which the extinct moas of New Zealand are not sister to the extant kiwis (2, 10, 11), as would be predicted by strict vicariance. Thus, fossil data confirm that simple vicariant models can be rejected. It may be possible to distinguish among evolutionary scenarios with time-calibrated trees, but the large differences in evolutionary rates evident in our analyses (e.g., Fig. 1) combined with the paucity of good fossil calibration points for paleognaths render this task difficult.

We amplified and sequenced 20 nuclear loci (see SI Methods for detailed methods). Sequences (accession numbers listed in Table S7) were aligned manually (see SI Methods); ambiguously aligned regions and sparsely sampled sites (those not present in at least four birds and three paleognaths) were excluded from analyses.

Four nonpasserine neognaths were used as outgroups in most analyses of the 14 taxon, 23,902-bp dataset (e.g., Fig. 1). Crocodilians and passerines were included only in analyses designed to test the position of the root of the avian tree using the subset of sequences for which crocodilian and avian sequences could be aligned (4,668 bp; Fig. 2).

Optimal MP and ML trees using unpartitioned data were identified by using PAUP* 4.0b10 (58). MP analyses used branch-and-bound searches with equally weighted characters and assessed support using 1,000 bootstrap replicates. The appropriate nucleotide substitution model (Table S8) for ML analyses was determined by using the Akaike information criterion and Modeltest 3.6 (59) or the models appropriate for RY-coded data (see SI Methods). ML analyses used heuristic searches with 10 random addition sequence replicates and tree bisection and reconnection (TBR) branch swapping; support was assessed by using 100 bootstrap replicates.

Three ML analyses with the data partitioned by locus were conducted. First, a partitioned RAxML (60) analysis, using the general time-reversible (GTR)+Γ model with distinct parameter values for each partition and linked branch length parameters, was conducted and support assessed by using 100 bootstrap replicates. The other two analyses used PAUP* and nhPhyML (61), allowing use of a more diverse set of nucleotide substitution models. Branch lengths were unlinked and monophyly of well established clades was assumed a priori, focusing on a “plausible set” of 315 trees (105 arrangements of the five major paleognath groups and three arrangements within tinamous; Fig. S2). The log likelihoods of all partitions for each tree were summed after being calculated using PAUP* (Table S9) or nhPhyML (not shown). PAUP* analyses used the best-fit model for each locus (Table S8) and nhPhyML analyses used the Galtier and Gouy (62) model, which allows the GC content to vary across the tree. In these analyses, support was assessed by using 1,000 resampling of estimated log-likelihoods (RELL) (63) bootstrap replicates.

Bayesian Markov chain Monte Carlo analyses were performed by using MrBayes 3.1.1 (64). We ran four chains for 10 million generations, sampling every 500 generations and discarding the first 500 trees, with other run parameters set at defaults. Partitioned Bayesian analyses used linked branch lengths.

The relationship between the probability that a clade is correct and bootstrap support is complex (65, 66), but the bootstrap is conservative under many circumstances. In contrast, Bayesian posterior probabilities can overestimate the probability that a clade is correct (e.g., 66). We also assessed branch (Bremer) support as well as hidden support and conflict for specific branches (31, 32).

We used two topology tests appropriate for comparisons of trees that were not specified a priori (67) to examine the position of the ostrich. The SH test (68) was performed by using the plausible set of 315 distinct trees (Fig. S2) because the test requires the inclusion of all trees that can be entertained as the true topology (67). The parametric bootstrap (SOWH) test (39) compares the difference in optimality scores for the empirical data on the optimal tree and a null hypothesis tree to a distribution generated by simulation. The null hypothesis topology was the most likely one with ratite monophyly (tree 19a of Fig. S2); 1,000 simulated datasets were analyzed by using MP and ML (SI Methods).

Base-compositional clustering was performed by using minimum evolution in PAUP* with a matrix of Euclidean distances between base-composition vectors for each taxon (see SI Methods). Topologies examined were limited to the plausible set (Fig. S2).

Individual gene trees may differ from the species tree. We tested whether the observed number of gene trees showing monophyly of non-ostrich paleognaths was unexpected given the null hypothesis of a completely polytomous species tree. We used binomial tests and two appropriate models of tree probability to assess this possibility (SI Methods).

To identify low-homoplasy indels, all gaps in a 19-gene, 171-taxon dataset (25) were coded by using the simple gap coding method (69). These indels were mapped by using MP on the ML tree for the 171-taxon dataset (25) and the same tree rearranged so ratites were monophyletic. Indels mapping unambiguously on the branch of interest with a consistency index ≥0.5 were examined further, and ambiguously aligned regions were removed. CLTCL1 was not included in the 171-taxon analysis (25), so it was examined independently by using a similar methodology.