A new time tree reveals Earth history’s imprint on the evolution of modern birds

A new time tree for birds

We identified avian clades that have a relatively old and well-characterized fossil whose affinities are well established by phylogenetic analysis and/or derived morphologies. These fossils set a minimum bound for the stem age of the corresponding clade and its sister group (Fig. 1A) (31–33). However, phylogenetic divergence time estimation also requires constraints on the maximum age of clades; otherwise, an arbitrarily old tree would be compatible with any minimum bound. Maximum age constraints can be based on the absence of fossils in earlier fossiliferous beds, but the absence may reflect fossilization and recovery biases rather than true absence. For example, even if a detailed fossil record suggests the absence of a clade in the Paleocene of Europe, the clade may have originated earlier in another less-studied region such as South America. Fossil recovery biases can be quantified and incorporated into age estimates (34, 35) but not without theoretical and practical difficulties (32, 33). This problem is exacerbated by methods that constrain clade age based on maximum bounds (36–38) because bounds do not reflect the uncertainty associated with maximum age constraints. A better approach is to use probability density functions as priors in Bayesian analyses (32, 39), but choosing among different functions and the parameters of those functions has been a subjective and contentious issue (33, 40).

Here, we overcame these problems by empirically generating calibration priors based on the fossil record. We modeled clade age uncertainty based on distributions of fossils finds. A set of fossil ages t₁…t_n represents a sample from a distribution that has a maximum bound θ, which is the age of the clade. The objective is to obtain a probability distribution for the age θ. Several methods have been developed to estimate the stratigraphic range of taxa from observed fossil ages either by assuming simple distributions of fossil ages (41, 42) or by modeling fossil recovery potentials (34). In the simplest case, when fossil ages t₁…t_n are uniformly distributed, the likelihood of a hypothesized age θ is simply 1/θⁿ for θ > t_n (41, 43), and in the absence of prior information, this likelihood is proportional to the probability of θ (43).

An important factor that should be taken into account is the geographic dimension of the fossil record, which includes not only fossilization and exhumation biases but also the vagaries of dispersal and colonization processes. For example, a clade may have a temporally dense fossil record in a region including the oldest fossil of the clade worldwide; in a naïve assessment, this would suggest that the clade originated shortly before the age of its oldest fossil. However, if the clade originated in a different geographic region, this age estimate would correspond to the age of dispersion into the region, not the age of the clade. We minimized this problem by analyzing the fossil record at a large geographic scale using only the first fossil occurrence of a clade in each continent. Therefore, for each avian clade with a well-characterized old fossil, we determined the oldest fossil on each continent through bibliographic and database searches and kept only those clades whose set of ages do not depart from a uniform distribution (data set 1; table S1). We then parameterized probability density functions such that they matched the shape of the likelihood functions (fig. S1) and used them as clade age priors in Bayesian divergence time estimation (table S1; further details and extensions are explained in Fig. 1 and Materials and Methods).

We used the empirical priors derived above in Bayesian divergence time estimation in the program BEAST 2 (44), using a genomic data set of 1156 clock-like exons for 48 species from Jarvis et al. (3) and a data set with a much denser taxonomic representation of recombination-activating genes (RAG-1 and RAG-2) for 230 species representing 202 families and all avian orders. Both data sets resulted in very similar estimates of the avian time tree. Both placed the most recent common ancestor of extant birds in the early Late Cretaceous, with the genomic data set resulting in slightly older ages [median, 96.6 Ma; 95% highest posterior density (HPD), 84.2 to 114.3 Ma; Fig. 2A] than the RAG data set (median, 91.5 Ma; 95% HPD, 79.8 to 106.8 Ma; Fig. 2B). These estimates for the crown age of Neornithes are younger than most previous estimates (5, 10–12), but higher posterior densities overlap with the recent analysis of bird genomes (3), in which a maximum age constraint was assumed for the age of Neornithes at the Early-Late Cretaceous boundary. Our analysis recovered a relatively young age for Neornithes without prior constraints on its age.

Radiation within the three major clades—Palaeognathae, Galloanseres, and Neoaves—began around the K-Pg transition (Fig. 2). Initial cladogenesis in Galloanseres and Neoaves occurred just before the K-Pg transition, whereas the timing of initial cladogenesis in Palaeognathae is poorly constrained and spans pre– and post–K-Pg times (Fig. 2). Overall, it is possible that initial radiation of these three major crown clades occurred more or less simultaneously around (and probably before) the K-Pg event, thus suggesting common causal factors. All avian orders originated during the Paleocene (median stem age; Fig. 2), except for Anseriformes and Galliformes, which originated in the latest Cretaceous, and some orders within the higher ratites, which may have originated in the late Paleocene or early Eocene, taking into account uncertainty.

Large-scale biogeographic history of birds

Next, we integrate phylogenetic relationships and the avian time tree with the spatial history of the major lineages of birds, set in the context of paleogeography and paleoclimatology. Because biogeographic history is far too complex to be inferred from ancestral reconstructions alone, our interpretations of the spatial history of birds are also derived from the joint geographic and temporal history among deep stem lineages (that is, avifaunas), the fossil record, phylogenetic and biogeographic patterns within clades, and how those data can be reconciled with Earth history.

We identified the geographic origin of Neornithes and of its three major subclades (Palaeognathae, Galloanseres, and Neoaves) to be West Gondwana, here taken to encompass a united continental South America, West Antarctica, and portions of East Antarctica (Fig. 3). This conclusion is robust to methodological assumptions (parsimony versus probabilistic modeling) and alternative tree topologies, with the exception that the likelihood reconstruction for the most recent common ancestor of Palaeognathae is ambiguous and includes both West Gondwana and the Palearctic (fig. S3). In addition to these basal groups, we reconstructed 22 major groups of Neornithes as being present in South America by the end of the middle Paleocene (~59.2 Ma) (Fig. 3 and table S2). Because of its taxonomic composition, this avifauna likely represented a broad sample of morphological and ecological diversity, from terrestrial lineages (ratites and stem-Gruiformes), freshwater and marine (stem-Charadriiformes, Sphenisciformes, Procellariiformes, and Phoenicopteriformes), predators (Cariamiformes and Falconiformes), nocturnal and aerial specialists (Caprimulgiformes), and arboreal groups (Coraciimorphae, Psittaciformes, and Passeriformes). The estimated mean ages indicate that this diverse assemblage was present very soon after the K-Pg, but uncertainties inherent in dating do not rule out the existence of some Neoavian lineages in the latest Cretaceous. The Paleocene avian radiation in South America was paralleled by a buildup of eutherian mammal diversity (45). West Gondwana would have been a large landmass in the latest Cretaceous and Paleocene, with predominantly warm, moist, equable climates (46) that would have facilitated the accumulation of diversity within a variety of environments, with a modern rainforest being established over a large portion of South America (45, 47, 48), and with warm to cool subtropical/temperate forests in the far south on Antarctica (46, 49).

We hypothesize that Neornithes spread from West Gondwana to the rest of the world via two principal routes. The first route is a trans-Antarctic corridor that linked West Gondwana to Australia and Zealandia during the Cretaceous and early Paleogene (Fig. 3) (44). This corridor enabled the establishment of a widespread southern biota with a diversity of organisms, including plants (50), dinosaurs (51), mammals (52), and birds (Fig. 3 and table S2). The exact configuration of the land connection between South America and West Antarctica is uncertain, but this terrestrial corridor was potentially disrupted as early as ~50 Ma to as late as ~40 to 37 Ma, when a deep opening of the Drake Passage is inferred (53–56). Continental blocks comprising the corridor would have remained in proximity until ~30 Ma (55, 57).

Before ~83 Ma, Zealandia was conjugate to the eastern Australian margin and to West Antarctica in the south (58, 59). At that time, sea-floor spreading began in the Tasman Basin and south of the Campbell Plateau, severing a direct connection between Zealandia and West Antarctica, thus establishing the Tasman Gateway as the primary corridor for a K-Pg biota to enter Australia and then Zealandia (59). The South Tasman Rise did not rift from East Antarctica until the middle to late Eocene ~40 to 33 Ma; before then, the continental fragments comprising the South Tasman Rise were intermittently subaerial and shallow marine, with terrigenous sedimentation persisting through the Eocene (60, 61). Our results suggest that at least seven to nine major clades involving palaeognaths, Galloanseres, caprimulgiforms, and psittacopasseres (parrots and passerines) were shared with Australia before its separation from Antarctica (table S2). Two additional lineages, tree swifts (Hemiprocnidae; 95% HPD, 21.4 to 43.1 Ma) and the Kagu (Rhynochetidae) (95% HPD, 16.9 to 45.2 Ma), may have also been part of this interchange, considering age uncertainties. Previously, biogeographers had assumed that the ~83-Ma split between Zealandia and West Antarctica and the initial formation of the Tasman Sea implied that younger elements of the biota must have arisen via long-distance dispersal (23), but this assumption of a deep age for the Zealandia biota is no longer necessary (59). Rather, biogeographic patterns of deeper avian lineages were likely established by multiple episodes of vicariance within an early Tertiary biota that occupied Australia and emergent portions of Zealandia (table S2). The presence in Zealandia of two ancient lineages of flightless palaeognaths (kiwis and moas) [the flightless Kagu on New Caledonia and its flightless sister group (Aptornis) on New Zealand], as well as ancient lineages of parrots and passerines, is consistent with this interpretation.

We propose that there was a second Paleogene gateway for avifaunal expansion out of West Gondwana involving North America (here termed the North American Gateway hypothesis). Although it has long been assumed that there was a wide ocean barrier between South and Middle America through the Late Cretaceous and Paleogene until the Plio-Pleistocene, we hypothesize a northward avifaunal expansion at this time coincident with the first appearances of Laurasian metatherian and placental mammals in South America (45, 62), thus implying the presence of a land bridge. Although the paleogeographic setting of lower Middle America is complex and controversial, this interpretation is consistent with a growing number of tectonic models that posit a southern land bridge (63–65) that resulted from subduction, arc magnetism, and collisionary processes involving movement of the Galapagos large igneous province into the Caribbean, spanning K-Pg times.

Our results suggest a central role for North America as a gateway for biotic dispersion into the Holarctic, especially during the Paleocene to early/middle Eocene. Northern Hemisphere Paleogene environments were broadly tropical and subtropical at lower latitudes and subtropical-warm temperate at higher latitudes, especially toward the Paleocene-Eocene Thermal Maximum (66–69). We infer that these avifaunas diversified within the Nearctic and spread into the Palearctic by the end of the Paleocene and earliest Eocene ~56 to 53 Ma (Fig. 3 and table S2).

The biogeographic analyses suggest that North America played a formative role in the origins of lineages that subsequently diversified in the Old World. Thus, the Green River Formation (53 to 51 Ma) (70) in North America and the Fur Formation (55 to 54 Ma) (71) and Messel oil shales (47 Ma) (29) in Europe record many of the stem lineages that were elements of a widely distributed Paleogene avifauna (table S2). The time tree is consistent with these lineages having initiated diversification before deposition of these fossils. For example, instead of the Coraciimorphae arising in Africa (72), our optimizations indicate that Leptosomatiformes, Coliiformes, and Coraciiformes colonized the Paleotropics independently from North American ancestors (Fig. 3). Given the presence of many of these lineages in the early Paleogene of Europe, we infer that they reached the western Palearctic through a North Atlantic corridor before ~52 Ma (table S2). The tectonic history of the North Atlantic in the Paleogene indicates land connections among northeastern North America, Greenland, and Europe (73, 74), and paleontological evidence records mammalian faunal interchange until sometime in the early Eocene (68, 75), when sea-floor spreading between northeastern Greenland and Europe would have created an ocean barrier ~53 to 52 Ma (73, 76).

Expansion into other Old World landmasses is more difficult to circumscribe temporally and geographically. The avian Paleogene record for Africa is sparse, and that for Madagascar is lacking entirely (29). Some old stem clades now characteristic of these two land masses are found in the North American or European fossil record (29). Both Malagasy endemic cuckoo rollers (Leptosomatiformes) and true rollers (Coracioidea) have Holarctic Paleogene records, but when they did arrive in Africa or Madagascar is unclear. Our reconstructions suggest that lineages arrived in Madagascar in the Paleocene (mesites, Mesitornithidae), the Eocene (ground rollers, Brachypteraciidae), and the Oligocene (vangas, Vangidae; Fig. 3), suggesting that the buildup of the Malagasy endemic avifauna may have been a protracted process that spanned at least the entire Paleogene.

A second route from North America to the Old World—the Bering land bridge—would have been emergent over much of the Tertiary and, depending on eustatic sea levels and high-latitude temperature gradients, would have facilitated frequent biotic interchange between Asia and North America (77, 78). However, a Beringian corridor for Paleogene avifaunas is not evident in our analysis, partly as a result of a poor fossil record in northern Asia. Cranes and their allies (Gruiformes) provide an exception, with Eocene-Oligocene fossil families Geranoididae in North America and Eogruidae in Asia showing close affinities (29). Several other stem lineages associated with the North Atlantic corridor may also have had Paleogene histories across Beringia, including stem-Struthioniformes (ostriches) and Caprimulgiformes such as frogmouths (Podargidae) and swifts (Apodidae) (Supplementary Materials).

A novel implication of the North American Gateway hypothesis is that the roots of the “Old World” suboscines (Eurylaimides) are in the New World. Four well-supported lineages are known, but their interrelationships are uncertain (79, 80). One of these lineages is Sapayoa aenigma of the rainforests of eastern Panama and northwestern Colombia. The other three lineages are tropical African–Southeast Asian lineages, with a few young species reaching Australasia. We estimate the crown radiation of Eurylaimides to be ~37 Ma; thus, biogeography and age suggest that one lineage or possibly two lineages independently colonized the Old World, with Sapayoa being relictual. If our estimate for the age of crown Eurylaimides at ~37 Ma is correct, then it implies that these lineages entered the paleotropics via Beringia, not across the North Atlantic.

A key event in avian evolution was the diversification of songbirds (Passeri). Our reconstructions indicate that dispersion into East Gondwana of the Acanthisitti probably occurred independently of Passeri (25, 81) (Fig. 3), and we place the most recent common ancestor of the latter on the Australian landmass at 47.3 Ma (HPD, 43.1 to 51.6), which is younger than current estimates (24, 82, 83). Over the next 5 to 7 Ma, oscines underwent substantial phenotypic and ecological evolution in greater Australia; by ~40 Ma, this had resulted in the endemic corvidan radiation and early lineages that reached Zealandia (Fig. 3). Contemporaneously, oscine lineages were transferred from Australia to Southeast Asia and Africa, including the Eupetidae (40.6 Ma) and the Promeropidae (39.5 Ma), the sister group of the Laurasian Passerides (Fig. 3). These lineages mark the beginnings of the great passeridan radiation across the paleotropics, Eurasia, and eventually the New World (25). This early Australian-Asian exchange also involved lineages of corvidans, including at least stem-vireonids (34.6 Ma), campephagids (34.4 Ma), and pachycephalids (31.7 Ma). These ages are coincident in time with the early northward drift of Australia and are estimates of the origin of stem lineages, not of crown clades, which are younger (for example, the crown age for Eupetidae is 26.8 Ma). By the early Oligocene (around 33.9 Ma), portions of northwestern continental Australia were adjacent to Asia, and collision followed in the Miocene (84, 85). In addition, the Miocene-Pliocene left-lateral motion along the Sorong Fault transported continental slivers to the west toward Asia. Thus, hypotheses of long-distance dispersal across wide ocean gaps (83, 86) may not be necessary to explain this complex avifaunal interchange. Once in tropical Asia by the late Eocene, each of the major songbird lineages radiated extensively into the Oligocene and early Miocene environments of the Afrotropics and Eurasia (Fig. 3) (25, 83, 86). Many lineages then entered North America at various times across an emergent Bering land bridge (25).

Temporal dynamics of avian diversification

To explore the effects of Earth history on avian diversification dynamics, we analyzed changes in diversification rates through time. We found that net diversification rates (speciation − extinction) track changes in global climate over the entire timespan of modern bird evolution: rates increase during periods of climate cooling and decrease during periods of climate warming (Fig. 4A). Relative extinction (extinction/speciation) was low, and periods of increased extinction coincide with periods of high net diversification (Fig. 4B), suggesting that the overall dynamics is dominated by changes in rates of speciation. Environment-dependent birth-death likelihood models, in which diversification rates are functions of a time-dependent environmental variable (87), also indicate that rates of speciation and extinction are negatively correlated with global temperature (table S3).

We propose that this macroevolutionary dynamics is largely the result of climate change inducing biome fragmentation and vicariance (that is, climate-induced vicariance). During major cooling trends of the Late Cretaceous (88) and the Cenozoic (89), tropical biomes would have experienced successive waves of retraction and fragmentation (49, 66, 67), which resulted in population fragmentation and increased risk of extinction across entire faunas, but also increased chances of speciation among the newly isolated populations (90). Thus, provided that speciation rates were higher than extinction rates, this scenario results in higher net diversification of entire avifaunas during periods of climate deterioration. Conversely, expansion of thermophilic forests during warm periods (49, 66, 67) promoted connectivity and homogenization of avifaunas (also suggested by a wave of worldwide dispersion; Fig. 3), a phenomenon also documented among austral fossil floras (91) and boreal mammalian faunas (68, 75). The final burst of diversification in our analysis occurred soon after the middle Miocene cooling episode, which is associated with global retraction of tropical and mesic biomes and expansion of steppes and deserts (92–95). Therefore, the middle Miocene cooling may have fragmented distributions and increased diversification of not only humid tropical but also subtropical and mesic taxa.

Although climate-induced vicariance predicts a final major burst of diversification during Pleistocene glaciations, this burst cannot be detected in our time tree because of sampling (that is, there are no branching events in the last 10 Ma). That most avian time trees show slowdowns in diversification toward the present (95) seems to falsify this prediction, but apparent slowdowns are likely to be artifacts of nucleotide model misspecification (96), deep gene-tree coalescence (97), and incomplete species sampling (98). Alternatively, climate-induced vicariance can also be compatible with low diversification during the Pleistocene if extinction rates increased more than speciation rates, thus dampening net diversification. Complete and calibrated species trees are needed to resolve this matter.

Other events in Earth history may also have influenced diversification rates in Neornithes. In particular, the rate increase at the end of the Cretaceous occurred near the K-Pg transition, suggesting a possible association with the K-Pg extinction event. Although the rate increase started before the K-Pg transition, coinciding with the latest Cretaceous cooling trend, the magnitude of the rate increase suggests a role for the K-Pg event in increasing rates further. The extent of this effect and its particular mechanisms requires further investigation. Although the K-Pg event marked the extinction of close relatives such as Enantiornithes and basal Ornithurae, the extent to which it affected modern birds directly is unclear from the fossil record because of the scarcity of Neornithes remains before and after the event (6, 17, 29, 96). Our biogeographic reconstructions suggest that Neornithes occupied West Gondwanan continents during the K-Pg transition and, therefore, like other animals and plants (91, 97, 98), may have escaped the catastrophic effects of the impact. On the other hand, if multiple early lineages of Neornithes faced extinction during the K-Pg event, that may have created a false signal of rate increase in the reconstructed phylogenies, in which case the apparent increase is merely a consequence of missing lineages before the event (99, 100). Finally, the big-bang model posits that surviving lineages of Neornithes would have experienced ecological opportunity that may have stimulated diversification (15, 18). In any case, catastrophic extinction events may not be general drivers of diversification in modern birds because they cannot explain changes in diversification rates over most of the Cenozoic (fig. S4). Thus, there are no major meteorite impacts, volcanism, or mass extinctions in the Miocene, when diversification rates increased dramatically (Fig. 4A and fig. S4). Our data suggest, in contrast, that global changes in climate are correlated with changes in diversification rates throughout the entire history of modern birds. Statistical comparison of diversification models also favors climate over extinction events as drivers of avian diversification (Supplementary Materials). Ultimately, we cannot reject a potential effect of the K-Pg mass extinction event on modern bird diversification, but further research on this matter should consider the potentially pervasive background effect of climate-induced vicariance on diversification.

Climate-induced vicariance may have a prominent role in controlling macroevolutionary rates (90), but it may be difficult to detect in paleontological data sets. Actually, except for a clear signal of mass extinctions associated with catastrophic events, general drivers of temporal macroevolutionary dynamics in terrestrial tetrapods have been difficult to identify (87, 101–103). Climate-induced vicariance may be obscured in paleontological data sets because it results in speciation in fragmented populations that are difficult to detect in a geographically sparse fossil record (104). Conversely, during benign climates, populations expand and have higher chances of being detected, resulting in a bias in which new species remain undetected until they expand during warm periods. This hinders the detection of climate-induced vicariance and creates a false signal of positive associations between diversification and warming in the fossil record.

Climate-induced vicariance has received considerable attention in microevolutionary studies of diversification. There is now widespread evidence of population differentiation and speciation associated with habitat fragmentation and biogeographic refugia (105–110). Climate-induced vicariance also explains how niche conservatism might stimulate speciation (111). Finally, analyses of geographic variation in diversification rates among extant birds also suggest higher speciation and extinction rates in regions with harsher climates (14, 112, 113). Therefore, climate-induced vicariance may explain avian diversification through space and time, providing a unified framework for understanding large-scale biodiversity dynamics in birds.

Estimation of clade age from the fossil record

If fossil occurrences belonging to a clade are uniformly distributed between the present and the time of origin of the clade θ, the likelihood of a hypothesized age θ given the observed fossil age occurrences t₁…t_n is simply 1/θⁿ for times older than the oldest fossil t_n (41, 43). This likelihood is proportional to the probability density of θ in the absence of prior information (43). Instead of using the present as the baseline, the age of the most recent fossil can be used as the baseline, in which case the likelihood becomes 1/(θ − t₁)^{n − 1} for θ > t_n (42). Therefore, the likelihood depends on the age of the oldest fossil t_n, the number of fossil occurrences n, and the timespan encompassing those occurrences (θ − t₁). The higher the number of fossil occurrences n and the narrower the timespan between the youngest fossil and clade age, the higher the concentration of the likelihood near the oldest fossil t_n. Then, a probability density function is parameterized in a way that replicates the shape of the likelihood. This probability density function can be implemented as a clade age prior in molecular Bayesian divergence time estimation. We also considered the case in which fossil ages have uncertainty, specifically when the oldest or youngest fossils are assigned to a geological time interval that spans more than 1 Ma. In such cases, we uniformly sampled pseudoreplicated sets of fossil ages from the time interval, estimated likelihood distributions for clade age, and computed an average probability distribution across pseudoreplicates. The methodology is explained with an example in Fig. 1. Functions for generating empirical calibration priors from the fossil record in the R language are available at https://github.com/evolucionario/cladeage.

Phylogenetic time-tree analysis

We conducted a Bayesian time-tree estimation in the program BEAST 2 (44) using two molecular data sets: one that emphasizes genomic coverage and another that emphasizes taxonomic coverage. The genome-scale data set consisted of the first and second codon positions of 1156 clock-like exons for 48 species from Jarvis et al. (3), which we filtered for missing data (positions with missing data or gaps for any taxa were deleted), and resulted in a final alignment of 124,196 bases. To maximize taxonomic coverage, we used sequences of the slowly evolving recombination-activating genes (RAG-1 and RAG-2) for 230 species representing 202 families and all avian orders (alignment length, 4092 bases). Amplification and sequencing protocols for new RAG sequences followed previous studies of avian RAGs (25). For the RAG alignment, we determined an optimal substitution model and partitioning scheme simultaneously using the PartitionFinder algorithm (114). The genomic data set was analyzed under a single GTR+ gamma substitution model. Bayesian inference was performed using BEAST 2 (44). Rate heterogeneity across lineages was modeled with a relaxed lognormal clock (39). Priors on rates of substitution and clocks were set to defaults (as in Beauti 2.1). We ran analyses using a birth-death tree prior with incomplete sampling that takes into account the fact that we sampled a small fraction of all avian diversity (Supplementary Materials). Priors on calibration clades were set using exponential or log-normal functions, with parameters chosen as to mimic the empirical density distributions for the origin of clades (fig. S1 and table S1). We ran analyses without sequence data to evaluate the behavior of calibration priors (115). The tree topology of the RAG data set was estimated together with divergence times, except for some clades (Fig. 2B) that were constrained to match relationships supported by recent multilocus and genomic analyses (3, 28). The tree of the genomic analysis was fixed to the ExaML Total Evidence Nucleotide Tree topology of Jarvis et al. (3). Posterior samples were obtained from six independent Markov chain Monte Carlo runs of 80 million generations, each of which sampled every 8000 generations for the RAG data set, and from four runs of 60 million generations sampled every 6000 generations for the genomic data set. Adequacy of sampling and convergence was evaluated by examining the traces and effective sample sizes of likelihoods and parameters.

Ancestral area reconstruction

We divided the globe into eight regions that reflect both current global distributional patterns and major continental plates to ensure that area definitions are meaningful throughout the Cenozoic (Supplementary Materials). We assigned regions to the tips of the tree based on the current distribution of the clades they represent (data set 2). Regions were coded as different states of a single character, and clades distributed across multiple regions were treated as “polymorphic.” We reconstructed ancestral areas using Fitch parsimony, a nonparametric method that implies an island model (116). This method assumes geographic transitions rather than agglomeration of ancestral areas and is appropriate when transitions are rare (116, 117)—two conditions that may be more appropriate for large biogeographic regions because most avian species are restricted to single continents among which dispersal events are rare. In addition, more complex methods require complete phylogenies in which the tips of the tree are individual lineages (116, 118). In the present case, because the tips of the tree are solely representatives of clades, some of which are large, statistical methods may result in biased parameter estimates, compromising ancestral area estimation. Parsimony optimizations were obtained using the function ancestral.pars in the phangorn library (119) with the option MPR to obtain all possible parsimony reconstructions. Nevertheless, to explore the effect of methodological assumptions, we also used a probabilistic dispersal-vicariance likelihood model implemented in the library BioGeoBEARS in R (118) (Supplementary Materials). We also explored the effect of using an alternative tree topology (4) (Supplementary Materials). The fossil record indicates that many taxa, which are today restricted to tropical latitudes, were once present in the northern continents (29); therefore, including only extant taxa can bias ancestral reconstructions toward tropical areas. To minimize this bias, we added 25 fossil taxa to the tree representing Holarctic distributions for clades now restricted to the tropics (data set 2). Fossil taxa were attached at the midpoint between their age and the stem age of the clade to which they belong, using a new R function (available at https://github.com/evolucionario/fossilgraft), and treated as terminals in biogeographic reconstructions.

Diversification dynamics analysis

To analyze variation in diversification rates through time, we estimated rates for 5-Ma intervals using the results of the RAG data set and the function bd.ME.optim in the TreePar library (120). We analytically accounted for missing taxa using the option “groups” (121) in which we specified the actual number of species in taxa represented by each terminal of the tree; the richness of each terminal taxon was estimated using current taxonomic knowledge (122, 123) and can be found in data set 2. We also fitted environment-dependent birth-death models in which diversification rates are functions of a time-dependent environmental variable (87): we modeled different scenarios in which speciation or extinction rates vary through time as functions of changes in paleotemperatures. Deep-sea paleotemperatures were derived from a global compilation of benthic foraminifera oxygen isotope (δ¹⁸O) data (124) and new estimating equations (125). Rate dependency on global temperature could be linear [r(t) = r₀ + αT(t)] or exponential [r(t) = r₀e^αT(t), in which r₀ and α are estimated parameters and T(t) is a function that describes changes in temperature over time, in this case, a smooth-spline function with 50 degrees of freedom of deep-ocean temperature data].