Multilocus, phenotypic, behavioral, and ecological niche analyses provide evidence for two species within Euphonia affinis (Aves, Fringillidae)

Abstract The integration of genetic, morphological, behavioral, and ecological information in the analysis of species boundaries has increased, allowing integrative systematics that better reflect the evolutionary history of biological groups. In this context, the goal of this study was to recognize independent evolutionary lineages within Euphonia affinis at the genetic, morphological, and ecological levels. Three subspecies have been described: E. affinis godmani, distributed in the Pacific slope from southern Sonora to Guerrero; E. affinis affinis, from Oaxaca, Chiapas and the Yucatan Peninsula to Costa Rica; and E. affinis olmecorum from Tamaulipas and San Luis Potosi east to northern Chiapas (not recognized by some authors). A multilocus analysis was performed using mitochondrial and nuclear genes. These analyses suggest two genetic lineages: E. godmani and E. affinis, which diverged between 1.34 and 4.3 My, a period in which the ice ages and global cooling fragmented the tropical forests throughout the Neotropics. To analyze morphometric variations, six morphometric measurements were taken, and the Wilcoxon Test was applied to look for sexual dimorphism and differences between the lineages. Behavioral information was included, by performing vocalization analysis which showed significant differences in the temporal characteristics of calls. Finally, Ecological Niche Models were estimated with MaxEnt, and then compared using the method of Broennimann. These analyses showed that the lineage distributed in western Mexico (E. godmani) has a more restricted niche than the eastern lineage (E. affinis) and thus we rejected the hypotheses of niche equivalence and similarity. Based on the combined evidence from genetic, morphological, behavioral, and ecological data, it is concluded that E. affinis (with E. olmecorum as its synonym) and E. godmani represent two independent evolutionary lineages.


Introduction
The integration of genetic, morphological, ecological, and behavioral data in systematic studies provides information on the evolutionary history of species and their populations, allowing a better assessment of species limits (Cadena and Cuervo 2010, Köhler et al. 2010, Pavlova et al. 2014, as well as understanding the role of geographical and ecological factors on population divergence within species . De Queiroz (2007) proposed that different types of data (e.g., morphological, ethological, ecological, molecular, etc.) are needed to determine operationally whether the lineages under study are evolving separately, and thus can be considered to represent different species. Species differentiation is affected by the time elapsed since the speciation event, a problem that should be considered in species delimitation studies. Padial et al. (2010) explained two ways to approach this problem: one is integration by congruence, and the other is integration by accumulation. In the first case (integration by congruence), taxonomists will consider two lineages as different species when there are concordant patterns of divergence among several taxonomic characters which result a full lineage separation. Meanwhile the integration by accumulation framework implies that divergences in any number of attributes (taxonomic characters) can provide evidence for the existence of a species, and in this case it is important to distinguish the group of characters (or even a single character) that promotes divergence and is reflected in the separation of lineages.
Species limits on birds have been studied using different approaches, including the use of morphological characters (Navarro-Sigüenza et al. 2001), coloration (Frith andFrith 1983, Rathbun et al. 2015), genetic variation (Olsson et al. 2013), songs (Sosa-López et al. 2013), and ecological niche modeling (Ruiz-Sánchez et al. 2015). In general, the objectives of these studies have been to resolve boundaries within species complexes and to evaluate subspecies recognized by taxonomic authorities such as the AOU (American Ornithology Union) and the IOC World Bird List (International Ornithology Committee).
DNA sequences have been useful to complement morphological and geographical information. Phylogenetics, molecular clocks, diversification rates, genetic populations and coalescence analyses have documented that geological complexity, heterogeneity of the environment, and climatic oscillations may have influenced patterns of genetic diversity, demography and divergence within species (de-Nova et al. 2012, Gu et al. 2013, Smith et al. 2014, Rodríguez-Gómez and Ornelas 2015. On the other hand, ecological niche modeling has provided information that supports the results of genetic studies on species delimitation (Raxworthy et al. 2007, Ruiz-Sánchez et al. 2015 and has helped discern whether speciation has been mediated by niche conservatism or ecological niche divergence (Brook et al. 2006, Wiens 2004, Wiens and Graham 2005. While vocal displays can be important prezygotic barriers to interspecific mating (Catchpole and Slater 2008), and many avian lineages have been discovered and, in part, diagnosed as distinct on the basis of differences in vocalization Ranft 2003, Halley et al. 2017). However, with the application of the integrative taxonomy it has become possible to incorporate diverse information such as multilocus, morphological and ecological data to hypothesize species limits (Mckay et al. 2014, Minoli et al. 2014, Perez and Borges-Martins 2019, Venkatraman et al. 2018. This approach has been very useful even for cryptic species (Ramos et al. 2019).
Euphonia affinis is a member of the family Fringillidae, subfamily Euphoniinae (Zuccon et al. 2012). Euphonia affinis appears to be phylogenetically related to E. chlorotica, E. luteicapilla, E. finschi, E. plumbea, E. concinna, and E. trinitiis (Isler andIsler 1987, Imfeld et al. 2020). Like all Euphoniinae, its distribution is restricted to the Neotropics. Specifically, E. affinis is a resident of the tropical lowlands from Mexico to Costa Rica (AOU 2003). In Mexico, the species is distributed along both slopes, from Sonora in the west and San Luis Potosi in the east, south to Central America (Howell and Webb 1995: fig. 1). Two subspecies are currently recognized based on geographical and morphological descriptions (AOU 1998). Euphonia affinis godmani Brewster, 1889 is endemic to western Mexico and is distributed from Sonora to central Guerrero. Euphonia affinis affinis Lesson, 1842 is distributed from eastern San Luis Potosi, southeastern Tamaulipas, Veracruz, Puebla and north-southwestern Oaxaca and the Yucatan Peninsula to Honduras, and on the Pacific Coast of Central America from Nicaragua to northwestern Costa Rica. The morphological characteristics that distinguish both subspecies are the subcaudal covert feathers, which are white in E. affinis godmani and yellow in E. affinis affinis in both males and females (Fig. 1). Dickerman (1981) described a third subspecies, Euphonia affinis olmecorum, based on differences in female coloration, which is paler than females of E. affinis affinis. Euphonia affinis olmecorum is distributed along the Gulf Coast of Mexico, from southeastern Tamaulipas and eastern San Luis Potosi to northern Chiapas (Hilty 2018); however, some taxonomic authorities do not recognize this subspecies, treating it as part of E. affinis affinis (Clement 2011). Currently, there are no studies of intraspecific limits for the species described in Euphoniinae. However, there is morphological and biogeographic evidence that the number of species is underestimated, with several species having wide ranges of distribution and more than one morphotype. Also, Imfeld et al. (2020) found genetic divergence between two subspecies of E. xanthogaster of similar magnitude to that between recognized species. Taken together, these observations indicate that species limits studies in Euphoniinae are needed.
In the present work, we applied integrative taxonomy to identify the independent evolutionary lineages within Euphonia affinis using four types of characters, multilocus genetic data, morphometric data, behavioral, and environmental niches. Based on the allopatric distribution of subspecies E. a. affinis (Eastern Mexico and Central America) and E. a. godmani (West of Mexico), as well as in the distinctive character of subcaudal feathers, we expect to recognize at least two independent evolutionary lineages that can be proposed to elevate at the species level.  Our goals were to: 1) obtain a phylogenetic hypothesis for Euphonia affinis subspecies using multilocus genetic data. 2) Associate the genetic variation and divergence times with historical geographic processes and barriers. 3) Describe the pattern of morphometric, behavioral, and environmental variation in Euphonia affinis, and associate it with genetic variation and phylogenetic relationships. Hence, our hypothesis is that multiple independent evolutionary lineages exist within the Euphonia affinis complex, and our objective is to define them with the integration of multilocus genetic data, morphometric, behavioral, and environmental data. Furthermore, we discuss a potential promotion of those lineages to species status.

Taxon sampling and sequencing procedures
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 S2).

Phylogenetic analysis
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  and discarded the first 25% of generations as burn-in.

Genetic diversity and structure
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.

Tests of divergence times
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 . 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).

Morphometrics
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.

Vocalization
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-affiniscomplex.

Ecological niche modeling and paleodistribution
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), 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. 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/ni-chetoolb2/), 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.

Genetic diversity and phylogenetic analyses
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.

Divergence times
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.

Vocalization
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).

Distribution modeling, paleodistribution, and ecological niche overlap
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(Warren et al. , 2010Broennimann et al. 2012).  (Lesson, 1842)  Females. Forehead olive-yellow and gray, olive-green back. The throat is olive-yellow, with a yellow belly, subcaudal coverts feathers also in yellow (Hilty 2018 (Hilty 2018).

Discussion
We provide molecular, morphological, behavioral, and environmental niche evidence supporting the existence of two evolutionary lineages within the Euphonia affinis complex (E. godmani and E. affinis). De Queiroz (2007) proposed that different types of data (morphological, ethological, ecological, molecular, etc.) can support species delimitation, if the lineages under study are evolving separately and can then be considered different species. Therefore, in our study we show these lineages have disjunct distributions, they have different plumages (white undertail coverts in E. godmani, and yellow undertail coverts in E. affinis, Hilty 2018), and we found significant differences in morphometric measurements. We did not find evidence of a third lineage corresponding to the subspecies E. affinis olmecorum proposed by Dickerman (1981); rather, it represents intraspecific variation within E. affinis.

Phylogenetics and genetic variation
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 , Ramírez-Barrera et al. 2018.

Morphometrics
Our analysis revealed significant differences between E.

Vocalization
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.

Ecological niche similarity
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.

Biogeographical history
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 2005de-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.

Conclusions
We incorporated different kinds of information to help us identify lineages within the Euphonia affinis species complex and understand the speciation process (De Queiroz 2005, 2007. We have demonstrated a sharp genetic split between E. a. affinis and E. a. godmani and we found a similar pattern in morphometrics, vocalizations, as well as in ecological niche data. So, we can conclude that our data support the consideration of E. affinis and E. godmani as two species. Zoology at Occidental College (MLZ), Museum of Natural Science (LSU), Museum of Comparative Zoology at Harvard University (MCZ), Academy of Natural Sciences of Drexel University (ANSP), the Delaware Museum of Natural History (DMNH) and Colección Nacional de Aves del Instituto de Biología, UNAM (CNAIB), for the facilities to access of skin specimens. We thank xeno-canto.org for making vocalization recordings available to the public. Alejandro Gordillo, Fabiola Ramírez, Raúl Ivan Martinez and Isabel Vargas for technical help, and to all the collectors at MZFC. We also appreciate the comments of David Prieto Torres, which were very important for the methods and the interpretation of results in ENM. Part of this research was done to obtain the masters degree of Alma Melisa Vázquez López, who was supported with a masters scholarship grant by CONACyT. The research was supported by PAPIIT/ DGAPA, UNAM through a grant to BEHB (IN204017).