Deep divergence of Red-crowned Ant Tanager (Habia rubica: Cardinalidae), a multilocus phylogenetic analysis with emphasis in Mesoamerica

Taxon sampling and laboratory procedures

We sequenced 125 individuals of H. rubica species that covered the species distribution, including samples from Mexico (using field collection on permit from Instituto Nacional de Ecología, SEMARNAT: FAUT-0169), Guatemala, El Salvador, Honduras, Nicaragua, Costa Rica, Panama, Peru, Bolivia, Paraguay, Brazil, Argentina and Venezuela (Fig. 1A). We also included 16 samples from the genus Chlorothraupis (C. olivacea, C. carmioli and C. stolzmanni), the sister group of H. rubica (Klicka, Burns & Spellman, 2007), and three congeneric samples that were used as outgroups (H. fuscicauda, H. atrimaxillaris, and H. gutturalis; see Table S1).

We isolated genomic DNA from tissues using Qiagen DNAeasy kit (Qiagen, Valencia, CA, USA) following the manufacturer’s protocol. We amplified five molecular markers via the polymerase chain reaction (PCR) in 12.5 µl volumes. The mitochondrial marker ND2 (1,041 bp), was generated for all individuals (n = 144) using the primers L5219 and H6313 (Sorenson et al., 1999). We then amplified four nuclear markers for a subset of individuals with the highest DNA quality (n = 37). These nuclear loci included two Z-linked introns, the aconitase I gene (ACOI-I9, 852 bp) amplified using the primers ACOI-I9F and ACOI-I9R (Kimball et al., 2009), and the receptor of the tyrosine kinase MUSK gene (MUSK, 576 bp) amplified using the primers MUSK-I3F and MUSK-I3R (Kimball et al., 2009). We also amplified two autosomal loci, a beta-fibrinogen gene (FGB-I5, 577), using the primers Fib5 and Fib6 (Kimball et al., 2009) plus the Ornithine Decarboxylase ODC gene (ODC, 711 bp) using the primers ODCF and ODCR (Primmer et al., 2002). We used the following annealing temperatures: 54 °C (ND2), 56 °C (MUSK), 60 °C (ACOI and FGB-I5) and 65 °C (ODC). PCR success was verified on 1% agarose gels.

DNA sequencing was performed in the High Throughput Genomic Unit Service of the University of Washington. We edited and aligned chromatograms with Sequencher v4.8 (GeneCodes Corporation, Ann Arbor, MI, USA). We corroborated the origin of all of our sequences by combining at least two of the following methods: amplifying overlapping gene segments, amplifying or sequencing one region with different primer sets, sequencing both DNA strands for all amplified fragments, or using multiple individuals of a single species. No evidence of NUMT contamination was found (Zhang & Hewitt, 1996; Sorenson & Quinn, 1998; Bensasson et al., 2001), and all nuclear sequences could be easily aligned with the corresponding genes from other avian species. The overall length of the concatenated mitochondrial and nuclear sequences was ∼3,800 bp. All sequences were deposited in FigShare (https://figshare.com/s/fe9f9f6fed1686782f62).

Phylogenetic analyses

Phylogenies were constructed from two different matrices; the first matrix included 144 sequences of the mitochondrial marker ND2 (hereafter referred to as the mtDNA dataset), and the second matrix included 37 concatenated sequences of the mitochondrial marker ND2 plus the four nuclear markers describe above (hereafter referred to as multilocus dataset). Phylogenetic analyses were performed for each matrix using a Bayesian Inference method (BI), with independent partitions and evolutionary rates defined for each gene region. The model parameters of nucleotide evolution for each DNA region were obtained with jModeltest v3.8 (Posada, 2008) using the Akaike Information Criterion (Akaike, 1973). Analyses were conducted in MrBayes v3.2 (Ronquist et al., 2012), and consisted of four MCMC chains of 50 million generations each, starting from random trees and using a uniform prior distribution of parameters. Trees were sampled every 250 generations, using Tracer v1.6.0 (Rambaut & Drummond, 2013), after discarding the initial 20% as burn-in. The remaining trees were used to construct a majority rule consensus tree with a posterior probability distribution. The final tree was visualized in FigTree v1.2.3 (http//tree.bio.ed.ac.uk/software/figtre/). Finally, we used MEGA v7 (Kumar, Stecher & Tamura, 2016) to estimate the pairwise genetic distances between lineages from mitochondrial and nuclear markers, applying a bootstrap of 1,000 iterations; this allowed us to indirectly determine the relative amounts of information provided by the mitochondrial and nuclear markers.

Diversity and genetic structure

A haplotype network was obtained for the mtDNA dataset with the Median-Joining algorithm available in Network v4.6 (Bandelt, Foster & Röhl, 1999). We used DnaSP v5.0 (Librado & Rozas, 2009) for estimating the number of haplotypes (H) and segregating sites (S), the haplotype (Hd) and nucleotide (π) diversities, and the Tajima’s D (Tajima, 1989) and Fu’s F (Fu, 1997) values for the genetic clusters recovered in the BI mitochondrial phylogeny mentioned above. Significance for the last two summary statistics was inferred from 1,000 replicates of the coalescent algorithm.

Species delimitation

The phylogeny obtained from the multilocus dataset was used to assess species delimitations under the Bayesian framework implemented in BP&P (Yang & Rannala, 2010). This method allows the estimation of posterior distributions of competing models with differing number of species. It incorporates a reversible-jump Markov chain Monte Carlo (rjMCMC) algorithm and a user-specified tree, where it assigns speciation probability values (Rannala & Yang, 2003; Yang & Rannala, 2010). BP&P assumes constant population sizes, no gene flow among lineages, and no population structure within lineages after speciation occurred. It furthers includes divergence time (τ) and mutation-scaled effective population size (θ) parameters. We first performed several preliminary runs with the two available algorithms (0: fixed tree, and 1: tree as guide tree) and different seed numbers to confirm consistency and fine-tuning the priors (ε and τ). After selecting algorithm 0, we conducted a first set of analyses allowing the program to make automatic adjustments to the priors (finetune = 1). Because different values of τ and θ can affect the model posterior probabilities for a same guide tree (Leaché & Fujita, 2010; McKay et al., 2013), and to explore a wide range of speciation histories, we performed a series of simulations using three different parameter combinations: 1) large population size/deep divergence [G (1, 10) for both θ and τ], 2) small population size/ shallow divergence [G (2, 2000) for both θ and τ], and 3) large population size/shallow divergence [G (1, 10) for θ and G (2, 2000) for τ] (Leaché & Fujita, 2010; Smith et al., 2013). We then performed a final set of analyses that incorporated the finetune parameters determined in the previous analysis (finetune = 0) with the same three different combinations for θ and τ parameters described above. Each set was composed of 100,000 generations, and samples were taken every 5 steps after discarding the initial 50,000 iterations as burn-in. All analyses were performed on the BPPX graphics user interface (http://abacus.gene.ucl.ac.uk/software) using the command line.

Tests of divergence times

Divergence times were inferred from the multilocus dataset. Analyses were performed in Beast v1.7 (Drummond et al., 2013) using the previously defined gene partitions with the corresponding evolutionary rates and models as estimated above (Ellegren, 2007; Smith & Klicka, 2010) and the seven phylogroups identified with the multilocus dataset. Clade divergence was simulated as a Yule’s first tree (Yule, 1924) with a lognormal relaxed molecular clock (Drummond et al., 2006). Analyses included 100 million iterations, with samples taken every 1,000 generations after eliminating 25% of the trees as burn-in. Mixing and likelihood stability were confirmed with TRACER v1.6 (Rambaut & Drummond, 2013) to make sure that the appropriate effective sample size (ESS) was over 200 units for each parameter (Drummond & Rambaut, 2007). The best supported tree was chosen with TreeAnnotator 1.8.0 (Drummond & Rambaut, 2007) through Maximum Clade Credibility.

Given the absence of fossils for H. rubica, or any close relatives, two calibrations based on geological data were used. First, the closing of the Isthmus of Panama, estimated between 3.1 and 4.0 Ma (Daza, Castoe & Parkinson, 2010; Smith & Klicka, 2010), was used for the separation of South and Central American populations; and second, the uprising of the Talamanca Cordillera near the Costa Rica–Panama border, estimated to have occurred between 2.5 to 3.9 Ma (Marshall & Liebherr, 2000; Daza, Castoe & Parkinson, 2010), was employed as the time of divergence between populations in Panama and those occurring elsewhere throughout Mesoamerica. In addition, as estimating divergence times based on external event that could be unrelated to a species biology can result in wrong inferences and circular reasoning, we estimated divergence dates a second time but using only the mtDNA dataset (which has a higher sample size) and no calibration points. This could be seen as an independent test of confidence for the estimates obtained for the multilocus analysis with geological constraints. We assumed a mutation rate of 1.25 ×10⁻² substitutions/site/lineage/Myr (Smith & Klicka, 2010), other parameters and priors were kept identical.

Phylogenetic analyses

We obtained mitochondrial (1,041 bp) and concatenated (3,757 bp) datasets. The best-fit models for each molecular markers were follows: ND2, TrN+I+G; ACOI-I9, TIM2+G; MUSK, TrN+G; FGB-I5, HKY+G and ODC, TrN.

The phylogenies obtained with mitochondrial and multilocus dataset revealed well-supported topologies, independently of the method employed. Both topologies are composed of three main clades with seven and eight subclades (considered here as “phylogroups”) for the multilocus and mitochondrial datasets, respectively. The first clade in the multilocus phylogeny (I) is comprised of two phylogroups with an apparent separation between populations occupying the western (A; Fig. 1B) and southern (B; Fig. 1B) Pacific coast of Mexico. The second clade (II) was composed of three well-supported phylogroups distributed from eastern Mexico to Panama. The first includes all individuals from the Mexican Gulf Coast, west of the Isthmus of Tehuantepec (C; Fig. 1B); the second one comprised samples collected east of this Isthmus and southward into Nicaragua (D; Fig. 1B); and the third one was composed by samples from Panama (E; Fig. 1B). The last clade (III) was composed of (at least) two phylogroups, one with individuals from the Amazon basin (F; Fig. 1B) and another one with birds from eastern South America (G; Fig. 1B). A single sample from Venezuela appears distinctive and is sister to this eastern clade. However, additional samples from Venezuela will be required for a proper systematic assessment.

The mitochondrial phylogeny, although based on many more samples, was almost identical to the one obtained with the multilocus dataset (Fig S1); it did however, show that phylogroup D was composed of two subgroups: one comprised of individuals from the Yucatan Peninsula, Honduras, Nicaragua and Costa Rica, and the other one of individuals from Guatemala and El Salvador.

The pairwise genetic distances showed significant genetic differentiation in the mitochondrial but not in the nuclear markers (Table 1). That is, the genetic distances were smaller for nuclear than for mitochondrial markers (Table 1). This implies that the mitochondrial marker provides larger amounts of relative information than the nuclear gene-regions (Table 2).

Diversity and genetic structure

The haplotype mtDNA network (Fig. 1C) revealed the same eight groups obtained with BI phylogeny (Fig S1). A high number of haplotypes had to be reconstructed to link some of these groups. For instance, 30 mutational steps separated the two groups from the Mexican Pacific coast from the other clusters, 45 more steps were needed to link the clades from South and Central America, while the two groups from South America were 41 steps apart; an additional 29 mutations were necessary to join the lone sample from Venezuela to the southeastern South America clade.

Overall, diversity estimates differed greatly between mtDNA phylogroups, with the South American and Panamanian ones being the more diverse, and the groups from the Mexican Pacific coast the least genetically rich. Tajimas D and Fu’s F values were negative in all phylogroups, excepting phylogroup D (eastern Mexico to northern Central America). However, these values were significant only for the two South American Clades (F & G; Fig. S1); thus suggesting recent population expansions (Table S2).

Species delimitation

Considering a hypothesis of seven species for this taxon, the Bayesian species delimitation supported the guide tree displaying speciation probabilities of 1.0 in most nodes, independently of the finetune parameter used (finetune = 1 or 0, Figs. 2A and 2B respectively). Seven species were supported with high speciation probabilities (0.97 to1.0) on the tree guide. These species corresponded to the lineages A, B, C, D, E, F and G. The different prior distributions for θ and τ had no effect on most nodes, the only exception being that of lineages A and B which were influenced by at least one of the prior combinations (from 0.37 to 0.97). However, while these and other analyses with similar topologies and highly supported clades suggest a multi-species complex within H. rubica, other types of data are still necessary to support this hypothesis (Sukumaran & Knowles, 2017; see below).

Tests of divergence times

According the multilocus dataset, the divergence between H. rubica and Chlorothraupis dated back to the late Neogene, between 6.65 and 3.91 Myr ago (Fig. 3A), and the ensuing diversification of H. rubica was rather rapid. The first split was estimated some 3.75 to 3.35 Myr, between the South American (clade III) and Central/Mesoamerican lineages (clades I and II); these two groups split again around 3 Myr ago. In Central America, the first divergence event (3.4–2.7 Myr) occurred between the lineages from the Mexican Pacific coast (Clade I) and the other phylogroups (Clade II). Within clade II, there was a gradual and northward divergence of lineages from Panama to eastern Mexico; these three diversification events during the last ∼3Myr or ∼2Myr. An independent analysis performed on the mtDNA dataset and with no calibration constraints, produced an identical topology with older mean divergence times (i.e., divergence within H. rubica occurred during the last 5.6 Myr). However, as confidence intervals were much wider (Fig. 3B) and overlapped with those obtained with the multilocus dataset, they provide support for the dates estimated based on geological times.

Phylogeny, diversity and genetic structure

The mitochondrial DNA topology indicates that populations of Mexico, Central and South America are profoundly differentiated. In this topology eight phylogroups were identified, six in Mesomerica and two in South America (Fig S1). The relationship between mitochondrial haplotypes could be determined by the action of various historical processes that promoted a deep genetic structure in Mesoamerica, resulting in a clear example of cryptic speciation. This complicated phylogenetic structure in Mesoamerica is consistent with the complex geological and biogeographic history of the region (Coates & Obando, 1996), and the emergence of geographic features such as the Rio Balsas depression, the Isthmus of Tehuantepec in southern Mexico, the Motagua-Polochic-Jocotán fault system in the North of Central America, and the Talamanca Cordillera in Costa Rica, which presumably played roles in shaping the observed geographic and genetic patterns (Daza, Castoe & Parkinson, 2010; Gutiérrez-García & Vázquez-Domínguez, 2012; Gutiérrez-García & Vázquez-Domínguez, 2013; Suárez-Atilano, Burbrink & Vázquez-Domínguez, 2014). Our results together with the study of Lavinia et al. (2015) provide a robust hypothesis of the evolutionary history of H. rubica. For instance, our study confirms the main results of Lavinia et al. (2015) for South America, and improves the resolution for Mesoamerica given our larger sample size and wider area covered.

In South America, our results, in combination with those of Lavinia et al. (2015), showed two distinct phylogroups that reflected population divergence according to the habitats preferred by H. rubica (Pennington et al., 2004; Pennington, Lavin & Oliveira-Filho, 2009; Lavin, 2006; Smith et al., 2013; Smith et al., 2014; Pennington & Lavin, 2016): one western phylogroup distributed in the mesic forest of the Amazon basin, and another one in the seasonal forests of the eastern and northern parts of South America. This last lineage included the single sample from Venezuela, instead of the western clade that is geographically more proximate (Prates et al., 2016; Prates et al., 2017). This could be indicating that the deep evolutionary history of the seasonal forests in South America (Pennington et al., 2004; Pennington, Lavin & Oliveira-Filho, 2009) are driving the divergence of other taxa. Additional sampling in this region is thus necessary to test this hypothesis.

The multilocus DNA topology identified seven of the eight phylogroups found with the mitochondrial DNA (Fig. 1B), in all cases the phylogroups are the same with the exception of an additional separation of a Yucatan Peninsula and northern of Central America clade identified using the more rapidly evolving mtDNA marker. All identified phylogroups (mitochondrial and multilocus dataset) are consistent with clear geographical discontinuities and are supported with high posterior probabilities.

Although, the haplotype network recovered the same groups identified in the mitochondrial DNA phylogeny (Fig. 1C), it was necessary to invoke five hypothetical haplotypes between Mexican and Central American groups. These results support the idea about the complex geological history in Mesoamerica resulting in a high degree of genetic structure in the northernmost distribution of this species. In South America, 12 hypothetical haplotypes were required, but unlike the previous case, this could be explained by the lack of sampling in this specific area.

The genetic distances between phylogroups were relatively high and supported the groups identified in our phylogenetic and BPP analyses. The largest difference found was between the group from the Western Pacific Coast of Mexico (group A) and the group from the Southeast of South America (group G) (Table 1). This result is not surprising because the geographical distance between them is also the largest and is congruent with the divergence date (Figs. 3A and 3B). On the other hand, the more similar phylogroups were the Western (A) and the Southern Pacific Coast of Mexco (B) groups, which represent the most recent divergence event (Figs. 3A and 3B). As expected the mitochondrial DNA provides more polymorphic sites because its rate of evolution is faster (Table 2) (Brown, George Jr & Wilson, 1979).

Species delimitation and divergence times

The deep divergence times among lineages estimated herein (Figs. 3A and 3B), most of which dated back to Pleistocene and Pliocene, support the criticisms raised regarding the allopatric speciation hypothesis that invoked the Pleistocene glacial cycles as the major drivers of neotropical bird diversity and closing of the Isthmus of Panama (e.g., Moritz et al., 2000; Cowling, Maslin & Sykes, 2001; Glor, Vitt & Larson, 2001; Mayle et al., 2004). According to the ND2 divergence analyses, the first split was during the Miocene (5.69 Myr), and this is the main discordance with respect to the calibrated multilocus divergence analyses (3.55 Myr). These dates are also consistent with the compelling evidence for speciation drawn from the BP&P analyses, which suggest that enough differences have accumulated between most of the H. rubica lineages to consider them a complex of cryptic species (Fig. 2).

At present, intraspecific delimitation within H. rubica is based on geographic variation in plumage color, with 17 subspecies currently recognized (https://avibase.bsc-eoc.org/avibase.jsp?lang=EN). Many of these subspecies are not supported by the lineages described herein (Fig. 2). For instance, the subspecies rosea and rubicoides are paraphyletic (lineages A&B, and C&D, from BP&P analysis, respectively). The first taxon is further intermingled with ssp. affinis individuals in lineage B, while the second one is combined with individuals of ssp. holobrunnea and nelsoni in the lineages C and D. Some subspecies also appear indistinguishable with the markers that we used (e.g., ssp. hesterna, peruviana and rhodinolaema within lineage F). Overall, our results support a taxonomic revision of H. rubica.

Our test for species delimitation suggests that H. rubica should minimally be split into the following species: (1) H. rosea, distributed along the Pacific coast of western Mexico (Jalisco, Nayarit and Colima; lineage A); (2) H. affinis, distributed along the Pacific coast of southern Mexico (Michoacan, Guerrero and Oaxaca; lineage B); (3) H. holobrunnea, distributed over the slope of Gulf of Mexico (lineage C); (4) H. rubicoides, distributed from eastern Mexico to southern Nicaragua (lineage D); (5) H. vinacea, allopatric in Panama (lineage E); (6) H. rhodinolaema, distributed in the Amazon basin (lineage F); and (7) H. rubica, distributed in eastern Brazil, Argentina, and Paraguay (lineage G). Further tests and sampling should be necessary for a more adequate description of the Venezuela’s individuals, currently classified as spp. perijana (Hilty, 2011).

There are criticisms to the BP&P analysis implemented herein. For instance, simulation data have shown that it tends to misidentify population structure as species boundaries (Sukumaran & Knowles, 2017). However, it must be noted that our delimitation not only detected genetics lineages, but lineages whose distribution coincides with the geographical barriers and ecological differences (habitat preference) of each taxon. That is, it highlights the role of other factors, like vicariance, dispersal and ecology, under a framework of ecological allopatric speciation, which is one the drivers proposed for generating bird species diversity in the neotropics (Hooghiemstra & Van der Hammen, 1998; Hewitt, 2000; Bennett, 2004; Smith et al., 2013). It also shows that this bird diversity is underestimated, and that species cannot be described or estimated based only on morphological criteria. From this point of view, diversifying the techniques for species delimitation, like incorporating bioacoustics, high-throughput sequencing, niche modelling and advanced statistical tools, will also help addressing consequential questions in taxonomy and biogeography. For instance, such multifactorial data should help determining operationally whether the H. rubica lineages inferred herein are indeed evolving separately (see De Queiroz, 2007).