Skip to main content

Phylogeography and morphometric variation in the Cinnamon Hummingbird complex: Amazilia rutila (Aves: Trochilidae)

Abstract

Background

The Mesoamerican dominion is a biogeographic area of great interest due to its complex topography and distinctive climatic history. This area has a large diversity of habitats, including tropical deciduous forests, which house a large number of endemic species. Here, we assess phylogeographic pattern, genetic and morphometric variation in the Cinnamon Hummingbird complex Amazilia rutila, which prefers habitats in this region. This resident species is distributed along the Pacific coast from Sinaloa—including the Tres Marías Islands in Mexico to Costa Rica, and from the coastal plain of the Yucatán Peninsula of Mexico south to Belize.

Methods

We obtained genetic data from 85 samples of A. rutila, using 4 different molecular markers (mtDNA: ND2, COI; nDNA: ODC, MUSK) on which we performed analyses of population structure (median-joining network, STRUCTURE, FST, AMOVA), Bayesian and Maximum Likelihood phylogenetic analyses, and divergence time estimates. In order to evaluate the historic suitability of environmental conditions, we constructed projection models using past scenarios (Pleistocene periods), and conducted Bayesian Skyline Plots (BSP) to visualize changes in population sizes over time. To analyze morphometric variation, we took measurements of 5 morphological traits from 210 study skins. We tested for differences between sexes, differences among geographic groups (defined based on genetic results), and used PCA to examine the variation in multivariate space.

Results

Using mtDNA, we recovered four main geographic groups: the Pacific coast, the Tres Marías Islands, the Chiapas region, and the Yucatán Peninsula together with Central America. These same groups were recovered by the phylogenetic results based on the multilocus dataset. Demography based on BSP results showed constant population size over time throughout the A. rutila complex and within each geographic group. Ecological niche model projections onto past scenarios revealed no drastic changes in suitable conditions, but revealed some possible refuges. Morphometric results showed minor sexual dimorphism in this species and statistically significant differences between geographic groups. The Tres Marías Islands population was the most differentiated, having larger body size than the remaining groups.

Conclusions

The best supported evolutionary hypothesis of diversification within this group corresponds to geographic isolation (limited gene flow), differences in current environmental conditions, and historical habitat fragmentation promoted by past events (Pleistocene refugia). Four well-defined clades comprise the A. rutila complex, and we assess the importance of a taxonomic reevaluation. Our data suggest that both of A. r. graysoni (Tres Marías Islands) and A. r. rutila (Pacific coast) should be considered full species. The other two strongly supported clades are: (a) the Chiapas group (southern Mexico), and (b) the populations from Yucatán Peninsula and Central America. These clades belong to the corallirostris taxon, which needs to be split and properly named.

Background

The Mesoamerican dominion is the area where the Neartic and Neotropics overlap, and includes most of the Mexican and Central American lowlands as well as the Mexican Transition Zone (Morrone 2014). This biogeographic dominion and surrounding areas are well known for possessing high levels of species richness, which are thought to have originated from its complex topography, wide range of environmental conditions, and climatic history (Savage 1966). One of the most distinctive habitats found in the Mesoamerican dominion is the tropical deciduous forest, which has a high number of endemic species (Gordon and Ornelas 2000).

Tropical deciduous forest has the driest conditions of the tropical forests, and precipitation is strongly seasonal (Rzedowski 1978; Meave et al. 2012). In addition, it is considered one of the most threatened habitats in the Neotropics due to land use change, habitat fragmentation, and human population growth (Lerdau et al. 1991; Koleff et al. 2012). The diversity of the tropical deciduous forest has been attributed to historical factors such as the climatic fluctuations of the Pleistocene, in which colder and drier periods promoted its fragmentation, followed by periods of expansion with warmer and more humid conditions (Werneck et al. 2011).

Some studies of bird species distributed in Mesoamerican deciduous forests have shown high levels of genetic and morphological variation associated with historic and demographic events (e.g. Arbeláez-Cortés et al. 2014), while others have focused on allopatric populations, where long-term isolation, moderate or low gene flow, presence of geographic barriers, and historical events explained the overall geographic structure (e.g. Smith et al. 2011). One clear example of allopatric differentiation is the taxa inhabiting the Tres Marías Islands (Nayarit, Mexico), where there are several endemic subspecies (Cortés-Rodríguez et al. 2008; Smith et al. 2011; Montaño-Rendón et al. 2015). These subspecies share the common feature of clear morphological differentiation, mainly in body size, from their continental counterparts, representing independent evolutionary lineages with defined genetic structure (Ortíz-Ramírez et al. 2018). In addition to the Tres Marías Islands, the Yucatán Peninsula and the Coastal Plain of southwestern Mexico have been considered to have high potential for diversification (García-Deras et al. 2008; González et al. 2011; Smith et al. 2011; Ramírez-Barrera et al. 2018; Hernández-Baños et al. 2020). Furthermore, a study based on molecular data showed that species distributed in the Neotropical lowlands have a low and constant rate of diversification over time; extrapolating those rates of diversification to the present leads to a greater number of lineages than is currently described (Weir 2006). The number of species present in the tropics is thus not properly reflected in the currently accepted taxonomy.

Hummingbirds belong to one of the most diverse bird families (Trochilidae: 361 spp.; Gill et al. 2021), with specialized nectarivorous feeding mode and high rates of speciation. This has resulted in a huge diversity of bill morphologies and color patterns among nine highly supported groups (Bleiweiss et al. 1997; McGuire et al. 2014). In this study, we examined the evolutionary history of the Cinnamon Hummingbird (Amazilia rutila), a medium-sized hummingbird that inhabits tropical deciduous forests ranging from 0 to 1600 m of elevation (Howell and Webb 1995). This species is mainly distributed among three biogeographic provinces: the Pacific lowlands, Yucatán Peninsula, and Mosquito, based on the regionalization of the Neotropical region (Smith 1941; Ryan 1963; West 1964; Morrone 2014). It is a resident species from Sinaloa in northwestern Mexico, south along the Pacific Coast (including the Tres Marías Islands) to Costa Rica, and along the coastal plain of the Yucatán Peninsula south to Belize (including “offshore cays”) and northeastern Nicaragua (Howell and Webb 1995). Based on morphology and allopatric distributions, four subspecies are recognized (Billerman et al. 2020; Gill et al. 2021). A. r. rutila (DeLattre 1843) is distributed in western and southern Mexico (from Jalisco to Oaxaca) and has green plumage on the upperparts, cinnamon underparts and a rufous tail. A. r. diluta (van Rossem 1938) inhabits northwestern Mexico (from Sinaloa to Nayarit); its upper parts are more golden-bronze and underparts are paler and less reddish than rutila. A. r. corallirostris (Bourcier and Mulsant 1846) is distributed from south and southeastern Mexico (from Chiapas to the Yucatán Peninsula) to Costa Rica and is much more deeply colored than rutila. Finally, located within the Pacific lowlands province, the subspecies A. r. graysoni (Lawrence 1866) inhabits the Tres Marías Islands; it is darker in color and larger than rutila (Ridgway 1901). In general, the systematics of A. rutila are unresolved and controversial. Weller (1999) considered that A. rutila forms a superspecies with A. yucatanensis (Buff-bellied Hummingbird), and in a multilocus molecular study, the genus Amazilia formed a polyphyletic group nested with the genera Hylocharis, Chrysuronia, Lepidopyga, and Damophila (McGuire et al. 2007).

The aims of this study were to: (1) analyze the genetic variation and morphometric differentiation across the geographic distribution of A. rutila, (2) reconstruct phylogenetic relationships and past scenarios throughout the geographic range to analyze divergence times and demographic changes over time and space, and (3) propose an evolutionary history and taxonomic hypothesis of intraspecific lineages. We expected to find high levels of genetic structure associated with the isolation of allopatric populations and environmental conditions over continuous ranges.

Methods

Sampling, PCR and sequencing

We obtained 84 tissue samples from the A. rutila complex and one sample from A. yucatanensis as the outgroup (see Additional file 1: Tables S1, S2). We supplemented our sampling with GenBank sequences from A. rutila, A. yucatanensis and A. tzacatl (McGuire et al. 2014). Tissues samples were provided by the collection of the Museo de Zoología (MZFC) “Alfonso L. Herrera” (Universidad Nacional Autónoma de México), the Burke Museum (University of Washington), and the Biodiversity Institute (University of Kansas).

DNA was extracted using the Qiagen DNAEasy kit, following the manufacturer’s protocols. We amplified and sequenced two mitochondrial markers: partial and complete NADH dehydrogenase subunit 2 (ND2: 605–1041 bp), and the full length of cytochrome C oxidase subunit I gen (COI: 568 bp). We also amplified two nuclear markers: the regions between exons 4 and 5 of the Muscle Skeletal Receptor Tyrosine Kinase gene (MUSK: 628 bp), and a segment comprising the end of exon 6 and the beginning of exon 8 of the Ornithine Decarboxylase gene (ODC: 571). We used primers L5219, H6313, L5758 and H5766 for the amplification of ND2 (Sorenson et al. 1999); primers F1-COI and R2-COI for COI (Chaves et al. 2008); ODC2-F and ODC2-R for ODC (McGuire et al. 2007); and MUSK R3 and MUSK F3 for the amplification of MUSK (McGuire et al. 2007). All PCR products were confirmed on 1% agarose gel, and sequencing reactions were performed by the Washington University Genomics Center. All new sequences have been deposited in GenBank (Accession numbers MZ998668-MZ998740, OK624605-OK624657, OK614044-OK614080, OK618334-OK618361).

Sequences were proofread using Sequencher v.4.8 (Gene Codes Corporation 2007), and aligned with the Clustal W function in BioEdit (Thompson et al. 1994; Hall 1999). We corroborated the mitochondrial origin of our sequences using at least two of the following methods: amplifying/sequencing overlapping gene segments, amplifying/sequencing one region with different primer sets, and sequencing both DNA strands. We found no evidence of contamination of mtDNA sequences.

Phylogenetic analyses, genetic structure and neutrality tests

We used the Bayesian approach implemented in PHASE v.2.1 (Stephens et al. 2001; Stephens and Scheet 2005) to reconstruct haplotypes and avoid heterozygotes in nuclear sequences, selecting the pairs of haplotypes with a posterior probability higher than 0.90. We obtained the models of evolution for each gene using PartitionFinder (Lanfear et al. 2017), using the following parameters: linked branched lengths, greedy search algorithm (Lanfear et al. 2012), and the Bayesian Information Criterion (BIC) for model selection. The phylogenetic analyses were conducted using each locus separately and using the concatenated matrix including nuclear and mitochondrial loci (multilocus). We conducted phylogenetic analysis under Bayesian Inference (BI), using MrBayes 2.0 (Huelsenbeck and Ronquist 2002) in CIPRES Science Gateway (Miller et al. 2010). Each analysis consisted of four independent chains, random starting trees, and uniform prior distribution of parameters. The chains were run for 30 million generations, sampling trees every 1000 generations. The asymptote was determined visually, 5000 burn-in trees were discarded, and the remaining trees from the plateau phase were used to estimate Bayesian posterior probabilities. We considered that clades were strongly supported if they were present in ~ 95% of the sampled trees (Huelsenbeck and Ronquist 2002; Wilcox et al. 2002). Also, a Maximum Likelihood analysis was conducted using RaxML v.8.0.0 (Stamakis 2014) in CIPRES Science Gateway (Miller et al. 2010), using separate partitions for each locus (see “Results”), with 100 bootstrap replicates.

We defined four different groups, corresponding to allopatric populations from sampled localities of A. rutila (see Fig. 1): (1) Tres Marías Islands, (2) Pacific coast, (3) Chiapas, and (4) Yucatán Peninsula and Central America (hereafter: Tres Marías, Pacific, Chiapas and Yucatán_CA). For the comparison between populations, we specified geographic groups based on the main clades recovered from phylogenetic analysis (concatenated-multilocus), and tested for substructure in Pacific and Yucatán_CA groups (see “Results”). A median-joining network was constructed to visualize the structure of populations, the number of haplotypes, their frequencies, and the relationships among individuals using the program Network 4.5.1.0 (Bandelt et al. 1999). We tested locus neutrality with Fu’s Fs (Fu 1997) and Tajima’s D (Tajima 1989) metrics in DnaSP v.5 (Librado and Rozas 2009), and obtained the summary statistics for each population. We tested the differentiation hypothesis using an analysis of molecular variance AMOVA and the FST fixation index using Arlequin 3.1 (Excoffier and Schneider 2006). These analyses are useful for observing population structure, estimating population differentiation directly from molecular data, and testing the hypothesis of differentiation (Excoffier et al. 1992). Additionally, we used STRUCTURE 2.3.2 (Pritchard et al. 2000) for population assignments under an admixed and LocPrior model (Hubisz et al. 2009), with 10,000 iterations of burn-in and 20,000 MCMC (Markov chain Monte Carlo), for 20 replicates of each K value (from 1 to 6). To evaluate STRUCTURE results, we used the Evano method to choose K value by ∆K (Evanno et al. 2005), as implemented on the Structure Harvester website (Earl and vonHoldt 2012). Genetic distances were obtained with MEGA under Maximum Likelihood model, with 100 bootstrap replicates (Stecher et al. 2020).

Fig. 1
figure1

Geographic location of Amazilia rutila samples used for genetic analyses (left). Haplotype networks of single ND2 and COI mtDNA markers (right). Different colors represent different groups based on haplotype networks and allopatric populations: Pacific: Pacific coast, Tres Marías: Tres Marías Islands, Yucatán_CA: Yucatán Peninsula and Central America

Divergence time estimates

We estimated divergence times with the concatenated matrix using Beast v.2.6.2 (Bouckaert et al. 2019). Different partitions by gene were defined based on Pacheco et al. (2011) for ND2 and COI, Lerner et al. 2011 for ODC, and Ellegren (2007) for MUSK. We employed a strict clock and a constant population model as a tree prior. We chose the strict clock model empirically based on data fit, and the constant population model based on our demographic reconstruction (see “Results”). The analysis was run for 200 million generations, sampling every 1000. We used Tracer v.1.7 (Rambaut et al. 2018) to ensure adequate Effective Sample Size (ESS > 200), which was reached for most statistics, with the exception of prior and probability statistics. With TreeAnnotator v. 2.6.2 (Rambaut and Drummond 2013) we discarded the first 25% or trees as burn-in and produced a maximum clade credibility tree with 95% highest probability density intervals. The final tree was visualized in FigTree v.1.4.2 (Rambaut 2014).

Morphometric analysis

We took linear measurements of five morphological traits from 210 study skins housed at the American Museum of Natural History (AMNH) and at the Museo de Zoología “Alfonso L. Herrera” (MZFC). We measured wing chord (WC) and tail length (TL) using an Avinet calibrated ruler to the nearest 0.5 mm; also, we measured culmen length (CL), beak width (BW), and beak height (BH) to the nearest 0.1 mm with a Mitutoyo digital caliper.

We assigned each specimen to geographic groups defined based on the haplotype network (see “Results”; Fig. 1). We performed Mann–Whitney U tests between sexes within each geographic group, with Bonferroni correction to evaluate sexual dimorphism. We also carried out Kruskall-Wallis test with “geographic group” as the grouping variable to assess geographic variation in the measured characters. Finally, we carried out a PCA to visually examine the structure of variation in multivariate space. All statistical analyses were performed in R (R Core Team 2020).

Ecological Niche Modelling under past scenarios and demography

To infer suitability of environmental conditions under past scenarios, we performed projections using Ecological Niche Models. Presence records of A. rutila were downloaded from GBIF (Global Biodiversity Information Facility) and accessed from R via rgbif (https://github.com/ropensci/rgbif; taxon key: 2476417; https://doi.org/10.15468/dl.7k98v7). These records were cleaned over multiple steps using the following R packages: CoordinateCleaner (Zizka et al. 2019), nichetoolbox (Osorio-Olvera et al. 2016), dplyr (Wickham et al. 2020), and raster (Hijmans 2020). We used the sets of bioclimatic layers (current and past) from WorldClim (Hijmans et al. 2005; Braconnot et al. 2007; Booth et al. 2014). Our models for the past were projected onto the Last InterGlacial (LIG: 140–120 kya), the Last Glacial Maximum (LGM: 21 kya), and the Mid Holocene (MH: 6 kya) scenarios. We performed these projections using the continental phylogroups (Pacific coast, Chiapas, and Yucatán and Central America); the Tres Marías group was not included in these analyses due to the low number of occurrence records (n = 3). The selection of bioclimatic variables was based primarily on a principal component analysis, where we considered the most important variables of the first four principal components (from one to three variables), ensuring that no intercorrelated variables were selected (Pearson correlation r < 0.75, Additional file 2: Fig. S1). To define the area of accessibility for the species, we used the ellipsenm R package (Cobos et al. 2020) to delimit the area that contained occurrence points across the distribution of the species with a polygon representing a buffer area of 75 km (“M” area). Models were performed in Maxent v.3.4.1 (Phillips et al. 2021) with regularization multiplier = 2 and with 10 replicates of cross-validation with no extrapolation. To binarize outputs, we took into account the logistic threshold from maximum training sensitivity plus specificity.

We obtained Bayesian Skyline Plots (BSP) using Beast v.2.6.3 (Bouckaert et al. 2019) to infer population dynamics and changes in demography over time. We used the ND2 dataset and performed four different runs: three of them based on geographic groups (Pacific, Chiapas and Yucatán_CA), and one based on the whole A. rutila complex. The Tres Marías group was not included in an independent run because of the lack of informative sites for ND2. We set a GTR + substitution model, a strict molecular clock, a Coalescent Bayesian Skyline as the tree prior, and a Piecewise-constant skyline model (Drummond et al. 2005) with five groups. We used a mutation rate of 0.0145 substitutions per site per lineage per million years for ND2, following Lerner et al. (2011), and ran each analysis through 20 million generations.

Results

Phylogenetic analysis

The model of sequence evolution selected for the mtDNA data set was GTR + G (nst = 6, rates = gamma). The nucleotide composition of mtDNA was as follows: T = 23.9%, C = 34.5%, A = 30.8%, G = 10.6% for ND2; and T = 25.9%, C = 15.7%, A = 26.4%, G = 31.8% for COI. The parsimony informative sites for ND2 were 72 (605 bp), and 38 for COI (568 bp). For nuclear markers, the models of sequence evolution were: F81 (nst = 1) for MUSK, and GTR + I + G (nst = 6, rates = invgamma) for ODC. The nucleotide composition was: T = 29.6%, C = 19.4%, A = 29.4%, G = 21.4% for MUSK; and T = 36.4%, C = 17.3%, A = 24.7%, G = 21.4% for ODC. The parsimony informative sites for MUSK were 19 (628 bp), and 83 for ODC (571 bp).

According to the topology from the concatenated data set (mitochondrial and nuclear genes) under Bayesian Inference and Maximum Likelihood, A. rutila forms a monophyletic group with respect to A. yucatanensis and A. tzacatl (see Fig. 2). Within A. rutila, three major lineages can be identified: Pacific coast, Chiapas, and Yucatán Peninsula and Central America. There was also a well-supported clade (posterior probability = 0.93) within the Pacific phylogroup, which contained all individuals from the Tres Marías Islands. Within the Yucatán_CA phylogroup, four of the five individuals from Nicaragua grouped together into a well-supported clade (pp = 0.99; bootstrap = 99). The individual phylogenies for mitochondrial loci (ND2 and COI) recovered three major lineages (Pacific, Chiapas, and Yucatán_CA) (see Additional file 2: Fig. S1). However, individual nuclear locus phylogenies did not recover any geographic structure (see Additional file 2: Fig. S1).

Fig. 2
figure2

Phylogenies with concatenated data (mtDNA: ND2, COI; nuclearDNA: MUSK, ODC) under Bayesian Inferences and Maximum Likelihood. Four main groups are represented in different colors (see geographic location in Fig. 1). Illustrations represent three different morphotypes based on differences in body size: Tres Marías–Islands (largest size, top), Pacific–continental (middle), Yucatán_CA (smallest size, bottom). Node supports are shown on main clades (posterior probability/bootstrap). Illustrations by Ana Hernández Ramírez

Genetic structure and neutrality tests

Deviations from neutrality were rejected for mitochondrial loci (ND2 and COI) and nuclear loci, except for the ODC locus (see Table 1). We found 21 haplotypes for ND2: 9 corresponded to the Pacific group, 9 to the Yucatán-Central America group, 2 to Chiapas, and 1 to the Tres Marías Islands. The Yucatán-Central America group was separated from the rest by 57 mutational steps (Fig. 1). For COI, we found 11 haplotypes in total: 4 corresponded to the Pacific group, 3 to Chiapas, 3 to Yucatán-Central America, and 1 to the Tres Marías Islands (shared with 2 individuals from the Pacific group from Guerrero). The Yucatán-Central America group was separated from the rest by 15 mutational steps (Fig. 1). The Pacific coast and Yucatán-Central America groups had the highest haplotype and nucleotide diversity (see Table 2).

Table 1 Neutrality test by locus
Table 2 mtDNA summary statistics

The mtDNA genetic distances among groups showed an overall genetic distance of 5.24% for the Pacific group, 5.69% for Chiapas, and 7.52% for the Yucatán_CA group (Table 3). The pairwise FST fixation indices were all significant and indicated strong population structure of geographic groups for mtDNA data (Table 4). We did not find strong differentiation between the Pacific and its subgroup from Guerrero. However, the Guerrero subgroup showed a more robust differentiation with respect to the Tres Marías group than to the rest of the Pacific subgroup. Substructure was detected within the Yucatán_CA group, with some Nicaraguan samples (UWBM-6911, UWBM-68991, UWBM-69261, and UWBM-69388) separated from the rest of the group, with an FST of 0.78549. AMOVA analyses revealed that most of the genetic variation was among geographic groups (96.26%), and there was little variation within populations (3.74%; Table 3). The FST value from general AMOVA results indicated strong genetic structure and differentiation among populations (FST = 0.96; Table 5).

Table 3 mtDNA pairwise genetic distances among geographic groups in the Amazilia rutila complex
Table 4 Pairwise FST indices (below diagonal, *P ≤ 0.05) among populations in the Amazilia rutila complex
Table 5 mtDNA AMOVA comparing geographic groups in the Amazilia rutila complex

The STRUCTURE analyses detected two genetic groups (K = 2; first-level analysis; See Additional file 3: Figs. S2, S3; Fig. 3). One group comprised all individuals from Sinaloa to Chiapas, including Tres Marías. The second group included the Yucatán Peninsula and Central American individuals. To evaluate potential substructure within these groups, we ran a hierarchically structured analysis. Individual structure analyses for the first group revealed two genetic groups (K = 2), with no ancestral mixing between them. One group was composed of individuals from Chiapas and the second group recovered all of the Pacific clade including Tres Marías individuals. The substructure within the Yucatán_CA group recovered two groups (K = 2). One group included samples from western Nicaragua and three samples from Yucatán. The second group included one sample from western Nicaragua, one sample from Guatemala and the rest of the Yucatán Peninsula individuals.

Fig. 3
figure3

STRUCTURE results. The best fit K = 2 with the entire data set (above) and hierarchically within each genetic cluster (below: K = 2 for Pacific and Chiapas groups, and K = 2 for Yucatán_CA group)

Divergence time estimates

Divergence times estimated that the A. rutila complex split from its sister clade around 9.4 million years before present (Myr BP) (95% highest posterior density [HPD] 8.04‒10.87; Fig. 4), while the A. rutila complex root was dated 7.05 Myr ago (95% HPD 5.88‒8.29). This estimate also corresponds to the split between the two main clades: (1) Pacific + Tres Marías and Chiapas, and (2) Yucatán_CA. The Pacific clade and Chiapas split around 3.99 Myr BP (95% HPD 3.16‒4.84). The Pacific Clade node origin was 1.65 Myr BP (95% HPD 0.81‒1.81 Ma). According to these estimates, A. rutila arrived at the Tres Marías Islands 0.2 Myr BP (95% HPD 0.6‒0.44 Ma). The age of the continental subclades within the Pacific clade suggests that Pleistocene climatic changes may have caused contractions in the continental populations of A. rutila and resulted in incipient substructure, however although the Tres Marías Islands is well supported, other clades are not supported. The Chiapas clade arose 0.76 Myr BP (95% HPD 0.43‒1.15 Ma).

Fig. 4
figure4

Ultrametric phylogenetic tree obtained with Beast using concatenated matrix and different locus rates. Four main groups are represented with a vertical bar in different colors (see geographic location in Fig. 1). Posterior probabilities in blue

The crown node for the Yucatán_CA clade emerged 1.98 Myr BP (95% HPD 1.42‒2.6 Ma). Internal clades in the Yucatán_CA group were detected by Beast analysis, but they had weak node support, with the exception of the Nicaragua clade, which had intermediate support (0.7 posterior probability) and a date of 0.12 Myr BP (95% HPD 0.03‒0.25 Ma).

Morphometrics

We found minor sexual dimorphism in this species. Females were slightly but statistically significantly smaller than males in two body size variables we measured (WC, W = 2196.5, p < 0.001; TL, W = 2194, p < 0.001). However, the effect sizes were very small, so we carried out subsequent tests of geographic variation with both sexes pooled together.

There were statistically significant differences among geographic groups (Kruskall-Wallis tests WC: \(\chi_{4}^{2}\) = 78.86, p < 0.001; TL: \(\chi_{4}^{2}\) = 73.09, p < 0.001; BH: \(\chi_{4}^{2}\) = 47.45, p < 0.001; CL: \(\chi_{4}^{2}\) = 62.61, p < 0.001; BW: \(\chi_{4}^{2}\) = 90.30, p < 0.001). The Tres Marías phylogroup had larger values than the three continental groups for all variables. This was consistent with the results of our PCA, which showed that the continental phylogroups largely overlapped in morphological measurements, but the Tres Marías group was significantly larger. Principal component 1 explained 74.75% of the cumulative variance and roughly corresponds to overall body size, with large contributions from Wing Chord (WC) and Tail Length (TL) (Fig. 5).

Fig. 5
figure5

Principal component analysis (PCA) of morphological data, showing the first two Principal Components, which explain 91.2% of the total variation (top); and box plots of each morphometric character measured. BW beak width, BH beak height, CL culmen length, TL tail length, WC wing chord

Ecological Niche modelling onto past scenarios and demography

We obtained a total of 445 presence records for A. rutila after filtering steps. The number of occurrence points for each group was: (1) Pacific: 144 records, (2) Chiapas: 20 records, and (3) Yucatán_CA: 281 records. The results of the PCA of climate variables for the whole complex showed that the first four principal components explained most of the variation (88%, see Additional file 4: Fig. S4). According to these results and correlation values, we selected the following bioclimatic variables: BIO5 (Max Temperature of Warmest Month), BIO7 (Temperature Annual Range), BIO11 (Mean Temperature of Coldest Quarter), BIO12 (Annual Precipitation), BIO15 (Precipitation Seasonality), and BIO18 (Precipitation of Warmest Quarter). The average evaluation metrics for the model results were: training AUC = 0.7598 (St.Dv ± 0.005), and a threshold value of 0.4178 (Maximum training sensitivity plus specificity Logistic threshold). Projections onto past scenarios for the whole complex are shown in Additional file 5: Figs. S5, S6. The suitability of environmental conditions has changed somewhat over time, with subtle reductions from LIG to LGM and posterior arrangements after LGM through the MH period. Projections during LGM showed discontinuous distribution, which contrasts with the current continuous distribution for the species (Pacific coast and Central America). Independent analyses were run for each geographic group, in which the first two principal components explained most of the variation (see Additional file 5: Figs. S5, S6). The bioclimatic variables selected for the Pacific group were BIO1, BIO3, BIO4, BIO7, BIO13, BIO17 and BIO19; for the Chiapas group were BIO4, BIO5, BIO12, BIO13, BIO14, BIO15, BIO16, and BIO18; and for the Yucatán_CA group were BIO3, BIO4, BIO7, BIO11, and BIO13. The resulting metrics were as follows: Pacific AUC = 0.8046 (St.Dv ± 0.0473), threshold = 0.358; Chiapas AUC = 0.8184 (St.Dv ± 0.1119), threshold = 0.4255; Yucatán_CA AUC = 0.6515 (St.Dv ± 0.0470), threshold = 0.4515. Projections for Pacific group predicted some areas where the Chiapas and Yucatán_CA groups are currently distributed (see Additional file 5: Figs. S5, S6). Also, for the projections in the Yucatán_CA group, suitable areas were predicted where the Pacific and Chiapas groups now occur. The most conservative projections for suitable habitats were recovered when the Chiapas group was modeled. Past projections revealed a reduction of suitability during LGM and an increase during MH periods in all cases. During the LIG period, the Pacific showed continuous suitable habitat along the Pacific coast to Central America. In LIG projections for Chiapas and Yucatán_CA, no suitable conditions were recovered (see Additional file 5: Figs. S5, S6).

The BSP results showed that population dynamics and demography in the A. rutila complex have been generally constant over time with a notable population reduction close to the present. Effective population size patterns showed a stationary trend in the geographic groups analyzed (Pacific, Chiapas, Yucatán_CA), with recent population reductions in the Pacific and Chiapas groups, and a reduction followed by a slight increase towards the present in the Yucatán_CA group (see Additional file 5: Figs. S5, S6).

Discussion

Genetic variation and phylogeographic pattern

We found both deep and shallow genetic differentiation throughout the geographic distribution of A. rutila. The lineages corresponding to the Pacific clade (pp = 0.99), Chiapas (pp = 0.99) and the Yucatán Peninsula-Central America (pp = 0.99) show deep structuring according to genetic distances, FST index, AMOVA, haplotype network, multilocus phylogeny and mtDNA phylogenies. Within the Pacific clade, we can distinguish clear structure for Tres Marías Islands individuals with a monophyletic clade (pp = 0.94; geographically isolated), FST index, one haplotype for ND2 locus, and one haplotype for COI (shared with two samples from Guerrero); however, the structure analysis shows mixed ancestries between the Tres Marías Islands and Continental samples. Similarly, a multilocus phylogeny suggests a shallow substructure in the Yucatán Peninsula-Central America clade for West Nicaraguan samples (pp = 0.99), which share ancestry with some samples from Yucatán.

Despite minor differences in median-joining networks based on the ND2 versus COI markers (Fig. 1), the geographic groups had no shared haplotypes, with the exception of one haplotype that was shared between the Tres Marías group and three individuals from the Pacific group (COI marker). The information from nuclear genes that was included to construct the multilocus phylogenetic tree (see Additional file 2: Fig. S1) did not weaken the pattern found using mtDNA; the multilocus phylogeny maintained the same geographic structure with high node support values (posterior probabilities and bootstrap values), even of the retained ancestral polymorphisms, favored by nuclear signal and its reduced substitution rates compared to mtDNA (Moore 1995). A. rutila belongs to the third youngest subclade of nine Trochilidae subfamilies, known as the Emeralds (Bleiweiss et al. 1997; McGuire et al. 2007, 2014). The tempo and mode of evolution in the hummingbird higher clades has been considered to be the product of a series of historical events including colonizations, extinctions, recolonizations, and the influence of mountain chain uplift, resulting in the radiation of highland species and the establishment of lowland taxa (Bleiweiss 1998; McGuire et al. 2014). Therefore, timing is crucial to understanding phylogeographic patterns in A. rutila. This complex split from its sister group 7.05 Myr BP (95% HPD 5.88‒8.29 Myr BP), and the four main clades reported here arose 3.99‒0.20 Myr BP. During the late Miocene, Pliocene and Pleistocene, Middle American biological diversity was influenced by several major processes. One is the uplift of major mountain chains, which generated geographic barriers and created open niches in high-altitude and low-altitude habitats that favored diversification (McCormack et al. 2008). Biodiversity was also influenced by the climatic changes of Pleistocene climatic oscillations, which created a mosaic of habitats that resulted in isolated populations.

The Tres Marías phylogroup is embedded in the larger Pacific clade with a node origin dated 0.20 Myr BP (95% HPD 0.6‒0.44 Myr BP); this incomplete separation is a pattern that has been reported in other birds at early stages of speciation (e.g. Cortés-Rodríguez et al. 2008). STRUCTURE analysis did not recover substructure for Tres Marías individuals or any Pacific samples. However, the ∆K selection method has been previously shown to underestimate population structure (Janes et al. 2017). Results of the AMOVA and FST index suggest incipient population structure, and haplotype networks showed a clear differentiation at the population level. Additionally, the Tres Marías group was the most morphologically distinct compared to the other groups. The Pacific clade is represented as a polytomy in the resulting phylogenetic tree, and some clades are recovered within it but with no node support (see Additional file 2: Fig. S1). As mentioned, the single valid clade is composed of the Tres Marías Islands individuals. There are several examples of bird taxa in this area (The Pacific coast of Mexico and Tres Marías Islands) which show similar patterns of geographic variation and population structure (e.g. Habia rubica: Ramírez-Barrera et al. 2018; Phaethornis longirostris‒P. mexicanus: Arbeláez-Cortés and Navarro-Sigüenza 2013; Vireo hypochryseus: Arbeláez-Cortés et al. 2014). Those studies found that the main forces acting as drivers of divergence were the presence of geographic and ecological barriers, as well as historical habitat fragmentation. This is not the case for A. rutila along the Pacific coast, since no structured pattern was found. There are two possible explanations for this lack of geographic substructure. First, small sample sizes may cause some haplotypes to be missing, since we found many intermediate haplotypes. Another explanation is an incipient substructure produced by short periods of isolation during the Pleistocene. The age of the Pacific clade is 1.65 Myr BP (95% HPD 0.81‒1.81), and the age ranges estimated for subclades are 1.32 to 0.20 Myr BP, which corresponds with Pleistocene cycles. Thus, these mixed haplotypes could be the result of short periods of isolation during the LGM with secondary contact at the LIG and the present. This idea is supported by the LGM Pacific projections (see Additional file 5: Figs. S5, S6) and BSP analysis (Additional file 6: Fig. S7). In contrast, STRUCTURE hierarchical analysis did not show different groups within the Pacific Clade (Fig. 2).

The Chiapas clade is well defined, supported by genetic distance (5.69%) and FST values (0.97‒1.00). This region east of the Isthmus of Tehuantepec is mainly influenced by historical events between the middle and late Pliocene and marine transgressions and regressions during the late Pliocene (Barber and Klicka 2010; Peterson et al. 2010). In general, the Isthmus of Tehuantepec represents a geographic barrier for highland bird species (e.g. Zamudio-Beltrán and Hernández-Baños 2018), but not for lowland species. However, some phylogeographic studies of lowland birds have described the influence of this low altitude zone (e.g. Campylorhynchus rufinucha: Vázquez-Miranda et al. 2009). According to our time-calibrated tree, the Chiapas clade separated during Pleistocene climatic fluctuations ~ 760 kyr BP (95% HPD 0.43‒1.15 Ma), as there was a disruption of suitable conditions during the Last Glacial Maxima across the Isthmus of Tehuantepec (see Fig. 5). The Yucatán Peninsula and Central America clade also represents an independent entity, separated by 57 and 15 mutational steps, respectively, in the ND2 and COI trees (see Fig. 1), assuming low levels of gene flow, with an origin 1.98 Myr BP (95% HPD 1.42‒2.6 Myr BP). This genetic signal of isolation from the other groups may be influenced by the geographic distance and differences in environmental spaces. The Pacific and Chiapas group are distributed along the tropical deciduous forest belt, which has marked dry and rainy seasons and has been repeatedly identified as a distinct biogeographical province (i.e. Gordon and Ornelas 2000), while the Yucatán Peninsula has moderately pronounced dry and rainy seasons with relatively constant average temperature throughout the year (Paynter 1955). Also, our ENM past projections suggest periods of habitat contractions, during which populations from the Yucatán Peninsula and Central America could have been closer from each other (LGM MIROC, see Additional file 5: Figs. S5, S6), followed by periods of isolation. Although A. rutila is allopatrically distributed in the Yucatán Peninsula and Central America, the studied individuals were placed in a single group. This link between individuals close to the Yucatán Peninsula and Central America was also found in A. tzacatl (Miller et al. 2011), in which populations in Tabasco (Mexico), Belize, Honduras, and Costa Rica shared several haplotypes despite the long distances between them, and the most differentiated populations were the southernmost populations in Panama. The periods of isolation mentioned above combined with demographic decreases (see BSP results) could explain the shallow substructure found in the West Nicaragua samples (pp = 0.99). However, our STRUCTURE analysis shows a mixed ancestry between those samples and some Yucatán individuals (Fig. 2).

The influence of Pleistocene climatic fluctuations on processes of diversification and incipient speciation have been reported for many birds (e.g. Euphonia affinis: Vázquez-López et al. 2020), including hummingbirds (e.g. Eugenes fulgens: Zamudio-Beltrán et al. 2020). These habitat fluctuations are well documented as factors that promote population contractions, limiting gene flow (refugia hypothesis; Arango et al. 2021), as well as expansions, resulting in zones of contact between distant populations, with current signals of shared haplotypes and polyphyletic patterns (e.g. Mesoamerican Amazilia complex: Jiménez and Ornelas 2016; White-chested Amazilia complex: Rodríguez-Gómez and Ornelas 2018). Diversifications occurring throughout the Miocene, when the estimates of divergence dates of the genus Amazilia are presumed ca. 19.85 Mya (i.e. Ornelas et al. 2014), may have depended less on these climatic conditions. However, for events traced after this period, the rapid changes in climate cycling were crucial (Bleiweiss 1998). Furthermore, phylogeographic patterns of species sharing (totally or partially) the geographic distribution of A. rutila show similarities (e.g. Turdus rufopalliatus: Montaño-Rendón et al. 2015; Ortíz-Ramírez et al. 2018; de Silva et al. 2020). In general, the diversification processes have resulted from a variety of events, sharing the explanations of an intricate orographic history, climatic fluctuations during the Pleistocene, and thus low signals of gene flow favored by the isolation of populations due to distance or historical refuges.

Morphometrics

Our results showed little variation among the continental phylogroups of this species in the morphological characters we measured, despite A. rutila’s large geographic range. However, the Tres Marías Islands phylogroup is larger than the continental phylogroups in four out of five measurements (Fig. 3). This is consistent with several other avian species from Tres Marías (Grant 1965) and is a clear example of the “Island Rule” (van Valen 1973; Lomolino 1985; Clegg and Owens 2002; Boyer and Jetz 2010, Benítez-López et al. 2021). This shift in body size has been attributed to changes in thermal ecology (Clegg and Owens 2002), as the Tres Marías Islands are significantly drier and warmer than the continental habitat of the species. However, dramatic divergence in body size can be the result of geographic isolation independent of environmental pressure (Seeholzer and Brumfield 2018). This divergence in body size with relatively similar beak sizes suggest that beak morphology might be more tightly constrained than body size due to this species’ feeding ecology.

Implications for taxonomy

We present a taxonomic proposal for the A. rutila complex based on the recovery of four independent clades with high support values corresponding to natural geographic groups.

Our results support the first differentiated clade, A. r. graysoni, distributed in the Tres Marías Islands (Nayarit, Mexico). This differentiation was also confirmed by our morphometric analyses, which showed that the Tres Marías Islands population is the largest morphotype. A. r. graysoni was first suggested as an independent species by Ridgway (1901), who described Amazilia graysoni (Grayson’s Hummingbird) as a separate species similar in coloration to A. rutila, but darker and much larger. Then following the alternative phylogenetic/evolutionary species taxonomy of the Mexican avifauna, Navarro-Sigüenza and Peterson (2004) proposed A. r. graysoni as an independent evolutionary species, which is morphologically distinct in size and coloration from the rest of the subspecies within A. rutila. The proposal of elevating this taxon to species level has been discussed recently (de Silva et al. 2020). The destruction of native habitat, the introduction of exotic species, and reduced population sizes increase the vulnerability of birds on the Pacific Islands of Mexico and must be immediately addressed in conservation policies (Hahn et al. 2012; de Silva et al. 2020).

The second differentiated clade corresponds to the group distributed along the Pacific coast. Within this second group, two subspecies have been described (A. r. rutila and A. r. diluta), which, according to our results, are not observable as separate groups. Therefore, we propose that the Pacific coast phylogroup be represented by A. r. rutila (the first named).

The third phylogroup is distributed in Chiapas, Mexico (east of the Isthmus of Tehuantepec). This group corresponds to A. r. corallirostris, which is distributed from south and southeastern Mexico to Costa Rica. The Chiapas phylogrouop is distributed in an important and well-studied main fauna region, considered an area of endemism in western Mexico with one of the highest values of avian species richness (Peterson and Navarro 2000; García-Trejo and Navarro-Sigüenza 2004).

Finally, the fourth independent group corresponds to the clade of the Yucatán Peninsula—Central America, the most genetically differentiated (7.52%), which also corresponds to A. r. corallirostris. According to our results, A. r. corallirostris is a polyphyletic subspecies, since it is composed of two independent groups. A. r. corallirostris was described by Bourcier and Mulsant (1846) with specimens from Escuintla, Guatemala. Therefore, the name A. corallirostris should be applied to the Yucatán-CA clade (Ridgway 1901). This left the Chiapas group without a name, but specimens from Huehuetán, Chiapas were described under A. cinnamomea saturata by Nelson in 1898, and this name was later considered a synonym of A. r. corallirostris (Ridgway 1901). As a consequence, the Chiapas lineage could be named A. saturata (Nelson 1898).

Conclusions

The evolutionary hypothesis of diversification within this group is of geographic isolation (limited gene flow), differences in current environmental conditions, and historical habitat fragmentation promoted by past scenarios (Pleistocene refuges). The A. rutila complex consists of four well-defined clades that, in our opinion, merit taxonomic reevaluation. Our data suggest that A. r. graysoni (Tres Marías Islands) as well as A. r. rutila (Pacific coast) should be considered full species. The other two strongly supported clades are: (a) the Chiapas group (southern Mexico), and (b) the populations from Yucatán Peninsula and Central America. These clades belong to the A. rutila corallirostris taxon as currently defined. We propose that this group be split and renamed as two species: A. saturata for the Chiapas group and A. corallirostris for the Yucatán Peninsula and Central America group.

Availability of data and materials

The sequences generated during the current study are available in GenBank (accession numbers: MZ998668-MZ998740, OK624605-OK624657, OK614044-OK614080, OK618334-OK618361).

References

  1. Arango A, Villalobos F, Prieto-Torres DA, Guevara R. The phylogenetic diversity and structure of the seasonally dry forest in the Neotropics. J Biogeogr. 2021;48:176–86.

    Google Scholar 

  2. Arbeláez-Cortés E, Navarro-Sigüenza AG. Molecular evidence of the taxonomic status of western Mexican populations of Phaethornis longirostris (Aves: Trochilidae). Zootaxa. 2013;3716:81–97.

    PubMed  Google Scholar 

  3. Arbeláez-Cortés E, Roldán-Piña D, Navarro-Sigüenza AG. Multilocus phylogeography and morphology give insights into the recent evolution of a Mexican endemic songbird: Vireo hypochryseus. J Avian Biol. 2014;45:253–63.

    Google Scholar 

  4. Bandelt HJ, Foster P, Röhl A. Median-joining networks for inferring intraspecific phylogenies. Mol Biol Evol. 1999;16:37–48.

    CAS  PubMed  Google Scholar 

  5. Barber BR, Klicka J. Two pulses of diversification across the Isthmus of Tehuantepec in a montane Mexican bird fauna. Proc R Soc B Biol Sci. 2010;277:2675–81.

    CAS  Google Scholar 

  6. Benítez-López A, Santini L, Gallego-Zamorano J, Milá B, Walkden P, Huijbregts MA, et al. The island rule explains consistent patterns of body size evolution in terrestrial vertebrates. Nat Ecol Evol. 2021;5:768–86.

    PubMed  Google Scholar 

  7. Billerman SM, Keeney BK, Rodewald PG, Schulenberg TS. Birds of the world. Ithaca: Cornell Laboratory of Ornithology; 2020. https://birdsoftheworld.org/bow/home.

  8. Bleiweiss R. Tempo and mode of hummingbird evolution. Biol J Linn Soc. 1998;65:63–76.

    Google Scholar 

  9. Bleiweiss R, Kirsch JA, Matheus JC. DNA hybridization evidence for the principal lineages of hummingbirds (Aves: Trochilidae). Mol Biol Evol. 1997;14:325–43.

    CAS  PubMed  Google Scholar 

  10. Booth TH, Nix HA, Busby JR, Hutchinson MF. BIOCLIM: the first species distribution modelling package, its early applications and relevance to most current MAXENT studies. Divers Distrib. 2014;20:1–9.

    Google Scholar 

  11. Bouckaert R, Vaughan TG, Barido-Sottani J, Duchêne S, Fourment M, Gavryushkina A, et al. BEAST 2.5: an advanced software platform for Bayesian evolutionary analysis. PLoS Comput Biol. 2019;15:e1006650.

    CAS  PubMed  PubMed Central  Google Scholar 

  12. Bourcier J, Mulsant E. Description de vingt espèces nouvelles d’oiseaux-mouches. Ann Sci Phys Nat Agric. 1846;9:312–32.

    Google Scholar 

  13. Boyer AG, Jetz W. Biogeography of body size in Pacific island birds. Ecography. 2010;33:369–79.

    Google Scholar 

  14. Braconnot P, Otto-Bliesner B, Harrison S, Joussaume S, Peterchmitt JY, Abe-Ouchi A, et al. Results of PMIP2 coupled simulations of the Mid-Holocene and Last Glacial Maximum-Part 1: experiments and large-scale features. Clim Past. 2007;3:261–77.

    Google Scholar 

  15. Cavender-Bares J, Gonzalez-Rodriguez A, Pahlich A, Koehler K, Deacon N. Phylogeography and climatic niche evolution in live oaks (Quercus series Virentes) from the tropics to the temperate zone. J Biogeogr. 2011;38:962–81.

    Google Scholar 

  16. Chaves AV, Clozato CL, Lacerda DR, Sari H, Santos FR. Molecular taxonomy of Brazilian tyrant-flycatchers (Passeriformes, Tyrannidae). Mol Ecol. 2008;8:1169–77.

    CAS  Google Scholar 

  17. Clegg SM, Owens PF. The “island rule” in birds: medium body size and its ecological explanation. Proc R Soc B Biol Sci. 2002;269:1359–65.

    Google Scholar 

  18. Cobos ME, Osorio-Olvera L, Soberón J, Peterson AT. ellipsenm: ecological niche characterizations using ellipsoids. R package; 2020. https://github.com/marlonecobos/ellipsenm.

  19. Cortés-Rodríguez N, Hernández-Baños BE, Navarro-Sigüenza AG, Omland KE. Geographic variation and genetic structure in the Streak-backed Oriole: low mitochondrial DNA differentiation reveals recent divergence. Condor. 2008;4:729–39.

    Google Scholar 

  20. de Silva HG, Pérez Villafaña MG, Cruz-Nieto J, Cruz-Nieto MÁ. Are some of the birds endemic to the Tres Marías Islands (Mexico) species? Bull Br Ornithol Club. 2020;140:7–37.

    Google Scholar 

  21. DeLattre PA. Oiseaux-mouches nouveaux au peu connus, découverts au Guatimala. L’echo Du Monde Savant. 1843;7:1068–70.

    Google Scholar 

  22. Drummond AJ, Rambaut A, Shapiro B, Pybus OG. Bayesian coalescent inference of past population dynamics from molecular sequences. Mol Biol Evol. 2005;22:1185–92.

    CAS  PubMed  Google Scholar 

  23. Earl DA, vonHoldt BM. STRUCTURE HARVESTER: a website and program for visualizing STRUCTURE output and implementing the Evanno method. Conserv Genet Res. 2012;4:359–61.

    Google Scholar 

  24. Ellegren H. Molecular evolutionary genomics of birds. Cytogenet Genome Res. 2007;117:120–30.

    CAS  PubMed  Google Scholar 

  25. Evanno G, Regnaut S, Goudet J. Detecting the number of clusters of individuals using the software STRUCTURE: a simulation study. Mol Ecol. 2005;14:2611–20.

    CAS  PubMed  PubMed Central  Google Scholar 

  26. Excoffier LG, Schneider S. Arlequin v. 3.1: an integrated software package for population genetics data analysis. Evol Bioinform Online. 2006;1:47–50.

    Google Scholar 

  27. Excoffier L, Smouse PE, Quattro JM. Analysis of molecular variance inferred from metric distances among DNA haplotypes: application to human mitochondrial restriction region. Genetics. 1992;131:479–91.

    CAS  PubMed  PubMed Central  Google Scholar 

  28. Fu Y. Statistical tests of neutrality of mutations against population growth, hitchhiking and background selection. Genetics. 1997;147:915–25.

    CAS  PubMed  PubMed Central  Google Scholar 

  29. García-Deras GM, Cortés-Rodríguez N, Honey M, Navarro-Sigüenza AG, García-Moreno J, Hernández-Baños BE. Phylogenetic relationships within the genus Cynanthus (Aves: Trochilidae), with emphasis on C. doubledayi. Zootaxa. 2008;1742:61–8.

    Google Scholar 

  30. García-Trejo EA, Navarro-Sigüenza AG. Patrones biogeográficos de la riqueza de especies y el endemismo de la avifauna en el oeste de México. Acta Zool Mex. 2004;20:167–85.

    Google Scholar 

  31. Gill F, Donsker D, Rasmussen P. IOC World Bird List (v.11.1). 2021. https://www.worldbirdnames.org/.

  32. González C, Ornelas JF, Gutiérrez-Rodríguez C. Selection and geographic isolation influence hummingbird speciation: genetic, acoustic and morphological divergence in the wedge-tailed sabrewing (Campylopterus curvipennins). BMC Evol Biol. 2011;11:38.

    PubMed  PubMed Central  Google Scholar 

  33. Gordon CE, Ornelas JF. Comparing endemism and habitat restriction in Mesoamerican tropical deciduous forest birds: implications for biodiversity conservation planning. Bird Conserv Int. 2000;10:289–303.

    Google Scholar 

  34. Grant PR. The adaptive significance of some size trends in island birds. Evolution. 1965;19:355–67.

    Google Scholar 

  35. Hahn IJ, Hogeback S, Römer U, Vergara PM. Biodiversity and biogeography of birds in Pacific Mexico along an isolation gradient from mainland Chamela via coastal Marías to oceanic Revillagigedo Islands. Vertebr Zool. 2012;62:123–44.

    Google Scholar 

  36. Hall TA. BioEdit: a user-friendly biological sequence alignment editor and analysis program for Windows 95/98/NT. Nucleic Acids Symp Ser. 1999;41:95–8.

    CAS  Google Scholar 

  37. Hernández-Baños BE, Zamudio-Beltrán LE, Milá B. Phylogenetic relationships and systematics of a subclade of Mesoamerican emerald hummingbirds (Aves: Trochilidae: Trochilini). Zootaxa. 2020;4748:581–91.

    Google Scholar 

  38. Hijmans RJ, Cameron SE, Parra JL, Jones PG, Jarvis A. Very high resolution interpolated climate surfaces for global land areas. Int J Climatol. 2005;25:1965–78.

    Google Scholar 

  39. Hijmans RJ. raster: geographic data analysis and modeling. R package. version 3.0-12. 2020. https://CRAN.R-project.org/package=raster.

  40. Howell SN, Webb S. A guide to the birds of Mexico and northern Central America. Oxford: Oxford University Press; 1995.

    Google Scholar 

  41. Hubisz M, Falush D, Stephens M, Pritchard J. Inferring weak population structure with the assistance of sample group information. Mol Ecol Resour. 2009;9:1322–32.

    PubMed  PubMed Central  Google Scholar 

  42. Huelsenbeck JP, Ronquist F. MrBayes: a program for the Bayesian inference of phylogeny. Bioinformatics. 2002;17:754–5.

    Google Scholar 

  43. Janes JK, Miller JM, Dupuis JR, Malenfant RM, Gorrell JC, Cullingham CI, et al. The K = 2 conundrum. Mol Ecol. 2017;26:3594–602.

    PubMed  Google Scholar 

  44. Jiménez A, Ornelas JF. Historical and current introgression in a Mesoamerican hummingbird species complex: a biogeographic perspective. PeerJ. 2016;4:e1556.

    PubMed  PubMed Central  Google Scholar 

  45. Koleff P, Urquiza-Haas T, Contreras B. Prioridades de conservación de los bosques tropicales en México: reflexiones sobre su estado de conservación y manejo. Ecosistemas. 2012;21:6–20.

    Google Scholar 

  46. Lanfear R, Calcott B, Ho SY, Guindon S. PartitionFinder: combined selection of partitioning schemes and substitution models for phylogenetic analyses. Mol Biol Evol. 2012;29:1695–701.

    CAS  PubMed  Google Scholar 

  47. Lanfear R, Frandsen PB, Wright AM, Senfeld T, Calcott B. PartitionFinder 2: new methods for selecting partitioned models of evolution for molecular and morphological phylogenetic analyses. Mol Biol Evol. 2017;34:772–3.

    CAS  PubMed  Google Scholar 

  48. Lawrence GN. Annals of the Lycaeum of natural history of New York. New York: Lyceum of Natural History; 1866.

    Google Scholar 

  49. Lerdau M, Whitbeck J, Holbrook NM. Tropical deciduous forest: death of a biome. Trends Ecol Evol. 1991;6:201–2.

    CAS  PubMed  Google Scholar 

  50. Lerner HRL, Meyer M, James HF, Hofreiter M, Fleischer RC. Multilocus resolution of phylogeny and timescale in the extant adaptive radiation of Hawaiian honeycreepers. Curr Biol. 2011;21:1838–44.

    CAS  PubMed  Google Scholar 

  51. Librado P, Rozas J. DnaSP v5: a software for comprehensive analysis of DNA polymorphism data. Bioinformatics. 2009;25:1451–2.

    CAS  PubMed  Google Scholar 

  52. Lomolino MV. Body size of mammals on islands: the island rule reexamined. Am Nat. 1985;125:310–6.

    Google Scholar 

  53. McCormack JE, Peterson AT, Bonaccorso E, Smith TB. Speciation in the highlands of Mexico: genetic and phenotypic divergence in the Mexican jay (Aphelocoma ultramarina). Mol Ecol. 2008;17:2505–21.

    CAS  PubMed  Google Scholar 

  54. McGuire JA, Witt CC, Altshuler DL, Remsen JV Jr. Phylogenetic systematics and biogeography of hummingbirds: Bayesian and maximum likelihood analyses of partitioned data and selection of an appropriate partitioning strategy. Syst Biol. 2007;56:837–56.

    CAS  PubMed  Google Scholar 

  55. McGuire JA, Witt CC, Remsen JV Jr, Corl A, Rabosky DL, Altshuler DL, et al. Molecular phylogenetics and the diversification of hummingbirds. Curr Biol. 2014;24:910–6.

    CAS  PubMed  Google Scholar 

  56. Meave JA, Romero-Romero MA, Salas-Morales SH, Pérez-García EA, Gallardo-Cruz JA. Diversidad, amenazas y oportunidades para la conservación del bosque tropical caducifolio en el estado de Oaxaca, México. Ecosistemas. 2012;21:85–100.

    Google Scholar 

  57. Miller MJ, Lelevier MJ, Bermingham E, Klicka JT, Escalante P, Winker K. Phylogeography of the Rufous-tailed Hummingbird (Amazilia tzacatl). Condor. 2011;113:806–16.

    Google Scholar 

  58. Miller MA, Pfeiffer W, Schwartz T. Creating the CIPRES Science Gateway for inference of large phylogenetic trees. Gateway computing environments workshop. New Orleans: IEEE; 2010.

  59. Montaño-Rendón M, Sánchez-González LA, Hernández-Alonso G, Navarro-Sigüenza AG. Genetic differentiation in the Mexican endemic Rufous-backed Robin, Turdus rufopalliatus (Passeriformes: Turdidae). Zootaxa. 2015;4034:495–514.

    PubMed  Google Scholar 

  60. Moore WS. Inferring phylogenies from mtDNA variation: mitochondrial-gene trees versus nuclear-gene trees. Evolution. 1995;49:718–26.

    PubMed  Google Scholar 

  61. Morrone JJ. Biogeographical regionalisation of the Neotropical region. Zootaxa. 2014;3782:1–110.

    PubMed  Google Scholar 

  62. Navarro-Sigüenza AG, Peterson AT. An alternative species taxonomy of the birds of Mexico. Biota Neotrop. 2004;4:1–32.

    Google Scholar 

  63. Nelson EW. Descriptions of new birds from the Tres Marias Islands, Western Mexico. Proc Biol Soc Washington. 1898;12:63-4.

    Google Scholar 

  64. Ornelas JF, González C, de los Monteros AE, Rodríguez-Gómez F, García-Feria LM. In and out of Mesoamerica: temporal divergence of Amazilia hummingbirds pre-dates the orthodox account of the completion of the Isthmus of Panama. J Biogeogr. 2014;41:168–81.

    Google Scholar 

  65. Ortiz-Ramírez MF, Sánchez-González LA, Castellanos-Morales G, Ornelas JF, Navarro-Sigüenza AG. Concerted Pleistocene dispersal and genetic differentiation in passerine birds from the Tres Marías Archipelago, Mexico. Auk. 2018;135:716–32.

    Google Scholar 

  66. Osorio-Olvera L, Barve V, Barve N, Soberón J. nichetoolbox: from getting biodiversity data to evaluating species distribution models in a friendly GUI environment. R package. version 0.2.0.0. 2016. http://shiny.conabio.gob.mx:3838/nichetoolb2/.

  67. Pacheco MA, Battistuzzi FU, Lentino M, Aguilar RF, Kumar S, Escalante AA. Evolution of modern birds revealed by mitogenomics: timing the radiation and origin of major orders. Mol Biol Evol. 2011;28:1927–42.

    CAS  PubMed  PubMed Central  Google Scholar 

  68. Ridgway R. Part V: Family Trochilidae in: The birds of North and Middle America: a descriptive catalogue. Washington: Bulletin of the United States National Museum No. 50; 1901.

  69. Paynter RA. The ornithogeography of the Yucatán Peninsula. New Haven: Yale University; 1955.

    Google Scholar 

  70. Peterson AT, Navarro AG. Western Mexico: a significant center of avian endemism and challenge for conservation action. Cotinga. 2000;14:42–6.

    Google Scholar 

  71. Peterson AT, Navarro-Sigüenza AG, Li X. Joint effects of marine intrusion and climate change on the Mexican avifauna. Ann Assoc Am Geogr. 2010;100:908–16.

    Google Scholar 

  72. Phillips SJ, Dudík M, Schapire RE. Maxent software for modeling species niches and distributions. version 3.4.1. 2021.

  73. Pritchard JK, Stephens M, Donnelly P. Inference of population structure using multilocus genotype data. Genetics. 2000;155:945–59.

    CAS  PubMed  PubMed Central  Google Scholar 

  74. R Core Team. R: a language and environment for statistical computing. Vienna: R Foundation for Statistical Computing; 2020. https://www.R-project.org/.

  75. Rambaut A, Drummond AJ, Xie D, Baele G, Suchard MA. Posterior summarization in Bayesian phylogenetics using Tracer 1.7. Syst Biol. 2018;67:901–4.

    CAS  PubMed  PubMed Central  Google Scholar 

  76. Rambaut A, Drummond AJ. TreeAnnotator v.2.6.2:MCMC Output analysis. 2013; http://beast.bio.ed.ac.uk.

  77. Rambaut A. FigTree v1.4.2: Tree figure drawing tool. 2014; http://tree.bio.ed.ac.uk.

  78. Ramírez-Barrera SM, Hernández-Baños BE, Jaramillo-Correa JP, Klicka J. Deep divergence of Red-crowned Ant Tanager (Habia rubica: Cardinalidae), a multilocus phylogenetic analysis with emphasis in Mesoamerica. PeerJ. 2018;6:e5496.

    PubMed  PubMed Central  Google Scholar 

  79. Rodríguez-Gómez F, Ornelas JF. Genetic structuring and secondary contact in the white-chested Amazilia hummingbird species complex. J Avian Biol. 2018;49:e01536.

    Google Scholar 

  80. Ryan RM. The biotic provinces of Central America. Acta Zool Mex. 1963;6:1–54.

    Google Scholar 

  81. Rzedowski J. Vegetación de México. México: Limusa; 1978.

    Google Scholar 

  82. Savage JM. The origins and history of the Central America herpetofauna. Copeia. 1966;4:719–66.

    Google Scholar 

  83. Seeholzer GF, Brumfield RT. Isolation by distance, not incipient ecological speciation, explains genetic differentiation in an Andean songbird (Aves: Furnariidae: Cranioleuca antisiensis, Line-cheeked Spinetail) despite near threefold body size change across an environmental gradient. Mol Ecol. 2018;27:279–96.

    PubMed  Google Scholar 

  84. Smith H. Las provincias bióticas de México, según la distribución geográfica de las lagartijas del género Sceloporus. Anales De La Escuela Nacional De Ciencias Biológicas. 1941;2:103–10.

    Google Scholar 

  85. Smith BT, Escalante P, Hernández-Baños BE, Navarro-Sigüenza AG, Rohwer S, Klicka J. The role of historical and contemporary processes on phylogeographic structure and genetic diversity in the Northern Cardinal, Cardinalis Cardinalis. BMC Evol Biol. 2011;11:136.

    PubMed  PubMed Central  Google Scholar 

  86. Sorenson MD, Ast JC, Dimcheff DE, Yuri T, Mindell DP. Primers for a PCR-based approach to mitochondrial genome sequencing in birds and other vertebrates. Mol Phylogenet Evol. 1999;12:105–14.

    CAS  PubMed  Google Scholar 

  87. Stamakis A. RaxML version 8: a tool for phylogenetic analysis and post-analysis of large phylogenies. Bioinformatics. 2014;30:1312–3.

    Google Scholar 

  88. Stecher G, Tamura K, Kumar S. Molecular evolutionary genetics analysis (MEGA) for macOS. Mol Biol Evol. 2020;37:1237–9.

    CAS  PubMed  PubMed Central  Google Scholar 

  89. Stephens M, Scheet P. Accounting for decay of linkage disequilibrium in haplotype inference and missing-data imputation. Am J Human Genet. 2005;76:449–62.

    CAS  Google Scholar 

  90. Stephens M, Smith N, Donelly P. A new statistical method for haplotype reconstruction from population data. Am J Human Genet. 2001;68:978–89.

    CAS  Google Scholar 

  91. Tajima F. Statistical method for testing the neutral mutation hypothesis by DNA polymorphism. Genetics. 1989;123:585–95.

    CAS  PubMed  PubMed Central  Google Scholar 

  92. Thompson JD, Higgins DG, Gibson TJ. CLUSTAL W: improving the sensitivity of progressive multiple sequence alignment through sequence weighting, position-specific gap penalties and weight matrix choice. Nucleic Acids Res. 1994;22:4673–80.

    CAS  PubMed  PubMed Central  Google Scholar 

  93. van Rossem AJ. A northwest race of the Cinnamon Hummingbird. Condor. 1938;40:226–7.

    Google Scholar 

  94. van Valen LM. Pattern and the balance of nature. Evol Theor. 1973;1:31–49.

    Google Scholar 

  95. Vázquez-López AM, Morrone JJ, Ramírez-Barrera SM, López-López A, Robles-Bello SM, Hernández-Baños BE. Multilocus, phenotypic, behavioral, and ecological niche analyses provide evidence for two species within Euphonia affinis (Aves, Fringilidae). ZooKeys. 2020;952:129–57.

    PubMed  PubMed Central  Google Scholar 

  96. Vázquez-Miranda H, Navarro-Sigüenza AG, Omland KE. Phylogeography of the Rufous-naped Wren (Campylorhynchus rufinucha): speciation and hybridization in Mesoamerica. Auk. 2009;126:765–78.

    Google Scholar 

  97. Weir JT. Divergent timing and patterns of species accumulation in lowland and highland neotropical birds. Evolution. 2006;60:842–55.

    PubMed  Google Scholar 

  98. Weller AA. Cinnamon Hummingbird. In: del Hoyo JJ, Elliott A, Sargatal J, editors. Handbook of the birds of the world. Barn-owls to hummingbirds. Barcelona: Lynx Editions; 1999. p. 596‒597.

  99. Werneck FP, Costa GC, Colli GR, Prado DE, Sites JW. Revisiting the historical distribution of Seasonally Dry Tropical Forests: new insights based on palaeodistribution modelling and palynological evidence. Global Ecol Biogeogr. 2011;20:272–88.

    Google Scholar 

  100. West RC. The natural regions of Middle America. In: West RC, editor. Handbook of Middle American Indians. Austin: University of Texas Press; 1964. p. 363–83.

    Google Scholar 

  101. Wickham H, Francois R, Henry L, Müller K. dplyr: A Grammar of Data Manipulation. R package. version 0.8.5. 2020. https://CRAN.R-project.org/package=dplyr.

  102. Wilcox TP, Zwickl DJ, Heath TA, Hillis DM. Phylogenetic relationships of the dwarf boas and a comparison of Bayesian and bootstrap measures of phylogenetic support. Mol Phylogenet Evol. 2002;25:361–71.

    CAS  PubMed  Google Scholar 

  103. Zamudio-Beltrán LE, Hernández-Baños BE. Genetic and morphometric divergence in the Garnet-Throated Hummingbird Lamprolaima rhami (Aves: Trochilidae). PeerJ. 2018;6:e5733.

    PubMed  PubMed Central  Google Scholar 

  104. Zamudio-Beltrán LE, Ornelas JF, Malpica A, Hernández-Baños BE. Genetic and morphological differentiation among populations of the Rivoli’s Hummingbird (Eugenes fulgens) species complex (Aves: Trochilidae). Auk. 2020;137:ukaa032.

    Google Scholar 

  105. Zizka A, Silvestro D, Andermann T, Azevedo J, Duarte-Ritter C, Edler D, et al. CoordinateCleaner: standardized cleaning of occurrence records from biological collection databases. Method Ecol Evol. 2019;10:744–51.

    Google Scholar 

Download references

Acknowledgements

We thank M. Robbins of the University of Kansas (Natural History Museum), and S. Birks of the University of Washington (Burke Museum of Natural History and Culture) for providing tissue samples used in this study. We also thank Alejandro Gordillo Martínez, Fabiola Ramírez, Isabel Vargas-Fernandez and Raúl Ivan Martínez Becerril for technical support. We acknowledge the General Direction of Wildlife (Instituto de Ecología, SEMARNAT, México) for providing collecting permits. L. Kiere reviewed the English.

Funding

This research was supported by PAPIIT/DGAPA, Universidad Nacional Autónoma de México (UNAM) through a grant to Blanca E. Hernández-Baños (IN220620). LEZ-B acknowledges the Postdoctoral scholarship provided by DGAPA-UNAM.

Author information

Affiliations

Authors

Contributions

BEH-B and MV-L designed the study. BEH-B secured financial support. MV-L carried out the laboratory work. MV-L, BEH-B, SMR-B and LEZ-B analyzed the data. BEH-B, MV-L, NC-R, AB-H, LEZ-B and KR contributed to the writing and improvement of the manuscript. All authors read and approved the final manuscript.

Corresponding author

Correspondence to Blanca E. Hernández-Baños.

Ethics declarations

Ethics approval and consent to participate

This study was performed in line with the principles of the Declaration of Helsinki.

Consent for publication

Not applicable.

Competing interests

The authors declare that they have no competing interests.

Supplementary Information

Additional file 1: Table S1.

Information of voucher numbers (ID), geographic locations (longitude, latitude), and biological collections of tissue samples used in this study. Table S2. Information of accession numbers (AN) and geographic locations (locality, longitude, latitude) of GenBank’s sequences used in this study.

Additional file 2: Figure S1.

The Bayesian Inference and Maximum Likelihood of the mitochondrial DNA and nuclear DNA.

Additional file 3: Figure S2.

STRUCTURE analysis of Amazilia rutila complex. Figure S3. Hierarchically structure analysis.

Additional file 4: Figure S4.

Statistics obtained for the selection of bioclimatic variables (19 variables from WorldClim) that were used on Ecological Niche Modelling projections, for the whole complex (Amazilia rutila), and each geographic group (Pacific, Chiapas, Yucatán_CA).

Additional file 5: Figure S5.

Ecological Niche Modelling for Amazilia rutila under different past scenarios: LIG (Last Inter Glacial), LGM (Last Glacial Maxima), and MH (Mid Holocene) periods. Figure S6. Ecological Niche Modelling for each geographic group of Amazilia rutila under different past scenarios: LIG (Last Inter Glacial), LGM (Last Glacial Maxima), and MH (Mid Holocene) periods.

Additional file 6: Figure S7.

Bayesian skyline plots (BSP) for geographic groups (Pacific, Chiapas, Yucatán_CA), and for the A. rutila complex.

Rights and permissions

Open Access This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if changes were made. The images or other third party material in this article are included in the article's Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article's Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this licence, visit http://creativecommons.org/licenses/by/4.0/. The Creative Commons Public Domain Dedication waiver (http://creativecommons.org/publicdomain/zero/1.0/) applies to the data made available in this article, unless otherwise stated in a credit line to the data.

Reprints and Permissions

About this article

Verify currency and authenticity via CrossMark

Cite this article

Vázquez-López, M., Córtes-Rodríguez, N., Robles-Bello, S.M. et al. Phylogeography and morphometric variation in the Cinnamon Hummingbird complex: Amazilia rutila (Aves: Trochilidae). Avian Res 12, 61 (2021). https://doi.org/10.1186/s40657-021-00295-0

Download citation

Keywords

  • Amazilia rutila
  • Genetic structure
  • Genetic variation
  • Phylogeography
  • Tres Marías Islands
  • Trochilidae
  • Tropical deciduous forest