Phylogeny

We used the maximum clade credibility tree for all Coerebinae as recovered by Burns et al. (2014), as this represents the most complete phylogeny for this clade to date. To incorporate phylogenetic uncertainty, analyses were also run across 100 trees from the posterior distribution of trees. Trees were absolute time calibrated using secondary node calibrations provided by Barker et al. (2015) placing the root node at ∼8.3 mya. To confirm the validity of this date, we compared our age of the ancestor of Darwin's finches to divergence date estimations from previous studies for this clade. Our age of 2.6 million years since its divergence with the sister Tiaris lineage is consistent with other hypothesized ranges of dates that places the Darwin's finch radiation within the last 3 million years (Grant 1999, Petren et al. 2005).

Biogeographic Analysis

Zoogeographic regions for Central and South America have previously been described using vegetation structure and endemic species ranges by Parker et al. (1996). We modified these regions for our biogeographic analyses because BioGeoBEARS cannot reasonably analyze a data set that contains all 23 regions that were identified by Parker et al. (1996). Thus, we combined some regions based on adjacency and similarity, resulting in 8 regions: Mexico and Central America (combined as Central America), the Caribbean Islands, the Pacific Coast of South America, the Andes, Northern South America, Eastern South America, Cocos Island, and the Galápagos Islands (Figure 1). In addition, to directly examine the Caribbean versus South American mainland hypotheses of origin, and to replicate analyses of Burns et al. (2002), we also ran all analyses using a simplified geographic regime of just 5 regions: Central America, the Caribbean Islands, South America, Cocos Island, and the Galápagos Islands. All regions were reconstructed as a shapefile using ecoregion polygons from Olson et al. (2001) and joined in ArcGIS for visualization and distance measurements. We used a model-based approach implementing maximum likelihood to estimate ancestral ranges using the supermodel provided by BioGeoBEARS (Matzke 2013a). We compared the 3 commonly used models: Dispersal–Extinction Cladogenesis (DEC), a likelihood version of the Dispersal–Vicariance model (DIVALIKE), and a likelihood version of the BayArea model (BAYAREALIKE). The DEC model emphasizes range changes at speciation events (cladogenesis). During these events, one daughter lineage will always occupy only a single region of the ancestral range, either sympatrically or through vicariance. The DIVALIKE model is similar to DEC, but allows a daughter lineage to inherit more than a single region of the ancestral range through a large vicariant event. The DIVALIKE model does not, however, allow a daughter lineage to inherit a small range that is sympatric with the range of the other daughter lineage. Finally, the BAYAREALIKE model does not emphasize range changes at speciation events, but instead models range changes between speciation events (anagenesis) through range expansion and contraction. A more thorough overview of these models appears in Matzke (2013a). Additionally, all models were modified to include the +J parameter. This parameter allows for cladogenetic jump-dispersal events, where a new lineage disperses to a region outside of the region(s) occupied by its ancestor during speciation. We compared nested models using likelihood ratio tests (LRT) and non-nested models using AIC scores. The maximum region occupancy was set to 6 for the 8-area regime and 3 for the 5-area regime to match the most widespread taxon in the clade, the Bananaquit (Coereba flaveola). BioGeoBEARS also has the ability to test models for the support of distance as a modifier on dispersal probability. These models include an additional parameter X, where the probability of a dispersal event between 2 regions is multiplied by the distance between regions taken to the X power. For example, if model optimization results in X = −1, the probability of dispersal between 2 regions is weighted by inverse the distance (i.e. the farther away 2 regions are, the lower the weight on a dispersal event). Van Dam and Matzke (2016) demonstrate the insight this approach can lend to biogeographical reconstructions, and suggest it may be particularly important in coarse-scale analyses that include oceanic islands, as is the case for Coerebinae. Thus, we also tested the above models including distance with the +X parameter. We incorporated distance by including a matrix of pairwise distances between regions and used the shortest distance between regions to represent the most conservative modifier on dispersal probability. Because distance modifies dispersal weight as an exponent, the model cannot incorporate a distance of zero. Thus, regions that are adjacent to each other were assigned a value of 1 m so that it was equivalent to a non-weighted dispersal probability. The distance measurements were made with the constructed shapefile using the Near Table function of ArcGIS with the geodesic setting to correct for the curved surface of the earth. X was incorporated as a free parameter and model comparisons were made using AIC scores and LRTs.

FIGURE 1.

Map of biogeographic areas used in BioGeoBEARS eight-area regime. Five-area regime combined the four regions of South America mainland that are shown here.

8-area Regime

Under the 8-area regime, model comparison of DEC, DIVALIKE, BAYAREALIKE, as well as all models with the +J parameter, revealed the DEC+J model received the highest likelihood score and the best AIC score of the 6 standard models (Table 1A). However, the inclusion of a distance matrix in the DEC+J model was significantly supported by our data, as determined by an LRT (P < 0.001), and suggests DEC+J+X was the overall best fit model. Of the 9 models compared, this model received 64% of model weight after correcting for sample size (AIC_c). Additionally, the DEC+J+X model performed best when the X parameter was allowed to be free and resulted in a negative dispersal modifier exponent, indicating that distance had a slight inverse effect on dispersal probability. The J parameter also performed best when allowed to be free, and indicates that allowing for jump-dispersal events at cladogenesis was important in estimating the ancestral ranges for this clade. Under this model, range changes can occur at speciation events (cladogenesis) or between speciation events (anagenesis). To model both types of changes, ancestral ranges are estimated just before speciation, depicted on the nodes of a phylogeny, and just after speciation, depicted on the corners of the branches immediately following a node. Thus, differences in range estimation between a node and a corner represent cladogenetic range changes, and differences in range estimation between a corner and a node represent anagenetic range changes.

TABLE 1.

Summary of model results for all 9 models compared under the (A) 8-area regime and the (B) 5-area regime. Models with +J indicate those allowing for jump dispersals. Models with +X incorporate distance as a dispersal modifier according to the pairwise distance matrix for all regions. Best-supported model shown in bold. K = number of parameters

The DEC+J+X model was able to reconstruct ancestral ranges with high probability for almost all nodes (Figure 2), except for the root node. However, the combined probability that this state is or at least includes the Caribbean is extremely high (99.5%). The cladogenetic event at this root node resulted in the single taxon lineage of what is now the Bananaquit, and a Caribbean-restricted ancestor to the rest of the Coerebinae clade. The subsequent speciation event resulted in the sister relationship between Tiaris olivaceus and the rest of Coerebinae. This speciation event likely occurred within the Caribbean, with T. olivaceus then dispersing anagenetically to mainland Central and South America after the divergence event.

FIGURE 2.

Ancestral range estimation for Coerebinae modeled under DEC+J+X using the 8-area regime. Node a corresponds to the ancestral range just prior to the cladogenetic event that led to Tiaris fuliginosus, Tiaris obscurus, and the Darwin's finch clade as sister to the Melanospiza richardsoni and Tiaris bicolor clade. Corner a corresponds to the ancestral range just after the above cladogenetic event. Node b corresponds to the ancestral range just prior to the cladogenetic event that led to the Darwin's finch clade, and the 2 Tiaris species sister to Darwin's finches. The 2 highest range probabilities are shown for each of the 3 ancestral ranges labeled in the tree. Dashed line corresponds to the oldest estimated age of the current formation of the Galápagos Islands. Clade label on right indicates species in the Darwin's finch radiation.

Many of the ancestral ranges after this event were also estimated as Caribbean with high probability. The next event not estimated as Caribbean came after the cladogenetic event that resulted in the Melanospiza richardsoni and Tiaris bicolor clade as sister to Tiaris fuliginosus, Tiaris obscurus, and the Darwin's finch clade. This event was most likely a jump-dispersal event from the Caribbean to the Galápagos (Figure 2 Node a and Corner a), providing support for the hypothesis of a Caribbean origin to the Darwin's finches. This event was followed by an additional jump dispersal at the next speciation event that resulted in the recolonization of the mainland by the ancestor of Tiaris fuliginosus and Tiaris obscurus, while the ancestor to Darwin's finches remained on the Galápagos. The ancestral state of these Tiaris is less certain, but was recovered as most likely Andean. One final event was recovered within the Darwin's finch clade which corresponds to the almost certain jump dispersal from the Galápagos to Cocos Island by Pinaroloxias inornata (Cocos Finch). Outside of the Darwin's finch clade, one additional range expansion was reconstructed, indicating Tiaris bicolor dispersed anagenetically to Northern South America.

Determining support for the origin of Darwin's finches comes from examination of states from 3 event locations on the tree: the range before and after the cladogenetic event leading to the Tiaris fuliginosus, Tiaris obscurus, and the Darwin's finch clade as sister to the Melanospiza richardsoni and Tiaris bicolor clade (Figure 2 Node a and Corner a, respectively), the states before and after the cladogenetic event corresponding to the split between the Tiaris fuliginosus and Tiaris obscurus clade and the Darwin's finch clade (Figure 2 Node b), and the anagenesis between these 2 events. Ancestral range estimation suggests with very high probability (>80%) a Caribbean-restricted ancestor to these Tiaris, Melanospiza, and Darwin's finches (Figure 2 Node a). Single-area states comprise 99.7% of the probable states post-cladogenesis at Node a (Figure 2 Corner a), indicating this ancestral range likely did not span more than one region. The Galápagos is recovered as most probable at almost 60%, the Andes at just over 20%, and all other states less than 10% probable. Range reconstruction prior to the next speciation event (Figure 2 Node b) also recovered the Galápagos as the most probable, suggesting no anagenetic range change between these 2 speciation events.

Running the same model across random trees from the posterior distribution produced largely similar results as the maximum clade credibility tree. Because ancestral ranges prior to the divergence of the Darwin's finch and Tiaris obscurus and T. fuliginosa clade were all inferred as Caribbean, changes in topology outside of this clade did not change the interpretation of the ancestral ranges. Only one topology was recovered that was different with respect to Darwin's finches and the 2 sister Tiaris species. These species were always recovered as a clade, but the Tiaris species were occasionally recovered within the Darwin's finch radiation (12% of trees examined). This change in topology resulted in very strong support for a Caribbean to Galápagos dispersal; however, this topology is unsupported overall and has not been recovered in other studies of Darwin's finches that have included these Tiaris species. Thus, it is unlikely this ancestral range reconstruction reflects a probable scenario.

5-area Regime

Running the same analyses under a 5-area regime produced different results than the 8-area regime. The best-supported model under the 5-area regime remained DEC+J+X (Table 1B); however, this model only received 32.72% of model support based on AIC_c score, with the DIVALIKE+J+X model receiving the next highest at 19.26%. Despite altering the geographic delineations, the data supported inclusion of distance as a modifier on dispersal, and best fit a model that allowed for jump-dispersal events. Under this regime, the J parameter was estimated at 0.16 and the X parameter was negative, again demonstrating that distance has an inverse effect on dispersal probability.

As with the 8-area regime, the ancestral range at the root node was unable to be reconstructed with any certainty, but indicates a strong probability it included the Caribbean (98%). After the divergence event at the root, all ancestral ranges were reconstructed as Caribbean restricted until the divergence between Tiaris bicolor and Melanospiza richardsoni, and the clade containing Darwin's finches and the 2 sister Tiaris species (Figure 3 Node c). Similar to the 8-area regime, the reconstruction under this model showed a high probability of a Caribbean-restricted ancestor to these Tiaris, Melanospiza, and Darwin's finches (Figure 3 Node c). Although the probability under the 5-area regime was ∼20% lower, reconstructions under both regimes estimate this probability to be over 50%. The ancestral range just after this event that led to the lineage containing Galápagos finches (Figure 3 Corner b) was estimated as either Galápagos or South American mainland with almost equal probability. Thus, ancestral range reconstruction under the 5-area regime suggests equal probability between the 2 hypothesized areas of origin for Darwin's finches. While the South American mainland was the most probable state at this corner (53%), the Galápagos received an almost equal probability (46%). No anagenetic change was reconstructed along this ancestral lineage with reasonable probability, with the Galápagos and the South American mainland again reconstructed with about equal probability at Node d (Figure 3). This final speciation event at Node d (Figure 3) separating Darwin's finches was reconstructed as an additional jump-dispersal event either back to mainland South America, or to the Galápagos from mainland South America. Determining the direction of this jump-dispersal would require a more definitive range estimation at Node d (Figure 3).

FIGURE 3.

Ancestral range estimation for Coerebinae modeled under DEC+J+X using the 5-area regime. Node c corresponds to the ancestral range just prior to the cladogenetic event that led to Tiaris fuliginosus, Tiaris obscurus, and the Darwin's finch clade as sister to the Melanospiza richardsoni and Tiaris bicolor clade. Corner b corresponds to the ancestral range just after the above cladogenetic event. Node d corresponds to the ancestral range just prior to the cladogenetic event that led to the Darwin's finch clade, and the 2 Tiaris species sister to Darwin's finches. The 2 highest range probabilities are shown for each of the 3 ancestral ranges labeled in the tree. Dashed line corresponds to the oldest estimated age of the current formation of the Galápagos Islands. Clade label on right indicates species in the Darwin's finch radiation.

Regime Comparisons

Our comparisons show strong conflict in model selection and ancestral range reconstruction based on 2 geographic regimes. These 2 sets of results confer opposing conclusions for the biogeographic origin to Darwin's finches. Under both geographic regimes, model selection analyses suggest range changes occur most frequently during speciation events, and often result in a narrowly distributed daughter lineage. Model selection also indicates that these range changes at speciation events involve a daughter lineage dispersing to a new region—one that is not occupied by the ancestor. This result appears to follow hypotheses that jump-dispersals and the incorporation of distance-informed dispersal probabilities are important for island systems (Matzke 2013b, Van Dam and Matzke 2016). Our 8-area regime represents a more detailed geographic division in an attempt to provide the model with meaningful geographic boundaries. Under this geographic delineation, the model estimated ∼60% probability of a Galápagos restricted ancestor to Darwin's finches and the 2 sister Tiaris species (range at Figure 2 Corner a), which was almost 40% higher than any other range estimated, and thus represents support for a direct Caribbean–Galápagos dispersal event. However, the reconstruction of this event under our 5-area regime produced different results. Unlike the 8-area regime, the dispersal of this ancestor was reconstructed as approximately equal probability of a Caribbean–Galápagos dispersal as a South American mainland–Galápagos dispersal (Figure 3, Corner b and Node d). Conservative interpretation of these data suggests, given the current phylogeny and biogeographic methods, there remains an equal probability of a Caribbean origin or a South American mainland origin to the Darwin's finch radiation. These results differ from Burns et al. (2002) who recovered strong support for a Caribbean origin.

The difference in results produced by the 2 geographic regimes used here highlights the potential problems of requiring a priori regional assignment. Regional delineation may not always be straightforward, and subjectivity can introduce bias into model estimations. While a simplified regime may be able to more explicitly examine one particular dispersal event, it may also produce unintended effects during parameter estimation. In general, simplifying a geographic area may mask meaningful biogeographic events that happen on a smaller scale. In our simplified regime, reducing South America from four regions to one has the potential to underestimate jump-dispersals by overemphasizing the prevalence of regionally sympatric speciation. This would then lower the weight of any given long-distance dispersal event. Conversely, over-dividing, or inaccurately dividing geographic regions, may produce the opposite effect and could overemphasize dispersal, anagenetically or cladogenetically. Other recent studies of South American biogeography have delineated geographic boundaries similar to our 8-area regime (e.g., Batalha-Filho et al. 2014, Lima et al. 2017). The results presented here suggest that models such as those run through BioGeoBEARS are sensitive to these delineations and should be considered during future investigation of biogeography in this area. The recognition of this problem is not novel, and is an explicit assumption of any ancestral range reconstruction, but this study demonstrates some of the potential problems that can result. One way to overcome these biases could be to take a more objective approach to regional delineation by incorporating explicit range data. These data could be analyzed using a statistical, model-based approach and use model comparison to select the best fit geographical regime. A similar approach implementing hierarchical clustering algorithms was used by Kennedy et al. (2017) in a biogeographic study of corvids, but was dependent on the inclusion of a large number of species.

Despite this uncertainty, the Caribbean remains a plausible biogeographic origin to the Darwin's finch radiation. The Galápagos Islands are located on the Nazca plate, and are the product of a volcanic hotspot (Christie et al. 1992). The young hypothesized age of the current islands (4–5 mya; Hickman and Lipps 1985), along with the Nazca plate's eastern movement, indicate the Galápagos Islands were never part of the mainland, and have never been located closer to mainland South America. Although the greater distance between Caribbean islands and the Galápagos might make a direct Caribbean-to-Galápagos dispersal seem unlikely, birds have previously demonstrated the ability to disperse long distances. For example, solitaires in the genus Myadestes have colonized both the Hawaiian Islands and Greater Antilles from Mesoamerica (Miller et al. 2007). Moreover, other Galápagos taxa exhibit affiliations with Caribbean taxa (Grehan 2001), indicating a Caribbean origin is plausible. Arbogast et al. (2006) used mitochondrial sequence data to investigate the colonization of the Galápagos Islands by mockingbirds. Despite relatives in nearby Ecuador, an alternative relationship was recovered linking Galápagos mockingbirds to Caribbean and North American species, suggesting dispersal propensity outweighed geographic proximity with regard to colonization probability. Another well-known bird of the Galápagos, the American Flamingo (Phoenicopterus ruber), exhibits a similar Caribbean-to-Galápagos connection, showing affinities to the birds of the West Indies rather than the geographically closer species of the Andes (Thornton 1971). This pattern has also been found in non-avian taxa such as snakes (Antillophus), isopods (Nesophiloscia), sponges (Rhabderemia), and moths (Oxydia) (Grehan 2001). Finally, neither of our analyses provided exclusive support for a South American mainland origin. Thus, the Caribbean remains as a likely source for the origin of Darwin's finches.

Between the Galápagos and the Caribbean lies the Panamanian Isthmus, with a complex geological history (Montes et al. 2012, Coates and Stallard 2013, Prevosti and Forasiepi 2018). While the closure of the Panamanian Isthmus may have facilitated dispersal through a number of mechanisms, including island-hopping or climate shifts, the date of the isthmus closure is still debated (O'Dea et al. 2016, Molnar 2017), and some dates indicate this closure occurred prior to the diversification of the Coerebinae clade. Thus, the potential role the formation of the isthmus played in the dispersal to the Galápagos remains unclear.

Early Diversification of Coerebinae

Due to the restriction of many early diverging taxa to the Caribbean, ancestral ranges early in the diversification of this clade were reconstructed as Caribbean-restricted regardless of the geographic regime or model used. Although the ancestral range at the root was not reconstructed with any certainty, there was a high probability this range at least included the Caribbean under both regimes. When considering this probability, the role of this region in the origin and the early diversification of this clade becomes apparent. These results support conclusions drawn by Burns et al. (2002) as well, who highlighted the importance of these islands in the evolution of Coerebinae.

While the ancestral range at the root may have been Caribbean restricted or widespread, examination of species' current distributions reveals the early diversification of this clade was likely driven by a series of specializations within small populations across the Caribbean. These repeated colonization events imply the biogeographic diversification may have been dependent on a high degree of “evolvability” within this clade. Previous examination of bill size and shape within Coerebinae species report a much higher degree of variation than in other close relatives (Burns et al. 2002, Tokita et al. 2016) and provide a likely candidate trait that may have been responsible for colonization success. Additional studies of bill development within Coerebinae suggest a large degree of flexibility in the genetic controls of size and shape variation (Abzhanov et al. 2004, Abzhanov et al. 2006, Mallarino et al. 2012). The large flexibility of this genetic component to bill variation may help explain how large changes in the bill can occur in a rapid radiation such as Darwin's finches. The evolvability of these traits likely played an integral role in the success of each colonization event and, as a result, heavily shaped the biogeographic history of this subfamily.

In addition to bill, the species of Coerebinae exhibit a propensity for dispersal much higher than other clades of birds in the tanager radiation. Of the 374 species of tanager, only 34 are native to isolated islands, with 27 of them being species of Coerebinae. This suggests a very strong phylogenetic component to island distributed species, and thus an underlying genetic predisposition for longer dispersal tendencies. A similar conclusion was drawn by Arbogast et al. (2006) in describing dispersal patterns seen for 3 species of mockingbirds. These species are distributed throughout North and South America, with one species colonizing Socorro Island off the west coast of Mexico. They suggest the ancestor to this lineage may have exhibited a tendency to wander much wider distances than other species, thus allowing it to come into contact with farther, more isolated islands. The same may be true for Coerebinae species, as only 2 species of this clade are not present on islands (though they exhibit large distributions across the mainland). An ancestor to this lineage that was more likely to wander would be much more likely to encounter oceanic islands such as those of the Caribbean, with the variance in bill shape providing the means for successful colonists to adapt to open ecological niches. An underlying genetic component would allow this dispersal propensity to persist throughout subsequent speciation events until the long-distance dispersal to, and proliferation throughout, the Galápagos Islands.