The geography of evolutionary divergence in the highly endemic avifauna from the Sierra Madre del Sur, Mexico

Background Mesoamerica is a remarkable region with a high geological and ecological complexity. Within northern Mesoamerica, the biotic province of the Sierra Madre del Sur (SMS) in southwestern Mexico harbors exceptionally high avian endemism and diversity. Herein, we searched for spatially and temporally concordant phylogeographic patterns, in four bird genera from three distinct avian orders co-distributed across Mesoamerica and investigated their causes through hypothesis testing regarding historical processes. Selected species include endemic and differentiated populations across the montane forests of Mesoamerica, and particularly within the SMS. Results We gathered mitochondrial DNA sequences for at least one locus from 177 individuals across all species. We assessed genetic structure, demographic history, and defined a framework for the coalescent simulations used in biogeographic hypothesis testing temporal and spatial co-variance. Our analyses suggested shared phylogeographic breaks in areas corresponding to the SMS populations, and between the main montane systems in Mesoamerica, with the Central Valley of Oaxaca and the Nicaragua Depression being the most frequently shared breaks among analyzed taxa. Nevertheless, dating analyses and divergence patterns observed were consistent with the hypothesis of broad vicariance across Mesoamerica derived from mechanisms operating at distinct times across taxa in the SMS. Conclusions Our study provides a framework for understanding the evolutionary origins and historical factors enhancing speciation in well-defined regions within Mesoamerica, indicating that the evolutionary history of extant biota inhabiting montane forests is complex and often idiosyncratic.


Background
Comparison of phylogeographic patterns across codistributed species provides valuable insights on the probable congruence of processes that may have driven intraspecific diversification in particular regions [1][2][3][4]. Several phylogeographic studies of individual bird species have shown that genetic structure is frequently associated to discontinuous ranges and geographical barriers that restrict admixture between populations, leading to genetic differentiation (see [5][6][7][8]) which may have proceeded independently or jointly with the divergence in morphological, behavioral or ecological traits (see [9][10][11][12][13]). However, a relatively lesser effort has been dedicated to provide robust biogeographic hypotheses seeking to explain shared distributional patterns, even when it may be expected that species sharing geographic areas should show congruent spatio-temporal patterns of differentiation [1,[14][15][16][17][18][19].
Situated at the northernmost Neotropics, Mesoamerica possess one of the highest levels of endemism and species diversity, but also one of the most globally threatened biota and high rates of deforestation [20,21]. The complex geological history of this region, as well as the cyclic changes in vegetation and climate, pose it as a challenging area for biogeographic and evolutionary studies; in addition, constant orogenic processes have promoted a highly broken topography characterized by highland isolated patches of humid montane forest between 600 and 3000 m, which includes both humid pine-oak forest and cloud forest [22][23][24]. This mosaiclike landscape has been associated to centers of diversification along elevational gradients, in which both a high species richness and endemism have evolved for the last 2 Myr, likely as a result of both environmental and geological complexity, as well as Pleistocene climatic fluctuations [6,15,18,22,[25][26][27][28][29], which is supported by relatively recent intraspecific differentiation processes in several groups of organisms (see [27,30]), therefore explaining the existence of numerous endemic species in different taxonomic groups, including birds [31][32][33][34][35].
The Sierra Madre del Sur (SMS) biogeographic province is included within the Mesoamerican highlands, in the Mexican Transition Zone (MTZ, [36]). This isolated mountain range is mostly surrounded by dry lowlands of the Isthmus of Tehuantepec (IT) to the east, the Pacific slope to the south and west, the Balsas Depression and the central Trans-Mexican Volcanic Belt to the north ( Fig. 1 [38][39][40][41];). The SMS with an overall extension of ca. 56,729 km 2 , is located along the Mexican Pacific slope (Jalisco to Oaxaca, 16-18°N, 95-102°W) [42]. The region hosts 622 avian species, 29 (4.6%) of which are classified as semi-endemic (species breeding elsewhere but for which the entire species' population spends the wintering season within the political limits of Mexico [43]), 15 (2.4%) as quasi-endemic (species whose geographic range extends no more than 35,000 km 2 to a neighboring country due to habitat continuity, [43]), and 54 (8.6%) are endemic to Mexico; in addition, 134 (21.5%) species are considered in some category of risk [34]. Furthermore, as in other Mesoamerican bird taxa [44][45][46][47][48], the SMS holds endemic and well-differentiated subspecies, suggesting the importance of isolation for the evolution of intraspecific variation [17,49].
Isolated montane forest patches are generally thought to act as promoters for divergence between populations of several species, which have served as models for the understanding of speciation processes in island-like settings [50][51][52][53][54]. These fragmented habitats influence the population dynamics and demography, given that resident populations typically show reduced effective population sizes [18,50,[55][56][57][58][59][60], they are prone to a reduced adaptive potential and increased extinction risk due to the stochastic effects of genetic drift [61][62][63][64][65][66][67]. Therefore, the study of evolutionary history of co-distributed species in these habitats is highly important for conservation purposes.
Herein, we compared phylogeographical structure patterns among four co-distributed and not closely-related Neotropical bird taxa with well-differentiated populations restricted to the SMS. We tested the following biogeographic hypotheses in an ABC framework: (a) a southern origin (southern Central America) followed by a northward expansion; (b) a Mexican origin and range expansion towards southern Central America; and (c) a northern Central American origin with north and southward dispersal. We also aimed to determine whether phylogeographic structure and divergence times in each taxon are congruent among taxa, therefore revealing patterns of simultaneous diversification. Specifically, we used mitochondrial datasets in a comparative phylogeography framework to test: (1) The presence of geographically-structured genetic variation; (2) the existence of spatially-congruent genetic variation across analyzed species; (3) if biogeographic history and similar diversification scenarios support vicariance, dispersal or population admixture. In addition, we also examined temporal patterns of diversification through hABC simulations, testing for either synchronous or multiple pulses for diversification in bird assemblages across the SMS region.

Phylogenetic and population structure analyses
All of our phylogenetic and genetic assignment analyses supported previous findings on the genetic structure for all of the taxa ( Fig. 2 and Fig. 3). Genetic structure is clearly delimited by landscape breaks corresponding to the main montane regions in Mexico and Central America. Nevertheless, we found that the detected genetic structure in these taxa is not coincident with the accepted taxonomy at the intraspecific level.
We found four distinctive lineages in Aulacorhynchus. Two of these, SMS (restricted to the SMS highlands) and EMx-NCA (widespread in Mesoamerica, including SMO, TM, Chiapas, and northern Central America) are differentiated by 27 mutational steps (Fig. 2c). The other two are distributed in southern Central America from Costa Rica to Panama. Phylogenetic analyses support this structure with ML bootstrap values and BI posterior probabilities above 80 and 0.9 respectively (Fig. 2e). We further observed two subclades in SMS samples, each representing Oaxaca and Guerrero populations.
For Chlorospingus, we observed six distinct mtDNA lineages, four of which included Mexican samples (Fig. 2d). Individuals representing the SMS (Sierra de Miahuatlán, Oaxaca and Guerrero) were recovered in a separated group from samples from northern Oaxaca, which clustered with samples from the eastern Sierra Madre Oriental highlands (SMO); both groups are separated by 26 mutational steps. A third group  corresponded to the Chimalapas (Oax.) and northern Chiapas region, while the fourth group corresponded exclusively to the TM, separated from the SMS by 32 mutations. This genetic structure is supported by BI, as well as by genetic assignment analyses with overall high posterior probabilities. Some nodes however, showed low bootstrap and posterior probability values (Fig. 2f). As observed in Aulacorhynchus, samples from the states of Oaxaca and Guerrero were recovered as distinct subclades within the SMS.
Cardellina versicolor, the sister species, is differentiated by more than 30 mutational steps from the C. rubra group. Within C. rubra, all recovered lineages coincide with previous findings, with one group found in the SMOc, another in the TMVB, and a third one in the SMS. We found 8 samples from a locality in central Guerrero (Carrizal de Bravo) from the SMS to be grouped in the TMVB clade, which also includes samples from localities of Jalisco and Michoacán. This paraphyletic relationship was recovered in all analyses, receiving strong support from both bootstrap (> 80) and Bayesian posterior probabilities (1; Fig. 3e). Nevertheless, Fig. 3 Geographical distribution of sampled haplotypes for (a) Cardellina and (b) Eupherusa species. c and d depict median-joining network of concatenated mtDNA loci for Cardellina and Eupherusa respectively, numbers on lines depict mutational steps between haplotypes, gray dots represent median vectors inferred for the data. e and f depict maximum clade credibility trees using BEAST with branch support (BI/bootstrap). Values between brackets indicate the 95% highest posterior densities (HDP) of the estimated times of divergence events (in Myr). Nodes that have no represented time frame of diversification are those whose lower and upper bounds of the HDP interval had posterior probabilities inferior to 0.5 genetic assignment analysis yielded a distinct result for these Guerrero samples, particularly for three of them (jk04-349, jk04-351, and jk04-353), which clustered with the SMS group; the remaining five clustered with the TMVB group.
Finally, in Eupherusa-Thalurania, we found five lineages corresponding to the currently recognized species. T. ridgwayi group is separated by 39 mutational steps from all of the Eupherusa samples; species from the SMS E. cyanoprhys and E. poliocerca haplotypes are separated by 22 mutational steps; and finally, species E. nigriventris and E. eximia are separated by 50 mutations (Fig. 3d). Topologies from BI, ML (Fig. 3f) recovered high posterior probabilities and bootstrap support values (> 0.8, > 75 respectively). All analyses recovered a monophyletic Eupherusa genus, where the SMS is represented by E. poliocerca and its sister species E. cyanophrys, and a clade from eastern Mexico to Central America, composed of E. eximia sister to E. nigriventris. Finally, for each of these species, we did not find evidence supporting further intraspecific genetic structure, even when their distribution, as in E. poliocerca (Oaxaca-Guerrero) and E. eximia (E. e. nelsoni, E. e. eximia, and E. e. egregia), has been described to be disjunct [28,68,69].
Hierarchical AMOVAs (Table 1) for the analyzed species showed that variation was explained by differences among groups. In this case, the fixation indices among groups (F CT ) were not significant, whereas fixation indices among populations within groups (F SC ) were. Genetic diversity varied among datasets. Overall, we obtained high haplotype diversity (Hd > 0.5) and low nucleotide diversity ( Table 2). The only low value of Hd was observed in E. cyanophrys. Fu's F S and Tajima's D values were not significant, excepting TMVB and SMS in Cardellina; however, we observed a tendency towards negative values (Table 2).
Overall, pairwise F ST values showed significant high genetic differentiation between analyzed groups. However, probably due to small sample sizes, some comparisons showed no significant values, as in the SP cluster of Aulacorhynchus (Table 3), as well as in comparisons between E. nigriventris vs T. ridgwayi, and E. nigriventris vs E. cyanophrys clusters ( Table 4). The lowest F ST value obtained occurred between Cardellina clusters from the TMVB vs SMS (Table 4).

Divergence time estimation and coalescent-based estimation of population history
Even when divergence dates varied, analyses showed that divergence between populations occurred during the Pleistocene, mainly within the last 2 Myr. The bestsupported scenario for four population groups in Aulacorhynchus was scenario 1 (Additional file 5: Figure S5), with a posterior probability (PP) of 0.57 and a 95%      Figure S8) received the better support (PP = 0.6, 95% CI: 0.5909-0.625; Fig. 4d), with high confidence (Additional file 3). Scenario 1 suggested diversification events similar to the latter, mean parameter estimates (  Fig. 4g), with high confidence (Additional file 3). Scenario 3 (Additional file 5: Figure S9) predicts an initial northward dispersal event from northern Central America (C. versicolor) population through the Isthmus of Tehuantepec at t4 (~1.3 Myr); followed by vicariance events driven by the effect of Mexican highlands at t3 (~426 kyr) between SMOc and an hypothetical ancestral population that at t2 (~260 kyr) is fragmented into the TMVB and SMS populations (Table 5). This scenario also supports a recent change in the TMVB effective population size during the Wisconsin Glacial stage at 40 kyr (95% CI: 3500 yr -122 kyr; Additional file 4 and Additional file 5).
The best scenario for Eupherusa including T. ridgwayi (Additional file 5: Figure S10 Table 5). When T. ridgwayi is excluded (Additional file 5: Figure S11), the best-supported was scenario 1 (PP = 0.868, 95% CI: 0.859-0.877; Fig. 4f), which is similar to the previous scenario. Confidence in scenario choice was high given low values of type I and II errors (Additional file 3). Mean parameter estimates of scenario 1 (Table 5)

Test of simultaneous divergence
The hABC analyses including four population pairs spanning along two identified putative barriers (Balsas Depression and Oaxaca's Valleys), showed a relatively strong posterior probability support for four independent divergence events (PP Ψ = 4 = 0.417) and the highest mode model-averaged posterior estimate of the hyperparameter Ω = 0.1415 (95% HDP interval = 0.0035-0.2998), in comparison with posterior probability of two (PP Ψ = 2 = 0.282) or three (PP Ψ = 3 = 0.275) divergence events (Fig. 5a). Subsequently, we tested subsets of the data set for simultaneous divergence. First we tested two taxon pairs within SMS (Guerrero and Oaxaca populations separated by the Río Verde drainage) of Aulacorhynchus and Chlorospingus, obtaining high support for a single divergence event (PP Ψ = 1 = 0.6749), along with a mode value of Ω consistent with synchronous divergence (Ω = 0.000475) and an Ω = 0, which is included within the 95% HDP intervals (Fig. 5b), whereas when including the Eupherusa taxon pair (cyanophrys-poliocerca), we found low support for a single divergence event, suggesting two divergence events across these taxa (PP Ψ = 2 = 0.5436) and a Ω mode of 0.05 consistent with asynchronous divergence (Fig. 5).
When evaluating southwestern SMS against eastern populations spanning Oaxaca's Central Valley barrier, there was a relatively strong support for a single divergence event between these populations (PP Ψ = 1 = 0.4536) and Ω = 0.0022, also indicating support for synchronous diversification between these taxa occurring at~1.06 Myr (Fig. 6d). Finally, the evaluation of population pairs distributed along Central America revealed a very high posterior probability of one divergence event (PP Ψ = 1 = 0.7228), and a Ω hyperparameter value = 0 (95% HDP interval = 0.0-0.1167) across the Nicaraguan Depression lowlands between Costa Rica and Nicaragua occurring 0.924 Myr (Fig. 6e).

Discussion
The aim of comparative phylogeography is to detect concordant patterns among co-distributed species, based on the idea that they may share common histories [1,14,16]. Thus, finding similar genetic patterns among different species may suggest similar historical processes influencing major biological components in a given region. Nevertheless, as shown in the results presented herein, complete concordance is commonly non-existent, which is probably a result of intrinsic characteristics of species (e.g., feeding habits) triggering different responses to particular events, therefore emphasizing the role of idiosyncratic events in the recent evolutionary history of genetic loci and bird taxa in the region [14,58].
Furthermore, altitudinal gradients may shape broad biogeographical patterns and potentially population genetic variation but its influence is complex and variable among species [70][71][72][73]. Models show that cold and humid conditions were extended along the Mexican highlands during the Pleistocene; nevertheless the dynamics under which each taxon responded to the climatic cycles depends on several factors (e. g., time of establishment in an area, environmental preferences, dispersal capabilities [74]). Our study also indicates that at least two unrelated bird genera (Aulacorhynchus & Chlorospingus) with distinct natural histories but which inhabit roughly the same regions and wide elevational range along the humid montane forests of Mesoamerica, present highly similar phylogeographic structure, including a significant split between the highlands in Guerrero and Oaxaca within the SMS, which may be expected due to discontinuous distribution ranges, as well as to their presumably reduced dispersal abilities. The pattern observed in Cardellina and Eupherusa do not fully agree with the expectation that the degree of isolation between montane ranges is responsible for shaping the genetic differentiation between populations, as has been reported in previous studies of montane birds (see [45,75]). The analyses in these two taxa evidenced exclusive haplotypes from Oaxaca, but genetic differentiation across the Río Verde drainage, as seen in Aulacorhynchus and Chlorospingus does not occur. Regarding Cardellina, despite the fact that the TMVB and the SMS montane ranges are separated by a significant geographic barrier such as the Balsas Depression (see [7]), incomplete lineage sorting between both regions is reported. These differences may also be due to the fact that both taxa inhabit different elevational ranges along the region: Eupherusa ranging to the lower limits of the cloudforest and Cardellina ranging in the upper reaches of the humid montane forest, suggesting that isolation and connectivity cycles in the vegetation types might have a differential impact on the genetic pool of birds depending on the elevational range, which might explain the observed patterns [74].

Phylogenetic and population structure analyses
Our results, as several others from the Neotropics (see [11,13,18,76,77]), found shared phylogeographic breaks at which genetic differentiation occurs among widespread Mesoamerican highland bird species. We observed major breaks, as seen in Aulacorhynchus and Chlorospingus, between Chiapas and Mexican-Central American lineages. This has been also supported by other bird species [75], plants [18], rodents [20,35], insects [51] and a snake [78]. Moreover, phylogeography shows isolation between populations inhabiting the mountains of southwestern Mexico (SMS) and northern Oaxaca (see Chlorospingus) evidencing that populations are clearly sorted within Oaxaca across its Central Valley, and more importantly morphological variation of these populations has also been documented as in E. eximia nelsoni [79]. Isolation between Central American highland populations (NCA-SCA) has been likely promoted by the Nicaraguan Depression (see [13,77,80]). The most relevant break in our study, occurred within the SMS region, between the states of Oaxaca and Guerrero highlands, which are separated by the Río Verde drainage, this break is supported by studies in plants [33], lizards [81,82], frogs [31] and even in lowland species, such as iguanas [83].
The patterns observed in the haplotype networks mostly reveal a consistent match between haplogroups and subspecies, as well as species boundaries (see Eupherusa). The overall structures in the haplotype networks resemble an expected pattern for ancient divergence among the sky-islands where reciprocally monophyletic groups occur on each sky-island with no mixing (see Aulacorhynchus and Chlorospingus), and is consistent with estimated divergence and the gradual cooling occurring in the Gelasian stage of late Plioceneearly Pleistocene [84]; whereas in Eupherusa it is difficult to discern specific events for each divergence estimate due to wide confidence intervals. Cardellina is an exception which resembles a star-shaped network of a post glacial sky-island divergence [53] which is consistent with an interglacial period during the Pre-Illinoian stage [84,85] and the estimated divergence in BEAST and DIYABC analyses.
AMOVA results highlighted the existence of genetic structure for all studied species. Most molecular variation was found among groups; nevertheless, the fixation indices (F SC ) among populations within groups were significant, as in Aulacorhynchus, Chlorospingus and Cardellina. This pattern suggests the existence of further structured populations in at least one of the defined clusters, as supported by haplotype networks (Fig. 2b  and Fig. 3d).
We obtained overall high values of haplotype diversity (Hd) and low nucleotide diversity (π) for the majority of analyzed groups within each taxon. These patterns have been attributed to population growth following demographic bottlenecks enhancing retention of novel mutations [86,87], which is congruent in sky-island populations supporting lower nucleotide variation given reduced gene flow to other populations [52,88,89]. Population expansions may have been promoted by climatic oscillations during the Pleistocene; nevertheless, we found significant support for population expansion only in Cardellina populations from the TMVB, and SMS. For the rest of the taxa, Tajima's D and Fu's F S values suggest that the observed nucleotide polymorphism is selectively neutral even when they tend towards negative figures, thus favoring demographic stability during the Pleistocene. Pairwise comparisons of F ST values for each taxon are high and mostly significant, which is expected given their disjunct geographic distribution. Overall, all populations concerned are allopatric and several lineages possess phenotypic and phylogenetic identity, thus the F ST fixation indices support the idea that these populations have followed their own evolutionary trajectory for a long time (Tables 3, 4). Despite the F ST pairwise comparison between SMS and TMVB populations was low; according to [90] these populations have great genetic differentiation. Moreover, some populations, despite allopatry, may be geographically close enough to show evidence of intermittent connections in response to environmental changes [91].

Divergence time estimation and coalescent-based estimation of population history
Our results highlight the importance of Pliocene-Pleistocene events in promoting the intra-specific genetic structure in the analyzed species, as reported in other studies. Divergence between SCA and northern populations is concordant with the formation of the Nicaragua Depression (ca. 4 Myr [92];), furthermore, divergence between NCA and Mexican lineages post-date the time frame of the formation of the volcanic arcs in Central America (Miocene-Pliocene or before [92];), the IT dry plains (ca. 6-3 Myr [38];), orogeny of SMS (ca. , and the formation of the TMVB as well as its period of volcanism in central Mexico (ca. 93,94];). Our results suggest that these geographic breaks may have acted indirectly as semipermeable barriers to dispersal, leading to an accumulation of processes involving isolation and expansion [95], as the Cardellina analyses revealed.
DIYABC scenarios allowed us to infer evolutionary processes affecting populations. Although estimated divergence dates varied among populations of the analyzed taxa, major splits occurred during the Pleistocene. Scenarios from Aulacorhynchus, Chlorospingus, and Cardellina, along with their respective pairwise F ST comparisons, support the hypothesis of sequential northward population range expansions from Central America with ancestral populations being split by vicariant events, such as fracture and changes in altitudinal ranges of montane forests linked to climatic oscillations, therefore promoting geographical isolation. Although our analyses support a southern origin, the pattern observed in Eupherusa + Thalurania ridgwayi suggest alternative scenarios as also probable. Based on previous evidence supporting T. ridgwayi as a member of the Eupherusa genus [96,97], we included samples of T. ridgwayi in the diversification scenario comparisons, and we were unable to discern if diversification according to the bestsupported scenario occurred from two potential regions that may have served as refugia (one in northwestern Mexico, another in southern Mexico-Central America), or whether SMS (E. poliocerca-E. cyanoprhys) and Central American lineages of Eupherusa come from a migration event from northwestern Mexico. Further studies should focus on the reconstruction of ancestral areas of distribution to clarify this issue (see [50]).
Overall, dates from the DIYABC analyses suggest more recent diversification times than our BEAST analyses, thus highlighting the effect that distinct assumed mutation rates have an effect of temporal frame estimation [98]; nevertheless, these values overlap their respective 95% HDP confidence intervals. Even in vicariant events, differences in time estimates can also be due to variance in the coalescent process and related to the demography of each species [18,50], however we acknowledge the inconveniences of single locus based analyses and that distinct calibration for divergence time estimation may yield potential errors, thus future surveys incorporating more data should be able to obtain stronger and more accurate estimates of timing and synchrony of divergence.

Test of simultaneous divergence
Even when species ranges overlap across their distributional areas, some divergence time estimation frames are not congruent with each other, therefore, asynchrony in diversification within Mesoamerica could be depicted as a result of Pleistocene climatic oscillations which promoted rapid pulses of diversification and speciation rates across multiple lineages [13,18,19]. Our analyses including four population pairs detected non-simultaneous divergence events within the Pleistocene, pointing to the existence of several temporal frames where diversification and fixation of genetic differences within the same montane ranges occurred among distinct lineages. Our results also highlight a high probability of simultaneous divergence occurrence between two taxa (Aulacorhynchus and Chlorospingus) in the SMS region across the Río Verde drainage lowland barrier, with a mean divergence date of 0.9339 Myr, which 95% HDP confidence intervals falls within the estimated divergence dates between the same populations. When Eupherusa population pair (cyanophrys-poliocerca) was added, estimations supported two diversification events, thus implying that even when the three taxa share their overall distributional ranges in the SMS, they do not fully share a community-wide historical pattern across the landscape.
Even when we do not fully out rule any orogenic effect over the analyzed species, major changes to population structure should be promoted by the onset of climatic oscillations in the Pleistocene, which have been invoked to explain patterns of genetic differentiation in the increasing rates of speciation and diversification in Neotropical biota, given the evidence of a downward altitudinal migration of the forest line which resulted in a montane forest-like dominated landscape in the lowlands with further retreat when temperatures raised [99][100][101]. Although for the C. rubra lineage the observed strong signature of expansion in TMVB and SMS populations due to Pleistocene conditions and forest migrations resulted in the admixture of two isolated lineages; our results reveal that demographic dynamics for other species persisted in isolated regions throughout the Pleistocene allowing the evolution of high endemism, which could be an overall result of high effective population sizes.
Genetic endemism as a result of isolation of intraspecific lineages, as seen within the SMS should be considered in conservation measures given that they preserve historical components and maintain the species capability of adaptation [3,66].

Conclusions
Our survey provides a framework for understanding the complex geographic history of Mesoamerica and its role in promoting lineage diversification in birds. From our analyses we recover multiple independently evolving lineages restricted to montane areas across the region, our results also highlight that congruence may be difficult to occur given idiosyncratic life histories of species; nevertheless, co-distributed taxa will share more common history than a random scenario could predict. Montane regions within Mesoamerica are strong drivers of lineage diversification, and as such our results point into major areas of great genetic diversity and endemism, such as isolated highlands of the SMS in the Mexican states of Oaxaca and Guerrero. Given ancient development of major montane systems, climatic oscillations and orogeny derived conditions were key drivers of diversification between these lineages as splits were calculated to occur mainly within the Pleistocene.

Phylogenetic and population structure analyses
For each taxon, we retrieved mitochondrial DNA (mtDNA) sequences deposited in GenBank (https:// www.ncbi.nlm.nih.gov/genbank/), and conducted analyses using alignments of concatenated mtDNA loci (Additional file 1). We estimated nucleotide substitution models and partition schemes for each species in PAR-TITIONFINDER [110], using the Bayesian Information Criterion (BIC) for model selection. Given the10 bp frameshift overlap in the ATPase 8 and 6 genes and the tRNA-Lys gene, we analyzed these as a single genetic region [111,112]. Resultant partition schemes and model substitution parameters were used for conducting a phylogenetic reconstruction using the Bayesian inference approach (BI) for each species, as implemented in MRBAYES 3.2 [113]. We ran two independent searches using four Markov-Chains MonteCarlo for 10 7 generations sampling every 1000 generations. Convergence across runs was evaluated using two methods: a) the examination of the standard deviation of split frequencies (with acceptance values < 0.01); and b) by verification of parameter estimates in TRACER v1.6 [114], based on acceptable effective sample sizes (ESS values > 200). After checking for convergence, the first 25% of the generated trees were discarded as burn-in and the remaining 75% were kept to calculate posterior probabilities. In addition, we also estimated phylogenetic trees using maximum likelihood (ML) criteria, as implemented in RAXMLGUI 1.5b1 [115][116][117], using the GTRCAT model, and estimated nodal support via 1000 bootstrap iterations using the selected partition. We selected closely-related taxa as outgroups, as suggested by previous published studies: Aulacorhynchus albivitta [108], Chlorospingus flavopectus phaeocephalus [103], Cardellina versicolor [5], and Thalurania ridgwayi [96,97].
We used NETWORK 4.6.1.1 (Fluxus Engineering, www.fluxus-engineering.com) to visualize the relationships among haplotypes by constructing networks using the median-joining algorithm [118,119], assigning equal weights to all variable sites and an epsilon parameter with default values (ε = 0).
To explicitly test for phylogeographic structure in each taxon avoiding a priori criteria to delineate populations, we used a Bayesian model-based clustering algorithm, as implemented in GENELAND 4.0.3 [120,121], which assigns samples to clusters (K) according to both geographical adjacency and genetic similarity through simulations with the Reversible Jump Markov Chain MonteCarlo (RJMCMC) algorithm. We varied the maximum number of expected clusters (K max ) from 2 to 10, following possible sub structuring of populations within the main montane ranges in Mesoamerica. We performed 10 independent runs of 10,000 iterations each, a thinning value of 1000 and a 10% burn-in. Best results were selected according to highest posterior probabilities.
Genetic structure was assessed through hierarchical analysis of molecular variance (AMOVA) using pairwise differences, based on the number of clusters (K) obtained by GENELAND for each taxon. We also assessed genetic diversity of each population group within each species through the estimation of haplotype diversity (Hd), and nucleotide diversity (π). Genetic divergence between groups was measured using pairwise F ST fixation index [122,123], interpreting the results following guides in [90]. Significance of F ST tests was assessed using 1000 permutations. To test for evidence of recent demographic changes in the selected species, we estimated historical demographic dynamics through the calculation of Fu's F S statistic [124] and Tajima's D statistic [125] neutrality tests. Significance of these tests (p < 0.02 in the case of F S statistic) was calculated through 1000 simulations. All analyses were conducted in ARLEQUIN 3.1 [126].

Divergence time estimation
We used a Bayesian MCMC-based approach to calculate divergence times among haplogroups within each species independently using BEAST v1.8.4 [127]. For each selected taxon, we first tested whether our dataset fits either to a strict clock model or to a relaxed clock model. We performed selection tests through the steppingstone method (SS [128];), as implemented in MRBAYES 3.2 [113]. Given our partitioning models, the mean marginal likelihood of a strict clock performed better than a relaxed clock model for all taxa (Additional file 2). Thus, chains were run using a UPGMA starting tree, under a strict clock with substitution models (according to results from PARTITIONFINDER) for 10 8 generations sampling every 1000 steps, using a Yule speciation process [129,130] with no topological constraints, and discarded the first 25% as burn-in. We used two approximations to convert branch lengths into time: (A) a body mass-corrected molecular clock rate of 0.0042 (min = 0.0011, max = 0.0158) subs/site/Myr [131], and (B) an uncorrelated lognormal relaxed clock substitution rate fixed at 0.01 average subs/site and SD = 0.003 [132][133][134] for each locus calculated in BEAST v1.8.4 [127]. Adjustment of the body mass-corrected and calculated clocks was evaluated using Bayes factors (BF), calculated from the marginal likelihoods from path sampling (PS) and stepping-stone (SS) methods in BEAST v1.8.4 [127,135,136]. Each marginal likelihood was estimated through 100 path steps with a Beta distribution (0.3, 1.0). We considered a 3 Log ml difference as strong evidence [137] against the null hypothesis (Additional file 2). We checked for stationarity using TRACER v1.6 [114]. Node ages are presented as mean heights and 95% credibility interval values with a posterior probability limit of 0.5 and resumed in a maximum clade credibility tree (MCCT). All of these were generated using TREEAN-NOTATOR v1.8.4 [127], and trees were visualized in FIGTREE v1.4.2 (http://tree.bio.ed.ac.uk/software/figtree/).

Coalescent-based estimation of population history
We tested for different hypotheses of population divergence and admixture using an approximate Bayesian computation (ABC) approach in DIYABC 2.1.0 [138]. We conducted initial simulations using 15 competing evolutionary scenarios per species (Additional file 3). Evolutionary scenarios for each species were built considering results of the phylogenetic and GENELAND analyses, biogeographic diversification scenarios proposed in previous studies (see [68,106,109]), and geographical breaks invoked as responsible of diversification events. We followed recommendations in [50,139,140] for the assessment of appropriate priors, as well as to select the highly-informative summary statistics from the simulations resembling datasets similar to the empirical ones. When phylogenies suggested the existence of two genetic groups within the SMS (p. e. distinctive Guerrero and Oaxaca lineages, as in Aulacorhynchus and Chlorospingus), we contrasted the same scenarios considering those clades as different populations in the analyses. All competing scenarios were eliminated in successive rounds, where the preferred hypotheses were the ones with highest posterior probabilities. From the selection process described above, we tested three final models for Aulacorhynchus, Chlorospingus, and Cardellina; and two for Eupherusa-Thalurania (Additional file 5). For each matrix analyzed, we simulated 1 × 10 6 datasets and obtained summary statistics per scenario in each simulation. Based on rates reported previously for birds [141][142][143][144][145], we used an HKY mutation model with a uniform prior distribution, a mean mutation rate with a gamma distribution set to 2.0 × 10 − 8 substitutions/site/year (min = 1.6 × 10 − 8 , max = 2.9 × 10 − 8 substitutions/site/year) The posterior probabilities of competing scenarios were computed using a logistic regression on the 1% of simulated datasets closest to the observed data. The selected scenario was that with highest probability value and a non-overlapping 95% confidence interval. For the best supported scenario, we performed a model-checking procedure using a Principal Components Analysis (PCA) on test statistics to visualize the fit between simulated and observed datasets. Confidence on the chosen scenario was assessed though the calculation of type I and type II error rate, from 500 simulated pseudo-observed datasets (PODs) generated with the data of the best-supported scenario [146,147]. Point estimates for demographic and temporal parameters were obtained by local linear regression on the 1% of simulations closest to the observed dataset for the best-supported scenario [138]. Divergence time obtained from DIYABC output was transformed assuming a conservative two-year generation time, which has been previously used in similar species groups [66,148].

Test of simultaneous divergence
We tested whether genetic differentiation occurred simultaneously between the main geographical barriers in the analyzed clusters by performing a hierarchical approximate Bayesian computation (hABC) analysis as implemented in MSBAYES [4]. Populations inhabiting the montane regions throughout the study area showed similar patterns of genetic structuring even when they are not closely related, we may therefore expect interspecific simultaneous divergence among regions if the processes driving diversification are common. We performed hABC analyses in 15 population pairs that spanned the same putative barriers to gene flow: between populations within the SMS in southwestern Mexico; between populations of the SMS and eastern populations across a lowland barrier in Oaxaca's Central Valleys, and between populations spanning the Nicaraguan Depression in Central America (Fig. 1).
We used jModelTest [149,150] to estimate a transition-transversion rate for each population pair, as required by MSBAYES. The analyses involved the estimation of a vector of summary statistics from the sequences utilized; afterwards 1 × 10 6 data sets were simulated under the specified multi-taxon model using prior distribution for the demographic parameters. The prior of the maximum possible number of divergence events (Ψ) was set to be equal to the number of lineage pairs tested. For the last stage of the MSBAYES analyses, we used the acceptance/rejection algorithm to approximate the posterior distribution for the hyper parameters that characterize the degree of variability of demographic and temporal parameters given the empirical data (e. g., τ, Ψ, Ω). We used Ω estimates to evaluate the support for each hypothesis: synchronous diversification is expected when Ω ≤ 0.01; asynchronous diversification if Ω > 0.01 [4,18,151,152]. To estimate mean divergence times, we converted model-averaged E(τ) estimates (provided in MSBAYES as coalescent units of 4N AVE generations) to absolute time (T div ), using the equation T div = E(τ) × (θ AVE /μ) × g, where g is the generation time in years, and θ AVE is the mean of the upper θ prior [151], and μ is the mutation rate. We used a mean mutation rate of 0.0042 substitutions/site/Myr prior [131], and an average generation time of 2 years [66,148].
Additional file 1: Genbank accession numbers. Samples obtained for the realization of analyses in this study.
Additional file 2: Marginal likelihoods for sequence evolution. Marginal likelihoods estimation for sequence fit to a strict or relaxed clock and comparison of evolutionary rate assumptions through Bayes factors (BFs).