For the ingroup we used 19 tissues from Euphonia affinis affinis and four from E. affinis godmani; for the outgroup we obtained one tissue sample from Chlorophonia occipitalis, two from Euphonia chlorotica, one from E. luteicapilla and one from Haemorhous mexicanus (Suppl. material 1, Table S1). We completed our outgroup dataset with sequences from Genbank of Euphonia chlorotica, E. finschi, E. lanirostris, E. minuta, E. musica, E. violacea, E. xanthogaster, Fringilla montifringilla, F. teydea, Carpodacus vinaceus and Haemorhous mexicanus (Suppl. material 1, Table S2).

Genomic DNA was isolated using the Qiagen DNeasy^TM kit (Qiagen Inc., Valencia, CA, USA) following the manufacturer’s protocol. We amplified five molecular markers: one mitochondrial ND2 (NADH Dehydrogenase Subunit 2, Sorenson et al. 1999) and four nuclear genes ODC (Ornithine Decarboxylase, Allen and Omland 2003), GAPDH intron 11 (Glyceraldehyde-3-phosphate dehydrogenase, Friesen et al. 1997) , BRM intron 15 (BRM transcription regulatory protein, Marthinsen et al. 2009), and MUSK (Muscle, skeletal receptor tyrosine-protein kinase, Kimball et al. 2009). Technique specifications are given in Suppl. material 1, Table S3. Amplifications were done via PCR in 12.5 μl reactions. PCR and products were visualized on a 1% agarose gel. DNA sequencing was performed by the High-Throughput Genomics Unit service (University of Washington). We edited and aligned chromatograms in Sequencher v4.8 (GeneCodes Corporation, Ann Arbor, MI). All sequences were deposited in GenBank; with the accession numbers ND2: MT452146–MT452170, ODC: MT452191–MT452215, GAPDH: MT452124–MT452145, BRM: MT452098–MT452123 and MUSK MT452171–MT452190. (Also included in https://github.com/almamelisa/Euphonia-affinis-complex).

We constructed the alignment of each gene using the CLUSTAL IW (Thompson et al. 1994) function in BIOEDIT (Hall 1999). Then, we used JModelTest 0.1.1 (Posada 2008) to estimate evolutionary models for each molecular marker. We conducted phylogenetic analyses in MrBayes v3.2.3 (Huelsenbeck and Ronquist 2003) for Bayesian Inference (BI), using a partitioned dataset including the four nucleotide genes ODC, MUSK, GAPDH, and BRM; ND2 was partitioned by codon position (see Table 2). For the BI, we ran the process for 10 million generations, sampling every 1,000 generations. We examined the convergence of the chains with Tracer v1.6 (Rambaut and Drummond 2013) and discarded the first 25% of generations as burn-in.

To resolve heterozygotes in nuclear sequences, we used a Bayesian approach in PHASE v 2.1 (Stephens et al. 2001), selecting the pairs of haplotypes with a posterior probability higher than 0.90. Then, we used DnaSP v5.0 (Librado and Rozas 2009) to estimate the number of haplotypes (H), haplotype (Hd) nucleotide diversities (π), and Fst. Genetic distances were obtained with MEGA 5 (Tamura et al. 2013). Finally, we obtained the haplotype network for each gene with the median-joining algorithm using Network v4.6 (Bandelt et al. 1999) through DnaSP.

Divergence times were estimated from the multilocus dataset with three genes, the mitochondrial gene ND2, and the two nuclear genes with sequences for all samples and outgroups (ODC and GAPDH) using Beast v1.8 (Drummond et al. 2013). As point calibration, we used a secondary dating based on the divergence between Fringillidae and the New World nine-primaried Oscines 17.1104 My with a 95% HPD of 14.7743–19.6278 calculated by Oliveros et al. (2019). We assigned a Normal distribution to the secondary dated point. Additionally, we defined partitions with different evolutionary rates corresponding to each gene fragment (ND2 = 0.029 s/s/My, ODC = 0.0014 s/s/My and GAPDH = 0.0012 s/s/My), based on Lerner et al. (2011). We used a Yule model as a tree prior (Gernhard 2008). For the molecular clock model, we selected the normal relaxed molecular clock following Drummond et al. (2006) and Li and Drummond (2012). We performed 20 million generations, sampling every 1,000 and corroborating the appropriate effective sample size (ESS > 200) with TRACER v1.6 (Rambaut and Drummond 2013). Finally, in TREE ANNOTATOR v 1.8.0 (Rambaut and Drummond 2013), we did a burn-in of 5,000 trees and produced the maximum clade credibility tree with 95% highest probability densities. The tree was visualized with FigTree v1.4.2 (Rambaut 2014).

Six morphometric measurements of 355 specimens (233 males and 122 females) were taken from the following collections (see Suppl. material 2, also included in https://github.com/almamelisa/Euphonia-affinis-complex): Museo de Zoología Alfonso L. Herrera (UNAM), Colección Nacional de Aves-Instituto de Biología (UNAM), Moore Laboratory of Vertebrate Zoology, American Museum of Natural History, Louisiana Museum of Natural History, Academy of Natural Sciences, Museum of Comparative Zoology, and Delaware Museum of Natural History. The morphometric measurements (following the recommendations of Baldwin et al. 1931) were: bill length (BL, from the upper base of the bill to the tip of the upper mandible), bill width (BW), bill depth (BD, from the upper mandible to the base of the bill at the distal edge of the nostrils), wing chord (WC, distance from the carpal joint the tip of the longest primary), tarsus length (TL), and tail length (TLE, distance from the uropygial gland to the tip of the longest rectrix). All measurements were taken only by the first author to avoid bias in the process using a dial caliper with a precision of 0.1 mm, except for tail length, which was taken with a millimeter ruler and in three independent events. To obtain our final data set we averaged the three independent events for every measure. Since both our molecular phylogenetic results and our analysis of the previously proposed plumage color differences do not provide evidence of E. affinis olmecorum as an independent evolutionary lineage, we decided to analyze morphometric variation, vocalization, and ecological niche only between E. affinis affinis and E. affinis godmani.

The normality was tested with the Shapiro-Wilk test of normality in R (R Core Team 2017). Since, the normality was rejected in all except one of the groups, we evaluated the sexual dimorphism with the Unpaired Two-Samples Wilcoxon Test, also with basic R functions. We obtained significance differences between male and females in three variables WC, TLE, and BD (see results), so we evaluated these variables differences between the lineages in a separated way for males and females with the Unpaired Two-Samples Wilcoxon Test. The rest of the variables were evaluated jointly for both sexes and also with the Wilcoxon test. With the previous arrangement, we also did a Principal Component Analysis based on the correlation matrix with the R package Factominer (Lê and Husson 2008). Graphs were generated in R package Factorextra (Kassambara and Mundt 2016). All the scripts and input data are in https://github.com/almamelisa/Euphonia-affinis-complex.

We obtained 19 recordings of Euphonia affinis calls from the Xeno-Canto (XC; http://www.xeno-canto.org) open access database. We used only call recordings in which the subspecies was identified and in which Euphonia was identified as the foreground species. We visualized and measured spectrograms of these recordings using the Raven Pro 1.6 software (Cornell University, Ithaca, NY). We visually inspected the spectrograms, and from each recording we selected one call section that did not overlap with any background vocalizations or other sounds. In recordings where more than one call variant occurred (for example, variants with differing number of notes), we selected one of the most frequent type. The most common call type for this species consists of a short series (2 to 4 notes) of whistled notes with decreasing pitch. Since recording conditions were not standardized, we only took frequency and duration measurements, which are not heavily affected by distance. We measured low and high frequencies (LowFreq and HiFreq), change in frequency (DeltaF), duration of call (DeltaT), number of notes (Notes) and emission rate (Speed; number of notes divided by duration). All measured variables were rescaled by log transforming them.

Unpaired Two-Samples Wilcoxon Test were carried out on individual variables to test for differences between the two groups. We also performed a principal component analysis (PCA) to explore the relation between the two groups in multivariate space. All the scripts and input data are in https://github.com/almamelisa/Euphonia-affinis-complex.

The georeferenced records were obtained from the specimens used in the morphometric and genetic analyses, 102 for E. affinis affinis and 29 for E. affinis godmani. To define the M area (accessibility area; sensu Soberón & Peterson, 2005) for each evolutionary lineage herein identified, we plotted the record points onto the biogeographic provinces of the Neotropical region (Morrone 2014) and chose the provinces that matched the record points for both lineages, using the shapefiles provided by Löwenberg-Neto (2014). Such considerations assumed that these regions may define the accessible historical area and specific restriction region for each lineage (Svenning and Skov 2004).

For the first explorative analysis, we used the 19 bioclimate layers from WorldClim and assessed which variables were the most important for the model, according to the Jackknife test calculated in MaxEnt (Royle et al. 2012). In a second modeling exercise, we generated the species models using those non-correlated (r < 0.8) environmental variables in combination with the most relevant environmental variables identified in the first approach. According to previous published works (Ortega-Andrade et al. 2015, Hernández et al. 2018), these additional steps allowed us to reduce overfitting of the generated distribution models, minimizing the collinearity problems among variables (Dormann et al. 2013). Pearson correlation test among bioclimatic variables was performed in R with the basic commands. Final models were performed considering only those 12 climatic variables: BIO3 = Isothermality (BIO2/BIO7) (* 100), BIO5 = Max Temperature of Warmest Month, BIO6 = Min Temperature of Coldest Month, BIO7 = Temperature Annual Range (BIO5-BIO6), BIO8 = Mean Temperature of Wettest Quarter, BIO9 = Mean Temperature of Driest Quarter, BIO10 = Mean Temperature of Warmest Quarter, BIO14 = Precipitation of Driest Month, BIO15 = Precipitation Seasonality (Coefficient of Variation), BIO16 = Precipitation of Wettest Quarter, BIO 18 = Precipitation of Warmest Quarter, and BIO 19 = Precipitation of Coldest Quarter. The climatic layers were used in ascii format and ~ 1 km resolution, and they were cut to the shape of the M area using the R package Raster (Hijmans 2019). We generated the final Ecological Niche Model (ENM) for each lineage using 75% of the record points as training data and 25% as testing data. We performed 25 replicates, 500 iterations, with 0.00001 as the convergence limit and 0.5 prevalence.

To evaluate the models we calculated the partial ROC (Receiver Operating Characteristic) in the web tool Niche Tool Box (https://shiny.conabio.gob.mx:3838/nichetoolb2/), the parameters were 0.05 proportions of omission, 50 random points percentage and 500 iterations. Also, in MaxEnt, we made four models projections one in the M area of each lineage and the remaining three were made to obtain the paleodistribution; we projected the ENM in the last maximum glacial period (~ 22,000 years ago) considering two general circular models: MIROC-ESM (Hasumim and Emori 2004) and CCSM (Collins et al. 2004). We also projected the ENM onto the last interglacial period (~120,000–140,000) (Otto-Bliesner et al. 2007). Finally, using the R packaged Ecospat (Di Cola et al. 2017) we evaluated the ecological overlap between the lineages, to quantify equivalence and similarity among groups using the Broennimann´s method (Broennimann et al. 2012). It consists of three steps: the first one is to calculate the density of occurrences and environmental factors across the axes of the environmental principal component analysis; the second is to evaluate the superposition niche along the gradient of multivariate analysis. Finally, the Schoener’s D observed (Schoener 1968) and the statistical similarity I observed (Warren et al. 2008) are compared with the 100 repetitions of randomly generated simulated values for D and I (Warren et al. 2008, Broennimann et al. 2012). This final step consists in testing two hypotheses, the equivalence and similarity among groups. The hypotheses of niche equivalence and similarity are rejected if the empirically observed D and I values are significantly different from the values expected from the pseudoreplicates. All the scripts and input information are in https://github.com/almamelisa/Euphonia-affinis-complex.

The multilocus dataset analyses revealed a well-supported monophyly for the Euphonia affinis complex and recovered two main phylogroups: one included the samples from western Mexico (E. affinis godmani) and the other comprised samples from eastern Mexico and Central America (E. affinis affinis and E. affinis olmecorum) (Fig. 1). The sister group of E. affinis complex was the clade including E. chlorotica, E. luteicapilla, and E. finschi. As shown in Table 2, the genetic distances between E. a. affinis and E. a. godmani have values similar to genetic distances between the other species.

The haplotype network obtained with ND2 sequences showed two geographically structured haplogroups: a western group and an eastern-CA group (Fig. 1), respectively, with two haplotypes of E. affinis godmani and 10 that included samples of E. affinis affinis and E. affinis olmecorum separated by 77 permutations. ND2 had the highest total haplotypic diversity (see Table 1). These same two groups were obtained for ODC, MUSK, and BRM genes. The only exception was the GAPDH network, which did not recover these geographically structured groups. Also, in Table 1 we show the results of nucleotide diversity, Tajima’s D, nucleotide composition, molecular evolutionary models, parsimony informative sites, monomorphic sites, and the alignment base pairs.

Euphonia affinis godmani and E. affinis affinis split 2.6 Mya (1.5–4.0 Mya, 95% HPD), during the Late Pliocene-Early Pleistocene (Fig. 2). According to our analyses, the family Fringillidae is divided in three subfamilies: the oldest, Fringillinae, originated 14.21 Mya (10.2–17.9 Mya, 95% HPD, Highest Posterior Density), the split between Carduelinae and Euphoniinae was 12.9 Mya (9.2–16.8, 95% HPD), the Euphoniinae origen was 8.5 Mya (5.9–11.2 Mya, 95% HPD), and Carduelinae diverged 8.1 Mya (4.9–11.4 Mya, 95% HPD. Our estimate for the split between Fringillidae and Plectrophenas nivalis was 16.38 Mya (13.3–19.5 Mya, 95% HPD) while our point of calibration was 17.1104 Mya with a 95% HPD (14.7743–19.6278; Oliveros et al. 2019). Also, our results for the Fringillinae node age and the split between Carduelinae and Euphoniinae are consistent with the ages calculated for Euphonia and Chlorophonia phylogeny in Imfeld et al. (2020). However, there are also some differences between the ages calculated by us and by Imfeld et al. (2020), since they calculated that the split between E. affinis and E. luteicapilla was less than 1 Mya, whereas we calculated that the split between E. affinis and the rest of Euphonias was 4.3 Mya ago, and E. affins seems to be a sister group of E. chlorotica, E. luteicapilla and E. finschi.

Figure 2. 

Ultrametric phylogenetic tree obtained by BEAST using ND2, ODC, and GAPDH concatenated matrix. The rhombus node represents the calibration point 17.1104 My with a 95% HPD of (14.7743, 19.6278) (see methods), dark gray circle node represents the E. affinis origin and light gray circle node represents the break between E. a. godmani and E. a. affinis. Above the branch the diversification dates (My) and in brackets the 95% HPD. Below branch the number indicated the posterior probability. The green area corresponds to the period when lowland dry forests had a greater expansion in Western Mexico.

A total of 355 specimens was analyzed, of which 180 were males and 97 females of E. affinis affinis, and 53 males and 25 females of E. affinis godmani. Morphometric sexual dimorphism was found in three variables: TLE (Tail Length), WC (Wing Chord), and BD (Bill Depth) (Fig. 3, Table 3). Males showed statistically significant differences between E. affinis godmani and E. affinis affinis for these three characters, whereas females differed significantly in only two characters (WC and BD, Fig. 3). For the variables TL (Tail Length), BL (Bill Length), and BW (Bill width) we found significant differences between both groups (Fig. 3, Table 3). PCA analyses for males showed an 83.5% proportion of variance explained for two principal components, while females showed an 80.2% proportion of explained variance explained. Both Component plots showed that E. affinis godmani was distributed in quadrant 2, while E. affinis affinis had a wide distribution. The dispersion plot showed a partial overlap between E. affinis affinis and E. affinis godmani, but wing chord and bill depth were clearly differentiated in boxplots. These were also the two most important variables included in the first principal component in both sexes, according to the eigenvectors (Fig. 3).

PCA analyses for both sexes showed that the first two principal components explained a large proportion of the variance (70.37%, Table 3). Contrary to the previous graphs, both sex component plots show a total overlap between the two lineages.

Figure 3. 

Morphometric analyses results. A) Females boxplots and PCA for WC, TLE, and BD morphometric characters. B) Males boxplot and PCA for WC, TLE, and BD morphometric characters. C) Boxplot and PCA for TL, BL, and BW. WC, TLE, and BD characters were analyzed by separated sex, because the analyses indicated sexual dimorphism (see results and Table 3). Bill length (BL, from the upper base of the bill to the tip of the upper mandible), bill width (BW), bill depth (BD, from the upper mandible to the base of the bill at the distal edge of the nostrils), wing chord (WC, distance from the carpal joint the tip of the longest primary), tarsus length (TL), and tail length (TLE, distance from the uropygial gland to the tip of the longest rectrix).

None of the frequency variables measured differed significantly between godmani and affinis groups (N =19). On the other hand, we found that the emission rate of E. a. affinis is much lower than that of E. a. godmani, on average 2.9 notes/s versus 5.09 notes/s. These differences are statistically significant (W = 0, p < 0.001). The first two Principal Components together explain 78.75% of variance. The first PC separates both groups unambiguously (Fig. 4), and has a highly positive correlation with call duration, as well as a highly negative correlation with emission rate (notes per second).

Figure 4. 

Vocalization analysis. Boxplot of note emission rate A and PCA of measured vocal characters B Calls differ between the two groups in temporal structure, but not in frequency or number of notes.

Our models obtained a high mean value for AUC (Area Under the Curve) ratio values and statistically significant, 1.68 for E. affinis affinis and 1.81 for E. affinis godmani (***P < 0.05), this indicates a good fit of ENM’s. According to the Jackknife test and contribution variables obtained by MaxEnt, the most important variable for E. affinis affinis model was BIO 15 (precipitation seasonality), and the variable BIO8 (mean temperature of wettest quarter) for E. affinis godmani. We present the ENM predictions in four levels in Fig. 4; the first ones are the present predictions for the two lineages and the overlap between them, in general, both models predicted the previously known area of distribution for both species (Fig. 4), with an overlap in the West Pacific coast. In the second level, we observed that E. affinis godmani has a limited ability to predict its ecological niche in the geographic areas where E. affinis affinis is distributed, while E. affinis affinis projected its ecological niche on a large geographic area of distribution for E. affinis godmani.

The third part is the projection of the models in Last Glacial Maximum conditions (LGM 21–18,000 years ago), for E. affinis affinis showing a reduction in their environmental suitability along the present distribution, with predictions in areas like the Yucatan Peninsula and the western coast of Mexico with a gap at the western coast of the Tehuantepec Isthmus, unlike Present predictions where Central America has only small patches with predictions for E. affinis affinis (Fig. 4). For E. affinis godmani we found predictions in the western coast of Mexico, with a large gap between the central Mexican Coast west and the western coast of the Tehuantepec Isthmus, it also has a small patch prediction in the Yucatan Peninsula. The fourth part is the Last Interglacial (~ 120,000–140,000 years ago), for both lineages, the areas with high environmental suitability increased with respect to LGM, for E. affinis affinis it including the western Yucatan Peninsula, the Central western Mexican coast, and the Tehuantepec Isthmus to the western coast of Central America. For E. affinis godmani, the prediction areas are the West Mexican coast and the western Yucatan Peninsula.

The results of ecological overlap for the environmental PCA exhibit a large niche of E. a. affinis, while E. a. godmani exhibits an ecological niche compaction. A total variance of 83.67% is explained for the three principal components, with 41.04% for PC1, 29.886% for PC2 and 12.733% for PC3 (Fig. 5). The D and I statistic observed values were close to zero 0.01174013 and 0.05643784 respectively. We can reject the niche equivalency because the D observed value does not fall within the density of 95% of the random simulated values, however we obtained a no significant p value (0.9901). In contrast, the niche similarity test between E. affinis affinis and E. affinis godmani shows that they are less similar than expected by chance (Fig. 5), since there is not significant climatic niche conservatism (p = 0.31683, p = 27723) between them. With this evidence we reject the niche conservatism hypothesis between E. a. affinis and E. a. godmani, so we can say that the niches are divergent (Warren et al. 2008, 2010; Broennimann et al. 2012).

Figure 5. 

Ecological niche modelling and its projection in the geographic areas for E. a. affinis (yellow) and E. a. godmani (blue). In all four panels (a-d), the contribution values of each environmental variable of ENM’s is illustrated in the left and the projection of the Ecological niche conditions in the geographic distribution area is shown in the maps. a Ecological Niche projected in the current geographic distribution area of E. affinis and E. a. godmani. b ENM’s projected into the geography for each lineage. c ENM of E. a. affinis and E. a. godmani projected in the Last Maximum Glacial ecological conditions. d ENM of E. a. affinis and E. a. godmani projected in the Last Inter Glacial ecological conditions.

Figure 6. 

Equivalence and similarity tests in environmental space for E. a. affinis and E. a. godmani. A PCA of Ecological niche for of E. affinis lineages and the variables contribution to the analyses. The gray gradient indicates the density of the occurrences of the lineages, and the dashed and solid line indicates the 50% and 100% of the environmental background B graphical results of the equivalency tests comparing the two lineages. For both tests (equivalence and similarity) we only presented values for the D metrics. For all graphs the D observed values of the overlap niche analyses are present with the black diamond. The p value is showing in each graphic, all of them not significant for these analyses C graphical results of the similarity test comparing the two lineages in both directions (E. a. affinis vs. E. a. godmani and vice versa), ns = Not significant, p > 0.05.

We found three lines of evidence in the molecular data to support the taxonomic split of the E. affinis complex at species level. The first is the reciprocal monophyly between E. affinis and E. godmani found in the multilocus analysis using mitochondrial and nuclear genes, where the samples of E. affinis olmecorum are included into E. affinis. These results agree with the proposals by Ridgway and Friedman (1901) and the taxonomic proposal of Navarro-Sigüenza and Peterson (2004). The second line of evidence is the genetic distance between the western and eastern groups, which is similar to other Euphonia species close to the E. affinis complex (Table 2). These genetic distance values are also similar to distances found in other bird complexes distributed in Mexico and Central America that have been recognized as distinct species (Puebla-Olivares et al. 2008, Arbeláez-Cortés and Navarro-Sigüenza 2013, Zamudio-Beltrán and Hernández-Baños 2015). The third line of evidence is the high index of genetic fixation FST (Table 2) and the haplogroups in the haplotype networks (Fig. 1). It is important to mention that the haplotype networks showed geographical correspondence, with E. godmani in western Mexico, and E. affinis in eastern Mexico and Central America. These results are consistent with studies of other bird species distributed in Mesoamerica (Smith et al. 2011, Ramírez-Barrera et al. 2018).

Our analysis revealed significant differences between E. a. godmani and E. a. affinis in six characters among lineages. Bill Depth and Wing Chord are bigger for E. a. godmani, while E. a. affinis has bigger dimensions on Tail Length. Even though the rest of characters have significant differences, in the PCA plots the Tail Length, Bill Length and Bill Width characters do not show dispersion between both lineages, so we can assign Wing Chord, Tail Lenth and Bill Depth as diagnostic characters for males, and Wing Chord and Bill Depth as diagnostic characters for females. These results are similar to Phaethornis mexicanus morphometric patterns, a species also distributed along the Atlantic and Pacific Slope, where the Pacific lineage also shows bigger dimensions vs. the Atlantic lineage (Arbeláez-Cortés and Navarro-Sigüenza, 2013).

Our results show that there are significant differences in the temporal characteristics of calls between E. godmani and E. affinis while we found that there is little divergence in spectral structure or frequency measurements. E. godmani emits call notes at a significantly faster rate than E. affinis. Many bird species are highly sensitive to temporal cues in recognizing conspecific vocalizations (Dooling and Prior, 2016), which suggests that while call structure and frequency in this complex has been conserved, variation in tempo could be an important cue in conspecific recognition.

Euphonia affinis and E. godmani represent two different lineages with no significant conservatism in their ecological niches (west vs. east). The env-PCA, also, showed a larger ENM for E affinis, respect to E. godmani, also the western lineage has a limited ability to predict its ENM in the geographic area of E. affinis. The western coast of Mexico is characterized by a highly contrasting dry season vs. a wet season over the year, this characteristic is unique with respect to the eastern tropical area, so E. godmani has become restricted to these conditions. These results are similar to other taxa with sister lineages distributed along the Pacific and Atlantic slopes in Mesoamerican (Hernández-Canchola and León-Paniagua, 2017). It is interesting that E. godmani shows a reduction in ecological niche, while E. affinis presents a broader ecological niche. That may suggest a scenario where E. godmani was able to invade the western area of Mexico, and, in the absence of ecological competition from other Euphonias, it adapted and specialized to the floristic resources, as well as to the temperature and precipitation conditions of the area. While E. affinis conserved a broader ecological niche, as reflected in its geographical distribution, allowed it to explore more regions and resources, even in the presence of different species of Euphoniinae.

Lineage divergence between E. godmani (western Mexico) and E. affinis (eastern Mexico and Central America) occurred ~ 2.6 Mya (1.5–4.0 Mya HPD 95%), a range between the Pliocene and Pleistocene epochs. During the Pliocene, the Sierra Madre Occidental and the Transmexican Volcanic Belt finished emerging, which made the Pacific Slope drier than the Atlantic slope, due to the hillside effect (Graham and Dilcher, 1995). Additionally, the drier conditions were favored by meteorological phenomena that made the Pacific coast warmer than the Atlantic coast in the northern hemisphere (Molnar and Cane 2007). These events were decisive for the conformation of the tropical deciduous forest that extended throughout the Pacific from western Mexico to western Panama. According to studies of paleontological and molecular evolution, botanical elements present in the dry forests today were already present in said area since the Miocene (Graham and Dilcher 1995, Becerra 2005 de-Nova et al. 2012) however, from the Middle Pliocene to the late Pliocene these elements were unified as a plant community, promoting the diversification of some botanical groups (Becerra et al. 2005, de-Nova et al. 2012). Also, significant isolated periods of dry forest have been attributed to diversification in the Pacific Slope area (Becerra 2005, de-Nova et al. 2012, Willis et al. 2015). As a consequence, this province is characterized by a pattern of a high number of endemic lineages and species (Zaldívar et al. 2004, García-Deras et al. 2007, Zarza et al. 2008, Ramírez-Barrera et al. 2018). We found two threads of evidence that support the relationship between divergence of lineages for E. affinis and the origin of dry forests. The first evidence is the age of 2.6 Mya when the West and East lineages diverged during the late Pliocene, which coincides with the establishment of dry forests in Western Mexico. The other evidence is the adaptation and restriction of the environmental niche of E. a. godmani to the environmental conditions of Western Mexico. Other biogeographic events of Mesoamerica that shaped the biota were the closure of the Isthmus of Panama during the late Pliocene and the orographic changes in the Atlantic slope by the last raise of Transmexican Volcanic Belt and the Sierra Madre Oriental. However, the Atlantic Slope shows a wide mosaic of environments and ecosystems (Graham and Dilcher 1995), in contrast to the dry forest-dominated West slope, which could explain the more extensive environmental niche of E. a. affinis.

In addition to the consequences of the orographic changes of the Pliocene, during the Late Pliocene, global and continuous cooling periods were frequent, and during the Pleistocene the climatic oscillations were defined by glacial and interglacial periods (Zachos et al. 2001). During the glaciations, the species inhabiting temperate zones expanded their distribution to lower altitudes (Moreno-Letelier et al. 2014), while the geographic distribution of tropical vegetation was reduced. Tropical forests were affected by periods of low humidity which favored the reduction of the distributional range of several species, thus probably promoting speciation in plant species (Gentry 1982) as well as in the fauna of these forests, including birds (Smith and Klicka 2010). The divergence between Mesoamerican lowlands species has been attributed to these climatic changes, for example, amphibians (Greenbaum et al. 2011), reptiles (Ruane et al. 2014), mammals (Castañeda-Rico et al. 2014), and birds (Arbeláez-Cortés and Navarro-Sigüenza 2013). This work shows that orographic and environmental changes promoted the divergence of two lineages within E. affinis, probably due to isolation events and environmental adaptations, which in turn could accentuate the present differences in morphological, genetic, behavioral, and ecological characteristics previously described.