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,55].
Furthermore, altitudinal gradients may shape broad biogeographical patterns and potentially population genetic variation but its influence is complex and variable among species [147,148, 149,150]. 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 [151]). 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 [42,121]). 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 cloud-forest 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 [151].
Phylogenetic and population structure analyses
Our results, as several others from the Neotropics (see [11,13,28,122,123]), 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 [121], plants [28], rodents [18,34], insects [48] and a snake [124]. 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 [125]. Isolation between Central American highland populations (NCA-SCA) has been likely promoted by the Nicaraguan Depression (see [13,123,126]). 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 [32], lizards [127,128], frogs [30] and even in lowland species, such as iguanas [129].
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 Pliocene-early Pleistocene [130]; 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 [50] which is consistent with an interglacial period during the Pre-Illinoian stage [130,131] 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 (FSC) 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 (Figures 2d, 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 [132,133], which is congruent in sky-island populations supporting lower nucleotide variation given reduced gene flow to other populations [49,134,135]. 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 FS 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 FST 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 FST fixation indices support the idea that these populations have followed their own evolutionary trajectory for a long time (Tables 3, 4). Despite the FST pairwise comparison between SMS and TMVB populations was low; according to [91] 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 [136].
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; [137]), 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; [137]), the IT dry plains (ca. 6-3 Myr; [36]), orogeny of SMS (ca. 35-20 Myr), and the formation of the TMVB as well as its period of volcanism in central Mexico (ca. 16-7 Myr; [37-39,138,139]). Our results suggest that these geographic breaks may have acted indirectly as semi-permeable barriers to dispersal, leading to an accumulation of processes involving isolation and expansion [140], 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 FST 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 [71,72], we included samples of T. ridgwayi in the diversification scenario comparisons, and we were unable to discern if diversification according to the best-supported 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 [47]).
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 [141]; 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 [28,47], 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,28,142]. 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 [143-145]. 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,63].