Phylogenomics and biogeography of the world's thrushes (Aves, Turdus): new evidence for a more parsimonious evolutionary history

(a) Phylogenomic and divergence time estimates

We performed unpartitioned concatenated maximum-likelihood (ML) analyses in RAxML 8.2.9 [12] for two datasets: (1) a matrix including UCEs (two different thresholds of missing data, 75% and 85%) and (2) a matrix including the additional non-UCE loci (intragenic and untranslated regions). We assessed support for the best ML topology by performing 150 nonparametric bootstrap (BS) replicates using the autoMRE option in RAxML with the GTR GAMMA site-rate substitution model.

To avoid loci that could bias species tree analyses, we calculated the number of parsimony-informative sites (PIS) for each UCE locus and created a subset of UCE loci having a number of PIS in the upper quartile of the range (every locus with equal or more than 65 PIS). For the additional non-UCE loci, we included all markers recovered. We performed nonparametric bootstrapping by sites for each locus using RAxML with the GTR-GAMMA model of site-rate substitution and sorted the bootstrapped gene trees following Seo [13], an approach for estimating phylogeny and calculating bootstrapping using multilocus sequence data in the context of a distance method. To infer the species tree and calculate nodal support from the frequency of splits in trees constructed, we inputted the results to ASTRAL-II 4.7.8 [14]. We set ASTRAL to run 150 BS replicates of species trees. We also performed a species tree analysis using SVDquartets [15] for the UCE loci, including all loci recovered. By analysing fewer loci for the species tree, we were able to compare the results using different approaches.

To estimate divergence times, we used two late Miocene fossils from Hungary, correlated to the Mammal Neogene zone 13 (MN13) (7.2–4.9 Ma). These two fossils were compared to T. viscivorus, T. pilaris, and T. torquatus, and T. merula and T. philomelos, respectively. We tentatively interpret this as that they represent forms contemporary with the stem node of the main global radiation of thrushes, and use them for setting the minimum age for the clade including all Turdus species minus T. philomelos and T. viscivorus (figure 3). To fix a maximum age for this same node, we used the minimum age for the split between Regulidae and Certhioidea (17.2 Ma), based on fossil data and assigned in [16] (figure 3). As an additional secondary calibration point, we added the split between Muscicapidae and Turdidae, with a confidence interval of 22–13.5 Ma gleaned from [16] (figure 3). A UCE ML tree (the topology recovering highest support values overall) was time-calibrated under penalized likelihood using treePL [17]. The rate smoothing parameter was set to 100. The smoothing value was established using the cross-validation parameter available in treePL, we left the cvstart = 1000 (default) and cvstop = 0.1 (default).

(b) Biogeographic analyses

To infer the spatial range evolution for Turdus, we used Biogeography with Bayesian (and likelihood) Evolutionary Analysis in R Scripts—BioGeoBEARS, version 0.2.1 [18]. We used the time-calibrated phylogeny for Turdus calculated in treePL and coded each species as present or absent, according to their ranges following del Hoyo [1]. Eight biogeographic areas were defined based on the zoogeographical realms estimated in [19], with adaptations (figure 3): 1. Afrotropical (Afr), 2. Western Palaearctic (WP), 3. Eastern Palaearctic (EP), 4. Oriental (O), 5. Nearctic (N), 6. Panamanian (C), 7. Antilles (Ant), 8. South America (SA).

(a) Phylogenetic relationships

Our results using different thresholds of concatenated UCE datasets (75% and 85%) recovered T. viscivorus (historical samples GNM8238 and GNM172, from Sweden) and T. philomelos (U4329, from Belarus) as incremental sisters to all other Turdus species (figure 1 and electronic supplementary material, figures S2 and S3). In the main clade, all taxa from the Afrotropical region, including T. pelios, previously contentiously placed as part of the South American radiation, form a strongly supported monophyletic radiation (Node A, figure 1 and electronic supplementary material, figure S2), that is sister to a clade composed of all other species. The Palaearctic clade contains T. iliacus and T. merula, with strong support. In previous studies, these two species have been found in varying phylogenetic positions that suggest different biogeographic histories from other WP species [2,3]. We consider the following cut-offs for BS values: 50–74% weak; 75–84% moderate and greater than 85% as strong support (electronic supplementary material, figure S2 and figure 1).

• Download figure
• Open in new tab
• Download PowerPoint

While support for both the African and the Eurasian clades is strong, support for a monophyletic group [Node D, figure 1] in the Americas is weak. The important first splits within this clade are weakly supported, but the topology suggests a first divergence between Antillean taxa (a strongly supported clade including T. aurantius from Jamaica, T. plumbeus schistaceus from Cuba, and T. plumbeus ardosiacus from Puerto Rico) and remaining taxa (figure 1). The following basal divergences are all between taxa now inhabiting the Nearctic, Panamanian, or Antillean regions [19]. Although the divergence pattern is mostly weakly supported, a biogeographic structure may be discernible. A strongly supported clade comprises mainly species found in Central America (T. rufitorques, T. nigrescens, T. infuscatus), but also T. migratorius, a species widely distributed in North America. An exclusively Panamanian clade (strongly supported) consists of two subspecies of T. plebejus (rafaelenses from Nicaragua and plebejus from Panama). A large and strongly supported monophyletic Neotropical radiation (Node E; figure 1) includes the early (but poorly supported) sequential divergence of two species from the Antilles (T. swalesi from Hispaniola and T. lherminieri from Dominica) and is further subdivided into two main subclades corresponding to highland and lowland species, respectively.

The topology based on 49 additional loci (electronic supplementary material, figure S4) includes a concatenated matrix of 65 574 bp, using a taxon set with 105 individuals. BS support for a monophyletic Turdus is 100% (electronic supplementary material, figure S4). Turdus viscivorus is sister to the remainder of the clade with strong support. The next lineage to diverge is T. philomelos, which in turn is sister to the remainder of the clade. A basal split separates two clades with low BS support values—one large group includes species from the Palaearctic, Africa, Nearctic, Panamanian, as well as the Antillean regions, and is sister to a group which includes species from the Antilles (T. swalesi and T. lherminieri) and all samples from the continental South America region (electronic supplementary material, figure S4, weakly supported). Support for a monophyletic African radiation is strong, and T. swalesi and T. lherminieri are recovered (although with weak support) as part of the Neotropical radiation, as for the UCE dataset (Node E, figure 1 and electronic supplementary material, figures S2 and S3).

The analysis of the cytb dataset (electronic supplementary material, figure S5) recovers a monophyletic Turdus clade (excluding Otocichla mupinensis and Psophocichla litsitsirupa, which were recovered as sisters to the Turdus clade and which have both previously been included in Turdus), where T. viscivorus and T. philomelos are again the first lineages to diverge, being sisters to the remainder of the clade with high support. Within the main clade, most early-diverging nodes are unsupported, but otherwise, there are many similarities to previous studies [2,3]. African species are divided into six clades, with T. pelios sister to Antillean species and a Eurasian clade with the exception of T. merula and T. iliacus. Unlike previous studies [2,3], there is, for the first time, high support for a monophyletic South American clade based on a single mitochondrial locus. A number of unexpectedly deep or shallow divergences of potential taxonomic interest are highlighted in electronic supplementary material, figure S5 and discussed in the electronic supplementary material.

The UCE species tree (electronic supplementary material, figure S6, with gene trees of 483 UCE loci as input in ASTRAL) recovers a similar topology to the concatenated UCEs and cytb regarding T. viscivorus and T. philomelos. The placement of these species was also recovered, although in reverse order compared to concatenated UCEs and cytb. Here, T. viscivorus and T. philomelos are sisters to a strongly supported clade, which consists of two major groups: one including Panamanian, WP, Antillean, African, Eastern Palaearctic; and one South American group (Node D, figure 1), with a few Panamanian and Antillean species, plus two Eastern Palaearctic species, T. eunomus and T. naumanni. Support for the diverging order is weak, but several of the same groups as in figure 1 are recovered with strong support, e.g. T. plebejus, the Central American clade (T. rufitorques, T. nigrescens, T. infuscatus, T. migratorius), the African clade, and the Eurasian clade. Also, the Neotropical clade is similar to the topology in figure 1, including the division between the highland and lowland clades.

The UCE species tree (electronic supplementary material, figure S7, with 1931 UCE loci as input in SVDquartets) shows T. lherminieri from the Antilles as sister species to all other species of the crown group. Again, this analysis recovers most of the same clades as the other analyses performed in this study, while there is no resolution regarding the phylogenetic relationships among these large groups. The similarities between analyses are the recovery of a monophyletic African group, a monophyletic Eurasian group, a monophyletic Panamanian group, and a generally similar Neotropical clade.

The species tree built from the 49 additional loci (electronic supplementary material, figure S8; with gene trees of each additional locus as input) recovers a tree topology where the species T. viscivorus and T. philomelos are again the first to diverge, followed by two main groups (electronic supplementary material, figure S8): a moderately supported Afrotropical clade sister to the remainder of the thrushes from all non-African regions. In this tree, the main topology is insufficiently supported, but major clades such as the Eurasian clade including T. iliacus and T. merula are recovered, although support is weak. The American radiation is recovered as monophyletic, but again with weak support. Although there are differences in detail, major patterns in the Panamanian and Neotropical clades proposed by previous analyses shown here are corroborated.

(b) Divergence time estimates

The crown age of Turdus is estimated to ca 9.3 million years (Myr) (figure 3 and electronic supplementary material, figure S9), in the late Miocene. The first event of the major global radiation starts ca 7.2 million years ago (Ma), corresponding to the divergence of the African clade, rapidly followed by the divergence at ca 5.7 Ma between the Palaearctic and Oriental clades from the clade comprising species with distributions in the New World. Around 5.3 Ma, the New World radiation begins. Panamanian and South American clades first begin to diverge at ca 4.4 Ma and 4.9 Ma, respectively.

(c) Geographical ancestral range reconstruction

Including the ‘jump dispersal’ (j) parameter in the DEC model (figure 3), implemented in BioGeoBEARS, suggests a similar scenario as the DEC reconstruction (not shown), but causes the reconstructions at the nodes to become less noisy. Independent dispersals originating in the WP and moving towards the Eastern Palaearctic continuing to the Oriental region are suggested. The reconstruction of the New World radiation (including N, C, GA, and SA, figure 3 for abbreviations) suggests that it originated from a single dispersal event from the WP to the Antilles, although weak clade support makes this conclusion tentative (figure 1). From the Antilles, one or two dispersals to the Central America region are inferred, but the interpretation is complicated. From the Panamanian region, one lineage (T. migratorius) was recovered to have dispersed to North America. The radiation in South America is inferred to have originated directly from an ancestor dispersing from the Great Antilles. Separate dispersals back to the Central America region from South America are implied for the ancestors of T. assimilis and T. grayi (figure 3).

(a) Data considerations

The different evolutionary scenarios that have been suggested by studies of Turdus [2,3] propose conflicting scenarios of how they spread across the world. Most likely, the main reasons for the discrepancies are the differences in the amount and source of genetic data, combined with the stochasticity introduced by the short internodes caused by the high rate of speciation during the radiation. The UCE topology (ML) with 1931 concatenated loci (75% completeness of taxa) is the best currently available overview of the Turdus radiation, due to the highest number of loci sampled and the better resolution recovered to date (figure 1 and electronic supplementary material, figure S2). The concatenated UCE dataset shows reduced incongruence and improved resolution of the topology, as already documented for this kind of approach [20]. The more limited genetic sampling used by Voelker et al. [2] and Nylander et al. [3] makes these analyses more vulnerable to disagreement between the genealogy of the loci analysed and the species tree, as well as to the effects of various evolutionary processes, such as incomplete lineage sorting, gene flow, and hybridization. By exploring the different kinds of datasets obtained by genomic sampling, we were able to estimate a novel phylogenetic hypothesis based on 1931 UCEs (figure 1 and electronic supplementary material, figures S2 and S3). These data suggest a much simpler and parsimonious dispersal scenario compared to previous studies.

(b) Divergence time estimates

The age of the split between Oenanthe and Turdidae was estimated to be 22 Ma, corresponding to the upper extent of the confidence limit given for this split in figure 2 in Oliveros et al. [16]. This places the radiation of Turdidae within a time frame compatible with several recent studies [16,21–24]. Furthermore, the age of divergences between species is also compatible with the standard 2.1% per million years of divergence rate in the cytb gene [25].

• Download figure
• Open in new tab
• Download PowerPoint

(c) Phylogenetic implications and the historical biogeography of Turdus

Our time-calibrated phylogeny based on 1931 UCEs (figure 3 and electronic supplementary material, figure S9) sheds new light on the relationships among the majority of species in the genus Turdus. Previous studies suggested a series of trans-Atlantic dispersal events, both from Africa to the New World and from the Caribbean Islands to the Old World. Our findings challenge these hypotheses and that the African thrushes acted as a source for dispersal to the New World [3,26] and suggest a different scenario for the evolution of the genus Turdus. Although we cannot completely rule out the possibility of more than one trans-Atlantic dispersal (see Phylogenetic incongruences leading to alternative phylogenetic histories in the electronic supplementary material), the topologies allowing for more than one trans-Atlantic dispersal are unsupported in that respect, and in some details most likely spurious. Our statistically most well-supported results based on 1931 UCEs differ most importantly from previous hypotheses in that (i) the African radiation is monophyletic, (ii) Eurasian taxa did not originate from a west to east trans-Atlantic dispersal, and (iii) that the American taxa originate from a single dispersal from the Old World (WP) to the Americas, a result recovered also with low support in the analysis of 49 independent loci.

• Download figure
• Open in new tab
• Download PowerPoint

Ancestral lineages are here inferred to having occurred in the WP, and from there, dispersals to the Afrotropical, Eastern Palaearctic, and South American regions took place. The Afrotropical radiation is inferred to have been entirely in situ, with no indication that it had any role in the colonization of South America (figure 3 and electronic supplementary material, figure S9). The radiation in the Old World seems to have originated in the WP, with a subsequent spread towards the Eastern Palaearctic and Oriental regions on at least two separate occasions (figure 3 and electronic supplementary material, figure S9). We note, however, that genomic sampling of Palaearctic species is incomplete, making these conclusions tentative. Contrary to previous hypotheses, T. pelios is recovered in the African clade and is no longer suggested to be more closely related to species from the Caribbean (figures 2 and 3), thus contradicting both a close relationship with other taxa from South America [2], as well as a reversed trans-Atlantic dispersal event, as inferred by both Voelker et al. [2] and Nylander et al. [3]. T. philomelos and T. viscivorus represent ancient lineages that seem to have been unaffected by the processes driving the global radiation from the Palaearctic region (figures 1 and 3 and electronic supplementary material, figures S2–S9). T. iliacus is inferred to be sister to T. merula and both are part of the Palaearctic clade, with strong support. By contrast, both Voelker et al. [2] and Nylander et al. [3] recovered T. iliacus within the New World clade, as sister to T. plebejus, a Panamanian species.

We find evidence for the Antilles as the location of the first colonization event (figure 3). In the ancestral range reconstruction shown in figure 3, an Antillean distribution receives the highest probability, but Central America is also a possibility. Turdus thus seems to be an example of a long-distance dispersal from the Palaearctic where the birds first settled in the Caribbean region and then proceeded to radiate and colonize the mainland. These new data receive higher phylogenetic support than previous studies and generate a both more parsimonious explanation and a more logical scenario for the colonization of the Americas by Turdus.

The initial events in the colonization of the Americas seem to be centred around the Caribbean region. The most likely scenario according to the ancestral range reconstruction (figure 3) is a dispersal from a founder population in the Greater Antilles to Central America, resulting in one limited radiation including the ancestors of the T. rufitorques clade and in the T. plebejus lineage, followed by a reverse dispersal back to the Antilles from the Panamanian region by the ancestors of T. swalesi and T. lherminieri. In an alternative scenario, including two dispersals from the Greater Antilles to the Panamanian region, the ancestor of the clade including T. rufitorques would have had to disperse first, followed by the ancestor of T. plebejus, while the ancestors of T. swalesi and T. lherminieri would emanate from dispersal events within the Antilles.

Thrushes probably entered the Panamanian region by the time a land bridge was already fully or partly established [27–29]. There is an ongoing debate on the emergence of the Panamanian land bridge and the timing and processes underlying the Great American Biotic Interchange [29], which is partly due to difficulties in the interpretation of the exposure of land and shallow waters [28]. The closure of the Panamanian land bridge has had an important impact on the diversification of a variety of New World birds, as shown by Smith and Klicka [30]. Weir et al. [31], Töpel et al. [32], Antonelli et al. [33], Oliveros et al. [34], among others, have shown a marked increase in the rates of interchange of South and Central American birds, after or in connection with land bridge completion. These studies suggest that the route of most dispersal events was primarily south to north.

While a common pattern for birds, thrushes are likely to only have been marginally affected by the emergence of the Panamanian uplift. The colonization of South America is instead inferred to be a dispersal event directly from the Antilles, with no evidence of any connection to the Panamanian land bridge. This is not surprising, since these birds evidently dispersed to and between the Antillean islands, and could thus have dispersed between the Panamanian or Antillean regions and continental South America independently of a land bridge. This result shows similarities to findings from studies of other taxa. For instance, Sturge et al. [35] showed, based on mitochondrial DNA (cytb and ND2), that the colonization of South America by New World orioles (Icterus) is likely to have involved ancestors from the Caribbean, and a recent biogeographic reconstruction recovered reverse colonization by short-faced bats from the islands (e.g. Antilles) to South America [36]. Other examples include the genus Megascops (Aves, Strigidae) endemic to the New World with a likely origin in the Caribbean Islands [37].

By improving the taxonomic and geographical sampling within South America, we recovered for the first time a single, monophyletic South American clade within Turdus. By contrast, Nylander et al. [3] inferred three South American clades, and thus suggested independent colonization events. Our results regarding the South American radiation are similar to those of Voelker et al. [2] in this respect, apart from the African species (T. pelios) that was part of the American clade in their analysis but is part of the African clade in this study (figures 1 and 2 and electronic supplementary material, figure S2).

According to our phylogenetic hypothesis, the South American clade is divided into two main groups. One comprises mainly highland species (mostly Andean and Tepuian, but one Atlantic Forest species) and corresponds to an Andean clade already identified in Voelker et al. [2]. The second main clade includes species from the lowlands in Amazonia, Chaco, Cerrado, and Pampas, as well as Central America (T. grayi and T. assimilis). Paraphyly of T. olivater was detected, with T. olivater kemptoni (U4322, Venezuela) being recovered as sister to T. fuscater fuscater (U4285, Bolivia) and the Tepuian T. olivater roraimae sister to these two. This pattern has been documented in a previous study using a multilocus dataset [37]. Similarly, T. fuscater appears as paraphyletic, with T. fuscater gigantoides (U4287, Peru) being sister to T. chiguanco conradi (U4279, Peru) [38]. From continental South America, a number of independent reversals—dispersal events to both the Antilles and the Panamanian region—are inferred, which in the latter case could have been aided by the land bridge.

Future studies could focus on clarifying the geographical range evolution of Turdus within South America and along altitudinal zones, as well as investigating the influence of landscape and climatic changes on their past and current distribution and genetic diversity.