Study Species

Cinnamon Hummingbird is distributed from Sinaloa south along the Pacific slope (including the Tres Marías Islands) to southern Chiapas, and along the coast of the Yucatán Peninsula in Mexico to Belize, with additional Central American populations extending to NW Costa Rica (Friedmann et al., 1950; Howell & Webb, 1995). Although sexes are similar (plumage), males have the bright red bill with black tip and females have it with mostly black above (Howell & Webb, 1995). In Mexico, the Cinnamon Hummingbird occupies the tropical dry forests, thorn forests and gallery forests from sea level up to 1600 m of altitude (Arizmendi et al., 2020; Howell & Webb, 1995). In the TDF of western Mexico and the Motagua Valley in Guatemala, Cinnamon Hummingbird is a resident and most abundant territorial species, and the core hummingbird species of the hummingbird-plant species network (Arizmendi & Ornelas, 1990; Bustamante Castillo et al., 2018, 2020; Díaz Infante et al., 2020). The clade composed of Amazilia rutila and sister species A. yucatanensis Cabot, 1845 and A. tzacatl De la Llave, 1833 is retrieved as monophyletic in molecular phylogenies (McGuire et al., 2014; Ornelas et al., 2014) and split from South American species c. 14–11 million years ago (Ornelas et al., 2014). The split rutila / tzacatl-yucatanensis is estimated to have occurred c. 12 million years ago, and the age of the A. rutila clade is estimated at 8 million years ago while the split A. tzacatl / A. yucatanensis occurring at c. 7–5 million years ago (Ornelas et al., 2014). These estimates coincide with high diversification rates of Bursera trees (Burseraceae, torchwood, copal) on the Pacific Coast and Balsas River Basin around 7.5 million years ago and the lower diversification in the Atlantic Coast and recent entry into Central American lowlands (Becerra, 2005), where A. yucatanensis and A. tzacatl are currently distributed.

Four subspecies are recognized on the basis of geographic variation in size and plumage colour (Dickinson & Remsen, 2013; Gill et al., 2020; Friedmann et al., 1950; Schuchmann, 1999; Weller, 1999): (1) A. r. diluta Van Rossem, 1938 from W Mexico (Sinaloa and Nayarit), sometimes considered a synonym of A. r. rutila (Friedmann et al., 1950; Van Rossem, 1938); (2) A. r. rutila (DeLattre, 1843) on the Pacific slope of Mexico, from Jalisco to Oaxaca, and south to Chiapas and the Yucatán Peninsula south to the arid valleys of interior Guatemala to NW Costa Rica; (3) A. r. graysoni Lawrence, 1866 off W Mexico on Tres Marías Islands; proposed to be elevated to species Amazilia graysoni (Gómez de Silva et al., 2020; Vázquez-López et al., 2021); and (4) A. r. corallirostris (Bourcier & Mulsant, 1846) SE Mexico on the Pacific slope in southern Chiapas, Mexico, south through western Guatemala to the Lempa River, El Salvador (Friedmann et al., 1950); or includes the populations from the Yucatán Peninsula south to Costa Rica, usually attributed to nominate rutila (Weller, 1999). Amazilia rutila graysoni of Tres Marías Islands are darker (more bronzy above and duskier below) and larger than mainland birds (Howell & Webb, 1995; Vázquez-López et al., 2021), A. r. diluta is similar to A. r. rutila but coloration below is paler and more pinkish cinnamon and upperparts slightly more golden bronze (less greenish), whereas A. r. corallirostris birds have darker rufous coloration underparts with chin more conspicuously spotted white (Howell & Webb, 1995; Schuchmann, 1999; Van Rossem, 1938). Populations on the Yucatán Peninsula are allopatrically distributed from those of A. r. rutila on the Pacific slope and from those of A. r. corallirostris from southern Chiapas, Mexico through western Guatemala to El Salvador and Costa Rica (Howell & Webb, 1995; Weller, 1999). However, there is no agreement with regard the phenotypic differences distinguishing these subspecies and their distribution limits and taxonomic status (Dickinson & Remsen, 2013; Friedmann et al., 1950; Gill et al., 2020; Schuchmann, 1999; Weller, 1999).

Sample Collection and DNA Sequencing

For molecular analysis, we sampled 71 individuals from 20 sites across the distribution of A. rutila in Mexico and Guatemala between almost sea level and 1095 m of altitude (Table 1,  Supplementary Figure S1). Hummingbirds were captured using mist nets, and two rectrices were collected from each hummingbird as a source of DNA for subsequent genetic analysis before the bird was released. Samples were collected under the required permits and using approved animal welfare protocols. According to current known range of subspecies, we sampled four individuals of A. r. diluta, 30 individuals of A. r. rutila, and 37 of A. r. corallirostris (Table 1,  Supplementary Figure S1). The sampling presented in this study practically covers the continental distribution limits of A. rutila in northern Mesoamerica; however, large extensions along the Pacific slope and Central America remain unsampled and the limits and precise distribution of subspecies are disputed (Weller, 1999).

Table 1.

Geographic Information and Sample Sizes of the 20 Amazilia rutila Localities Sampled in this Study.

DNA was extracted from tail feathers with the DNeasy Blood and Tissue extraction kit (Qiagen, Valencia, CA, USA), following the protocol recommended by the manufacturer. We amplified and sequenced three gene regions: NADH nicotinamide dehydrogenase subunit 2 (ND2, 398 bp), and the mitochondrial adenosine triphosphatase synthase 6 and 8 genes (ATPase 6 and ATPase 8, hereafter ATPase, 682 bp). Amplification of ND2 was conducted with primers L5215–H5578 (Hackett, 1996), whereas for the ATPase we used primers L8929 (Sorenson et al., 1999) and H9855 (Eberhard & Bermingham, 2004). Protocols for DNA extraction, PCR, and for sequencing the PCR products are described elsewhere (Licona-Vera & Ornelas, 2014). All newly acquired sequences have been submitted to GenBank (Accession nos. ND2: OP837824–OP837894; ATPase 6–8: OP837895–OP837953).

Relationships among Haplotypes

Statistical parsimony networks for the single (ND2, ATPase 6–8) and combined (ND2+ATPase 6–8) mtDNA datasets were constructed to infer relationships among haplotypes using TCS v1.2.1 (Clement et al., 2000), with the 95 % probability connection limit and gaps treated as single evolutionary events, and visualized using POPART (Leigh & Bryant, 2015). Loops were resolved following the criteria given by Pfenninger & Posada (2002). The aligned sequences of ND2 and ATPase for haplotype networks (Figure 1(a) and (b)) were 398 (71 sequences) and 682 (59 sequences) base pairs (bp) long, respectively. The combined ND2+ATPase sequences (59 sequences) were 1080 bp. An additional statistical parsimony network was constructed as described using a composite ND2 dataset that included our newly generated sequences and those in Vázquez-López et al. (2021) downloaded from the GenBank (Accession nos. MZ998668–MZ998740). All populations are newly sampled in this study, particularly those from most of the locations from Jalisco (locations 3–8), Santa Efigenia, Oaxaca (location 12), and Esquipulas, Guatemala (location 14) were not sampled in the previous work ( Supplementary Table S1; see also  Table S1 in Vázquez-López et al., 2021). However, the resulting matrix (143 sequences, 1079 bp after alignment) has a lot of missing data after alignment. Our ND2 sequences had a length of 398 bp and those in Vázquez-López et al. study had 605 to 1041 bp. Haplotype networks can be misleading in the presence of missing data because produce biased network relationships (collapsing of different sequences into a single haplotype) and/or fail to indicate alternative positions for sequences (Joly et al., 2007). Thus, we cut the sequences at the same length (235 bp) to avoid bias. The composite ND2 dataset (143 sequences) retrieved a single network and yielded 15 haplotypes ( Supplementary Figure S2); however, we observed that in some sequences from the Vázquez-López et al. study insertions were presumably introduced during alignment and some sequences with missing data were coded incorrectly. These problems can affect results of analyses that rely on network topology, such as haplotype diversity and population differentiation indices or past population history inferences, and could bias estimates of population structure and migration rates (Joly et al., 2007). For these reasons, we do not use the composite ND2 dataset in further analysis.

Figure 1.

Statistical parsimony networks of (a) single ND2 and (b) ATPase data and (c) combined ND2+ATPase data of Amazilia rutila overlaid on a map of North America. Numbers correspond to sampling localities. Refer to Table 1 and  Supplementary Table S1 (Supporting information) for sampling locality information and haplotype distribution. Haplotypes are coded with a different colour according to geographic regions, Pacific slope west of the Isthmus of Tehuantepec (PAC, yellow), Pacific slope of Chiapas and Oaxaca east of the Isthmus of Tehuantepec (dark red), and Yucatán Peninsula (YUC, turquoise), and a number codes each of the haplotypes. The size of sections of the pie charts corresponds to the number of individuals with that haplotype. In the inset (d), a Bayesian analysis of population genetic structure (BAPS) is presented based on the ND2+ATPase sequences. Colours indicate different genetic clusters (K = 3), as explained above.

Then, to estimate the most likely number of genetically differentiated clusters without making a priori assumptions about the partitioning of genetic diversity, we performed a nonspatial genetic mixture analysis with a Bayesian model-based approach employing the program BAPS v5.3 (Corander et al., 2008) and the combined ND2+ATPase 6–8 dataset. BAPS assesses the most likely number of genetically different clusters using the module for linked molecular data. The codon linkage model, appropriate for our sequence data, was applied. We surveyed the probability of a different number of genetic clusters (PP) under the independent loci model in two independent runs with the number of proposed clusters (K) ranging from 2 to 15, with 10 runs for each K.

Population Indices and Phylogeographic Structuring

To describe intraspecific genetic variation of A. rutila, molecular diversity indices (h, gene diversity; π, nucleotide diversity) and pairwise comparisons of F _ST values between genetic groups were calculated using ARLEQUIN v3.5 (Excoffier & Lischer, 2010) with 20,000 permutations. Haplotype diversity indices for each population (h _S, v _S) and at the species level (h _T, v _T), and coefficients of population differentiation (G _ST, N _ST) were estimated using PERMUT v2.0 (Pons & Petit, 1996). We further compared the G _ST and N _ST values and tested for phylogeographic structure using PERMUT with 10,000 permutations and the U-statistic. An N _ST value significantly higher than the G _ST value provides evidence of phylogeographic structure (Pons & Petit, 1996). Note that ‘populations’ are sampling localities, whereas ‘groups’ are sets of pooled populations, as specified in Table 1. Locations with one or two samples were lumped with closest location (seven samples from locations 2, 4 and 7 into population 2, four samples from locations 3, 5 and 8 into population 3, and five samples from locations 15, 16 and 17 into population 11).

Geographic Structure of Populations

To test for the presence of hierarchical population structure, analyses of molecular variance (AMOVA) were run in ARLEQUIN v3.5. Populations were treated (a) as one group to determine how much variation is explained by differences between populations sampled, and grouped into (b) two groups of populations separated by the Isthmus of Tehuantepec, from Sinaloa and Durango to Oaxaca on the Pacific slope of Mexico (west of the Isthmus of Tehuantepec; locations 1–11, Table 1) and from southern Chiapas, Mexico through western Guatemala to El Salvador and the Yucatán Peninsula (east of the Isthmus of Tehuantepec; locations 12–20, Table 1), or (c) three groups (PAC = locations 1–11, CHIS_OAX = locations 12–13, YUC = locations 14–20; Table 1). Significance of AMOVAs was determined with 20,000 permutations each. These groupings allow us to test several models for the presence of hierarchical population structure and genetic variability found among geographic regions.

We also infer the optimal number of geographically homogeneous and maximally differentiated groups (K) using a spatial analysis of molecular variance (SAMOVA) implemented in SAMOVA v1.0 (Dupanloup et al., 2002) and testing values of K from 2 to 5, with 10 replicates per each K value.

Historical Demography

Signatures of demographic expansion in A. rutila were addressed by means of neutrality tests, Fu’s Fs (Fu, 1997), Tajima’s D (Tajima, 1989) and Ramos-Onsins and Rozas’ R2 (Ramos-Onsins & Rozas, 2002) statistics of neutrality. R2 is the most powerful tests used to detect population growth (Ramos-Onsins & Rozas, 2002). Significance was evaluated by comparing observed values with null distributions generated by 10,000 replicates, using the empirical population sample size and the observed number of segregating sites in the pegas package of R v4.1.1 (Paradis, 2010; R Core Team, 2020). Also, we compared the distributions of pairwise nucleotide differences among haplotypes (mismatch distribution) to expectations of a sudden-expansion model (Rogers, 1995) for the genetic groups. Mismatch distributions were generated with ARLEQUIN v3.5 using the sum of squared deviations test (SSD) and Harpending’s raggedness index (Harpending, 1994), both of which are higher in stable, nonexpanding populations (Rogers & Harpending, 1992). Deviations from the null model of population expansion were tested with 1,000 parametric bootstrapping (Schneider & Excoffier, 1999) in ARLEQUIN. Significant negative values of D and Fs and small positive values of R2 indicate an excess of low frequency mutations relative to expectations under the standard neutral model (i.e. strict selective neutrality of variants, constant population size, and lack of subdivision and gene flow). The mismatch distribution analysis was carried out using the sudden demographic expansion model of Schneider & Excoffier (1999) in unsubdivided populations and the spatial expansion model in a subdivided population (Excoffier, 2004). We employed 20,000 replicates to test the goodness-of-fit of the observed mismatch distribution to that expected under the spatial and sudden demographic expansion model using the sum of squares differences (SSD) and Harpending’s raggedness index (Hri index; Harpending, 1994) according to Rogers & Harpending (1992).

The time at which the spatial expansion event took place was dated following the expression, t = τ/2 µk, where τ is the estimated number of generations since the expansion, µ is the mutation rate per site per generation, and k is the sequence length (Schneider & Excoffier, 1999). The expansion parameter tau (τ) was estimated using ARLEQUIN in genetic lineages in which signs of sudden demographic expansion were evident. To convert the time since expansion (t), we used a 2.1-years generation time based on the observation that the age of maturity begins 1 year after hatching, and an assumed low annual adult survival rate of 0.52 reported for Bassilina (Hylocharis) leucotis (Ruiz-Gutiérrez et al., 2012). The approximate average generation time (T) is calculated according to T = a + [s / (1–s)] (Lande et al., 2003), where a is the time to maturity and s is the adult annual survival rate. Based on this, the estimate for T was 2.1 years.

Lastly, changes in effective population size (N _e) through time were estimated using Bayesian skyline plots (BSP) in BEAST v2.4.2 (Bouckaert et al., 2014). We chose the HKY+I (for A. rutila sensu lato) and HKY (for PAC and YUC groups) substitution models with empirical base frequencies, a strict clock model, and a piecewise-linear coalescent Bayesian skyline tree prior with five starting groups. One independent run of 30 million generations each was run, with trees and parameters sampled every 1000-iterations, with a burn-in of 10 %. Results of each run were visualized using TRACER to ensure that stationarity and convergence had been reached (ESS > 200). The time axis was scaled using the mtDNA geometric mean substitution rate of 0.001214 s/s/l/Myr (Lerner et al., 2011). The BSP for the CHIS_OAX group was not estimated due to small sample size.

Ecological Niche Modeling (ENM)

We constructed an ENM model in MAXENT v3.3.3k (Phillips et al., 2006) to predict the current distribution of suitable habitat occupied by Amazilia rutila. Coordinates of occurrence data obtained from the Global Biodiversity Information Facility ( GBIF.org (13 October 2021) GBIF Occurrence Download,  https://doi.org/10.15468/dl.ps6zn8) were assembled and supplemented with our geo-referenced records from field collection. We debugged these data by selecting the records from museum collections and iNaturalist research-grade observations but eliminating those outside the known distribution of the species. After careful verification of every data location, we excluded duplicate occurrence records or in close proximity to each other (c. 1 km²) to reduce the effects of spatial autocorrelation using GridSample package in R (Thomson et al., 2017). The dataset was restricted to 312 unique presence records for the analysis.

We used 19 climatic variables summarizing data of precipitation and temperature (BIO1–BIO19 variables) as climate layers from WorldClim v1.4 (Booth et al., 2014; Hijmans et al., 2005) at c. 1 km² (30 arc-sec) spatial resolution to predict the current distribution of suitable habitat occupied by A. rutila. An important step in ENM is to delineate a realistic calibration region (‘M’, BAM diagram; Soberón & Peterson, 2005); that is, the set of sites accessible to a species over which models are calibrated (Atauchi et al., 2020; Barve et al., 2011; Freeman et al., 2019; Soberón & Peterson, 2005). In this study, we calibrated the ENM model for A. rutila (in practice a mask or GIS polygon) using a geographical clipping to the species range based on the ecoregions of Mexico and Central America proposed by Olson et al. (2001) as the ‘M’ accessible areas considering potential boundaries on the landscape to dispersal and altitude range limits ( Supplementary Table S2 and  Figure S3; Barve et al., 2011). Climate layers were clipped with the ‘M’ region for use in MAXENT analysis. We used the variance inflation factor (VIF) in the usdm package (Naimi, 2015) from the R v4.1.1 (R Core Team, 2020) to exclude highly correlated variables and avoid multicollinearity. We retained bioclimatic variables with VIF values lower than 10 for all species records in each run, until all remaining VIF values were less than 10, as VIF values higher than 10 indicate strong collinearity (e.g., Ranjitkar et al., 2016). After removing the highly correlated variables, eight variables were used in the final analysis (BIO2 = Mean Diurnal Range, BIO3 = Isothermality, BIO8 = Mean Temperature of Wettest Quarter, BIO13 = Precipitation of Wettest Month, BIO14 = Precipitation of Driest Month, BIO15 = Precipitation Seasonality, BIO18 = Precipitation of Warmest Quarter, and BIO19 = Precipitation of Coldest Quarter). Final model was constructed in MAXENT with ten cross-validation replicates without extrapolation and considering the average output grids as the final predictive model. We applied a binary transformation (absence or presence, zero or one) using the 10th percentile training presence logistic threshold (T10LT). The geographic representations of the climatically suitable areas were constructed using a continuous representation of environmental suitability values, which were spatially projected using ArcMap. The area under the receiver operating characteristic (AUC) curve was used to evaluate the prediction performance of the model, in which values around 0.5 represent distribution models no better than random and those around 1 represent a perfect fit between the observed and the predicted species distribution; acceptable models are those with > 0.7 AUC values (Phillips et al., 2006). However, several criticisms have been associated with this approach (e.g., Cobos et al., 2019; Lobo et al., 2008; Merow et al., 2014). For this reason, we also evaluated model using the partial-ROC test (Peterson et al., 2008). Within a value range from 0 to 2, values over 1 suggest a better performance than chance, by analysing the presences versus the absence against the total area predicted by MaxEnt (Osorio-Olvera et al., 2020). Lastly, we also used the true skill statistic (TSS), which is a threshold-dependent measure of model performance, to evaluate the accuracy of predictive maps generated by presence-only data (Allouche et al., 2006; Liu et al., 2013), where TSS values ranging between 0.4 and 0.8 are considered useful (Fielding & Bell, 1997; Landis & Koch, 1977). For each replicate, TSS was calculated using the T10LT and then TSS values were averaged among replicates using the sp package (Pebesma & Bivand, 2005) in R v4.1.1.

Niche Divergence

We quantified the differentiation (or overlap) between climatic niches of the Amazilia rutila genetic groups (PAC, CHIS_OAX, YUC) based on their ENMs. The databases were worked individually to avoid identification errors at the limits of the distribution areas. Climate niche overlaps among groups were estimated using the PCA-env method proposed by Broennimann et al. (2012) with variables selected in the ENM by assessing their variance inflation factors (VIF) and then eliminating the multicollinear predictors. After that, we reduced the 19 variables to the same eight variables used in our ENM. Principal component analysis (PCA) was used to transform the environmental space of the investigated or selected environmental variables into a two-dimensional space defined by the first and second principal components (Strubbe et al., 2015). The PCA-env is carried out to transform the climate layers into a reduced number of linearly uncorrelated variables, i.e. principal components (Broennimann et al., 2012). The first component represents the largest possible amount of variability in the original variables, and each subsequent component represents the largest part of the remaining variability. This test compares the available environmental conditions for each species within a defined study extension (background) with their observed occurrences and calculates the available environmental space defined by the first two PCA axes. The differences in the position of the species along the principal components reflect their environmental differences.

Subsequently, the overlapping of the niches by pairs of the groups (PAC vs. YUC, PAC vs. CHIS_OAX, YUC vs. CHIS_OAX) was calculated using the Schoener’s D metric (Schoener, 1970). The values ​​of this metric range from 0 (meaning that the niches are completely different) to 1 (meaning that the niches completely overlap) (Broennimann et al., 2012) and the graphs were made to observe the surfaces density of occurrences for each group. To facilitate the interpretation of results, the outputs were condensed into five classes as suggested by Rödder & Engler (2011): 0–0.2=no or very limited overlap, 0.2–0.4=low overlap, 0.4–0.6=moderate overlap, 0.6–0.8=high overlap and 0.8–1.0=very high overlap. Finally, two different randomization tests were used to test the niche evolution hypotheses. An equivalency test, which determines whether the niches of two entities in two geographic ranges are equivalent (that is, whether the niche overlap is constant by randomly reallocating the occurrences of both entities between the two ranges), and a similarity test which compares the niche overlap of one randomly distributed range on its background while keeping the other unchanged, and then performs reciprocal comparison. Each randomization process is repeated 100 times (to ensure that the null hypothesis can be rejected with a high level of confidence), producing a null distribution of overlapping values ​​against which the observed score was compared. If the observed value of D is located outside 95 % of the density of the simulated values, the null hypothesis is rejected (H0= the niches are similar or equivalent), which implies that the groups occupy different segments of the environmental space. The stochastically simulated values were generated using the ENM of each Cinnamon Hummingbird genetic group and an ENM created with random points drawn from the minimum convex polygon surrounding the original occurrence records of the other genetic groups of Cinnamon Hummingbird. The geographic ranges of the genetic groups were used as backgrounds individually. For the background test, we used the minimum convex polygon surrounding the original occurrence records and the entire Mexico and Guatemala as the species range. All analyses were computed with the ecospat package (Di Cola et al., 2017) in R.

Haplotype network

The combined ND2+ATPase of A. rutila from 20 localities yielded 23 haplotypes (Figure 1(c), Table 1,  Supplementary Table S3). The statistical parsimony analysis retrieved a single network, in which three haplogroups (mtDNA groups) were revealed: populations 1–11 west of the Isthmus of Tehuantepec (PAC, A. r. diluta+A. r. rutila), populations 12–13 from Oaxaca and Chiapas east of the Isthmus of Tehuantepec (CHIS–OAX, A. r. rutila or A. r. corallirostris), and populations 14–20 from the Yucatán Peninsula and Guatemala (YUC, A. r. rutila or A. r. corallirostris) (Figure 1(c), Table 1,  Supplementary Table S3).

The most widespread haplotype was H2, which forms the core of the first haplogroup composed of ten haplotypes (YUC), eight haplotypes exclusively found in populations from the Yucatán Peninsula and two haplotypes exclusively found in Guatemala (Figure 1(c),  Supplementary Table S2). The second most frequent haplotype (H1) recovered for the second haplogroup composed of eight haplotypes was distributed exclusively in PAC populations west of the Isthmus of Tehuantepec (Figure 1(c),  Supplementary Table S3). Forming the third haplogroup, haplotypes H7, H10 and H21–H23 were shared among individuals from populations east of the Isthmus of Tehuantepec, from Chiapas and Oaxaca (CHIS_OAX). No haplotypes separated by several mutational steps were shared between populations of the three-mtDNA groups (Figure 1(c),  Supplementary Table S3).

BAPS analysis with mtDNA sequences and spatial clustering of groups of individuals resulted in three congruent clusters (K = 3) as the best partition (log marginal likelihood = –686.98, PP = 1.0; Figure 1(d)). In agreement with the haplotype network, populations west of the IT fell into a single cluster (A. r. diluta + A. r. rutila). Populations east of the IT from Oaxaca and Chiapas clustered together (A. r. rutila or A. r. corallirostris), and populations from the Yucatán Peninsula and Guatemala (A. r. rutila or A. r. corallirostris) formed the third cluster in the BAPS plot (Figure 1(d)).

Population Indices and Phylogeographic Structuring

Number of haplotypes varied among groups, from five in CHIS_OAX with six samples to ten in YUC with 31 samples (Table 2). Gene diversity and nucleotide diversity values were highest for CHIS_OAX (0.87, 0.002), respectively, followed by those for PAC (0.73, 0.001) and YUC (0.66, 0.001) (Table 2).

Table 2.

Results of Neutrality Tests and Mismatch Distribution of Amazilia rutila Samples (ND2+ATPase) by Genetic Group to Infer Demographic Expansion.

Observed genetic differentiation among populations based on ND2+ATPase (G _ST = 0.437, SE = 0.0838) indicated that A. rutila is genetically subdivided. Genetic diversity across all populations (h _T = 0.993, SE = 0.0407; v _T = 0.950, SE = 0.0490) was higher than the average within-population value (h _S = 0.525, SE = 0.0947; v _S = 0.256, SE = 0.1093). PERMUT analysis showed that N _ST (0.730, SE = 0.1177) and G _ST values were statistically different (p < 0.05), indicating phylogeographical structuring. Pairwise comparisons of F _ST values were high and significant when groups of populations were compared (PAC vs. YUC, F _ST = 0.9777, p < 0.001; PAC vs. CHIS_OAX, F _ST = 0.9616, p < 0.001; YUC vs. CHIS_OAX, F _ST = 0.9779, p < 0.001).

Geographic Structure of Populations

The AMOVAs showed significant genetic differentiation at each hierarchical level (Table 3). When groups were not defined, the AMOVA showed that the highest percentage of variation (97%) was explained by differences among populations and only 2.2 % by differences within populations (Table 3). When grouping populations as separated by the Isthmus of Tehuantepec, the highest percentage of genetic variance (73.4 %) was explained by differences between groups, 1.5 % by differences within populations, and 25.3 % of the variance accounted for differences among populations within groups. The F _CT value was high and significant (0.73, p < 0.01), indicating genetic differentiation between populations separated by the Isthmus. When AMOVA was performed for three groups (PAC, CHIS_OAX, YUC), the genetic differentiation was higher and significant (F _CT = 0.97, p < 0.001, Table 3).

Table 3.

Results of AMOVA and SAMOVA Models on Amazilia rutila Populations for Mitochondrial DNA (ND2 and ATPase).

The SAMOVA detected strong geographical structure with same three groups of populations inferred as the optimal number of geographical clusters (Table 3).

Historical Demography

For ND2+ATPase, Tajima’s D and Fu’s Fs values were negative and significantly different from zero, except Tajima’s D in PAC and CHIS_OAX (Table 2). Also, the low and non-significant values of SSD and Hri indicate a good fit to the demographic expansion model and consistent with a scenario of a sudden demographic expansion, so the expansion model was not rejected (Table 2). R2 statistic showed positive, small and highly significant values for PAC and YUC groups, indicating that these groups presented past demographic expansion (Table 2). Based on our estimated values of τ, the average time since the demographic expansion was 52.37 ka BP for PAC, 84.45 ka BP for CHIS_OAX, and 39.40 ka BP for YUC (Table 2).

The Bayesian skyline plots suggest that the effective population size was stable over time in YUC and PAC, except a marginal decrease in Amazilia rutila as a whole (100,000–50,000 years ago; Figure 2).

Figure 2.

Bayesian skyline plots showing historical demographic trends for YUC and PAC genetic groups and for Amazilia rutila as a whole using mitochondrial sequences (ND2 and ATPase data).

Ecological Niche Modelling

The current distribution of A. rutila was supported by high predictive power (AUC, mean ± SD, 0.796 ± 0.33) indicating adequate model performance. The partial ROC test (1.651 ± 0.03) showed that models were statistically significant (p < 0.01). Thus, performance values for the model assessment approach indicated that the distribution model was statistically accurate. Determination of the threshold probability for predicted presence using TSS resulted in a mean proportion of correctly classified training observations (TSS, mean ± SD, 0.476 ± 0.056). The projections of estimated current distribution suggest that areas of suitable habitat for A. rutila are continuous along the Pacific slope and that the predicted distribution on the Yucatán Peninsula is separated from the distribution predicted along the Pacific slope (Figure 3(a)).

Figure 3.

(a) Predicted distribution of Amazilia rutila at present. The output of MAXENT consists of grid maps with each cell having an index of suitability indicated between 0 and 1. Biogeographic barriers Isthmus of Tehuantepec (IT) and Motagua-Polochic-Jocotán (MPJ) fault system shown as dashed lines. The approximate geographic range separation between A. r. diluta and A. r. rutila subspecies is shown on the map by solid horizontal line. (b) PCA plot showing the first two components of PCA-env in A. rutila. (c) PAC, (d) CHIS_OAX and (e) YUC population niche displayed on the same multi-dimensional scale represented by the first two axes of a principal components analysis (PCA) summarizing the entire study area. Grey shadings show the density of the occurrences of the species by cell. The solid and dashed contour lines illustrate 100 % and 50 % of the available (background) environment, respectively. Abbreviations: PAC, populations west of the Isthmus of Tehuantepec (A. r. diluta + A. r. rutila); CHIS_OAX, populations from Oaxaca and Chiapas east of the Isthmus of Tehuantepec (Amazilia ‘saturata’ according to Vázquez-López et al., 2021); and YUC, populations from the Yucatán Peninsula and Guatemala (A. r. corallirostris).

Niche divergence

The PCA indicated three niche axes that together explain 64.38 % of the total environmental variation along climatic gradients in A. rutila (PC1 = 49.95 %, PC2 = 14.43 %; Figure 3(b)). PC1 is positively associated with mean diurnal range of temperature (BIO2) and precipitation seasonality (BIO15) and negatively associated with precipitation variables (BIO13, BIO14, BIO18, BIO19), and the second niche axis is positively associated with precipitation of warmest quarter (BIO18) and negatively associated with temperature variables (BIO3, BIO8) (Table 4). The occurrence density surfaces in environmental space, as determined by PCA-env, showed that the position in environmental space varied among genetic groups (Figure 3(c)–(e)). The contribution of the climatic variables on the two axes of the PCA-env and the percentage of inertia explained by the two axes is presented in Figure 4(a). Each genetic group differed in their position in environmental space. In general, the Schoener’s D niche overlap scores, which ranged from 0.17 to 0.39, indicate a low to very limited overlap in the fundamental climatic niche dimensions of all groups analyzed, particularly the very limited niche overlap between YUC and PAC (Table 5). The occupied niches by the genetic groups were not identical (p < 0.009); the null hypothesis of niche equivalency was rejected for all comparisons between genetic groups (Figure 4(b)–(d), Table 5). For niche similarity tests, however, all comparisons were rejected (Figure 4(e)–(j)), indicating that their environmental space is more similar to each other than expected by chance (Table 5).

Table 4.

Contributions of the First Two PCA-env Axes to Environmental Space.

Figure 4.

(a) The contribution of the climatic variables on the two axes of the PCA-env and the percentage of inertia explained by the two axes (PC1 and PC2). Histograms show simulated niche overlaps (grey bars) and the observed niche overlap between groups (bars with a circle) on which tests of niche equivalency (b–d) and niche similarity (e–j) were calculated from 100 iterations. The significance of the tests is shown. Abbreviations: PAC, populations west of the Isthmus of Tehuantepec (A. r. diluta + A. r. rutila); CHIS_OAX, populations from Oaxaca and Chiapas east of the Isthmus of Tehuantepec (Amazilia ‘saturata’ according to Vázquez-López et al., 2021); and YUC, populations from the Yucatán Peninsula and Guatemala (A. r. corallirostris).

Table 5.

Ecological Niche Comparisons for the Amazilia rutila Genetic Groups.

Genetic Differentiation among Amazilia rutila Populations

In this study, we elucidate the geographic structure of genetic variation of A. rutila populations in Mexico and Guatemala using mitochondrial DNA sequence data, and determine the effects of major geographic barriers, geographic distribution of suitable habitat and environmental variability on population divergence. The haplotype network, F _ST statistics, PERMUT, BAPS, AMOVA, and SAMOVA estimates revealed that individuals from the Yucatán Peninsula and Guatemala are genetically distinct from those genetically differentiated in the Pacific slope separated by the Isthmus of Tehuantepec (west of the Isthmus from Sinaloa to Oaxaca and east of the Isthmus in Oaxaca and Chiapas), supporting the hypothesis that geographical or ecological barriers have limited gene flow and promoted isolation between populations in the three geographic areas. In a recent study, Vázquez-López et al. (2021) evaluated the genetic differentiation of Amazilia rutila populations, recognizing the existence of three genetic groups: Mexican Pacific slope, Chiapas, and Yucatán Peninsula. In their phylogeny, individuals from the Tres Marías Islands were nested within the Mexican Pacific slope group. Despite this, they proposed elevating each group to species status, Amazilia rutila (Mexican Pacific slope), Amazilia graysoni (Tres Marías Islands), Amazilia saturata (Chiapas), and Amazilia corallirostris (Yucatán Peninsula and Central America). Using samples from a different set of individuals and the implementation of a new mitochondrial marker (ATPase 6–8), we recovered three genetic groups in continental populations and these results partially coincide with the groups that present Vázquez-López et al. (2021). In our study the CHIS_OAX genetic group (Chiapas lineage or A. saturata according to Vázquez-López et al., 2021) also includes populations of Oaxaca located east of the Isthmus of Tehuantepec (Figure 1(c)).

The geographical break between PAC and CHIS_OAX at the Isthmus of Tehuantepec favours the hypothesis of genetic divergence due to isolation between populations separated by this barrier (allopatric fragmentation). This break is spatially congruent with those of other hummingbird species (Cortés-Rodríguez et al., 2008a; González et al., 2011; Jiménez & Ornelas, 2016; Malpica & Ornelas, 2014; Ornelas et al., 2016; Rodríguez-Gómez et al., 2013; Rodríguez-Gómez & Ornelas, 2018), and those of other bird (e.g., Álvarez et al., 2016; Barrera-Guzman et al., 2012; Cortés-Rodríguez et al., 2013; Maldonado-Sánchez et al., 2016) and vertebrate taxa (e.g., León-Paniagua et al., 2007; Castoe et al., 2009; Mendoza et al., 2019; Mulcahy et al., 2006) inhabiting montane woodlands in Mesoamerica, including plants (e.g., Gutiérrez-Rodríguez et al., 2011; Martínez de León et al., 2022; Ornelas & Rodríguez-Gómez, 2015). However, these vicariant events seemly occurred at different times (Barber & Klicka, 2010; Ornelas et al., 2013, 2015), indicating that genetic divergence occurred variously during different temporal windows. In addition, phylogeographic studies have found low levels of genetic differentiation and gene flow between populations separated by the Isthmus, for both highland and lowland hummingbird species (Hernández-Soto et al., 2018; Rodríguez-Gómez et al., 2013, 2021, Zamudio-Beltrán et al., 2020a, 2020b) and other taxa (Arbeláez-Cortés et al., 2010; Cortés-Rodríguez et al., 2008b; Mendoza et al., 2019; Navarro-Sigüenza et al., 2008; Ornelas et al., 2010; Vázquez-Miranda et al., 2009), leading these authors to suggest that the Isthmus is a semipermeable barrier to gene flow. Divergence time estimates for the split between populations of A. rutila separated by the Isthmus (4.84–3.16 million years ago; Vázquez-López et al., 2021) support the hypothesis that the split between populations by the Isthmus pre-date the Pleistocene glacial periods. Although the existence of several mutational steps between populations separated by the Isthmus is indeed indicative of ancient allopatric fragmentation, it could also result from non-sampled haplotypes from intermediate populations for this lowland hummingbird species. Thus, increased population sampling along the coast of Oaxaca (between localities 11 and 12; Figure 1(c)) would be needed to further assess the hypothesis of allopatric fragmentation and the evolutionary distinctiveness of the PAC and CHIS_OAX mtDNA groups of A. rutila separated by the Isthmus.

Our survey of molecular sequence data in populations of A. rutila identified a third genetic group in the Yucatán Peninsula (YUC). According to Vázquez-López et al. (2021), the split between populations of A. rutila on the Yucatán Peninsula and those along the Pacific slope occurred 8.29–5.88 million years ago, with an age of 2.6–1.42 million years ago for the crown node of the YUC clade. The existence of unique haplotypes on the Yucatán Peninsula is consistent with a hypothesis of isolation (restricted northward gene flow). The ecological features and biogeography of the Yucatán Peninsula show the climatic effects of its 65-Myr history, submerging under warm tropical waters by numerous marine transgressions to have occurred since to the Pliocene-Miocene that resulted in basically a large limestone slab slowly emerging from south to north with a rain gradient, from a very humid in the south-southeast to dry in the north-northwest pattern (Espadas-Manrique et al., 2003; Licona-Vera et al., 2018b; Vázquez-Domínguez & Arita, 2010). The vegetation also follows the SE-NW rain gradient, from tropical rainforests in the southern Petén region to tropical scrubland in the extreme NW portion of the Yucatán peninsula, with extensive areas covered with deciduous or semideciduous tropical forests between these two extremes (Ramírez-Barahona et al., 2009; Vázquez-Domínguez & Arita, 2010). Because of this distinctiveness, the Yucatán Peninsula is biogeographically divided into two provinces, the Petén province in the south and the dry Yucatán province in the north (e.g., Espadas-Manrique et al., 2003; Ramírez-Barahona et al., 2009; Vázquez-Domínguez & Arita, 2010).

The biogeography of the region is not only influenced by the effects of the formation and emergence of the Yucatán Peninsula but also by those related to the geology of and tectonic activity and mountain ranges on the Maya block, between the biogeographic breaks Isthmus of Tehuantepec and Motagua-Polochic-Jocotán fault system (Gutiérrez-García & Vázquez-Domínguez, 2012). Although these geological and topographical features and isolation by mountain ranges in the Maya block surely influenced the evolutionary history of Mesoamerican species (e.g., González et al., 2011; Guevara-Chumacero et al., 2010; Gutiérrez-García & Vázquez-Domínguez, 2012; Jiménez & Ornelas, 2016; Rodríguez-Gómez & Ornelas, 2014; Williford et al., 2016), few studies have included phylogeographies of species with patterns of historical divergence and population genetic differentiation between groups of populations from Chiapas and from the Yucatán Peninsula (Gutiérrez-García & Vázquez-Domínguez, 2012; Licona-Vera et al., 2018b; Ortiz-Rodriguez et al., 2020; Oyama et al., 2016). In these studies, the observed patterns of genetic divergence are consistent with the hypothesis that isolation of the dry Yucatán province by semideciduous tropical rain forest along the Petén region and Chiapas restricted northward gene flow from locations of the Pacific slope in Mexico into the TDF in the extreme NW portion of the Yucatán peninsula (an “ecological barrier”; Licona-Vera et al., 2018b and references therein). Together, our results are consistent with both a model of allopatric fragmentation and reduced gene flow by both physical and environmental barriers. Although individuals from Guatemala had non-shared haplotypes with those in the Yucatán Peninsula, in support of the hypothesis that the Motagua-Polochic-Jocotán fault system is a geographic barrier, only one mutational step separated haplotype H9 from most frequent haplotype in the Yucatán Peninsula (H2). Increased molecular markers and population sampling along the coast of Chiapas (between localities 12 and 13; Figure 1(c)) and within Guatemala would be needed to further assess the hypothesis of genetic divergence by isolation and the evolutionary distinctiveness of A. rutila on the Yucatán Peninsula. Given that this region is poorly studied at the phylogeography level, the phylogeography of A. rutila in combination with a niche-modelling approach and nuclear DNA markers with higher resolution (microsatellites) is particularly important to explore the genetic structure and to understand the effects of the Motagua-Polochic-Jocotán fault system and evolution of the TDF biodiversity on the Yucatán Peninsula.

Historical Demography

Our results reveal a high level of gene diversity across the 20 populations (Table 2). The high levels of genetic differentiation (pairwise comparisons of F _ST values) and the significant phylogeographic structure may result from low levels of historical gene flow among the studied populations and implies that isolation is important in shaping the genetic differentiation. The genetic divergence between groups of populations highlights the importance of geographical barriers in driving isolation followed by intraspecific differentiation and genetic structuring of A. rutila populations (N _ST > G _ST). However, low frequency haplotypes within each of the mtDNA genetic groups (PAC, CHIS_OAX, YUC) differ from each other by only one mutation in most cases and from high frequency central haplotypes within each haplogroup, suggesting its recent formation. Also, the combination of high genetic diversity, low nucleotide diversity values, and the presence of many low frequency single haplotypes separated by few mutational steps, indicate rapid population growth from ancestral populations with small effective population size (Avise, 2000). Indeed, the negative and significant values of neutrality tests and BSPs suggest past demographic expansion without effective population size changes over time for the PAC and YUC genetic groups, and the average time since the demographic expansion was between 52.37 ka BP and 39.40 ka BP, respectively. Projections on the distribution of suitable habitat under past conditions in Vázquez-López et al. (2021) study revealed that suitable habitat for A. rutila was continuous along the Pacific slope but more restricted than its current distribution during the Last Inter Glacial (140–120 kyr BP), and isolated and restricted to the extreme NW portion of the Yucatán peninsula. They proposed that Pleistocene climatic changes were a crucial factor in the divergence and emergence of the current lineages within A. rutila. However, estimated divergence times between groups contradict this hypothesis. Although recent studies using genomic sequencing of hummingbirds have demonstrated inconsistencies of dated phylogenies analysed using nuclear and mitochondrial DNA (Andermann et al., 2019), time since the demographic expansion based on mtDNA sequences in our study suggest that population expansion within mtDNA lineages is more recent than the origin and formation of geographical barriers.

Our results highlight a pattern of low intraspecific genetic divergence in a hummingbird species with a continuous and widespread distribution in the TDF, which is unique when contrasted to higher genetic divergence observed in other co-distributed bird species along the Mexican Pacific slope. For instance, Arbeláez-Cortés et al. (2010) showed marked phylogeographic structure in three co-distributed bird species along the Mexican Pacific slope. All three species showed marked phylogeographic structure, with breaks found in roughly similar areas reported in other taxa, such as the border between the Mexican states of Nayarit and Jalisco, southern Jalisco and Michoacán, Guerrero and Oaxaca, and between Oaxaca and Chiapas, genetic breaks partially compatible with climatically stable areas.

Instead, our results are similar to the ones of other widespread hummingbird species and other bird taxa with ranges divided by geographic barriers such as the Isthmus of Tehuantepec and the Motagua-Polochic-Jocotán fault system. The climatic changes of the Pleistocene would explain the current distribution of the genetic groups, but would not provide an explanation for the marked genetic differentiation among mtDNA lineages without haplotype sharing and the distribution of the CHIS_OAX group between the eastern portion of the Isthmus of Tehuantepec and the Motagua-Polochic-Jocotán fault system, indicating that geographical barriers played a role on the differentiation of these mtDNA genetic groups. Based on past distribution models in Vázquez-López et al. (2021), the ancestral population of the CHIS_OAX and PAC groups probably presented a continuous distribution along the Pacific slope during the Pliocene. After the rise in sea level and the flooding of the Isthmus of Tehuantepec, this distribution was interrupted, causing populations to the east and west of the Isthmus to remain separated. Therefore, it is likely that the ancestral population that remained west to the Isthmus expanded its distribution along the Pacific slope and to the Tres Marías Islands (Vázquez-López et al., 2021).

Niche Divergence

The genetic structuring within A. rutila suggests isolation and low levels of gene flow between mtDNA groups of populations across their geographic ranges. This scenario is supported by niche overlap and niche equivalency tests, which indicated that the three-mtDNA groups have low-to-very limited niche’s overlap (D = 0.39–0.17) and similar but not identical niches (Table 5). The three groups of A. rutila hummingbirds occupy similar environmental space (fundamental niche), probably due to shared ancestral habitat preferences (niche conservatism; Soberón & Nakamura, 2009; Wiens & Graham, 2005). Subspecies differences in environmental space (niche divergence) have been shown in Buff-Bellied Hummingbird (Amazilia yucatanensis), sister to A. rutila (Vásquez-Aguilar et al., 2021), suggesting that the distribution of one subspecies or genetic group cannot be implied by the distribution of another one. If allopatric fragmentation occurred in A. rutila, one would expect genetically divergent groups of populations to retain certain aspects of their fundamental niche. However, significant differences between mtDNA groups in niche equivalency (niches spaces are not identical) suggest niche divergence; while allopatric fragmentation may have been important for population genetic differentiation of A. rutila, our data support an environmental differentiation scenario.

The environmental niches of the A. rutila mtDNA groups are similar in terms of temperature and precipitation based on the PCA-env results. The similarity of environmental conditions between regions might be directly related to the availability of floral resources, which in turn constitutes the limiting factor in the distribution of genetic groups (Abrahamczyk & Kessler, 2015). However, the environmental niches for the three genetic groups were not equivalent. For the PC1 (50 % of total variation), the PCA-env showed that the greater differences among the three mtDNA groups correspond to differences in temperature and precipitation variables: PC1 is positively associated with mean diurnal range in temperature (BIO2) and precipitation seasonality (BIO15), and negatively associated with precipitation of wettest and driest month (BIO13, BIO14) and precipitation of warmest and coldest quarter (BIO18, BIO19). These results indicate that precipitation variables may be more important than temperature in determining the limits in the distribution of all three mtDNA genetic groups. While these data alone do not represent a rigorous test of the causes of population divergence or adaptive evolution, the tests of niche equivalency and niche similarity show that the climatic niches of the mtDNA lineages of A. rutila are more different than expected based on random predictions, and that the differences among genetic groups exceed those predicted by the background environments of the regions they inhabit. Thus, our findings based on the non-equivalency of the fundamental niches support the hypothesis that the mtDNA lineages of A. rutila have undergone coarse-scale niche divergence and are constrained by a set of climatic and macro-environmental conditions that might determine the distribution and availability of floral resources to which A. rutila interact with within each of the regions.

Implications for Population Conservation

Significant genetic differentiation in mtDNA sequences and genetic structure among A. rutila populations suggest that isolation by geographical barriers played a role on population differentiation. We further provide evidence for rejection of the niche equivalency hypothesis that confirmed that the retrieved genetic groups exist in distinct environmental niche spaces. Given the high population genetic diversity and identification of possible processes influencing the genetic structure of the species, our study highlights recognition of evolutionary significant units for decision-making in preserving this widespread hummingbird species. Additional work is needed to understand factors other than simply variability in the background environment maintaining population separation and intraspecific divergence within A. rutila.