Taxonomic Sampling

We obtained loans of frozen or ethanol-preserved tissue samples for 4 Chlorophonia species and for 24 of 27 Euphonia species currently recognized by most taxonomic authorities (Gill and Donsker 2019, Chesser et al. 2018, Remsen et al. 2018) from the following natural history collections in the United States: the Louisiana State University Museum of Natural Science (LSUMZ), the Field Museum of Natural History (FMNH), the Yale Peabody Museum (YPM), and the United States National Museum of Natural History (USNM). We obtained 2 samples for E. xanthogaster that represent allopatric and phenotypically disparate subspecies, quitensis and ruficeps, of this widespread South American species. For the 3 Euphonia species (chalybea, concinna, and trinitatis) and 1 Chlorophonia (occipitalis) species where fresh tissues were not available, we obtained toe pad subsamples from study skin specimens in the following collections: FMNH, the Bell Museum at the University of Minnesota (MMNH), and the Museum of Vertebrate Zoology at the University of California, Berkeley. Additionally, we utilized tissues from Coccothraustes verpertinus and Spinus tristis from the MMNH collection, plus existing UCE phylogenomic sequence data for Fringilla montifringilla and Emberiza citrinella (Moyle et al. 2016, Oliveros et al. 2019), to include these taxa as outgroups spanning relationships within, and sister to, the family Fringillidae, respectively. A complete list of the museum specimens included in this study and their collection localities is shown in Table 1.

DNA Extraction and Illumina Sequencing

For all but 3 of the tissue specimens utilized in this study, we extracted total DNA from a tissue subsample using a DNeasy Blood & Tissue Kit (Qiagen, Valencia, California, USA) following the manufacturer's instructions. We extracted total DNA from the remaining 3 tissue specimens using a phenol-chloroform extraction and Phase Lock Gel following standard protocols from QuantaBio (Beverly, Massachusetts, USA). To extract DNA from the 3 toe pad subsamples, we washed the samples twice in 100% EtOH, soaked them in 100% EtOH for 15 min, and repeated this washing and soaking procedure for 70% EtOH and for ddH₂O, but with a 30-min soak in ddH₂O. We then followed a modified Qiagen DNeasy procedure. Specifically, we added an additional 20 µL of Proteinase K and 30 µL of dithiothreitol at 100 ng µL^–1 concentration during the digestion, and these samples were incubated on a shaker at 56°C set to 30 rpm overnight. Afterwards, we eluted the samples twice in 50 µL of elution buffer instead of 200 µL, and these separate elutions were pooled and reduced by vacuum drying to a volume of 60 µL. We assessed the concentration of all DNA extracts using a QuBit dsDNA high-sensitivity assay and the quality and range of fragment sizes were assessed by running 2 µL of genomic DNA on a 1% Tris-borate EDTA agarose gel.

Target capture and dual indexing for each specimen, except for C. occipitalis, was performed by Arbor Biosciences (Ann Arbor, Michigan, USA) using the Tetrapod 5k 1.8 MYbaits probe kit (ultraconserved.org) to target ∼5,000 UCE loci from each sample of genomic DNA (Faircloth et al. 2012). Assembly of individual baited libraries for each specimen, enrichment reactions for 7 pooled libraries (with 5 specimens each), and library equilibration, was also performed by Arbor Biosciences following their MYreads protocols. The resulting equilibrated library was sequenced in one half lane of an Illumina NextSeq 550 (San Diego, California, USA) run to obtain 150 base pairs (bp) paired-end reads at the University of Minnesota Genomics Center (St. Paul, Minnesota, USA). For C. occipitalis, we targeted and sequenced ∼1 kb of cytb by polymerase chain reaction (PCR). We amplified this locus in 5 segments due to the degraded state of the DNA extracted from the toepad subsample for this specimen, and each reaction had a volume of 25 µL and otherwise followed the procedures described by Barker et al. (2008). After performing an additional round of PCR, they were sequenced at the University of Minnesota Genomics Center. Raw read data from this study were accessioned in the National Center for Biotechnology Information (NCBI) Short Read Archive (BioProject ID: PRJNA604822).

Contig Assembly and Sequence Alignment

The resulting Illumina reads were trimmed to remove adapters and low-quality bases using illumiprocessor 2.0.9 (Faircloth 2013). We then used phyluce 1.5 to assemble sequence contigs and for all remaining steps prior to constructing the actual phylogenies (Faircloth 2016). We performed 2 contig assemblies of our read data using both Trinity 2.1.1 and SPAdes 3.12.0 to compare assembly quality and assess efficiency in mitogenome assembly. These contigs were then screened against the Tetrapod 5K UCE probe set to identify and then extract contigs that mapped to UCE element loci. We extracted these loci from the FASTA data into an incomplete sequence matrix and aligned these sequences using MAFFT 7.407. Using scripts in R 3.5.1 (R Core Team 2013), we then quantified these data to determine the number of UCE-matching contigs per specimen, the distribution of sequence lengths from these loci, and the total number of bp to be included in subsequent analyses. This incomplete sequence matrix was pared down to obtain an additional complete sequence matrix that only included UCE loci sequenced for every specimen in the dataset.

ML Phylogeny Inference

We converted the incomplete sequence matrix to a concatenated PHYLIP file to perform concatenated ML inference of the phylogeny in RAxML 8.2.9 (Stamatakis 2014). These analyses were preceded by running a UCE-specific partitioning algorithm, sliding window site characteristics of entropy (SWSC-EN; Tagliacollo and Lanfear 2018). Given the distinct structure of UCE loci, in which the core sequence is strongly conserved between species and the flanking regions vary considerably, this algorithm divides a concatenated alignment of UCE loci by splitting the sequence of each locus into 3 partitions representing the core and 2 flanking regions, and then partitioning the alignment given this heterogeneity in sequence variation. We modified this method slightly for the complete UCE alignment by increasing the window size from 50 to 100 bp to increase the minimum potential size of partitions in the final analysis. The output of this method for both the incomplete and complete sequence matrices was then used as input to optimize the partitioning scheme and assign models to the resulting partitions in PartitionFinder2 (Lanfear et al. 2017). Optimizing the partitioning scheme for the incomplete sequence matrix computationally failed due to the sheer number of possible partitioning schemes, so we manually set a scheme that partitioned the total aligned sequence into flanking regions, cores, and unsplit UCE loci. We then ran RAxML using the manual partitioning scheme for the incomplete sequence matrix and the best-supported partitioning scheme for complete sequence matrix and bootstrapped both analyses 100 times with the “-f a” option.

TABLE 1.

Taxonomic sampling of this study. Collection abbreviations follow the main text, with the addition of the University of Kansas Biodiversity Institute (KU). An asterisk (*) indicates the sample used was a toe pad instead of a tissue subsample. A caret (^) indicates a female specimen. A tilde (∼) indicates that previously generated UCE sequence data were used from this specimen. Group # corresponds to the Euphonia species groups following Isler and Isler (1999)

We performed a ML coalescent-based analysis using the incomplete sequence matrix in 2 steps. First, we generated individual PHYLIP files for each UCE locus in the incomplete sequence matrix, removed any loci that contained sequences for fewer than 4 taxa, and inferred gene trees for every individual locus using RAxML. Second, we used these gene trees as input to estimate a species tree for these data using ASTRAL-III 5.6.3 (Zhang et al. 2018). The nodal support and quantification of gene tree-species tree discordance for this coalescent-based topology were assessed using several metrics, including the localized posterior probability (LPP) computed by ASTRAL. In addition, we calculated internode certainty (IC) scores for each branch in the species tree as implemented in RAxML, which calculates the frequency of gene tree topologies that differ from that of the species tree with regard to individual branches and whether these different topologies occur with equal frequencies as expected in the absence of gene flow. IC scores near 1 indicate high concordance between the species tree and gene trees, scores near 0 indicate equal frequencies of discordance with the species tree among alternative gene tree topologies, and IC scores near –1 indicate strong discordance with the species tree and highly skewed representation among the conflicting gene trees (Salichos et al. 2014). Further, we calculated 2 node-specific concordance factors as implemented in IQ-TREE (1.7-beta9), gene concordance factors (gCF) and site concordance factors (sCF) (Nguyen et al. 2015, Minh et al. 2018). These values quantify the proportion of informative gene trees and sites in the aligned sequences for a given node in the species tree that are concordant with that node, with both values ranging from 0 (high discordance of gene trees and informative sites) to 1 (high concordance of gene trees and informative sites). Next, we also assessed support for this topology by bootstrapping the gene tree topologies. We generated 100 bootstrap replicates by sampling all the gene tree topologies with replacement and, for each replicate, inferred a species tree of those sampled gene trees using ASTRAL. We then utilized this distribution of 100 bootstrapped species trees as the distribution of “gene trees” to compare with the true species tree topology and to compute a set of bootstrap values using the “f -b” option in RAxML. Finally, to further test the robustness of the species tree topology obtained from the incomplete sequence, we performed ASTRAL analyses on the complete sequence matrix and a 50% complete sequence matrix (containing at least 19 taxa), and we calculated concordance factors and internode certainty values as nodal support for these species trees. We then visually compared the resulting topologies to assess any conflicting nodes.

Mitochondrial Phylogenomics

In addition to analysis of the UCE data, we assembled and analyzed complete or near-complete mitochondrial genomes for the taxa that were sampled by target enrichment. Off-target mitochondrial reads were assembled in 1 of 2 ways. First, Trinity and SPAdes assemblies were inspected for single large contigs of ∼16 kb in length that were good candidate mitogenomes. A set of mitochondrial genes from a related passerine (Spinus spinus, GenBank accession  NC015198) were mapped to these contigs in Geneious 10.2.4 (Biomatters, Auckland, New Zealand), and gene presence and length and open reading frames were confirmed. Where possible, these contigs were circularized and their coordinate systems reset to start at the first base of tRNA-Phe. For samples without single large contigs in their corresponding assembly, we used Geneious to map reads from those samples to 1 of 3 complete, high coverage euphoniine mitogenomes (C. cyanea, Euphonia musica, or Euphonia finschi) most closely related to that sample based on UCE results. Assemblies with <90% estimated mitogenome coverage were excluded from further analysis. The consensus of the remaining assemblies was extracted and examined as described above, and these mitogenomes were annotated and submitted to GenBank (accession IDs: MT063155–MT063184).

Assembled euphoniine and outgroup mitogenomes were aligned using MUSCLE 3.8.425. An additional outgroup (Agelaius phoeniceus, GenBank accession  JX516062) was included in this analysis to facilitate mapping of 48 standardized character sets (RNA stem and loop, protein coding gene codon positions, etc.; see Powell et al. 2013) onto the alignment. Model fitting for this alignment was performed recognizing these character sets, using PartitionFinder2, allowing only the GTR and GTR+G models with proportional branch lengths, using greedy optimization with AIC_c as an optimality criterion (Yang 1994). The mitogenome data were analyzed under the best fit partitioning and models using RAxML 8.2.11, calculating 1,000 bootstrap replicates.

Inference of a Time-calibrated Tree

While we consider the species tree topology inferred by ASTRAL to be the best estimate for the entirety of the UCE sequence data, ASTRAL is unable to estimate terminal branch lengths and only estimates internal branch lengths in coalescent units, rather than in units of absolute time. Another issue to overcome was the fact that C. occipitalis was not included in this tree. To estimate ultrametric, time-scaled branch lengths and node ages across this topology and to include all species, we inferred a tree using the partitioned, concatenated alignment of the complete UCE sequence matrix and the cytb sequences as one additional partition in BEAST 2.5.2 (Bouckaert et al. 2014). We extracted a cytb sequence from each complete mitochondrial genome and aligned these sequences with that from C. occipitalis using MUSCLE 3.8.425. Prior to the BEAST2 analysis, we performed an additional RAxML analysis on this new partitioned concatenated alignment using the “-g” option, which establishes a fixed starting tree, the ASTRAL species tree in this case, and optimizes the location of previously unsampled tips on this topology. Doing so allowed us to place C. occipitalis on the species tree topology. The tree resulting from this analysis was made ultrametric using the chronos function in the R package ape (Paradis et al. 2004), and we then used this resulting ultrametric phylogeny as the starting tree constraint in the BEAST2 analysis. We constrained the tree search in BEAST2 such that it would optimize branch lengths for this topology given the complete-partitioned sequence alignment with cytb and not change the topology. We calibrated the absolute age of this tree by assigning a log-normal prior to the root node using information from a recent UCE-based phylogenomic study of all families in the order Passeriformes using 13 fossil calibrations and overlapping sequence data for 2 species shared between that study and the present one, Fringilla montifringilla and Emberiza citrinella (Oliveros et al. 2019). This prior separates Emberiza from the family Fringillidae and was set to have a mean of 2.8365 and a standard deviation of 0.0745 to match the highest posterior density interval of the node age and uncertainty inferred from the family-level maximum-clade credibility (MCC) tree. We configured the run to have a strict clock model linked across all partitions with a uniform prior set with the default values, and we set identical HKY+G substitution models for each individual partition in the alignment with the default distributions set by BEAUti 2.5.2 (Hasegawa et al. 1985). We originally specified a birth-death tree model with gamma-distributed priors of α = 0.001 and β = 1,000 for both the birth and death rate and ran this analysis for 50,000,000 generations on the CIPRES Science Gateway (Miller et al. 2010).

This analysis reached convergence for every parameter in the model with high effective sample sizes except for the prior and the posterior, both of which had low effective sample size (ESS) values of only 8, which is a symptom of the model prior not having enough time to thoroughly mix and sample the entirety of the possible parameter space. Therefore, we modified our control file for the analysis by changing the birth-death tree model to a Yule model with a default uniform prior and by increasing the run length to 100,000,000 generations. We sampled the trace of this analysis every 10,000 generations, and sampled trees every 100,000 generations to generate a final posterior sample of 10,000 observations and a posterior distribution of 1,000 trees. Lastly, we determined whether this analysis had achieved convergence after its completion by analyzing its log in Tracer 1.7 (Rambaut et al. 2018).

Biogeographic Reconstruction

We obtained a MCC tree from the posterior distribution of trees sampled by the BEAST2 run in TreeAnnotator 2.5.2, specifying a 10% burn-in and obtaining median node heights, to reconstruct the biogeographic history of this lineage. We scored the biogeographic range of each species on this MCC tree as presence or absence in the following 4 biogeographic regions: North America (Mexico south to Darién, Panama); South America (Darién, Panama south to Brazil, Paraguay, and Argentina; including Trinidad); the Caribbean (any island in the Greater or Lesser Antilles, excluding Trinidad); and the Eastern Hemisphere (Eurasia, given the scope of sampling in this study). Only 3 species' ranges included more than 1 of these biogeographic regions given this scoring strategy (E. fulvicrissa, E. laniirostris, and E. minuta), so they were assigned as present in 2 regions, North and South America.

We compared the fit of 3 likelihood-based models of biogeographic range evolution, DEC, DIVALIKE, and BAYAREALIKE, using the BioGeoBEARS package in R (Matzke 2013a, b). All 3 models optimize parameters for anagenetic rates of range expansion or contraction (d or e) along the branches of a phylogeny but differ in their parameterization of cladogenetic range changes through vicariant or sympatric events. The BAYAREALIKE model only parameterizes simple sympatric cladogenetic events with the y parameter, whereas both the DEC and DIVALIKE models account for both sympatric and vicariant range change events with the parameters y and v. These 2 are further differentiated in that the DEC model contains an additional parameter, s, which defines sympatric speciation in a subset of the lineage's total range, whereas the DIVALIKE model allows for both narrow-range and widespread vicariance events. Additionally, BioGeoBEARS can modify these 3 models by including an additional parameter, j, which models the probability of founder-event speciation, but recent work has debated whether the likelihood of such models is directly comparable with those of alternative models that lack this parameter, making assessments of model fit problematic (Ree and Sanmartín 2018). Therefore, we only report results from fitting the DEC, DIVALIKE, and BAYAREALIKE models without j to these biogeographic range data.

Sequencing and UCE Contig Assembly

Illumina sequencing of the equilibrated library with the 35 specimens in this study yielded 136,882,528 reads with an average read length of 302 bp, altogether yielding a total of 41.4 Gb of sequence data. This includes an average of 3,910,929 reads per specimen, although there was a notable range from the lowest to the highest number of reads for an individual (142,031 for Chlorophonia pyrrhophrys to 6,501,897 for Euphonia cyanocephala). Using the Trinity assembly algorithm to analyze these read data, we obtained an average of 4,289 contigs matching UCE loci from each specimen, although the number of UCE-matching contigs ranged from a minimum of 3,779 in C. pyrrhophrys to a maximum of 4,486 in Euphonia musica. A total of 4,956 unique UCE loci were recovered across all samples, and the mean and median sequence lengths of these UCE contigs were 928 and 989 bp, respectively. A complete matrix including only UCE loci with sequences from all 37 specimens included in the study contained only 150 loci; however, a 95% complete matrix contained 1,278 loci, a 90% complete matrix contained 2,870 loci, and the 50% complete matrix contained 4,829 loci. The SPAdes assembly performed comparably well in the recovery of UCE-matching contigs. Screening the results of this assembly yielded an average of 4,307 UCE contigs per specimen, with the minimum and maximum ranging from 3,776 to 4,486 in C. pyrrhophrys and E. musica, respectively. The mean sequence length for these UCE-matching contigs was 786, but the SPAdes assembly also included many additional contigs that, when screened against existing avian mitochondrial genome sequences, were found to comprise most, if not all, of the mitochondrial genome. These mitochondrial genomes were successfully assembled for every tissue specimen and were retained for separate phylogenetic analyses as a complement to the analyses performed with the Trinityassembled nuclear UCE sequence data. The number of reads and the number of UCE-matching contigs for each specimen in the separate Trinity and SPAdes assemblies is provided in Appendix Table 3.

Maximum Likelihood UCE Phylogenies

The concatenated incomplete sequence alignment of all 4,956 loci used to infer a phylogeny in RAxML was 4,746,327 bp in length and was split into 3 partitions defining the flanking regions, the core regions, and the UCE loci that were not split by SWSC-EN. The tree resulting from this analysis was, overall, strongly resolved as interpreted from bootstrap percentages and only 3 nodes received bootstrap values less than 100% (Figure 1A). These nodes involved the placement of E. saturata (91%), E. chlorotica (88%), and the separation of (E. anneae, E. xanthogaster) from (E. mesochrysa, [E. cayennensis, {E. rufiventris, E. pectoralis}]) (93%). The best-supported partitioning scheme for the concatenated alignment of the complete sequence matrix contained 119 partitions, and the total length of this 150-locus alignment was 159,894 bp in length. The RAxML phylogeny inferred with this complete, partitioned sequence alignment was also strongly resolved overall with only 7 nodes receiving <100% bootstrap support and this tree possessed only 2 topological differences from the incomplete phylogeny (Figure 1B). The incomplete topology inferred a sister relationship between 2 small clades of Euphonia species (anneae, xanthogaster) sister to (mesochrysa, [cayennensis, {rufiventris, pectoralis}]) to the exclusion of another small Euphonia clade containing (imitans, [gouldi, fulvicrissa]), whereas the complete topology inferred the (anneae, xanthogaster) clade being 1 node lower and sister to the other 2 clades. Neither arrangement of these clades had full bootstrap support in their respective trees, although both are >90%. The incomplete tree had E. jamaica as the most divergent branch in its respective clade, whereas E. saturata is the most divergent in the complete topology. Both of these placements have bootstrap support >90%. A notable feature of these phylogenies was that the 3 blue-hooded Euphonia species (cyanocephala, elegantissima, and musica) form a monophyletic group that is sister to the monophyletic Chlorophonia, a result that is consistent with previous results that found E. musica and C. cyanea as sister taxa (Zuccon et al. (2012)). This relationship was recovered in all subsequent phylogenetic analyses of both the nuclear UCE and mitochondrial genomic data and its taxonomic implications for these genera are discussed further below.

A total of 4,944 UCE gene trees containing >3 species were utilized to estimate the species tree for this lineage in ASTRAL. The topology of this species tree had maximum support values at all but 2 nodes as indicated by the LPP scores obtained via ASTRAL, with only the nodes defining the placement of E. chlorotica and E. saturata having scores less than 1 (Table 2; 0.98 and 0.76, respectively). There were no negative IC scores for any relationship in the species tree, but these scores did range from 0 to 0.643 and were generally low across the tree (median = 0.147). Both sets of concordance factors had fairly broad ranges as well, with gCF values ranging from 10.4 to 79.2 (median = 32.5) and sCF ranging from 20.8 to 96.3 (median = 58.95). gCF and sCF values were both positively correlated with IC values (Figure 2), with sCF values generally higher as expected given the relatively low levels of homoplasy in these slowly evolving loci. The support values from our manual bootstrapping were largely concordant with the LPP scores for the tree obtained directly from ASTRAL; the only 2 nodes in the entire topology with bootstrap scores less than 100% were those placing E. chlorotica (90%) and E. saturata (87%). Taken together, these metrics illustrate that the species tree inferred by ASTRAL does not have any hard, consistent conflict with any alternative topologies in the full distribution of gene trees, despite the variation in topologies that were inferred among individual gene trees.

This 4,944-locus species tree possessed 1 distinct relationship not found in either of the concatenated RAxML trees, with C. pyrrhophrys as the deepest split in the genus Chlorophonia with relatively high support (gCF = 34.2, sCF = 77.9, IC = 0.091, bootstrap = 100%) rather than being sister to C. cyanea. Interestingly, the species tree's placement of the 3 Euphonia clades mentioned above is concordant with that of the incomplete RAxML phylogeny, although the support values for this node are fairly low (gCF = 10.9, sCF = 35.0, IC = 0.0, bootstrap = 100%), and the placement of E. jamaica and saturata in this species tree also follows that of the incomplete RAxML phylogeny albeit with relatively low certainty in the placement of saturata as noted above.

The species tree generated from analyzing the 50% complete sequence matrix was topologically identical to the incomplete matrix topology and is not discussed further here, but the complete matrix species tree did possess a few topological differences when compared with the incomplete matrix tree (Appendix Figure 5, Appendix Table 5). The order of E. jamaica and E. saturata flip, in that jamaica is sister to all other species within its clade in the incomplete matrix topology whereas saturata is sister to all others in this same clade in the complete matrix topology. A similar flip of branching order occurs with E. concinna and E. finschi in this same clade, with the incomplete matrix topology having finschi diverge before concinna and the opposite relationship being supported in the complete matrix topology. Lastly, the complete matrix topology yielded a unique ordering of the 3 small Euphonia clades described above in that the (mesochrysa, [cayennensis, {rufiventris, pectoralis}]) clade was sister to the other 2, and the root node of this arrangement was well supported (gCF = 32.7, sCF = 75.3, IC = 0.555), more so than the root node of these 3 clades in the incomplete matrix topology.

FIGURE 1.

Concatenated phylogenies inferred from UCE sequences. The left tree (A) was generated from a partitioned analysis of 4,956 incompletely sampled loci and the right tree (B) from a partitioned analysis of 150 completely sampled loci. Black circles on nodes indicate 100% bootstrap support for that topology given its dataset and only nodes receiving <100% support have their bootstrap values reported on the tree. The white arrows point to the 2 topological differences between the trees: the arrangement of the 3 small Euphonia clades and of E. jamaica and saturata. The Euphonia clade marked with an asterisk designates the 3 blue-hooded species in this figure and in all following figures.

Mitochondrial Phylogeny

The phylogeny inferred from only the mitochondrial genome assemblies was largely concordant with our phylogenies inferred from the nuclear UCE loci, regardless of whether these loci were concatenated or analyzed with maximum likelihood or Bayesian methods (Figure 3). One sister species pair was inferred to have short branch lengths on this topology, with a total sequence divergence between Euphonia affinis and luteicapilla of 0.273%, and this shallow divergence was also reflected in the total UCE sequence alignment with 0.113% sequence divergence. Overall, the mitochondrial phylogeny was fairly well resolved, with all but 6 nodes receiving 100% bootstrap support. The other nodes lacking full support had bootstrap support ranging from 45 to 89%. Two points in this topology differed from either the concatenated phylogenies or the species tree. First, the relationship of the 3 small Euphonia clades matched that of the incomplete, concatenated phylogeny and the species tree with 100% bootstrap support, and only differs from the complete, concatenated phylogeny. Second, the mitochondrial data support a sister relationship of C. pyrrhophrys and cyanea with 100% bootstrap support, as is the case in the concatenated phylogenies with comparably high bootstrap support (97–100%). By contrast, the species tree reconstructed pyrrhophrys as the sister to all other Chlorophonia with reasonably high concordance factor support but low support from the IC index (gCF = 34.2, sCF = 77.9, IC = 0.091).

TABLE 2.

Node support values and node ages from the ASTRAL species tree topology of Euphoniinae. The full time-tree including the outgroups is provided in Appendix Figure 6 and the corresponding values for the nodes of those outgroups are provided in Appendix Table 4

Three additional points in the mitogenome tree support relationships that were not found in any of the UCE-based phylogenies. First, E. chrysopasta was found to be sister to the set of 3 small Euphonia clades. This placement has only 74% bootstrap support in the mitochondrial phylogeny, but its placement in the 2 concatenated trees as the sister to (minuta, [chalybea, {violacea, (hirundinacea, laniirostris)}]) had 100% bootstrap support. Its placement in the ASTRAL species tree matched that of the 2 other UCE-based phylogenies, but had moderate to low concordance factor and IC support (gCF = 18.5, sCF = 47.7, IC = 0.093). We searched the individual gene tree topologies containing this relationship to find the frequency at which the UCE sequences supported this relationship and found that only 6.9% (107/1,558) of gene trees did. Second, E. jamaica and saturata are sister species in the mitochondrial tree with reasonably high bootstrap support (84%), whereas jamaica is sister to all other species in its lineage with high support in both the ASTRAL species tree (gCF = 34.9, sCF = 70.0, IC = 0.322) and the incomplete RAxML phylogeny (bootstrap = 100%), and saturata is sister to all other species within its clade in the complete, concatenated RAxML phylogeny with full bootstrap support. Only 18.1% (701/3,867) of individual UCE gene trees supported a sister relationship between the 2 species. Third, E. chlorotica diverges prior to plumbea within its respective clade, whereas the opposite is true in every other analysis of UCE loci. This arrangement received 100% bootstrap support, although the placement of E. plumbea received only 59% support, and the placement of chlorotica within this clade consistently received only moderate support in the concatenated (88 and 73%) and species tree (LPP = 0.98, gCF = 10.4, sCF = 33.9, IC = 0.000) analyses. Only 2.3% (61/2,597) of gene trees supported this arrangement. Overall, we found evidence for only 3 instances of phylogenetic discordance between the mitochondrial and the nuclear UCE sequence data, with only moderate support.

FIGURE 2.

Internode certainty scores and concordance factors for the 4,944-locus ASTRAL species tree. Each point represents a node on the species tree, and the plot contains 2 full sets of points showing the gene and site concordance factors for each node in black and white, respectively. The dashed vertical line divides the plot in half: any points to the left of this line would indicate strong discordance of the gene trees and species tree for a particular node. While many points are on or directly adjacent to this line, no nodes in the species tree have negative internode certainties.

Time-calibrated Bayesian Tree and Biogeographic Modeling

We extracted the full sequence of cytb from the mitochondrial genomes and used the sequence from C. cyanea to assemble the fragments of cytb that were sequenced from C. occipitalis into a single consensus sequence. The cytb alignment for this locus was 1,145 bp in length and was concatenated as the 120th partition in the sequence alignment for the BEAST2 analysis. After modifying the BEAST2 analysis, we drastically improved convergence on the prior and posteriors (ESS ≥ 3,000), so we generated our MCC tree from this posterior tree distribution to obtain the first species-level, time-calibrated phylogeny for this subfamily (Figure 4; see Appendix Figure 6 for the complete tree with outgroups). As expected given our prior, the stem node for the entire phylogeny representing the divergence of Emberizidae and Fringillidae was estimated to be 16.8 million years ago (mya) (95% HPD = 14.5–19.3 mya), congruent with the age estimate and HPD originally obtained from Oliveros et al. (2019). Above this node, we recover the split between Fringillinae and the rest of the fringillids ∼14.6 mya (95% HPD = 12.4–16.5 mya), and the split between Euphoniinae and Carduelinae ∼13.8 mya (95% HPD = 11.8–15.8 mya). A long branch spanning nearly 7 myr was recovered connecting the stem and crown nodes for Euphoniinae, with the crown age for this subfamily being much younger at only 7.1 mya (95% HPD = 6.0–8.0 mya). The crown node for the genus Euphonia (minus the blue-hooded species) occurs shortly after at 6.5 mya (95% HPD = 5.4–7.3 mya), and the crown nodes for Chlorophonia and the blue-hooded Euphonia are even younger at 3.8 mya (95% HPD = 3.3–4.4 mya) and 2.0 mya (95% HPD = 1.6–2.3 mya), respectively. Following these crown ages, speciation events appear quite evenly spaced through time, with the youngest speciation events occurring less than 1 mya, as in the case of E. affinis and luteicapilla (0.68 mya, 95% HPD = 0.56–0.82 mya).

FIGURE 3.

RAxML phylogeny inferred from partitioned mitochondrial genomes. Black circles indicate 100% bootstrap support for nodes along the tree, and only bootstrap supports <100% are shown adjacent to their respective nodes. White arrows point to 3 points on this phylogeny that differ from all of the UCE-based phylogenies: the sister relationship of E. jamaica and saturata, the placement of E. chrysopasta, and the placement of E. chlorotica. As in Figure 1, the blue-hooded Euphonia clade is marked with an asterisk.

Of the 3 biogeographic models fit to this tree, the DIVALIKE model was best-supported, although the DEC model fit these data comparably well (ΔAIC = 1.437). The actual biogeographic histories between these 2 models differed only marginally in that 3 nodes were reconstructed as being widespread North and South American in the DEC model and only South American in the DIVALIKE model, so here we focus on the DIVALIKE model results. In this model, the most likely ancestral range of crown Euphoniinae was South America, which was also the most likely ancestral range for the node uniting the blue-hooded euphonias and Chlorophonia and for the ancestral node of the remaining Euphonia. Beyond this initial period of diversification, all but 1 diversification event appeared to have occurred in South America until ∼4 Ma, after which a number of nodes have inferred ranges occupying both North and South America. A total of 6 range expansion events into North America from South America were inferred under this model: (1) the crown node of blue-hooded Euphonia; (2) the ancestor of C. flavirostris, occipitalis, and callophrys; (3) the ancestor of E. trinitatis, luteicapilla, and affinis; (4) the ancestor of E. laniirostris and hirundinacea; (5) the ancestral node for the 3 small Euphonia clades; and (6) within the small Euphonia clades, the ancestor of E. anneae and xanthogaster. Additionally, 2 nodes were inferred to have ranges occupying both South America and some region of the Caribbean, and both of these nodes give rise to a single species currently found in the Caribbean on 1 side of the split and at least 1 South American species or node on the other side. However, these nodes leading to species occurring in the Caribbean were not temporally clustered like the South to North America range shifts and occur at 4.8 and 1.4 mya for E. jamaica and E. musica, respectively. Taken together, the DIVALIKE model of biogeographic range evolution showed a consistent South American origin for the euphoniine lineage that repeatedly expanded northward and dispersed twice to islands of the Caribbean.

FIGURE 4.

Biogeographic history and Bayesian time-scaled phylogeny of Euphoniinae. Node numbers correspond to the values of node support in Table 2, and the pie charts at each node show the proportional likelihood of possible ancestral ranges inferred by the DIVALIKE model. A legend describing these regions is given at the bottom with the following abbreviations: NA = North America, SA = South America, and C = Caribbean. Gray shaded regions along the phylogeny indicate 2 geologically significant time periods: the duration and end of Andean orogeny and the recently estimated age at which the Isthmus of Panama formed following O'Dea et al. (2016). Black lines delimit the euphoniine genera and gray bars adjacent to Euphonia delimit the clades that we describe in the Discussion. As previously, the blue-hooded euphonia clade is distinguished by the asterisk.

The results of the lesser supported DEC model suggest that 3 ancestors: (1) the ancestral node for Euphoniinae; (2) the node uniting Chlorophonia and blue-hooded euphonias; and (3) the node uniting Euphonia anneae + xanthogaster with E. mesochrysa, cayennensis, rufiventris, and pectoralis, occupied both North and South America, instead of only South America as in the DIVALIKE model. This biogeographic scenario implies only 4 range expansions from South America into North America instead of 6, but the same South America to Caribbean range expansions were inferred and the remainder of the biogeographic history inferred from this model is nearly identical to that of the DIVALIKE model.

Congruence of Phylogenomic Analyses and Previous Classifications

In this study, we applied molecular systematic methods to resolve the taxonomic relationships among the genera and species in the passerine subfamily Euphoniinae and used this phylogenetic framework to model the biogeographic history and timing of diversification within this unusual group of birds. Combining a well-developed target enrichment and sequencing pipeline with species-level sampling allowed us to infer the first species-level phylogeny for this subfamily with orders of magnitude more sequence data than previous studies. The phylogenetic hypotheses inferred from nearly 5,000 nuclear UCE loci and complete mitochondrial genomes using both concatenated and coalescent-based methods are largely congruent and effectively resolve the genus- and species-level relationships of Euphonia and Chlorophonia.

We found that all recognized species in the subfamily form a young, monophyletic lineage sister to outgroup representatives (Spinus and Coccothraustes) from the sister subfamily Carduelinae, confirming that Euphoniinae is a sound taxon. Immediately after the basal divergence of this lineage, we recover a divergence event that separates blue-hooded Euphonia species and Chlorophonia from the remaining Euphonia species, a relationship originally suggested by Zuccon et al. (2012) with a smaller sample of taxa. By including all species in this phylogeny, we confirmed that all 3 blue-hooded Euphonia species form a monophyletic lineage that is sister to the 5 Chlorophonia species with uniformly high support in every phylogenetic analysis, demonstrating that the genus Euphonia, as currently delimited, is paraphyletic with regard to Chlorophonia. We address the need for revising taxonomic limits of these genera and offer a resolution to this issue below. Within the blue-hooded Euphonia clade, every analysis recovered E. elegantissima as sister to musica and cyanocephala, and relationships within Chlorophonia varied only in the placement of pyrrhophrys. Every concatenation-based analysis of both the UCE and mitochondrial sequences infer a sister relationship between C. pyrrhophrys and cyanea to the exclusion of the other 3 species with high bootstrapping support, whereas the species tree analysis places C. pyrrhophrys as the sister to all other species within the genus, albeit with a short branch and with moderate to low IC support (0.091). These results suggest that the early diversification events within this genus may have happened rapidly.

Similar to Isler and Isler's (1999) treatment of Euphonia, our results support 5 clades within the genus that, with the exception of 1 species in the mitochondrial tree (E. chrysopasta), are found and well-supported in each phylogeny we inferred in this study. The first and largest clade we describe here contains the species E. jamaica, saturata, plumbea, chlorotica, finschi, concinna, trinitatis, luteicapilla, and affinis, and is an amalgamation of groups 1 and 2 sensu Isler and Isler. The only variable placements of any species in this clade are the oldest divergences defining the placement of E. jamaica, saturata, and chlorotica. The second clade we define includes the species E. chrysopasta, minuta, chalybea, violacea, laniirostris, and hirundinacea, and includes species from Isler and Isler's species groups 3, 5, and 6. Relationships within this clade were uniformly resolved across analyses with the exception of E. chrysopasta's placement in the mitochondrial genome tree as the sister taxon to the next 3 Euphonia clades. But the placement of this taxon at this point or within our second clade never received maximum support values in any phylogenetic analysis.

The last 3 clades we identified across our phylogenies form a monophyletic group in each tree, although the relationships between these 3 clades are some of the most variable points in the phylogeny of the entire subfamily. The concatenated phylogenies inferred from the 4,956 concatenated UCE loci and from the mitochondrial genomes plus the species tree inferred from the incomplete UCE matrix placed our third clade, containing the species E. imitans, gouldi, and fulvicrissa, as sister to the remaining 2 clades at this point in the tree. Excepting the exclusion of E. chrysopasta and mesochrysa, this clade is identical to Isler and Isler's group 5. Only the RAxML phylogeny inferred from the 150 concatenated and partitioned UCE loci proposed a different relationship among these 3 small clades in which our fourth clade, which is identical to Isler and Isler's group 7 and contains only E. anneae and xanthogaster, was sister to the other 2. Within this clade, we found surprisingly deep divergences between the 2 subspecies of E. xanthogaster in this study, and the significance of this result, along with its implications of potentially undescribed species-level diversity, are discussed below. Lastly, our fifth clade, which includes 4 species, E. mesochrysa, cayennensis, rufiventris, and pectoralis, has uniformly resolved relationships in every phylogeny, and is congruent with Isler and Isler's species group 8 except with the addition of mesochrysa.

Altogether, our phylogenetic hypotheses inferring the relationships within the subfamily Euphoniinae are largely concordant with one another and with previous taxonomic treatments of this group. With few exceptions, we recover the species groups defined by Isler and Isler in their classification of the genera with the only differing relationships inferred in our study involving species that were not concretely placed in their classification. At a higher level, the obvious exception to this concordance is Isler and Isler's species group 4, containing the blue-hooded euphonias, which we recover as sister to Chlorophonia and outside of Euphonia. With complete sampling and the analytical power afforded by our data, it is now possible and necessary to revise the generic taxonomy of these birds with 3 plausible alternative solutions outlined below.

Paraphyly of Euphonia and Historical Precedent for Three Genera of Euphoniine Finches

Here we present 3 potential taxonomic revisions that would resolve the paraphyly of Euphonia: (1) lump all species in the subfamily Euphoniinae into Euphonia, (2) add the blue-hooded euphonias to Chlorophonia, or (3) resurrect the previously described genus Cyanophonia for the blue-hooded euphonias. The genus Euphonia was established by Anselme-Gaëton Desmarest in 1806 in his Histoire naturelle des Tangaras, des Manakins et des Todiers and, thus, has priority over the genus Chlorophonia, which was described by Charles Lucien Bonaparte in 1851 in the second series and third volume of the Revue et Magasin de Zoologie Pure et Appliquée. While this priority lends support to the first proposed solution, the dissolution of Chlorophonia and lumping all species under Euphonia would remove a comparably long-standing genus and overturn 168 yr of recognition of the former genus. The second proposed solution would require the revision of 3 species binomials, rather than 5 in the first solution, but given that only females of the blue-hooded euphonia species and males in some populations of E. musica in the Lesser Antilles possess any green plumage as adults, the name Chlorophonia does not describe the phenotypes of these species particularly well. The third proposal stems from the fact that, in the same work in which he created the name Chlorophonia, Bonaparte also designated a separate genus for E. musica and cyanocephala (then described as aureata): “Nous donnons plus particulièrement le nom de CYANOPHONIA aux Euphones à tête bleue [More specifically we give the name CYANOPHONIA to the blue-headed Euphonias]” (Bonaparte, 1851). Because the name Cyanophonia has equal priority to Chlorophonia and was specifically given to describe the blue-hooded species as distinct from other euphonias and the chlorophonias, we recommend resurrecting this genus and assigning the blue-hooded euphonias to this genus. Additionally, this solution requires revising only 3 taxon names and will result in a descriptive generic name for them. Despite the description of this genus in 1851, to our knowledge, a type species has never been formally designated for the genus Cyanophonia. Of the 3 recognized species in this genus, musica was described by Johann Friedrich Gmelin in part 2 of the 13th edition of Systema Naturae (1789) and is older than both cyanocephala (1819) and elegantissima (1838). Therefore, we designate the Antillean Euphonia, Cyanophonia musica (Linnaeus and Gmelin 1789), as the type species for the genus.

Biogeographic History of Euphoniinae

Given the present-day distribution of the family Fringillidae and closely related passerine families, in addition to the biogeographic analyses performed in this study, we conclude that the euphoniine lineage dispersed to the Western Hemisphere from Eurasia during the mid- to late-Miocene between 13.8–7.1 mya. As described above, modern members of this lineage are distinct ecologically when compared with most other finches, so it is likely that this lineage dispersed from the Eastern Hemisphere and evolved these unique adaptations from an ancestral nomadic, granivorous phenotype during this 7 myr time frame. But by what route did this ancestral population arrive in the Americas? The vast majority of passerine dispersal events from the Eastern into the Western Hemisphere appear to have taken place through Beringia (Cracraft 1973, Barker et al. 2004), which facilitated overland dispersal into North America throughout most of the Cenozoic until the latest Miocene and Pliocene between 7.4 and 4.8 mya and then sporadically throughout the Pleistocene (Hopkins 1967, Marincovich and Gladenkov 1999). Although rarer, transoceanic dispersal events from Eurasia or Africa into the Americas, especially South America, have repeatedly occurred in a number of passerine lineages (Anthus, Voelker 1999; Turdus, Voelker et al. 2009; Donacobius, Fregin et al. 2012), in pelecaniform birds (Bubulcus ibis, Crosby 1972; Plegadis falcinellus, Oswald et al. 2019), and even in 2 orders of mammals (Primates and Rodentia, Simpson 1980). Many extant finch species are long-distance migrants and many others (euphonias included) are nomadic throughout their annual cycles, which has led to their colonization of isolated oceanic islands like the Azores, Canaries, and Madeira by Fringilla (Marshall and Baker 1999), and the Hawaiian archipelago by the ancestor of the honeycreepers (Lerner et al. 2011). Considering these confamilial examples of long-distance oceanic dispersal with the BioGeoBEARS results from this study, transoceanic dispersal from the Eastern Hemisphere to South America is indeed a possible, and perhaps even likely, means of arrival for Euphoniinae into the Americas. An alternative scenario, supported by the DEC model inferred from these data, proposes that the ancestral population giving rise to Euphoniinae was widespread in both North and South America. It is possible that the ancestral population could have dispersed to either North or South America over water or over land, expanded its range to occupy both continents, and underwent range contractions or extinctions within both Euphonia and Chlorophonia, only leaving descendants in South America. Given the variability in age estimates at which North and South America became one contiguous landmass (discussed below), we cannot rule out this scenario despite it not being selected as the one best-fitting the phylogeny and contemporary range data.

Following dispersal into the Western Hemisphere, both Euphonia and Chlorophonia apparently diversified in South America, while we infer that Cyanophonia remained widespread throughout North and South America until 1.8 mya, around which time it diversified. Species in this genus typically inhabit forests and other habitats where mistletoes (Santalales) are common, and their speciation during the Pleistocene may have been vicariant—driven by climate cycles that repeatedly shifted and isolated patches of such forests—with elegantissima diverging first in North America, followed by musica and cyanocephala in the Caribbean and South America, respectively. Interestingly, within Euphonia and in Chlorophonia, we found that all lineages presently occurring in North America evolved from South American ancestors. Assuming these lineages had evolved the strongly frugivorous and nonmigratory habits of modern species by their crown age, all these dispersal events likely required continuous expanses of tropical forest connecting both continents in order to take place. Such expanses of continuous habitat could only exist on a completely formed Isthmus of Panama. A substantial research effort spanning the fields of geology, paleobiology, and phylogenetics has focused on resolving the time frame over which the Central American Seaway closed and the Isthmus of Panama fully formed, with estimates ranging from older than 20 mya (Bacon et al. 2015), to the mid-Miocene around 13–10 mya (Montes et al. 2015), and as young as 2.8 mya (Coates and Obando 1996, O'Dea et al. 2016). The divergence times for each of the species and clades found in North America ranges from 4.3 to 1.8 mya, with the upper bound of the 95% HPD for any of these nodes extending to 4.9 mya. Assuming that these lineages dispersed northward as soon as adequate habitat was available and the route was traversable, which was likely the case in both North and South America by this time frame (Graham 1987a,b, Behrensmeyer et al. 1992, Jaramillo et al. 2014), our data support a younger age for the formation of the isthmus and corroborate previous findings that avian dispersal between North and South America reached its maximum frequency during the Pliocene (Bacon et al. 2015, Barker et al. 2015).

Additionally, our biogeographic modeling results clearly show 2 spatially and temporally independent dispersal events of a single species into the Caribbean from a South American ancestor, Cyanophonia musica and Euphonia jamaica. A large body of literature has investigated the biogeographic history of the Caribbean and the avenues by which this region has accumulated its biodiversity. Recent work describing the historical biogeography of both terrestrial mammals and birds has shown that dispersal to the islands of the Caribbean from South America is a prevalent and likely scenario for these broad taxonomic groups (Dávalos 2004, Ricklefs and Bermingham 2008). The geographic landscape of the Caribbean strongly resembled its modern state by ∼5 mya, at which time substantial tectonic uplift had led to the emergence of the majority of islands and, thus, could facilitate shorter distance dispersal events in a stepping-stone manner between landmasses (Perfit and Williams 1989). Given the absence of land bridges connecting South America to any of the Caribbean islands since the late Oligocene (Iturralde-Vinent and MacPhee 1999) and that both diversification events producing these species occurred after 5 mya, island-hopping and overwater dispersal appear to be the only possible route by which these 2 species could have come to occupy their current ranges. This explanation is especially convincing for Cyanophonia musica, which is found on most of the Lesser Antilles, Puerto Rico, and Hispaniola, and whose sister species is the widespread South American species C. cyanocephala, implying that this species may have expanded northward via successive northward colonization of islands from the mainland.

Potentially Unrecognized Species-level Diversity

Another notable outcome of this study was the finding that the magnitude of sequence and time divergence between the 2 subspecies of E. xanthogaster examined here is as great as, or even greater than, that observed between some other sister species in the tree. E. xanthogaster exhibits striking phenotypic variation throughout its extensive range in South America (Olson 1981), where populations in the northern Maritime Andes (subspecies exsul and badissima) and in La Paz and Cochabamba Bolivia (subspecies ruficeps) have rufous crowns and orange ventral plumage, and, throughout the remainder of this species' range, the crown and ventral plumage is yellow (including subspecies quitensis and the allopatric subspecies xanthogaster in the Atlantic Forest; Isler and Isler 1999). Although this species is sexually dimorphic, these phenotypic differences in crown and ventral plumage are found to some extent in both males and females. Considering the genetic and phenotypic distinctiveness of the subspecies included in this study, this suggests that a thorough population-level study of genetic and phenotypic divergence of this species may reveal undescribed species-level diversity in E. xanthogaster across the continent of South America.

Similar circumstances of allopatric and phenotypically distinctive subspecies are documented in several other species, including C. cyanea, E. affinis (Dickerman 1981), E. cyanocephala, E. laniirostris, and E. musica (Isler and Isler 1999, del Hoyo et al. 2019). A lineages-through-time plot generated from our ultrametric MCC tree shows a plateauing of species diversification toward the present and a significant, negative γ value (–3.758, P < 0.001) that can be interpreted in 2 ways (Appendix Figure 7). First, it is possible that Euphoniinae underwent rapid, early diversification that has slowed down toward the present in a density-dependent manner (Pybus and Harvey 2000, Rabosky 2009). Or second, and perhaps more likely given our findings with E. xanthogaster and our observations of other wide-ranging phenotypically disparate taxa, this plateauing of species diversification toward the present is an artifact of failing to account for undescribed species-level diversity within this lineage. Investigating potentially undescribed species-level diversity within this lineage and testing which scenario best explains the accumulation of lineages through time for this group will require much denser sampling and is beyond the scope of the present study.