Get 20M+ Full-Text Papers For Less Than $1.50/day. Start a 14-Day Trial for You or Your Team.

Learn More →

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

Phylogeography and morphometric variation in the Cinnamon Hummingbird complex: Amazilia rutila... 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, STRU CTU RE, F , AMOVA), Bayesian and Maximum Likelihood phylogenetic analyses, and divergence time estimates. In order to ST 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 projec‑ tions onto past scenarios revealed no drastic changes in suitable conditions, but revealed some possible refuges. Mor‑ phometric 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 geo ‑ graphic isolation (limited gene flow), differences in current environmental conditions, and historical habitat fragmen‑ tation promoted by past events (Pleistocene refugia). Four well‑ defined clades comprise the A. rutila complex, and we *Correspondence: behb@ciencias.unam.mx Museo de Zoología, Facultad de Ciencias, Universidad Nacional Autónoma de México, Apartado Postal, 04510 Mexico City, Mexico Full list of author information is available at the end of the article © The Author(s) 2021. 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:// creat iveco mmons. org/ licen ses/ by/4. 0/. The Creative Commons Public Domain Dedication waiver (http:// creat iveco mmons. org/ publi cdoma in/ zero/1. 0/) applies to the data made available in this article, unless otherwise stated in a credit line to the data. Vázquez‑López et al. Avian Res (2021) 12:61 Page 2 of 17 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. Keywords: Amazilia rutila, Genetic structure, Genetic variation, Phylogeography, Tres Marías Islands, Trochilidae, Tropical deciduous forest Background Mexico have been considered to have high potential for The Mesoamerican dominion is the area where the diversification (García-Deras et  al. 2008; González et  al. Neartic and Neotropics overlap, and includes most of the 2011; Smith et  al. 2011; Ramírez-Barrera et  al. 2018; Mexican and Central American lowlands as well as the Hernández-Baños et  al. 2020). Furthermore, a study Mexican Transition Zone (Morrone 2014). This biogeo - based on molecular data showed that species distributed graphic dominion and surrounding areas are well known in the Neotropical lowlands have a low and constant rate for possessing high levels of species richness, which are of diversification over time; extrapolating those rates of thought to have originated from its complex topography, diversification to the present leads to a greater number wide range of environmental conditions, and climatic of lineages than is currently described (Weir 2006). The history (Savage 1966). One of the most distinctive habi- number of species present in the tropics is thus not prop- tats found in the Mesoamerican dominion is the tropical erly reflected in the currently accepted taxonomy. deciduous forest, which has a high number of endemic Hummingbirds belong to one of the most diverse bird species (Gordon and Ornelas 2000). families (Trochilidae: 361 spp.; Gill et al. 2021), with spe- Tropical deciduous forest has the driest conditions cialized nectarivorous feeding mode and high rates of of the tropical forests, and precipitation is strongly sea- speciation. This has resulted in a huge diversity of bill sonal (Rzedowski 1978; Meave et al. 2012). In addition, it morphologies and color patterns among nine highly sup- is considered one of the most threatened habitats in the ported groups (Bleiweiss et al. 1997; McGuire et al. 2014). Neotropics due to land use change, habitat fragmenta- In this study, we examined the evolutionary history of the tion, and human population growth (Lerdau et  al. 1991; Cinnamon Hummingbird (Amazilia rutila), a medium- Koleff et al. 2012). The diversity of the tropical deciduous sized hummingbird that inhabits tropical deciduous for- forest has been attributed to historical factors such as the ests ranging from 0 to 1600  m of elevation (Howell and climatic fluctuations of the Pleistocene, in which colder Webb 1995). This species is mainly distributed among and drier periods promoted its fragmentation, followed three biogeographic provinces: the Pacific lowlands, by periods of expansion with warmer and more humid Yucatán Peninsula, and Mosquito, based on the region- conditions (Werneck et al. 2011). alization of the Neotropical region (Smith 1941; Ryan Some studies of bird species distributed in Mesoameri- 1963; West 1964; Morrone 2014). It is a resident species can deciduous forests have shown high levels of genetic from Sinaloa in northwestern Mexico, south along the and morphological variation associated with historic Pacific Coast (including the Tres Marías Islands) to Costa and demographic events (e.g. Arbeláez-Cortés et  al. Rica, and along the coastal plain of the Yucatán Peninsula 2014), while others have focused on allopatric popula- south to Belize (including “offshore cays”) and north - tions, where long-term isolation, moderate or low gene eastern Nicaragua (Howell and Webb 1995). Based on flow, presence of geographic barriers, and historical morphology and allopatric distributions, four subspecies events explained the overall geographic structure (e.g. are recognized (Billerman et  al. 2020; Gill et  al. 2021). Smith et  al. 2011). One clear example of allopatric dif- A. r. rutila (DeLattre 1843) is distributed in western and ferentiation is the taxa inhabiting the Tres Marías Islands southern Mexico (from Jalisco to Oaxaca) and has green (Nayarit, Mexico), where there are several endemic plumage on the upperparts, cinnamon underparts and a subspecies (Cortés-Rodríguez et  al. 2008; Smith et  al. rufous tail. A. r. diluta (van Rossem 1938) inhabits north- 2011; Montaño-Rendón et  al. 2015). These subspecies western Mexico (from Sinaloa to Nayarit); its upper parts share the common feature of clear morphological dif- are more golden-bronze and underparts are paler and ferentiation, mainly in body size, from their continental less reddish than rutila. A. r. corallirostris (Bourcier and counterparts, representing independent evolutionary Mulsant 1846) is distributed from south and southeast- lineages with defined genetic structure (Ortíz-Ramírez ern Mexico (from Chiapas to the Yucatán Peninsula) to et  al. 2018). In addition to the Tres Marías Islands, the Costa Rica and is much more deeply colored than rutila. Yucatán Peninsula and the Coastal Plain of southwestern Finally, located within the Pacific lowlands province, the V ázquez‑López et al. Avian Res (2021) 12:61 Page 3 of 17 subspecies A. r. graysoni (Lawrence 1866) inhabits the the Washington University Genomics Center. All new Tres Marías Islands; it is darker in color and larger than sequences have been deposited in GenBank (Accession rutila (Ridgway 1901). In general, the systematics of A. numbers MZ998668-MZ998740, OK624605-OK624657, rutila are unresolved and controversial. Weller (1999) OK614044-OK614080, OK618334-OK618361). considered that A. rutila forms a superspecies with A. Sequences were proofread using Sequencher v.4.8 yucatanensis (Buff-bellied Hummingbird), and in a mul - (Gene Codes Corporation 2007), and aligned with the tilocus molecular study, the genus Amazilia formed a Clustal W function in BioEdit (Thompson et  al. 1994; polyphyletic group nested with the genera Hylocharis, Hall 1999). We corroborated the mitochondrial origin of Chrysuronia, Lepidopyga, and Damophila (McGuire our sequences using at least two of the following meth- et al. 2007). ods: amplifying/sequencing overlapping gene segments, The aims of this study were to: (1) analyze the genetic amplifying/sequencing one region with different primer variation and morphometric differentiation across the sets, and sequencing both DNA strands. We found no geographic distribution of A. rutila, (2) reconstruct phy- evidence of contamination of mtDNA sequences. logenetic relationships and past scenarios throughout the geographic range to analyze divergence times and demo- Phylogenetic analyses, genetic structure and neutrality graphic changes over time and space, and (3) propose tests an evolutionary history and taxonomic hypothesis of We used the Bayesian approach implemented in PHASE intraspecific lineages. We expected to find high levels of v.2.1 (Stephens et  al. 2001; Stephens and Scheet 2005) genetic structure associated with the isolation of allopat- to reconstruct haplotypes and avoid heterozygotes in ric populations and environmental conditions over con- nuclear sequences, selecting the pairs of haplotypes with tinuous ranges. a posterior probability higher than 0.90. We obtained the models of evolution for each gene using Partition- Methods Finder (Lanfear et  al. 2017), using the following param- Sampling, PCR and sequencing eters: linked branched lengths, greedy search algorithm We obtained 84 tissue samples from the A. rutila com- (Lanfear et  al. 2012), and the Bayesian Information Cri- plex and one sample from A. yucatanensis as the terion (BIC) for model selection. The phylogenetic anal - outgroup (see Additional file  1: Tables S1, S2). We sup- yses were conducted using each locus separately and plemented our sampling with GenBank sequences from using the concatenated matrix including nuclear and A. rutila, A. yucatanensis and A. tzacatl (McGuire et  al. mitochondrial loci (multilocus). We conducted phylo- 2014). Tissues samples were provided by the collection genetic analysis under Bayesian Inference (BI), using of the Museo de Zoología (MZFC) “Alfonso L. Herrera” MrBayes 2.0 (Huelsenbeck and Ronquist 2002) in CIP- (Universidad Nacional Autónoma de México), the Burke RES Science Gateway (Miller et  al. 2010). Each analysis Museum (University of Washington), and the Biodiver- consisted of four independent chains, random starting sity Institute (University of Kansas). trees, and uniform prior distribution of parameters. The DNA was extracted using the Qiagen DNAEasy kit, chains were run for 30 million generations, sampling following the manufacturer’s protocols. We ampli- trees every 1000 generations. The asymptote was deter - fied and sequenced two mitochondrial markers: par - mined visually, 5000 burn-in trees were discarded, and tial and complete NADH dehydrogenase subunit 2 the remaining trees from the plateau phase were used to (ND2: 605–1041  bp), and the full length of cytochrome estimate Bayesian posterior probabilities. We considered C oxidase subunit I gen (COI: 568  bp). We also ampli- that clades were strongly supported if they were present fied two nuclear markers: the regions between exons 4 in ~ 95% of the sampled trees (Huelsenbeck and Ronquist and 5 of the Muscle Skeletal Receptor Tyrosine Kinase 2002; Wilcox et  al. 2002). Also, a Maximum Likelihood gene (MUSK: 628  bp), and a segment comprising the analysis was conducted using RaxML v.8.0.0 (Stamakis end of exon 6 and the beginning of exon 8 of the Orni- 2014) in CIPRES Science Gateway (Miller et  al. 2010), thine Decarboxylase gene (ODC: 571). We used primers using separate partitions for each locus (see “Results”), L5219, H6313, L5758 and H5766 for the amplification of with 100 bootstrap replicates. ND2 (Sorenson et al. 1999); primers F1-COI and R2-COI We defined four different groups, corresponding to for COI (Chaves et  al. 2008); ODC2-F and ODC2-R for allopatric populations from sampled localities of A. ODC (McGuire et  al. 2007); and MUSK R3 and MUSK rutila (see Fig.  1): (1) Tres Marías Islands, (2) Pacific F3 for the amplification of MUSK (McGuire et  al. coast, (3) Chiapas, and (4) Yucatán Peninsula and Cen- 2007). All PCR products were confirmed on 1% aga - tral America (hereafter: Tres Marías, Pacific, Chia - rose gel, and sequencing reactions were performed by pas and Yucatán_CA). For the comparison between Vázquez‑López et al. Avian Res (2021) 12:61 Page 4 of 17 Fig. 1 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 populations, we specified geographic groups based on Maximum Likelihood model, with 100 bootstrap repli- the main clades recovered from phylogenetic analysis cates (Stecher et al. 2020). (concatenated-multilocus), and tested for substructure in Pacific and Yucatán_CA groups (see “Results”). A Divergence time estimates median-joining network was constructed to visualize We estimated divergence times with the concatenated the structure of populations, the number of haplotypes, matrix using Beast v.2.6.2 (Bouckaert et  al. 2019). Dif- their frequencies, and the relationships among individ- ferent partitions by gene were defined based on Pacheco uals using the program Network 4.5.1.0 (Bandelt et  al. et  al. (2011) for ND2 and COI, Lerner et  al. 2011 for 1999). We tested locus neutrality with Fu’s Fs (Fu 1997) ODC, and Ellegren (2007) for MUSK. We employed a and Tajima’s D (Tajima 1989) metrics in DnaSP v.5 strict clock and a constant population model as a tree (Librado and Rozas 2009), and obtained the summary prior. We chose the strict clock model empirically based statistics for each population. We tested the differen - on data fit, and the constant population model based tiation hypothesis using an analysis of molecular vari- on our demographic reconstruction (see “Results”). The ance AMOVA and the F fixation index using Arlequin analysis was run for 200 million generations, sampling ST 3.1 (Excoffier and Schneider 2006). These analyses are every 1000. We used Tracer v.1.7 (Rambaut et  al. 2018) useful for observing population structure, estimating to ensure adequate Effective Sample Size (ESS > 200), population differentiation directly from molecular data, which was reached for most statistics, with the exception and testing the hypothesis of differentiation (Excoffier of prior and probability statistics. With TreeAnnotator v. et  al. 1992). Additionally, we used STRU CTU RE 2.3.2 2.6.2 (Rambaut and Drummond 2013) we discarded the (Pritchard et  al. 2000) for population assignments first 25% or trees as burn-in and produced a maximum under an admixed and LocPrior model (Hubisz et  al. clade credibility tree with 95% highest probability density 2009), with 10,000 iterations of burn-in and 20,000 intervals. The final tree was visualized in FigTree v.1.4.2 MCMC (Markov chain Monte Carlo), for 20 replicates (Rambaut 2014). of each K value (from 1 to 6). To evaluate STRU CTU RE results, we used the Evano method to choose K value Morphometric analysis by ∆K (Evanno et  al. 2005), as implemented on the We took linear measurements of five morphological traits Structure Harvester website (Earl and vonHoldt 2012). from 210 study skins housed at the American Museum of Genetic distances were obtained with MEGA under Natural History (AMNH) and at the Museo de Zoología V ázquez‑López et al. Avian Res (2021) 12:61 Page 5 of 17 “Alfonso L. Herrera” (MZFC). We measured wing chord account the logistic threshold from maximum training (WC) and tail length (TL) using an Avinet calibrated sensitivity plus specificity. ruler to the nearest 0.5  mm; also, we measured culmen We obtained Bayesian Skyline Plots (BSP) using Beast length (CL), beak width (BW), and beak height (BH) to v.2.6.3 (Bouckaert et al. 2019) to infer population dynam- the nearest 0.1 mm with a Mitutoyo digital caliper. ics and changes in demography over time. We used the We assigned each specimen to geographic groups ND2 dataset and performed four different runs: three of defined based on the haplotype network (see “Results”; them based on geographic groups (Pacific, Chiapas and Fig.  1). We performed Mann–Whitney U tests between Yucatán_CA), and one based on the whole A. rutila com- sexes within each geographic group, with Bonferroni cor- plex. The Tres Marías group was not included in an inde - rection to evaluate sexual dimorphism. We also carried pendent run because of the lack of informative sites for out Kruskall-Wallis test with “geographic group” as the ND2. We set a GTR + substitution model, a strict molec- grouping variable to assess geographic variation in the ular clock, a Coalescent Bayesian Skyline as the tree measured characters. Finally, we carried out a PCA to prior, and a Piecewise-constant skyline model (Drum- visually examine the structure of variation in multivari- mond et  al. 2005) with five groups. We used a mutation ate space. All statistical analyses were performed in R (R rate of 0.0145 substitutions per site per lineage per mil- Core Team 2020). lion years for ND2, following Lerner et al. (2011), and ran each analysis through 20 million generations. Ecological Niche Modelling under past scenarios and demography Results To infer suitability of environmental conditions under Phylogenetic analysis past scenarios, we performed projections using Eco- The model of sequence evolution selected for the logical Niche Models. Presence records of A. rutila were mtDNA data set was GTR + G (nst = 6, rates = gamma). downloaded from GBIF (Global Biodiversity Information The nucleotide composition of mtDNA was as follows: Facility) and accessed from R via rgbif (https:// github. T = 23.9%, C = 34.5%, A = 30.8%, G = 10.6% for ND2; and com/ ropen sci/ rgbif; taxon key: 2476417; https:// doi. org/ T = 25.9%, C = 15.7%, A = 26.4%, G = 31.8% for COI. The 10. 15468/ dl. 7k98v7). These records were cleaned over parsimony informative sites for ND2 were 72 (605  bp), multiple steps using the following R packages: Coordi- and 38 for COI (568 bp). For nuclear markers, the mod- nateCleaner (Zizka et  al. 2019), nichetoolbox (Osorio- els of sequence evolution were: F81 (nst = 1) for MUSK, Olvera et  al. 2016), dplyr (Wickham et  al. 2020), and and GTR + I + G (nst = 6, rates = invgamma) for ODC. raster (Hijmans 2020). We used the sets of bioclimatic The nucleotide composition was: T = 29.6%, C = 19.4%, layers (current and past) from WorldClim (Hijmans et al. A = 29.4%, G = 21.4% for MUSK; and T = 36.4%, 2005; Braconnot et al. 2007; Booth et al. 2014). Our mod- C = 17.3%, A = 24.7%, G = 21.4% for ODC. The parsi - els for the past were projected onto the Last InterGlacial mony informative sites for MUSK were 19 (628 bp), and (LIG: 140–120 kya), the Last Glacial Maximum (LGM: 21 83 for ODC (571 bp). kya), and the Mid Holocene (MH: 6 kya) scenarios. We According to the topology from the concatenated data performed these projections using the continental phy- set (mitochondrial and nuclear genes) under Bayesian logroups (Pacific coast, Chiapas, and Yucatán and Cen - Inference and Maximum Likelihood, A. rutila forms a tral America); the Tres Marías group was not included monophyletic group with respect to A. yucatanensis and in these analyses due to the low number of occurrence A. tzacatl (see Fig. 2). Within A. rutila, three major line- records (n = 3). The selection of bioclimatic variables was ages can be identified: Pacific coast, Chiapas, and Yuca - based primarily on a principal component analysis, where tán Peninsula and Central America. There was also a we considered the most important variables of the first well-supported clade (posterior probability = 0.93) within four principal components (from one to three variables), the Pacific phylogroup, which contained all individuals ensuring that no intercorrelated variables were selected from the Tres Marías Islands. Within the Yucatán_CA (Pearson correlation r < 0.75, Additional file  2: Fig. S1). phylogroup, four of the five individuals from Nicaragua To define the area of accessibility for the species, we used grouped together into a well-supported clade (pp = 0.99; the ellipsenm R package (Cobos et  al. 2020) to delimit bootstrap = 99). The individual phylogenies for mito - the area that contained occurrence points across the dis- chondrial loci (ND2 and COI) recovered three major lin- tribution of the species with a polygon representing a eages (Pacific, Chiapas, and Yucatán_CA) (see Additional buffer area of 75 km (“M” area). Models were performed file  2: Fig. S1). However, individual nuclear locus phylog- in Maxent v.3.4.1 (Phillips et al. 2021) with regularization enies did not recover any geographic structure (see Addi- multiplier = 2 and with 10 replicates of cross-validation tional file 2: Fig. S1). with no extrapolation. To binarize outputs, we took into Vázquez‑López et al. Avian Res (2021) 12:61 Page 6 of 17 Fig. 2 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 group was separated from the rest by 57 mutational steps Deviations from neutrality were rejected for mitochon- (Fig. 1). For COI, we found 11 haplotypes in total: 4 cor- drial loci (ND2 and COI) and nuclear loci, except for responded to the Pacific group, 3 to Chiapas, 3 to Yuca - the ODC locus (see Table  1). We found 21 haplotypes tán-Central America, and 1 to the Tres Marías Islands for ND2: 9 corresponded to the Pacific group, 9 to the (shared with 2 individuals from the Pacific group from Yucatán-Central America group, 2 to Chiapas, and 1 to Guerrero). The Yucatán-Central America group was the Tres Marías Islands. The Yucatán-Central America V ázquez‑López et al. Avian Res (2021) 12:61 Page 7 of 17 Table 1 Neutrality test by locus Group Tajimas’ D Fu and Li’s D ND2 COI MUSK ODC ND2 COI MUSK ODC ‡ ‡ ‡ ** ‡ ‡ ‡ * All − 0.14 − 0.57 − 1.45 − 2.3 − 1.47 − 2.39 0.07 − 3.91 ‡ ‡ ‡ ‡ ‡ ‡ ‡ * Pacific − 1.28 2.38 − 0.62 − 2.13 0.15 1.06 − 0.82 − 4.69 ‡ ‡ ‡ ‡ Tres Marías 0.0 0.0 − 1.13 − 1.23 0.0 0.0 − 1.15 − 1.27 ‡ ‡ ‡ ‡ ‡ ‡ Chiapas 0.0 − 0.17 − 0.31 0.24 0.0 1.44 0.13 0.71 ‡ ‡ ‡ * ‡ ‡ ‡ * Yucatán_CA − 0.57 − 1.08 0.70 − 1.93 − 0.43 − 1.31 1.149 − 2.56 ‡ ‡ Nicaragua − 0.75 ND ND ND − 0.75 ND ND ND ‡ * ** *** P > 0.5, P < 0.05, P < 0.01, P < 0.001 The mtDNA genetic distances among groups showed Table 2 mtDNA summary statistics an overall genetic distance of 5.24% for the Pacific group, Geographic group N H Hd Pi 5.69% for Chiapas, and 7.52% for the Yucatán_CA group Pacific 25 5 0.51 0.00196 (Table  3). The pairwise F fixation indices were all sig - ST Tres Marías 6 1 0.0 0.0 nificant and indicated strong population structure of Chiapas 5 1 0.0 0.0 geographic groups for mtDNA data (Table  4). We did Yucatán_CA 14 5 0.59 0.00212 not find strong differentiation between the Pacific and Nicaragua 4 3 0.83 0.00145 its subgroup from Guerrero. However, the Guerrero sub- All 50 12 0.83 0.03037 group showed a more robust differentiation with respect to the Tres Marías group than to the rest of the Pacific N number of sequences, H number of haplotypes, Hd haplotype diversity, Pi nucleotide diversity subgroup. Substructure was detected within the Yucatán_ CA group, with some Nicaraguan samples (UWBM-6911, UWBM-68991, UWBM-69261, and UWBM-69388) sep- Table 3 mtDNA pairwise genetic distances among geographic arated from the rest of the group, with an F of 0.78549. ST groups in the Amazilia rutila complex AMOVA analyses revealed that most of the genetic vari- Pacific Tres Marías Chiapas Yucatán_CA ation was among geographic groups (96.26%), and there was little variation within populations (3.74%; Table  3). Pacific The F value from general AMOVA results indicated ST Tres Marías 0.00629 strong genetic structure and differentiation among popu - Chiapas 0.03870 0.03766 lations (F = 0.96; Table 5). ST Yucatán_CA 0.07524 0.07389 0.07487 The STRU CTU RE analyses detected two genetic groups (K = 2; first-level analysis; See Additional file  3: Figs. S2, S3; Fig.  3). One group comprised all individuals separated from the rest by 15 mutational steps (Fig.  1). from Sinaloa to Chiapas, including Tres Marías. The sec - The Pacific coast and Yucatán-Central America groups ond group included the Yucatán Peninsula and Central had the highest haplotype and nucleotide diversity (see American individuals. To evaluate potential substructure Table 2). within these groups, we ran a hierarchically structured Table 4 Pairwise F indices (below diagonal, P ≤ 0.05) among populations in the Amazilia rutila complex ST Pacific Guerrero Tres Marías Chiapas Yucatán_CA Nicaragua Pacific ‒ Guerrero 0.33425 ‒ Tres Marías 0.55131 1 ‒ Chiapas 0.93300 1 1 ‒ Yucatán_CA 0.96654 0.99167 0.99209 0.99237 ‒ Nicaragua 0.96132 0.98693 0.98831 0.98791 0.78549 ‒ Geographic groups taking into account for this analysis: Pacific, Guerrero (clade within Pacific: CHIS529, CHIS533, CHIS534, CHIS542, CHIS549), Tres Marías, Chiapas, Yucatán_CA, Nicaragua (clade within Yucatán_CA: UWBM6911, UWBM68991, UWBM69261, UWBM69388) Vázquez‑López et al. Avian Res (2021) 12:61 Page 8 of 17 analysis. Individual structure analyses for the first group before present (Myr BP) (95% highest posterior density revealed two genetic groups (K = 2), with no ancestral [HPD] 8.04‒10.87; Fig.  4), while the A. rutila complex mixing between them. One group was composed of indi- root was dated 7.05 Myr ago (95% HPD 5.88‒8.29). viduals from Chiapas and the second group recovered This estimate also corresponds to the split between the all of the Pacific clade including Tres Marías individuals. two main clades: (1) Pacific + Tres Marías and Chia- The substructure within the Yucatán_CA group recov - pas, and (2) Yucatán_CA. The Pacific clade and Chi- ered two groups (K = 2). One group included samples apas split around 3.99 Myr BP (95% HPD 3.16‒4.84). from western Nicaragua and three samples from Yuca- The Pacific Clade node origin was 1.65 Myr BP (95% tán. The second group included one sample from western HPD 0.81‒1.81  Ma). According to these estimates, A. Nicaragua, one sample from Guatemala and the rest of rutila arrived at the Tres Marías Islands 0.2 Myr BP the Yucatán Peninsula individuals. (95% HPD 0.6‒0.44  Ma). The age of the continental subclades within the Pacific clade suggests that Pleis- Divergence time estimates tocene climatic changes may have caused contractions Divergence times estimated that the A. rutila com- in the continental populations of A. rutila and resulted plex split from its sister clade around 9.4 million years in incipient substructure, however although the Tres Table 5 mtDNA AMOVA comparing geographic groups in the Amazilia rutila complex Source of variation df Sum of squares Variance components Percentage of variation (%) Among groups 5 564.198 14.31 96.26 Within groups 47 26.141 0.55 3.74 Total 52 Fixation index F : 0.96 ST Significance test P < 0.0001 See geographic groups in Table 4 Fig. 3 STRU CTU RE 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) V ázquez‑López et al. Avian Res (2021) 12:61 Page 9 of 17 Fig. 4 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 Vázquez‑López et al. Avian Res (2021) 12:61 Page 10 of 17 Marías Islands is well supported, other clades are not (Annual Precipitation), BIO15 (Precipitation Season- supported. The Chiapas clade arose 0.76 Myr BP (95% ality), and BIO18 (Precipitation of Warmest Quarter). HPD 0.43‒1.15 Ma). The average evaluation metrics for the model results The crown node for the Yucatán_CA clade emerged were: training AUC = 0.7598 (St.D v ± 0.005), and a 1.98 Myr BP (95% HPD 1.42‒2.6  Ma). Internal clades threshold value of 0.4178 (Maximum training sensi- in the Yucatán_CA group were detected by Beast anal- tivity plus specificity Logistic threshold). Projections ysis, but they had weak node support, with the excep- onto past scenarios for the whole complex are shown tion of the Nicaragua clade, which had intermediate in Additional file  5: Figs. S5, S6. The suitability of envi - support (0.7 posterior probability) and a date of 0.12 ronmental conditions has changed somewhat over Myr BP (95% HPD 0.03‒0.25 Ma). time, with subtle reductions from LIG to LGM and posterior arrangements after LGM through the MH period. Projections during LGM showed discontinu- Morphometrics ous distribution, which contrasts with the current con- We found minor sexual dimorphism in this species. tinuous distribution for the species (Pacific coast and Females were slightly but statistically significantly Central America). Independent analyses were run for smaller than males in two body size variables we each geographic group, in which the first two princi - measured (WC, W = 2196.5, p < 0.001; TL , W = 2194, pal components explained most of the variation (see p < 0.001). However, the effect sizes were very small, so Additional file  5: Figs. S5, S6). The bioclimatic vari - we carried out subsequent tests of geographic varia- ables selected for the Pacific group were BIO1, BIO3, tion with both sexes pooled together. BIO4, BIO7, BIO13, BIO17 and BIO19; for the Chiapas There were statistically significant differences group were BIO4, BIO5, BIO12, BIO13, BIO14, BIO15, among geographic groups (Kruskall-Wallis tests WC: BIO16, and BIO18; and for the Yucatán_CA group 2 2 χ = 78.86, p < 0.001; TL: χ = 73.09, p < 0.001; BH: were BIO3, BIO4, BIO7, BIO11, and BIO13. The result - 4 4 2 2 χ = 47.45, p < 0.001; C L: χ = 62.61, p < 0.001; BW: ing metrics were as follows: Pacific AUC = 0.8046 (St. 4 4 χ = 90.30, p < 0.001). The Tres Marías phylogroup had Dv ± 0.0473), threshold = 0.358; Chiapas AUC = 0.8184 larger values than the three continental groups for all (St.Dv ± 0.1119), threshold = 0.4255; Yucatán_C A variables. This was consistent with the results of our AUC = 0.6515 (St.D v ± 0.0470), threshold = 0.4515. PCA, which showed that the continental phylogroups Projections for Pacific group predicted some areas largely overlapped in morphological measurements, where the Chiapas and Yucatán_CA groups are cur- but the Tres Marías group was significantly larger. rently distributed (see Additional file  5: Figs. S5, S6). Principal component 1 explained 74.75% of the cumu- Also, for the projections in the Yucatán_CA group, lative variance and roughly corresponds to overall suitable areas were predicted where the Pacific and body size, with large contributions from Wing Chord Chiapas groups now occur. The most conservative pro - (WC) and Tail Length (TL) (Fig. 5). jections for suitable habitats were recovered when the Chiapas group was modeled. Past projections revealed Ecological Niche modelling onto past scenarios a reduction of suitability during LGM and an increase and demography during MH periods in all cases. During the LIG period, We obtained a total of 445 presence records for A. the Pacific showed continuous suitable habitat along rutila after filtering steps. The number of occurrence the Pacific coast to Central America. In LIG projections points for each group was: (1) Pacific: 144 records, for Chiapas and Yucatán_CA, no suitable conditions (2) Chiapas: 20 records, and (3) Yucatán_CA: 281 were recovered (see Additional file 5 : Figs. S5, S6). records. The results of the PCA of climate variables The BSP results showed that population dynamics and for the whole complex showed that the first four prin - demography in the A. rutila complex have been generally cipal components explained most of the variation constant over time with a notable population reduction (88%, see Additional file  4: Fig. S4). According to these close to the present. Effective population size patterns results and correlation values, we selected the follow- showed a stationary trend in the geographic groups ing bioclimatic variables: BIO5 (Max Temperature of analyzed (Pacific, Chiapas, Yucatán_CA), with recent Warmest Month), BIO7 (Temperature Annual Range), BIO11 (Mean Temperature of Coldest Quarter), BIO12 (See figure on next page.) Fig. 5 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 V ázquez‑López et al. Avian Res (2021) 12:61 Page 11 of 17 Fig. 5 (See legend on previous page.) Vázquez‑López et al. Avian Res (2021) 12:61 Page 12 of 17 population reductions in the Pacific and Chiapas groups, Pleistocene, Middle American biological diversity was and a reduction followed by a slight increase towards the influenced by several major processes. One is the uplift of present in the Yucatán_CA group (see Additional file  5: major mountain chains, which generated geographic bar- Figs. S5, S6). riers and created open niches in high-altitude and low- altitude habitats that favored diversification (McCormack Discussion et  al. 2008). Biodiversity was also influenced by the cli - Genetic variation and phylogeographic pattern matic changes of Pleistocene climatic oscillations, which We found both deep and shallow genetic differentiation created a mosaic of habitats that resulted in isolated throughout the geographic distribution of A. rutila. The populations. lineages corresponding to the Pacific clade (pp = 0.99), The Tres Marías phylogroup is embedded in the larger Chiapas (pp = 0.99) and the Yucatán Peninsula-Central Pacific clade with a node origin dated 0.20 Myr BP (95% America (pp = 0.99) show deep structuring according to HPD 0.6‒0.44 Myr BP); this incomplete separation is genetic distances, F index, AMOVA, haplotype net- a pattern that has been reported in other birds at early ST work, multilocus phylogeny and mtDNA phylogenies. stages of speciation (e.g. Cortés-Rodríguez et  al. 2008). Within the Pacific clade, we can distinguish clear struc - STRU CTU RE analysis did not recover substructure for ture for Tres Marías Islands individuals with a monophy- Tres Marías individuals or any Pacific samples. However, letic clade (pp = 0.94; geographically isolated), F index, the ∆K selection method has been previously shown to ST one haplotype for ND2 locus, and one haplotype for COI underestimate population structure (Janes et  al. 2017). (shared with two samples from Guerrero); however, the Results of the AMOVA and F index suggest incipient ST structure analysis shows mixed ancestries between the population structure, and haplotype networks showed Tres Marías Islands and Continental samples. Similarly, a clear differentiation at the population level. Addition - a multilocus phylogeny suggests a shallow substructure ally, the Tres Marías group was the most morphologi- in the Yucatán Peninsula-Central America clade for West cally distinct compared to the other groups. The Pacific Nicaraguan samples (pp = 0.99), which share ancestry clade is represented as a polytomy in the resulting phy- with some samples from Yucatán. logenetic tree, and some clades are recovered within Despite minor differences in median-joining networks it but with no node support (see Additional file  2: Fig. based on the ND2 versus COI markers (Fig.  1), the geo- S1). As mentioned, the single valid clade is composed graphic groups had no shared haplotypes, with the of the Tres Marías Islands individuals. There are sev - exception of one haplotype that was shared between the eral examples of bird taxa in this area (The Pacific coast Tres Marías group and three individuals from the Pacific of Mexico and Tres Marías Islands) which show similar group (COI marker). The information from nuclear genes patterns of geographic variation and population struc- that was included to construct the multilocus phyloge- ture (e.g. Habia rubica: Ramírez-Barrera et al. 2018; Pha- netic tree (see Additional file  2: Fig. S1) did not weaken ethornis longirostris‒P. mexicanus: Arbeláez-Cortés and the pattern found using mtDNA; the multilocus phy- Navarro-Sigüenza 2013; Vireo hypochryseus: Arbeláez- logeny maintained the same geographic structure with Cortés et  al. 2014). Those studies found that the main high node support values (posterior probabilities and forces acting as drivers of divergence were the presence bootstrap values), even of the retained ancestral poly- of geographic and ecological barriers, as well as historical morphisms, favored by nuclear signal and its reduced habitat fragmentation. This is not the case for A. rutila substitution rates compared to mtDNA (Moore 1995). along the Pacific coast, since no structured pattern was A. rutila belongs to the third youngest subclade of nine found. There are two possible explanations for this lack Trochilidae subfamilies, known as the Emeralds (Blei- of geographic substructure. First, small sample sizes may weiss et al. 1997; McGuire et al. 2007, 2014). The tempo cause some haplotypes to be missing, since we found and mode of evolution in the hummingbird higher clades many intermediate haplotypes. Another explanation is has been considered to be the product of a series of his- an incipient substructure produced by short periods of torical events including colonizations, extinctions, recol- isolation during the Pleistocene. The age of the Pacific onizations, and the influence of mountain chain uplift, clade is 1.65 Myr BP (95% HPD 0.81‒1.81), and the age resulting in the radiation of highland species and the ranges estimated for subclades are 1.32 to 0.20 Myr BP, establishment of lowland taxa (Bleiweiss 1998; McGuire which corresponds with Pleistocene cycles. Thus, these et al. 2014). Therefore, timing is crucial to understanding mixed haplotypes could be the result of short periods of phylogeographic patterns in A. rutila. This complex split isolation during the LGM with secondary contact at the from its sister group 7.05 Myr BP (95% HPD 5.88‒8.29 LIG and the present. This idea is supported by the LGM Myr BP), and the four main clades reported here arose Pacific projections (see Additional file  5: Figs. S5, S6) 3.99‒0.20 Myr BP. During the late Miocene, Pliocene and and BSP analysis (Additional file  6: Fig. S7). In contrast, V ázquez‑López et al. Avian Res (2021) 12:61 Page 13 of 17 STRU CTU RE hierarchical analysis did not show different shows a mixed ancestry between those samples and some groups within the Pacific Clade (Fig. 2). Yucatán individuals (Fig. 2). The Chiapas clade is well defined, supported by genetic The influence of Pleistocene climatic fluctuations distance (5.69%) and F values (0.97‒1.00). This region on processes of diversification and incipient specia - ST east of the Isthmus of Tehuantepec is mainly influenced tion have been reported for many birds (e.g. Euphonia by historical events between the middle and late Plio- affinis : Vázquez-López et  al. 2020), including humming- cene and marine transgressions and regressions during birds (e.g. Eugenes fulgens: Zamudio-Beltrán et al. 2020). the late Pliocene (Barber and Klicka 2010; Peterson et al. These habitat fluctuations are well documented as fac - 2010). In general, the Isthmus of Tehuantepec represents tors that promote population contractions, limiting gene a geographic barrier for highland bird species (e.g. Zamu- flow (refugia hypothesis; Arango et  al. 2021), as well as dio-Beltrán and Hernández-Baños 2018), but not for low- expansions, resulting in zones of contact between distant land species. However, some phylogeographic studies of populations, with current signals of shared haplotypes lowland birds have described the influence of this low and polyphyletic patterns (e.g. Mesoamerican Amazilia altitude zone (e.g. Campylorhynchus rufinucha: Vázquez- complex: Jiménez and Ornelas 2016; White-chested Miranda et  al. 2009). According to our time-calibrated Amazilia complex: Rodríguez-Gómez and Ornelas tree, the Chiapas clade separated during Pleistocene cli- 2018). Diversifications occurring throughout the Mio - matic fluctuations ~ 760 kyr BP (95% HPD 0.43‒1.15 Ma), cene, when the estimates of divergence dates of the genus as there was a disruption of suitable conditions during the Amazilia are presumed ca. 19.85 Mya (i.e. Ornelas et al. Last Glacial Maxima across the Isthmus of Tehuantepec 2014), may have depended less on these climatic condi- (see Fig.  5). The Yucatán Peninsula and Central America tions. However, for events traced after this period, the clade also represents an independent entity, separated by rapid changes in climate cycling were crucial (Bleiweiss 57 and 15 mutational steps, respectively, in the ND2 and 1998). Furthermore, phylogeographic patterns of species COI trees (see Fig.  1), assuming low levels of gene flow, sharing (totally or partially) the geographic distribution with an origin 1.98 Myr BP (95% HPD 1.42‒2.6 Myr BP). of A. rutila show similarities (e.g. Turdus rufopalliatus: This genetic signal of isolation from the other groups Montaño-Rendón et  al. 2015; Ortíz-Ramírez et  al. 2018; may be influenced by the geographic distance and dif - de Silva et  al. 2020). In general, the diversification pro - ferences in environmental spaces. The Pacific and Chia - cesses have resulted from a variety of events, sharing the pas group are distributed along the tropical deciduous explanations of an intricate orographic history, climatic forest belt, which has marked dry and rainy seasons and fluctuations during the Pleistocene, and thus low signals has been repeatedly identified as a distinct biogeographi - of gene flow favored by the isolation of populations due cal province (i.e. Gordon and Ornelas 2000), while the to distance or historical refuges. Yucatán Peninsula has moderately pronounced dry and rainy seasons with relatively constant average tempera- Morphometrics ture throughout the year (Paynter 1955). Also, our ENM Our results showed little variation among the continental past projections suggest periods of habitat contractions, phylogroups of this species in the morphological char- during which populations from the Yucatán Peninsula acters we measured, despite A. rutila’s large geographic and Central America could have been closer from each range. However, the Tres Marías Islands phylogroup is other (LGM MIROC, see Additional file  5: Figs. S5, S6), larger than the continental phylogroups in four out of followed by periods of isolation. Although A. rutila is five measurements (Fig.  3). This is consistent with several allopatrically distributed in the Yucatán Peninsula and other avian species from Tres Marías (Grant 1965) and Central America, the studied individuals were placed is a clear example of the “Island Rule” (van Valen 1973; in a single group. This link between individuals close to Lomolino 1985; Clegg and Owens 2002; Boyer and Jetz the Yucatán Peninsula and Central America was also 2010, Benítez-López et  al. 2021). This shift in body size found in A. tzacatl (Miller et al. 2011), in which popula- has been attributed to changes in thermal ecology (Clegg tions in Tabasco (Mexico), Belize, Honduras, and Costa and Owens 2002), as the Tres Marías Islands are signifi - Rica shared several haplotypes despite the long distances cantly drier and warmer than the continental habitat of between them, and the most differentiated populations the species. However, dramatic divergence in body size were the southernmost populations in Panama. The can be the result of geographic isolation independent of periods of isolation mentioned above combined with environmental pressure (Seeholzer and Brumfield 2018). demographic decreases (see BSP results) could explain This divergence in body size with relatively similar beak the shallow substructure found in the West Nicaragua sizes suggest that beak morphology might be more tightly samples (pp = 0.99). However, our STRU CTU RE analysis constrained than body size due to this species’ feeding ecology. Vázquez‑López et al. Avian Res (2021) 12:61 Page 14 of 17 Implications for taxonomy name, but specimens from Huehuetán, Chiapas were We present a taxonomic proposal for the A. rutila com- described under A. cinnamomea saturata by Nelson plex based on the recovery of four independent clades in 1898, and this name was later considered a synonym with high support values corresponding to natural geo- of A. r. corallirostris (Ridgway 1901). As a consequence, graphic groups. the Chiapas lineage could be named A. saturata (Nelson Our results support the first differentiated clade, A. r. 1898). graysoni, distributed in the Tres Marías Islands (Nayarit, Mexico). This differentiation was also confirmed by our morphometric analyses, which showed that the Tres Conclusions Marías Islands population is the largest morphotype. A. The evolutionary hypothesis of diversification within this r. graysoni was first suggested as an independent species group is of geographic isolation (limited gene flow), dif - by Ridgway (1901), who described Amazilia graysoni ferences in current environmental conditions, and his- (Grayson’s Hummingbird) as a separate species similar torical habitat fragmentation promoted by past scenarios in coloration to A. rutila, but darker and much larger. (Pleistocene refuges). The A. rutila complex consists of Then following the alternative phylogenetic/evolution - four well-defined clades that, in our opinion, merit taxo - ary species taxonomy of the Mexican avifauna, Navarro- nomic reevaluation. Our data suggest that A. r. graysoni Sigüenza and Peterson (2004) proposed A. r. graysoni as (Tres Marías Islands) as well as A. r. rutila (Pacific coast) an independent evolutionary species, which is morpho- should be considered full species. The other two strongly logically distinct in size and coloration from the rest of supported clades are: (a) the Chiapas group (southern the subspecies within A. rutila. The proposal of elevating Mexico), and (b) the populations from Yucatán Peninsula this taxon to species level has been discussed recently and Central America. These clades belong to the A. rutila (de Silva et  al. 2020). The destruction of native habitat, corallirostris taxon as currently defined. We propose that the introduction of exotic species, and reduced popula- this group be split and renamed as two species: A. satu- tion sizes increase the vulnerability of birds on the Pacific rata for the Chiapas group and A. corallirostris for the Islands of Mexico and must be immediately addressed Yucatán Peninsula and Central America group. in conservation policies (Hahn et  al. 2012; de Silva et  al. 2020). Supplementary Information The second differentiated clade corresponds to the The online version contains supplementary material available at https:// doi. org/ 10. 1186/ s40657‑ 021‑ 00295‑0. group distributed along the Pacific coast. Within this second group, two subspecies have been described (A. r. Additional file 1: Table S1. Information of voucher numbers (ID), geo ‑ rutila and A. r. diluta), which, according to our results, graphic locations (longitude, latitude), and biological collections of tissue are not observable as separate groups. Therefore, we pro - samples used in this study. Table S2. Information of accession numbers pose that the Pacific coast phylogroup be represented by (AN) and geographic locations (locality, longitude, latitude) of GenBank’s sequences used in this study. A. r. rutila (the first named). The third phylogroup is distributed in Chiapas, Mex - Additional file 2: Figure S1. The Bayesian Inference and Maximum Likeli‑ hood of the mitochondrial DNA and nuclear DNA. ico (east of the Isthmus of Tehuantepec). This group Additional file 3: Figure S2. STRU CTU RE analysis of Amazilia rutila com‑ corresponds to A. r. corallirostris, which is distributed plex. Figure S3. Hierarchically structure analysis. from south and southeastern Mexico to Costa Rica. The Additional file 4: Figure S4. Statistics obtained for the selection of biocli‑ Chiapas phylogrouop is distributed in an important and matic variables (19 variables from WorldClim) that were used on Ecological well-studied main fauna region, considered an area of Niche Modelling projections, for the whole complex (Amazilia rutila), and each geographic group (Pacific, Chiapas, Yucatán_CA). endemism in western Mexico with one of the highest val- ues of avian species richness (Peterson and Navarro 2000; Additional file 5: Figure S5. Ecological Niche Modelling for Amazilia rutila under different past scenarios: LIG (Last Inter Glacial), LGM (Last Gla‑ García-Trejo and Navarro-Sigüenza 2004). cial Maxima), and MH (Mid Holocene) periods. Figure S6. Ecological Niche Finally, the fourth independent group corresponds to Modelling for each geographic group of Amazilia rutila under different the clade of the Yucatán Peninsula—Central America, the past scenarios: LIG (Last Inter Glacial), LGM (Last Glacial Maxima), and MH (Mid Holocene) periods. most genetically differentiated (7.52%), which also cor - Additional file 6: Figure S7. Bayesian skyline plots (BSP) for geographic responds to A. r. corallirostris. According to our results, groups (Pacific, Chiapas, Yucatán_CA), and for the A. rutila complex. A. r. corallirostris is a polyphyletic subspecies, since it is composed of two independent groups. A. r. corallirostris Acknowledgements was described by Bourcier and Mulsant (1846) with spec- We thank M. Robbins of the University of Kansas (Natural History Museum), imens from Escuintla, Guatemala. Therefore, the name A. and S. Birks of the University of Washington (Burke Museum of Natural History corallirostris should be applied to the Yucatán-CA clade and Culture) for providing tissue samples used in this study. We also thank Alejandro Gordillo Martínez, Fabiola Ramírez, Isabel Vargas‑Fernandez and (Ridgway 1901). This left the Chiapas group without a V ázquez‑López et al. Avian Res (2021) 12:61 Page 15 of 17 Bleiweiss R. Tempo and mode of hummingbird evolution. Biol J Linn Soc. Raúl Ivan Martínez Becerril for technical support. We acknowledge the General 1998;65:63–76. Direction of Wildlife (Instituto de Ecología, SEMARNAT, México) for providing Bleiweiss R, Kirsch JA, Matheus JC. DNA hybridization evidence for the collecting permits. L. Kiere reviewed the English. principal lineages of hummingbirds (Aves: Trochilidae). Mol Biol Evol. 1997;14:325–43. Authors’ contributions Booth TH, Nix HA, Busby JR, Hutchinson MF. BIOCLIM: the first species distribu‑ BEH‑B and MV ‑L designed the study. BEH‑B secured financial support. MV ‑L tion modelling package, its early applications and relevance to most carried out the laboratory work. MV‑L, BEH‑B, SMR‑B and LEZ ‑B analyzed the current MAXENT studies. Divers Distrib. 2014;20:1–9. data. BEH‑B, MV ‑L, NC‑R, AB‑H, LEZ ‑B and KR contributed to the writing and Bouckaert R, Vaughan TG, Barido‑Sottani J, Duchêne S, Fourment M, Gavry‑ improvement of the manuscript. All authors read and approved the final ushkina A, et al. BEAST 2.5: an advanced software platform for Bayesian manuscript. evolutionary analysis. PLoS Comput Biol. 2019;15:e1006650. Bourcier J, Mulsant E. Description de vingt espèces nouvelles d’oiseaux‑ Funding mouches. Ann Sci Phys Nat Agric. 1846;9:312–32. This research was supported by PAPIIT/DGAPA, Universidad Nacional Boyer AG, Jetz W. Biogeography of body size in Pacific island birds. Ecography. Autónoma de México (UNAM) through a grant to Blanca E. Hernández‑Baños 2010;33:369–79. (IN220620). LEZ‑B acknowledges the Postdoctoral scholarship provided by Braconnot P, Otto‑Bliesner B, Harrison S, Joussaume S, Peterchmitt JY, Abe ‑ DGAPA‑UNAM. Ouchi A, et al. Results of PMIP2 coupled simulations of the Mid‑Holo ‑ cene and Last Glacial Maximum‑Part 1: experiments and large ‑scale Availability of data and materials features. Clim Past. 2007;3:261–77. The sequences generated during the current study are available in GenBank Cavender‑Bares J, Gonzalez‑Rodriguez A, Pahlich A, Koehler K, Deacon N. (accession numbers: MZ998668‑MZ998740, OK624605‑ OK624657, OK614044‑ Phylogeography and climatic niche evolution in live oaks (Quercus OK614080, OK618334‑ OK618361). series Virentes) from the tropics to the temperate zone. J Biogeogr. 2011;38:962–81. Chaves AV, Clozato CL, Lacerda DR, Sari H, Santos FR. Molecular taxonomy Declarations of Brazilian tyrant‑flycatchers (Passeriformes, Tyrannidae). Mol Ecol. 2008;8:1169–77. Ethics approval and consent to participate Clegg SM, Owens PF. The “island rule” in birds: medium body size and its eco‑ This study was performed in line with the principles of the Declaration of logical explanation. Proc R Soc B Biol Sci. 2002;269:1359–65. Helsinki. Cobos ME, Osorio‑ Olvera L, Soberón J, Peterson AT. ellipsenm: ecological niche characterizations using ellipsoids. R package; 2020. https:// github. com/ Consent for publication marlo necob os/ ellip senm. Not applicable. 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: Competing interests low mitochondrial DNA differentiation reveals recent divergence. The authors declare that they have no competing interests. Condor. 2008;4:729–39. de Silva HG, Pérez Villafaña MG, Cruz‑Nieto J, Cruz‑Nieto MÁ. Are some of Author details 1 the birds endemic to the Tres Marías Islands (Mexico) species? Bull Br Museo de Zoología, Facultad de Ciencias, Universidad Nacional Autónoma de 2 Ornithol Club. 2020;140:7–37. México, Apartado Postal, 04510 Mexico City, Mexico. Department of Biol‑ 3 DeLattre PA. Oiseaux‑mouches nouveaux au peu connus, découverts au ogy, Ithaca College, Ithaca, NY 14850, USA. Museo de Zoología, Facultad de Guatimala. L’echo Du Monde Savant. 1843;7:1068–70. Estudios Superiores Zaragoza, Universidad Nacional Autónoma de México, 4 Drummond AJ, Rambaut A, Shapiro B, Pybus OG. Bayesian coalescent infer‑ Mexico City, Mexico. Department of Biology, Colorado State University, Fort ence of past population dynamics from molecular sequences. Mol Biol Collins, CO 80523, USA. Evol. 2005;22:1185–92. Earl DA, vonHoldt BM. STRU CTU RE HARVESTER: a website and program for Received: 19 February 2021 Accepted: 31 October 2021 visualizing STRU CTU RE output and implementing the Evanno method. Conserv Genet Res. 2012;4:359–61. Ellegren H. Molecular evolutionary genomics of birds. Cytogenet Genome Res. 2007;117:120–30. Evanno G, Regnaut S, Goudet J. Detecting the number of clusters of individu‑ References als using the software STRU CTU RE: a simulation study. Mol Ecol. Arango A, Villalobos F, Prieto‑ Torres DA, Guevara R. The phylogenetic diversity 2005;14:2611–20. and structure of the seasonally dry forest in the Neotropics. J Biogeogr. Excoffier LG, Schneider S. Arlequin v. 3.1: an integrated software package for 2021;48:176–86. population genetics data analysis. Evol Bioinform Online. 2006;1:47–50. Arbeláez‑ Cortés E, Navarro‑Sigüenza AG. Molecular evidence of the taxonomic Excoffier L, Smouse PE, Quattro JM. Analysis of molecular variance inferred status of western Mexican populations of Phaethornis longirostris (Aves: from metric distances among DNA haplotypes: application to human Trochilidae). Zootaxa. 2013;3716:81–97. mitochondrial restriction region. Genetics. 1992;131:479–91. Arbeláez‑ Cortés E, Roldán‑Piña D, Navarro ‑Sigüenza AG. Multilocus phylo ‑ Fu Y. Statistical tests of neutrality of mutations against population growth, geography and morphology give insights into the recent evolution hitchhiking and background selection. Genetics. 1997;147:915–25. of a Mexican endemic songbird: Vireo hypochryseus. J Avian Biol. García‑Deras GM, Cortés‑Rodríguez N, Honey M, Navarro ‑Sigüenza AG, García‑ 2014;45:253–63. Moreno J, Hernández‑Baños BE. Phylogenetic relationships within the Bandelt HJ, Foster P, Röhl A. Median‑joining networks for inferring intraspecific genus Cynanthus (Aves: Trochilidae), with emphasis on C. doubledayi. phylogenies. Mol Biol Evol. 1999;16:37–48. Zootaxa. 2008;1742:61–8. Barber BR, Klicka J. Two pulses of diversification across the Isthmus of García‑ Trejo EA, Navarro‑Sigüenza AG. Patrones biogeográficos de la riqueza Tehuantepec in a montane Mexican bird fauna. Proc R Soc B Biol Sci. de especies y el endemismo de la avifauna en el oeste de México. Acta 2010;277:2675–81. Zool Mex. 2004;20:167–85. Benítez‑López A, Santini L, Gallego ‑Zamorano J, Milá B, Walkden P, Huijbregts Gill F, Donsker D, Rasmussen P. IOC World Bird List (v.11.1). 2021. https:// www. MA, et al. The island rule explains consistent patterns of body size evo‑ world birdn ames. org/. lution in terrestrial vertebrates. Nat Ecol Evol. 2021;5:768–86. González C, Ornelas JF, Gutiérrez‑Rodríguez C. Selection and geographic Billerman SM, Keeney BK, Rodewald PG, Schulenberg TS. Birds of the world. isolation influence hummingbird speciation: genetic, acoustic and Ithaca: Cornell Laboratory of Ornithology; 2020. https:// birds ofthe world. org/ bow/ home. Vázquez‑López et al. Avian Res (2021) 12:61 Page 16 of 17 morphological divergence in the wedge‑tailed sabrewing (Campylop - Miller MA, Pfeiffer W, Schwartz T. Creating the CIPRES Science Gateway for terus curvipennins). BMC Evol Biol. 2011;11:38. inference of large phylogenetic trees. Gateway computing environ‑ Gordon CE, Ornelas JF. Comparing endemism and habitat restriction in Mesoa‑ ments workshop. New Orleans: IEEE; 2010. merican tropical deciduous forest birds: implications for biodiversity Montaño‑Rendón M, Sánchez‑ González LA, Hernández‑Alonso G, Navarro ‑ conservation planning. Bird Conserv Int. 2000;10:289–303. Sigüenza AG. Genetic differentiation in the Mexican endemic Rufous‑ Grant PR. The adaptive significance of some size trends in island birds. Evolu‑ backed Robin, Turdus rufopalliatus (Passeriformes: Turdidae). Zootaxa. tion. 1965;19:355–67. 2015;4034:495–514. Hahn IJ, Hogeback S, Römer U, Vergara PM. Biodiversity and biogeography Moore WS. Inferring phylogenies from mtDNA variation: mitochondrial‑ gene of birds in Pacific Mexico along an isolation gradient from mainland trees versus nuclear‑ gene trees. Evolution. 1995;49:718–26. Chamela via coastal Marías to oceanic Revillagigedo Islands. Vertebr Morrone JJ. Biogeographical regionalisation of the Neotropical region. Zool. 2012;62:123–44. Zootaxa. 2014;3782:1–110. Hall TA. BioEdit: a user‑friendly biological sequence alignment editor and Navarro‑Sigüenza AG, Peterson AT. An alternative species taxonomy of the analysis program for Windows 95/98/NT. Nucleic Acids Symp Ser. birds of Mexico. Biota Neotrop. 2004;4:1–32. 1999;41:95–8. Nelson EW. Descriptions of new birds from the Tres Marias Islands, Western Hernández‑Baños BE, Zamudio ‑Beltrán LE, Milá B. Phylogenetic relationships Mexico. Proc Biol Soc Washington. 1898;12:63‑4. and systematics of a subclade of Mesoamerican emerald humming‑ Ornelas JF, González C, de los Monteros AE, Rodríguez‑ Gómez F, García‑Feria birds (Aves: Trochilidae: Trochilini). Zootaxa. 2020;4748:581–91. LM. In and out of Mesoamerica: temporal divergence of Amazilia hum‑ Hijmans RJ, Cameron SE, Parra JL, Jones PG, Jarvis A. Very high resolution mingbirds pre‑ dates the orthodox account of the completion of the interpolated climate surfaces for global land areas. Int J Climatol. Isthmus of Panama. J Biogeogr. 2014;41:168–81. 2005;25:1965–78. Ortiz‑Ramírez MF, Sánchez‑ González LA, Castellanos‑Morales G, Ornelas JF, Hijmans RJ. raster: geographic data analysis and modeling. R package. version Navarro‑Sigüenza AG. Concerted Pleistocene dispersal and genetic dif‑ 3.0‑12. 2020. https:// CRAN.R‑ proje ct. org/ packa ge= raster . ferentiation in passerine birds from the Tres Marías Archipelago, Mexico. Howell SN, Webb S. A guide to the birds of Mexico and northern Central Auk. 2018;135:716–32. America. Oxford: Oxford University Press; 1995. Osorio‑ Olvera L, Barve V, Barve N, Soberón J. nichetoolbox: from getting Hubisz M, Falush D, Stephens M, Pritchard J. Inferring weak population struc‑ biodiversity data to evaluating species distribution models in a friendly ture with the assistance of sample group information. Mol Ecol Resour. GUI environment. R package. version 0.2.0.0. 2016. http:// shiny. conab io. 2009;9:1322–32.gob. mx: 3838/ niche toolb2/. Huelsenbeck JP, Ronquist F. MrBayes: a program for the Bayesian inference of Pacheco MA, Battistuzzi FU, Lentino M, Aguilar RF, Kumar S, Escalante AA. Evo‑ phylogeny. Bioinformatics. 2002;17:754–5. lution of modern birds revealed by mitogenomics: timing the radiation Janes JK, Miller JM, Dupuis JR, Malenfant RM, Gorrell JC, Cullingham CI, et al. and origin of major orders. Mol Biol Evol. 2011;28:1927–42. The K = 2 conundrum. Mol Ecol. 2017;26:3594–602. Ridgway R. Part V: Family Trochilidae in: The birds of North and Middle Jiménez A, Ornelas JF. Historical and current introgression in a Mesoamerican America: a descriptive catalogue. Washington: Bulletin of the United hummingbird species complex: a biogeographic perspective. PeerJ. States National Museum No. 50; 1901. 2016;4:e1556. Paynter RA. The ornithogeography of the Yucatán Peninsula. New Haven: Yale Koleff P, Urquiza‑Haas T, Contreras B. Prioridades de conservación de los University; 1955. bosques tropicales en México: reflexiones sobre su estado de conser ‑ Peterson AT, Navarro AG. Western Mexico: a significant center of avian end‑ vación y manejo. Ecosistemas. 2012;21:6–20. emism and challenge for conservation action. Cotinga. 2000;14:42–6. Lanfear R, Calcott B, Ho SY, Guindon S. PartitionFinder: combined selection Peterson AT, Navarro‑Sigüenza AG, Li X. Joint effects of marine intrusion of partitioning schemes and substitution models for phylogenetic and climate change on the Mexican avifauna. Ann Assoc Am Geogr. analyses. Mol Biol Evol. 2012;29:1695–701. 2010;100:908–16. Lanfear R, Frandsen PB, Wright AM, Senfeld T, Calcott B. PartitionFinder 2: new Phillips SJ, Dudík M, Schapire RE. Maxent software for modeling species niches methods for selecting partitioned models of evolution for molecular and distributions. version 3.4.1. 2021. and morphological phylogenetic analyses. Mol Biol Evol. 2017;34:772–3. Pritchard JK, Stephens M, Donnelly P. Inference of population structure using Lawrence GN. Annals of the Lycaeum of natural history of New York. New York: multilocus genotype data. Genetics. 2000;155:945–59. Lyceum of Natural History; 1866. R Core Team. R: a language and environment for statistical computing. Vienna: Lerdau M, Whitbeck J, Holbrook NM. Tropical deciduous forest: death of a R Foundation for Statistical Computing; 2020. https:// www.R‑ proje ct. biome. Trends Ecol Evol. 1991;6:201–2. org/. Lerner HRL, Meyer M, James HF, Hofreiter M, Fleischer RC. Multilocus resolution Rambaut A, Drummond AJ, Xie D, Baele G, Suchard MA. Posterior summariza‑ of phylogeny and timescale in the extant adaptive radiation of Hawai‑ tion in Bayesian phylogenetics using Tracer 1.7. Syst Biol. 2018;67:901–4. ian honeycreepers. Curr Biol. 2011;21:1838–44. Rambaut A, Drummond AJ. TreeAnnotator v.2.6.2:MCMC Output analysis. Librado P, Rozas J. DnaSP v5: a software for comprehensive analysis of DNA 2013; http:// beast. bio. ed. ac. uk. polymorphism data. Bioinformatics. 2009;25:1451–2. Rambaut A. FigTree v1.4.2: Tree figure drawing tool. 2014; http:// tree. bio. ed. Lomolino MV. Body size of mammals on islands: the island rule reexamined. ac. uk. Am Nat. 1985;125:310–6. Ramírez‑Barrera SM, Hernández‑Baños BE, Jaramillo ‑ Correa JP, Klicka J. Deep McCormack JE, Peterson AT, Bonaccorso E, Smith TB. Speciation in the high‑ divergence of Red‑ crowned Ant Tanager (Habia rubica: Cardinalidae), a lands of Mexico: genetic and phenotypic divergence in the Mexican jay multilocus phylogenetic analysis with emphasis in Mesoamerica. PeerJ. (Aphelocoma ultramarina). Mol Ecol. 2008;17:2505–21. 2018;6:e5496. McGuire JA, Witt CC, Altshuler DL, Remsen JV Jr. Phylogenetic systematics and Rodríguez‑ Gómez F, Ornelas JF. Genetic structuring and secondary contact biogeography of hummingbirds: Bayesian and maximum likelihood in the white‑ chested Amazilia hummingbird species complex. J Avian analyses of partitioned data and selection of an appropriate partition‑ Biol. 2018;49:e01536. ing strategy. Syst Biol. 2007;56:837–56. Ryan RM. The biotic provinces of Central America. Acta Zool Mex. 1963;6:1–54. McGuire JA, Witt CC, Remsen JV Jr, Corl A, Rabosky DL, Altshuler DL, et al. Rzedowski J. Vegetación de México. México: Limusa; 1978. Molecular phylogenetics and the diversification of hummingbirds. Curr Savage JM. The origins and history of the Central America herpetofauna. Biol. 2014;24:910–6. Copeia. 1966;4:719–66. Meave JA, Romero‑Romero MA, Salas‑Morales SH, Pérez‑ García EA, Gallardo‑ Seeholzer GF, Brumfield RT. Isolation by distance, not incipient ecological spe ‑ Cruz JA. Diversidad, amenazas y oportunidades para la conservación ciation, explains genetic differentiation in an Andean songbird (Aves: del bosque tropical caducifolio en el estado de Oaxaca, México. Ecosis‑ Furnariidae: Cranioleuca antisiensis, Line‑ cheeked Spinetail) despite near temas. 2012;21:85–100. threefold body size change across an environmental gradient. Mol Ecol. Miller MJ, Lelevier MJ, Bermingham E, Klicka JT, Escalante P, Winker K. Phyloge‑ 2018;27:279–96. ography of the Rufous‑tailed Hummingbird (Amazilia tzacatl). Condor. 2011;113:806–16. V ázquez‑López et al. Avian Res (2021) 12:61 Page 17 of 17 Smith H. Las provincias bióticas de México, según la distribución geográfica de Weller AA. Cinnamon Hummingbird. In: del Hoyo JJ, Elliott A, Sargatal J, edi‑ las lagartijas del género Sceloporus. Anales De La Escuela Nacional De tors. Handbook of the birds of the world. Barn‑ owls to hummingbirds. Ciencias Biológicas. 1941;2:103–10. Barcelona: Lynx Editions; 1999. p. 596‒597. Smith BT, Escalante P, Hernández‑Baños BE, Navarro ‑Sigüenza AG, Rohwer S, Werneck FP, Costa GC, Colli GR, Prado DE, Sites JW. Revisiting the historical Klicka J. The role of historical and contemporary processes on phylo‑ distribution of Seasonally Dry Tropical Forests: new insights based on geographic structure and genetic diversity in the Northern Cardinal, palaeodistribution modelling and palynological evidence. Global Ecol Cardinalis Cardinalis. BMC Evol Biol. 2011;11:136. Biogeogr. 2011;20:272–88. Sorenson MD, Ast JC, Dimcheff DE, Yuri T, Mindell DP. Primers for a PCR‑based West RC. The natural regions of Middle America. In: West RC, editor. Handbook approach to mitochondrial genome sequencing in birds and other of Middle American Indians. Austin: University of Texas Press; 1964. p. vertebrates. Mol Phylogenet Evol. 1999;12:105–14. 363–83. Stamakis A. RaxML version 8: a tool for phylogenetic analysis and post‑analysis Wickham H, Francois R, Henry L, Müller K. dplyr: A Grammar of Data Manipula‑ of large phylogenies. Bioinformatics. 2014;30:1312–3. tion. R package. version 0.8.5. 2020. https:// CRAN.R‑ proje ct. org/ packa Stecher G, Tamura K, Kumar S. Molecular evolutionary genetics analysis ge= dplyr. (MEGA) for macOS. Mol Biol Evol. 2020;37:1237–9. Wilcox TP, Zwickl DJ, Heath TA, Hillis DM. Phylogenetic relationships of the Stephens M, Scheet P. Accounting for decay of linkage disequilibrium in dwarf boas and a comparison of Bayesian and bootstrap measures of haplotype inference and missing‑ data imputation. Am J Human Genet. phylogenetic support. Mol Phylogenet Evol. 2002;25:361–71. 2005;76:449–62. Zamudio‑Beltrán LE, Hernández‑Baños BE. Genetic and morphometric diver ‑ Stephens M, Smith N, Donelly P. A new statistical method for haplotype recon‑ gence in the Garnet‑ Throated Hummingbird Lamprolaima rhami (Aves: struction from population data. Am J Human Genet. 2001;68:978–89. Trochilidae). PeerJ. 2018;6:e5733. Tajima F. Statistical method for testing the neutral mutation hypothesis by Zamudio‑Beltrán LE, Ornelas JF, Malpica A, Hernández‑Baños BE. Genetic and DNA polymorphism. Genetics. 1989;123:585–95. morphological differentiation among populations of the Rivoli’s Hum‑ Thompson JD, Higgins DG, Gibson TJ. CLUSTAL W: improving the sensitivity of mingbird (Eugenes fulgens) species complex (Aves: Trochilidae). Auk. progressive multiple sequence alignment through sequence weight‑ 2020;137:ukaa032. ing, position‑specific gap penalties and weight matrix choice. Nucleic Zizka A, Silvestro D, Andermann T, Azevedo J, Duarte‑Ritter C, Edler D, et al. Acids Res. 1994;22:4673–80. CoordinateCleaner: standardized cleaning of occurrence records from van Rossem AJ. A northwest race of the Cinnamon Hummingbird. Condor. biological collection databases. Method Ecol Evol. 2019;10:744–51. 1938;40:226–7. van Valen LM. Pattern and the balance of nature. Evol Theor. 1973;1:31–49. 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. 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. Weir JT. Divergent timing and patterns of species accumulation in lowland and highland neotropical birds. Evolution. 2006;60:842–55. Re Read ady y to to submit y submit your our re researc search h ? Choose BMC and benefit fr ? Choose BMC and benefit from om: : fast, convenient online submission thorough peer review by experienced researchers in your field rapid publication on acceptance support for research data, including large and complex data types • gold Open Access which fosters wider collaboration and increased citations maximum visibility for your research: over 100M website views per year At BMC, research is always in progress. Learn more biomedcentral.com/submissions http://www.deepdyve.com/assets/images/DeepDyve-Logo-lg.png Avian Research Springer Journals

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

Loading next page...
 
/lp/springer-journals/phylogeography-and-morphometric-variation-in-the-cinnamon-hummingbird-cL0w7Fkizv

References (94)

Publisher
Springer Journals
Copyright
Copyright © The Author(s) 2021
eISSN
2053-7166
DOI
10.1186/s40657-021-00295-0
Publisher site
See Article on Publisher Site

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, STRU CTU RE, F , AMOVA), Bayesian and Maximum Likelihood phylogenetic analyses, and divergence time estimates. In order to ST 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 projec‑ tions onto past scenarios revealed no drastic changes in suitable conditions, but revealed some possible refuges. Mor‑ phometric 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 geo ‑ graphic isolation (limited gene flow), differences in current environmental conditions, and historical habitat fragmen‑ tation promoted by past events (Pleistocene refugia). Four well‑ defined clades comprise the A. rutila complex, and we *Correspondence: behb@ciencias.unam.mx Museo de Zoología, Facultad de Ciencias, Universidad Nacional Autónoma de México, Apartado Postal, 04510 Mexico City, Mexico Full list of author information is available at the end of the article © The Author(s) 2021. 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:// creat iveco mmons. org/ licen ses/ by/4. 0/. The Creative Commons Public Domain Dedication waiver (http:// creat iveco mmons. org/ publi cdoma in/ zero/1. 0/) applies to the data made available in this article, unless otherwise stated in a credit line to the data. Vázquez‑López et al. Avian Res (2021) 12:61 Page 2 of 17 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. Keywords: Amazilia rutila, Genetic structure, Genetic variation, Phylogeography, Tres Marías Islands, Trochilidae, Tropical deciduous forest Background Mexico have been considered to have high potential for The Mesoamerican dominion is the area where the diversification (García-Deras et  al. 2008; González et  al. Neartic and Neotropics overlap, and includes most of the 2011; Smith et  al. 2011; Ramírez-Barrera et  al. 2018; Mexican and Central American lowlands as well as the Hernández-Baños et  al. 2020). Furthermore, a study Mexican Transition Zone (Morrone 2014). This biogeo - based on molecular data showed that species distributed graphic dominion and surrounding areas are well known in the Neotropical lowlands have a low and constant rate for possessing high levels of species richness, which are of diversification over time; extrapolating those rates of thought to have originated from its complex topography, diversification to the present leads to a greater number wide range of environmental conditions, and climatic of lineages than is currently described (Weir 2006). The history (Savage 1966). One of the most distinctive habi- number of species present in the tropics is thus not prop- tats found in the Mesoamerican dominion is the tropical erly reflected in the currently accepted taxonomy. deciduous forest, which has a high number of endemic Hummingbirds belong to one of the most diverse bird species (Gordon and Ornelas 2000). families (Trochilidae: 361 spp.; Gill et al. 2021), with spe- Tropical deciduous forest has the driest conditions cialized nectarivorous feeding mode and high rates of of the tropical forests, and precipitation is strongly sea- speciation. This has resulted in a huge diversity of bill sonal (Rzedowski 1978; Meave et al. 2012). In addition, it morphologies and color patterns among nine highly sup- is considered one of the most threatened habitats in the ported groups (Bleiweiss et al. 1997; McGuire et al. 2014). Neotropics due to land use change, habitat fragmenta- In this study, we examined the evolutionary history of the tion, and human population growth (Lerdau et  al. 1991; Cinnamon Hummingbird (Amazilia rutila), a medium- Koleff et al. 2012). The diversity of the tropical deciduous sized hummingbird that inhabits tropical deciduous for- forest has been attributed to historical factors such as the ests ranging from 0 to 1600  m of elevation (Howell and climatic fluctuations of the Pleistocene, in which colder Webb 1995). This species is mainly distributed among and drier periods promoted its fragmentation, followed three biogeographic provinces: the Pacific lowlands, by periods of expansion with warmer and more humid Yucatán Peninsula, and Mosquito, based on the region- conditions (Werneck et al. 2011). alization of the Neotropical region (Smith 1941; Ryan Some studies of bird species distributed in Mesoameri- 1963; West 1964; Morrone 2014). It is a resident species can deciduous forests have shown high levels of genetic from Sinaloa in northwestern Mexico, south along the and morphological variation associated with historic Pacific Coast (including the Tres Marías Islands) to Costa and demographic events (e.g. Arbeláez-Cortés et  al. Rica, and along the coastal plain of the Yucatán Peninsula 2014), while others have focused on allopatric popula- south to Belize (including “offshore cays”) and north - tions, where long-term isolation, moderate or low gene eastern Nicaragua (Howell and Webb 1995). Based on flow, presence of geographic barriers, and historical morphology and allopatric distributions, four subspecies events explained the overall geographic structure (e.g. are recognized (Billerman et  al. 2020; Gill et  al. 2021). Smith et  al. 2011). One clear example of allopatric dif- A. r. rutila (DeLattre 1843) is distributed in western and ferentiation is the taxa inhabiting the Tres Marías Islands southern Mexico (from Jalisco to Oaxaca) and has green (Nayarit, Mexico), where there are several endemic plumage on the upperparts, cinnamon underparts and a subspecies (Cortés-Rodríguez et  al. 2008; Smith et  al. rufous tail. A. r. diluta (van Rossem 1938) inhabits north- 2011; Montaño-Rendón et  al. 2015). These subspecies western Mexico (from Sinaloa to Nayarit); its upper parts share the common feature of clear morphological dif- are more golden-bronze and underparts are paler and ferentiation, mainly in body size, from their continental less reddish than rutila. A. r. corallirostris (Bourcier and counterparts, representing independent evolutionary Mulsant 1846) is distributed from south and southeast- lineages with defined genetic structure (Ortíz-Ramírez ern Mexico (from Chiapas to the Yucatán Peninsula) to et  al. 2018). In addition to the Tres Marías Islands, the Costa Rica and is much more deeply colored than rutila. Yucatán Peninsula and the Coastal Plain of southwestern Finally, located within the Pacific lowlands province, the V ázquez‑López et al. Avian Res (2021) 12:61 Page 3 of 17 subspecies A. r. graysoni (Lawrence 1866) inhabits the the Washington University Genomics Center. All new Tres Marías Islands; it is darker in color and larger than sequences have been deposited in GenBank (Accession rutila (Ridgway 1901). In general, the systematics of A. numbers MZ998668-MZ998740, OK624605-OK624657, rutila are unresolved and controversial. Weller (1999) OK614044-OK614080, OK618334-OK618361). considered that A. rutila forms a superspecies with A. Sequences were proofread using Sequencher v.4.8 yucatanensis (Buff-bellied Hummingbird), and in a mul - (Gene Codes Corporation 2007), and aligned with the tilocus molecular study, the genus Amazilia formed a Clustal W function in BioEdit (Thompson et  al. 1994; polyphyletic group nested with the genera Hylocharis, Hall 1999). We corroborated the mitochondrial origin of Chrysuronia, Lepidopyga, and Damophila (McGuire our sequences using at least two of the following meth- et al. 2007). ods: amplifying/sequencing overlapping gene segments, The aims of this study were to: (1) analyze the genetic amplifying/sequencing one region with different primer variation and morphometric differentiation across the sets, and sequencing both DNA strands. We found no geographic distribution of A. rutila, (2) reconstruct phy- evidence of contamination of mtDNA sequences. logenetic relationships and past scenarios throughout the geographic range to analyze divergence times and demo- Phylogenetic analyses, genetic structure and neutrality graphic changes over time and space, and (3) propose tests an evolutionary history and taxonomic hypothesis of We used the Bayesian approach implemented in PHASE intraspecific lineages. We expected to find high levels of v.2.1 (Stephens et  al. 2001; Stephens and Scheet 2005) genetic structure associated with the isolation of allopat- to reconstruct haplotypes and avoid heterozygotes in ric populations and environmental conditions over con- nuclear sequences, selecting the pairs of haplotypes with tinuous ranges. a posterior probability higher than 0.90. We obtained the models of evolution for each gene using Partition- Methods Finder (Lanfear et  al. 2017), using the following param- Sampling, PCR and sequencing eters: linked branched lengths, greedy search algorithm We obtained 84 tissue samples from the A. rutila com- (Lanfear et  al. 2012), and the Bayesian Information Cri- plex and one sample from A. yucatanensis as the terion (BIC) for model selection. The phylogenetic anal - outgroup (see Additional file  1: Tables S1, S2). We sup- yses were conducted using each locus separately and plemented our sampling with GenBank sequences from using the concatenated matrix including nuclear and A. rutila, A. yucatanensis and A. tzacatl (McGuire et  al. mitochondrial loci (multilocus). We conducted phylo- 2014). Tissues samples were provided by the collection genetic analysis under Bayesian Inference (BI), using of the Museo de Zoología (MZFC) “Alfonso L. Herrera” MrBayes 2.0 (Huelsenbeck and Ronquist 2002) in CIP- (Universidad Nacional Autónoma de México), the Burke RES Science Gateway (Miller et  al. 2010). Each analysis Museum (University of Washington), and the Biodiver- consisted of four independent chains, random starting sity Institute (University of Kansas). trees, and uniform prior distribution of parameters. The DNA was extracted using the Qiagen DNAEasy kit, chains were run for 30 million generations, sampling following the manufacturer’s protocols. We ampli- trees every 1000 generations. The asymptote was deter - fied and sequenced two mitochondrial markers: par - mined visually, 5000 burn-in trees were discarded, and tial and complete NADH dehydrogenase subunit 2 the remaining trees from the plateau phase were used to (ND2: 605–1041  bp), and the full length of cytochrome estimate Bayesian posterior probabilities. We considered C oxidase subunit I gen (COI: 568  bp). We also ampli- that clades were strongly supported if they were present fied two nuclear markers: the regions between exons 4 in ~ 95% of the sampled trees (Huelsenbeck and Ronquist and 5 of the Muscle Skeletal Receptor Tyrosine Kinase 2002; Wilcox et  al. 2002). Also, a Maximum Likelihood gene (MUSK: 628  bp), and a segment comprising the analysis was conducted using RaxML v.8.0.0 (Stamakis end of exon 6 and the beginning of exon 8 of the Orni- 2014) in CIPRES Science Gateway (Miller et  al. 2010), thine Decarboxylase gene (ODC: 571). We used primers using separate partitions for each locus (see “Results”), L5219, H6313, L5758 and H5766 for the amplification of with 100 bootstrap replicates. ND2 (Sorenson et al. 1999); primers F1-COI and R2-COI We defined four different groups, corresponding to for COI (Chaves et  al. 2008); ODC2-F and ODC2-R for allopatric populations from sampled localities of A. ODC (McGuire et  al. 2007); and MUSK R3 and MUSK rutila (see Fig.  1): (1) Tres Marías Islands, (2) Pacific F3 for the amplification of MUSK (McGuire et  al. coast, (3) Chiapas, and (4) Yucatán Peninsula and Cen- 2007). All PCR products were confirmed on 1% aga - tral America (hereafter: Tres Marías, Pacific, Chia - rose gel, and sequencing reactions were performed by pas and Yucatán_CA). For the comparison between Vázquez‑López et al. Avian Res (2021) 12:61 Page 4 of 17 Fig. 1 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 populations, we specified geographic groups based on Maximum Likelihood model, with 100 bootstrap repli- the main clades recovered from phylogenetic analysis cates (Stecher et al. 2020). (concatenated-multilocus), and tested for substructure in Pacific and Yucatán_CA groups (see “Results”). A Divergence time estimates median-joining network was constructed to visualize We estimated divergence times with the concatenated the structure of populations, the number of haplotypes, matrix using Beast v.2.6.2 (Bouckaert et  al. 2019). Dif- their frequencies, and the relationships among individ- ferent partitions by gene were defined based on Pacheco uals using the program Network 4.5.1.0 (Bandelt et  al. et  al. (2011) for ND2 and COI, Lerner et  al. 2011 for 1999). We tested locus neutrality with Fu’s Fs (Fu 1997) ODC, and Ellegren (2007) for MUSK. We employed a and Tajima’s D (Tajima 1989) metrics in DnaSP v.5 strict clock and a constant population model as a tree (Librado and Rozas 2009), and obtained the summary prior. We chose the strict clock model empirically based statistics for each population. We tested the differen - on data fit, and the constant population model based tiation hypothesis using an analysis of molecular vari- on our demographic reconstruction (see “Results”). The ance AMOVA and the F fixation index using Arlequin analysis was run for 200 million generations, sampling ST 3.1 (Excoffier and Schneider 2006). These analyses are every 1000. We used Tracer v.1.7 (Rambaut et  al. 2018) useful for observing population structure, estimating to ensure adequate Effective Sample Size (ESS > 200), population differentiation directly from molecular data, which was reached for most statistics, with the exception and testing the hypothesis of differentiation (Excoffier of prior and probability statistics. With TreeAnnotator v. et  al. 1992). Additionally, we used STRU CTU RE 2.3.2 2.6.2 (Rambaut and Drummond 2013) we discarded the (Pritchard et  al. 2000) for population assignments first 25% or trees as burn-in and produced a maximum under an admixed and LocPrior model (Hubisz et  al. clade credibility tree with 95% highest probability density 2009), with 10,000 iterations of burn-in and 20,000 intervals. The final tree was visualized in FigTree v.1.4.2 MCMC (Markov chain Monte Carlo), for 20 replicates (Rambaut 2014). of each K value (from 1 to 6). To evaluate STRU CTU RE results, we used the Evano method to choose K value Morphometric analysis by ∆K (Evanno et  al. 2005), as implemented on the We took linear measurements of five morphological traits Structure Harvester website (Earl and vonHoldt 2012). from 210 study skins housed at the American Museum of Genetic distances were obtained with MEGA under Natural History (AMNH) and at the Museo de Zoología V ázquez‑López et al. Avian Res (2021) 12:61 Page 5 of 17 “Alfonso L. Herrera” (MZFC). We measured wing chord account the logistic threshold from maximum training (WC) and tail length (TL) using an Avinet calibrated sensitivity plus specificity. ruler to the nearest 0.5  mm; also, we measured culmen We obtained Bayesian Skyline Plots (BSP) using Beast length (CL), beak width (BW), and beak height (BH) to v.2.6.3 (Bouckaert et al. 2019) to infer population dynam- the nearest 0.1 mm with a Mitutoyo digital caliper. ics and changes in demography over time. We used the We assigned each specimen to geographic groups ND2 dataset and performed four different runs: three of defined based on the haplotype network (see “Results”; them based on geographic groups (Pacific, Chiapas and Fig.  1). We performed Mann–Whitney U tests between Yucatán_CA), and one based on the whole A. rutila com- sexes within each geographic group, with Bonferroni cor- plex. The Tres Marías group was not included in an inde - rection to evaluate sexual dimorphism. We also carried pendent run because of the lack of informative sites for out Kruskall-Wallis test with “geographic group” as the ND2. We set a GTR + substitution model, a strict molec- grouping variable to assess geographic variation in the ular clock, a Coalescent Bayesian Skyline as the tree measured characters. Finally, we carried out a PCA to prior, and a Piecewise-constant skyline model (Drum- visually examine the structure of variation in multivari- mond et  al. 2005) with five groups. We used a mutation ate space. All statistical analyses were performed in R (R rate of 0.0145 substitutions per site per lineage per mil- Core Team 2020). lion years for ND2, following Lerner et al. (2011), and ran each analysis through 20 million generations. Ecological Niche Modelling under past scenarios and demography Results To infer suitability of environmental conditions under Phylogenetic analysis past scenarios, we performed projections using Eco- The model of sequence evolution selected for the logical Niche Models. Presence records of A. rutila were mtDNA data set was GTR + G (nst = 6, rates = gamma). downloaded from GBIF (Global Biodiversity Information The nucleotide composition of mtDNA was as follows: Facility) and accessed from R via rgbif (https:// github. T = 23.9%, C = 34.5%, A = 30.8%, G = 10.6% for ND2; and com/ ropen sci/ rgbif; taxon key: 2476417; https:// doi. org/ T = 25.9%, C = 15.7%, A = 26.4%, G = 31.8% for COI. The 10. 15468/ dl. 7k98v7). These records were cleaned over parsimony informative sites for ND2 were 72 (605  bp), multiple steps using the following R packages: Coordi- and 38 for COI (568 bp). For nuclear markers, the mod- nateCleaner (Zizka et  al. 2019), nichetoolbox (Osorio- els of sequence evolution were: F81 (nst = 1) for MUSK, Olvera et  al. 2016), dplyr (Wickham et  al. 2020), and and GTR + I + G (nst = 6, rates = invgamma) for ODC. raster (Hijmans 2020). We used the sets of bioclimatic The nucleotide composition was: T = 29.6%, C = 19.4%, layers (current and past) from WorldClim (Hijmans et al. A = 29.4%, G = 21.4% for MUSK; and T = 36.4%, 2005; Braconnot et al. 2007; Booth et al. 2014). Our mod- C = 17.3%, A = 24.7%, G = 21.4% for ODC. The parsi - els for the past were projected onto the Last InterGlacial mony informative sites for MUSK were 19 (628 bp), and (LIG: 140–120 kya), the Last Glacial Maximum (LGM: 21 83 for ODC (571 bp). kya), and the Mid Holocene (MH: 6 kya) scenarios. We According to the topology from the concatenated data performed these projections using the continental phy- set (mitochondrial and nuclear genes) under Bayesian logroups (Pacific coast, Chiapas, and Yucatán and Cen - Inference and Maximum Likelihood, A. rutila forms a tral America); the Tres Marías group was not included monophyletic group with respect to A. yucatanensis and in these analyses due to the low number of occurrence A. tzacatl (see Fig. 2). Within A. rutila, three major line- records (n = 3). The selection of bioclimatic variables was ages can be identified: Pacific coast, Chiapas, and Yuca - based primarily on a principal component analysis, where tán Peninsula and Central America. There was also a we considered the most important variables of the first well-supported clade (posterior probability = 0.93) within four principal components (from one to three variables), the Pacific phylogroup, which contained all individuals ensuring that no intercorrelated variables were selected from the Tres Marías Islands. Within the Yucatán_CA (Pearson correlation r < 0.75, Additional file  2: Fig. S1). phylogroup, four of the five individuals from Nicaragua To define the area of accessibility for the species, we used grouped together into a well-supported clade (pp = 0.99; the ellipsenm R package (Cobos et  al. 2020) to delimit bootstrap = 99). The individual phylogenies for mito - the area that contained occurrence points across the dis- chondrial loci (ND2 and COI) recovered three major lin- tribution of the species with a polygon representing a eages (Pacific, Chiapas, and Yucatán_CA) (see Additional buffer area of 75 km (“M” area). Models were performed file  2: Fig. S1). However, individual nuclear locus phylog- in Maxent v.3.4.1 (Phillips et al. 2021) with regularization enies did not recover any geographic structure (see Addi- multiplier = 2 and with 10 replicates of cross-validation tional file 2: Fig. S1). with no extrapolation. To binarize outputs, we took into Vázquez‑López et al. Avian Res (2021) 12:61 Page 6 of 17 Fig. 2 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 group was separated from the rest by 57 mutational steps Deviations from neutrality were rejected for mitochon- (Fig. 1). For COI, we found 11 haplotypes in total: 4 cor- drial loci (ND2 and COI) and nuclear loci, except for responded to the Pacific group, 3 to Chiapas, 3 to Yuca - the ODC locus (see Table  1). We found 21 haplotypes tán-Central America, and 1 to the Tres Marías Islands for ND2: 9 corresponded to the Pacific group, 9 to the (shared with 2 individuals from the Pacific group from Yucatán-Central America group, 2 to Chiapas, and 1 to Guerrero). The Yucatán-Central America group was the Tres Marías Islands. The Yucatán-Central America V ázquez‑López et al. Avian Res (2021) 12:61 Page 7 of 17 Table 1 Neutrality test by locus Group Tajimas’ D Fu and Li’s D ND2 COI MUSK ODC ND2 COI MUSK ODC ‡ ‡ ‡ ** ‡ ‡ ‡ * All − 0.14 − 0.57 − 1.45 − 2.3 − 1.47 − 2.39 0.07 − 3.91 ‡ ‡ ‡ ‡ ‡ ‡ ‡ * Pacific − 1.28 2.38 − 0.62 − 2.13 0.15 1.06 − 0.82 − 4.69 ‡ ‡ ‡ ‡ Tres Marías 0.0 0.0 − 1.13 − 1.23 0.0 0.0 − 1.15 − 1.27 ‡ ‡ ‡ ‡ ‡ ‡ Chiapas 0.0 − 0.17 − 0.31 0.24 0.0 1.44 0.13 0.71 ‡ ‡ ‡ * ‡ ‡ ‡ * Yucatán_CA − 0.57 − 1.08 0.70 − 1.93 − 0.43 − 1.31 1.149 − 2.56 ‡ ‡ Nicaragua − 0.75 ND ND ND − 0.75 ND ND ND ‡ * ** *** P > 0.5, P < 0.05, P < 0.01, P < 0.001 The mtDNA genetic distances among groups showed Table 2 mtDNA summary statistics an overall genetic distance of 5.24% for the Pacific group, Geographic group N H Hd Pi 5.69% for Chiapas, and 7.52% for the Yucatán_CA group Pacific 25 5 0.51 0.00196 (Table  3). The pairwise F fixation indices were all sig - ST Tres Marías 6 1 0.0 0.0 nificant and indicated strong population structure of Chiapas 5 1 0.0 0.0 geographic groups for mtDNA data (Table  4). We did Yucatán_CA 14 5 0.59 0.00212 not find strong differentiation between the Pacific and Nicaragua 4 3 0.83 0.00145 its subgroup from Guerrero. However, the Guerrero sub- All 50 12 0.83 0.03037 group showed a more robust differentiation with respect to the Tres Marías group than to the rest of the Pacific N number of sequences, H number of haplotypes, Hd haplotype diversity, Pi nucleotide diversity subgroup. Substructure was detected within the Yucatán_ CA group, with some Nicaraguan samples (UWBM-6911, UWBM-68991, UWBM-69261, and UWBM-69388) sep- Table 3 mtDNA pairwise genetic distances among geographic arated from the rest of the group, with an F of 0.78549. ST groups in the Amazilia rutila complex AMOVA analyses revealed that most of the genetic vari- Pacific Tres Marías Chiapas Yucatán_CA ation was among geographic groups (96.26%), and there was little variation within populations (3.74%; Table  3). Pacific The F value from general AMOVA results indicated ST Tres Marías 0.00629 strong genetic structure and differentiation among popu - Chiapas 0.03870 0.03766 lations (F = 0.96; Table 5). ST Yucatán_CA 0.07524 0.07389 0.07487 The STRU CTU RE analyses detected two genetic groups (K = 2; first-level analysis; See Additional file  3: Figs. S2, S3; Fig.  3). One group comprised all individuals separated from the rest by 15 mutational steps (Fig.  1). from Sinaloa to Chiapas, including Tres Marías. The sec - The Pacific coast and Yucatán-Central America groups ond group included the Yucatán Peninsula and Central had the highest haplotype and nucleotide diversity (see American individuals. To evaluate potential substructure Table 2). within these groups, we ran a hierarchically structured Table 4 Pairwise F indices (below diagonal, P ≤ 0.05) among populations in the Amazilia rutila complex ST Pacific Guerrero Tres Marías Chiapas Yucatán_CA Nicaragua Pacific ‒ Guerrero 0.33425 ‒ Tres Marías 0.55131 1 ‒ Chiapas 0.93300 1 1 ‒ Yucatán_CA 0.96654 0.99167 0.99209 0.99237 ‒ Nicaragua 0.96132 0.98693 0.98831 0.98791 0.78549 ‒ Geographic groups taking into account for this analysis: Pacific, Guerrero (clade within Pacific: CHIS529, CHIS533, CHIS534, CHIS542, CHIS549), Tres Marías, Chiapas, Yucatán_CA, Nicaragua (clade within Yucatán_CA: UWBM6911, UWBM68991, UWBM69261, UWBM69388) Vázquez‑López et al. Avian Res (2021) 12:61 Page 8 of 17 analysis. Individual structure analyses for the first group before present (Myr BP) (95% highest posterior density revealed two genetic groups (K = 2), with no ancestral [HPD] 8.04‒10.87; Fig.  4), while the A. rutila complex mixing between them. One group was composed of indi- root was dated 7.05 Myr ago (95% HPD 5.88‒8.29). viduals from Chiapas and the second group recovered This estimate also corresponds to the split between the all of the Pacific clade including Tres Marías individuals. two main clades: (1) Pacific + Tres Marías and Chia- The substructure within the Yucatán_CA group recov - pas, and (2) Yucatán_CA. The Pacific clade and Chi- ered two groups (K = 2). One group included samples apas split around 3.99 Myr BP (95% HPD 3.16‒4.84). from western Nicaragua and three samples from Yuca- The Pacific Clade node origin was 1.65 Myr BP (95% tán. The second group included one sample from western HPD 0.81‒1.81  Ma). According to these estimates, A. Nicaragua, one sample from Guatemala and the rest of rutila arrived at the Tres Marías Islands 0.2 Myr BP the Yucatán Peninsula individuals. (95% HPD 0.6‒0.44  Ma). The age of the continental subclades within the Pacific clade suggests that Pleis- Divergence time estimates tocene climatic changes may have caused contractions Divergence times estimated that the A. rutila com- in the continental populations of A. rutila and resulted plex split from its sister clade around 9.4 million years in incipient substructure, however although the Tres Table 5 mtDNA AMOVA comparing geographic groups in the Amazilia rutila complex Source of variation df Sum of squares Variance components Percentage of variation (%) Among groups 5 564.198 14.31 96.26 Within groups 47 26.141 0.55 3.74 Total 52 Fixation index F : 0.96 ST Significance test P < 0.0001 See geographic groups in Table 4 Fig. 3 STRU CTU RE 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) V ázquez‑López et al. Avian Res (2021) 12:61 Page 9 of 17 Fig. 4 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 Vázquez‑López et al. Avian Res (2021) 12:61 Page 10 of 17 Marías Islands is well supported, other clades are not (Annual Precipitation), BIO15 (Precipitation Season- supported. The Chiapas clade arose 0.76 Myr BP (95% ality), and BIO18 (Precipitation of Warmest Quarter). HPD 0.43‒1.15 Ma). The average evaluation metrics for the model results The crown node for the Yucatán_CA clade emerged were: training AUC = 0.7598 (St.D v ± 0.005), and a 1.98 Myr BP (95% HPD 1.42‒2.6  Ma). Internal clades threshold value of 0.4178 (Maximum training sensi- in the Yucatán_CA group were detected by Beast anal- tivity plus specificity Logistic threshold). Projections ysis, but they had weak node support, with the excep- onto past scenarios for the whole complex are shown tion of the Nicaragua clade, which had intermediate in Additional file  5: Figs. S5, S6. The suitability of envi - support (0.7 posterior probability) and a date of 0.12 ronmental conditions has changed somewhat over Myr BP (95% HPD 0.03‒0.25 Ma). time, with subtle reductions from LIG to LGM and posterior arrangements after LGM through the MH period. Projections during LGM showed discontinu- Morphometrics ous distribution, which contrasts with the current con- We found minor sexual dimorphism in this species. tinuous distribution for the species (Pacific coast and Females were slightly but statistically significantly Central America). Independent analyses were run for smaller than males in two body size variables we each geographic group, in which the first two princi - measured (WC, W = 2196.5, p < 0.001; TL , W = 2194, pal components explained most of the variation (see p < 0.001). However, the effect sizes were very small, so Additional file  5: Figs. S5, S6). The bioclimatic vari - we carried out subsequent tests of geographic varia- ables selected for the Pacific group were BIO1, BIO3, tion with both sexes pooled together. BIO4, BIO7, BIO13, BIO17 and BIO19; for the Chiapas There were statistically significant differences group were BIO4, BIO5, BIO12, BIO13, BIO14, BIO15, among geographic groups (Kruskall-Wallis tests WC: BIO16, and BIO18; and for the Yucatán_CA group 2 2 χ = 78.86, p < 0.001; TL: χ = 73.09, p < 0.001; BH: were BIO3, BIO4, BIO7, BIO11, and BIO13. The result - 4 4 2 2 χ = 47.45, p < 0.001; C L: χ = 62.61, p < 0.001; BW: ing metrics were as follows: Pacific AUC = 0.8046 (St. 4 4 χ = 90.30, p < 0.001). The Tres Marías phylogroup had Dv ± 0.0473), threshold = 0.358; Chiapas AUC = 0.8184 larger values than the three continental groups for all (St.Dv ± 0.1119), threshold = 0.4255; Yucatán_C A variables. This was consistent with the results of our AUC = 0.6515 (St.D v ± 0.0470), threshold = 0.4515. PCA, which showed that the continental phylogroups Projections for Pacific group predicted some areas largely overlapped in morphological measurements, where the Chiapas and Yucatán_CA groups are cur- but the Tres Marías group was significantly larger. rently distributed (see Additional file  5: Figs. S5, S6). Principal component 1 explained 74.75% of the cumu- Also, for the projections in the Yucatán_CA group, lative variance and roughly corresponds to overall suitable areas were predicted where the Pacific and body size, with large contributions from Wing Chord Chiapas groups now occur. The most conservative pro - (WC) and Tail Length (TL) (Fig. 5). jections for suitable habitats were recovered when the Chiapas group was modeled. Past projections revealed Ecological Niche modelling onto past scenarios a reduction of suitability during LGM and an increase and demography during MH periods in all cases. During the LIG period, We obtained a total of 445 presence records for A. the Pacific showed continuous suitable habitat along rutila after filtering steps. The number of occurrence the Pacific coast to Central America. In LIG projections points for each group was: (1) Pacific: 144 records, for Chiapas and Yucatán_CA, no suitable conditions (2) Chiapas: 20 records, and (3) Yucatán_CA: 281 were recovered (see Additional file 5 : Figs. S5, S6). records. The results of the PCA of climate variables The BSP results showed that population dynamics and for the whole complex showed that the first four prin - demography in the A. rutila complex have been generally cipal components explained most of the variation constant over time with a notable population reduction (88%, see Additional file  4: Fig. S4). According to these close to the present. Effective population size patterns results and correlation values, we selected the follow- showed a stationary trend in the geographic groups ing bioclimatic variables: BIO5 (Max Temperature of analyzed (Pacific, Chiapas, Yucatán_CA), with recent Warmest Month), BIO7 (Temperature Annual Range), BIO11 (Mean Temperature of Coldest Quarter), BIO12 (See figure on next page.) Fig. 5 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 V ázquez‑López et al. Avian Res (2021) 12:61 Page 11 of 17 Fig. 5 (See legend on previous page.) Vázquez‑López et al. Avian Res (2021) 12:61 Page 12 of 17 population reductions in the Pacific and Chiapas groups, Pleistocene, Middle American biological diversity was and a reduction followed by a slight increase towards the influenced by several major processes. One is the uplift of present in the Yucatán_CA group (see Additional file  5: major mountain chains, which generated geographic bar- Figs. S5, S6). riers and created open niches in high-altitude and low- altitude habitats that favored diversification (McCormack Discussion et  al. 2008). Biodiversity was also influenced by the cli - Genetic variation and phylogeographic pattern matic changes of Pleistocene climatic oscillations, which We found both deep and shallow genetic differentiation created a mosaic of habitats that resulted in isolated throughout the geographic distribution of A. rutila. The populations. lineages corresponding to the Pacific clade (pp = 0.99), The Tres Marías phylogroup is embedded in the larger Chiapas (pp = 0.99) and the Yucatán Peninsula-Central Pacific clade with a node origin dated 0.20 Myr BP (95% America (pp = 0.99) show deep structuring according to HPD 0.6‒0.44 Myr BP); this incomplete separation is genetic distances, F index, AMOVA, haplotype net- a pattern that has been reported in other birds at early ST work, multilocus phylogeny and mtDNA phylogenies. stages of speciation (e.g. Cortés-Rodríguez et  al. 2008). Within the Pacific clade, we can distinguish clear struc - STRU CTU RE analysis did not recover substructure for ture for Tres Marías Islands individuals with a monophy- Tres Marías individuals or any Pacific samples. However, letic clade (pp = 0.94; geographically isolated), F index, the ∆K selection method has been previously shown to ST one haplotype for ND2 locus, and one haplotype for COI underestimate population structure (Janes et  al. 2017). (shared with two samples from Guerrero); however, the Results of the AMOVA and F index suggest incipient ST structure analysis shows mixed ancestries between the population structure, and haplotype networks showed Tres Marías Islands and Continental samples. Similarly, a clear differentiation at the population level. Addition - a multilocus phylogeny suggests a shallow substructure ally, the Tres Marías group was the most morphologi- in the Yucatán Peninsula-Central America clade for West cally distinct compared to the other groups. The Pacific Nicaraguan samples (pp = 0.99), which share ancestry clade is represented as a polytomy in the resulting phy- with some samples from Yucatán. logenetic tree, and some clades are recovered within Despite minor differences in median-joining networks it but with no node support (see Additional file  2: Fig. based on the ND2 versus COI markers (Fig.  1), the geo- S1). As mentioned, the single valid clade is composed graphic groups had no shared haplotypes, with the of the Tres Marías Islands individuals. There are sev - exception of one haplotype that was shared between the eral examples of bird taxa in this area (The Pacific coast Tres Marías group and three individuals from the Pacific of Mexico and Tres Marías Islands) which show similar group (COI marker). The information from nuclear genes patterns of geographic variation and population struc- that was included to construct the multilocus phyloge- ture (e.g. Habia rubica: Ramírez-Barrera et al. 2018; Pha- netic tree (see Additional file  2: Fig. S1) did not weaken ethornis longirostris‒P. mexicanus: Arbeláez-Cortés and the pattern found using mtDNA; the multilocus phy- Navarro-Sigüenza 2013; Vireo hypochryseus: Arbeláez- logeny maintained the same geographic structure with Cortés et  al. 2014). Those studies found that the main high node support values (posterior probabilities and forces acting as drivers of divergence were the presence bootstrap values), even of the retained ancestral poly- of geographic and ecological barriers, as well as historical morphisms, favored by nuclear signal and its reduced habitat fragmentation. This is not the case for A. rutila substitution rates compared to mtDNA (Moore 1995). along the Pacific coast, since no structured pattern was A. rutila belongs to the third youngest subclade of nine found. There are two possible explanations for this lack Trochilidae subfamilies, known as the Emeralds (Blei- of geographic substructure. First, small sample sizes may weiss et al. 1997; McGuire et al. 2007, 2014). The tempo cause some haplotypes to be missing, since we found and mode of evolution in the hummingbird higher clades many intermediate haplotypes. Another explanation is has been considered to be the product of a series of his- an incipient substructure produced by short periods of torical events including colonizations, extinctions, recol- isolation during the Pleistocene. The age of the Pacific onizations, and the influence of mountain chain uplift, clade is 1.65 Myr BP (95% HPD 0.81‒1.81), and the age resulting in the radiation of highland species and the ranges estimated for subclades are 1.32 to 0.20 Myr BP, establishment of lowland taxa (Bleiweiss 1998; McGuire which corresponds with Pleistocene cycles. Thus, these et al. 2014). Therefore, timing is crucial to understanding mixed haplotypes could be the result of short periods of phylogeographic patterns in A. rutila. This complex split isolation during the LGM with secondary contact at the from its sister group 7.05 Myr BP (95% HPD 5.88‒8.29 LIG and the present. This idea is supported by the LGM Myr BP), and the four main clades reported here arose Pacific projections (see Additional file  5: Figs. S5, S6) 3.99‒0.20 Myr BP. During the late Miocene, Pliocene and and BSP analysis (Additional file  6: Fig. S7). In contrast, V ázquez‑López et al. Avian Res (2021) 12:61 Page 13 of 17 STRU CTU RE hierarchical analysis did not show different shows a mixed ancestry between those samples and some groups within the Pacific Clade (Fig. 2). Yucatán individuals (Fig. 2). The Chiapas clade is well defined, supported by genetic The influence of Pleistocene climatic fluctuations distance (5.69%) and F values (0.97‒1.00). This region on processes of diversification and incipient specia - ST east of the Isthmus of Tehuantepec is mainly influenced tion have been reported for many birds (e.g. Euphonia by historical events between the middle and late Plio- affinis : Vázquez-López et  al. 2020), including humming- cene and marine transgressions and regressions during birds (e.g. Eugenes fulgens: Zamudio-Beltrán et al. 2020). the late Pliocene (Barber and Klicka 2010; Peterson et al. These habitat fluctuations are well documented as fac - 2010). In general, the Isthmus of Tehuantepec represents tors that promote population contractions, limiting gene a geographic barrier for highland bird species (e.g. Zamu- flow (refugia hypothesis; Arango et  al. 2021), as well as dio-Beltrán and Hernández-Baños 2018), but not for low- expansions, resulting in zones of contact between distant land species. However, some phylogeographic studies of populations, with current signals of shared haplotypes lowland birds have described the influence of this low and polyphyletic patterns (e.g. Mesoamerican Amazilia altitude zone (e.g. Campylorhynchus rufinucha: Vázquez- complex: Jiménez and Ornelas 2016; White-chested Miranda et  al. 2009). According to our time-calibrated Amazilia complex: Rodríguez-Gómez and Ornelas tree, the Chiapas clade separated during Pleistocene cli- 2018). Diversifications occurring throughout the Mio - matic fluctuations ~ 760 kyr BP (95% HPD 0.43‒1.15 Ma), cene, when the estimates of divergence dates of the genus as there was a disruption of suitable conditions during the Amazilia are presumed ca. 19.85 Mya (i.e. Ornelas et al. Last Glacial Maxima across the Isthmus of Tehuantepec 2014), may have depended less on these climatic condi- (see Fig.  5). The Yucatán Peninsula and Central America tions. However, for events traced after this period, the clade also represents an independent entity, separated by rapid changes in climate cycling were crucial (Bleiweiss 57 and 15 mutational steps, respectively, in the ND2 and 1998). Furthermore, phylogeographic patterns of species COI trees (see Fig.  1), assuming low levels of gene flow, sharing (totally or partially) the geographic distribution with an origin 1.98 Myr BP (95% HPD 1.42‒2.6 Myr BP). of A. rutila show similarities (e.g. Turdus rufopalliatus: This genetic signal of isolation from the other groups Montaño-Rendón et  al. 2015; Ortíz-Ramírez et  al. 2018; may be influenced by the geographic distance and dif - de Silva et  al. 2020). In general, the diversification pro - ferences in environmental spaces. The Pacific and Chia - cesses have resulted from a variety of events, sharing the pas group are distributed along the tropical deciduous explanations of an intricate orographic history, climatic forest belt, which has marked dry and rainy seasons and fluctuations during the Pleistocene, and thus low signals has been repeatedly identified as a distinct biogeographi - of gene flow favored by the isolation of populations due cal province (i.e. Gordon and Ornelas 2000), while the to distance or historical refuges. Yucatán Peninsula has moderately pronounced dry and rainy seasons with relatively constant average tempera- Morphometrics ture throughout the year (Paynter 1955). Also, our ENM Our results showed little variation among the continental past projections suggest periods of habitat contractions, phylogroups of this species in the morphological char- during which populations from the Yucatán Peninsula acters we measured, despite A. rutila’s large geographic and Central America could have been closer from each range. However, the Tres Marías Islands phylogroup is other (LGM MIROC, see Additional file  5: Figs. S5, S6), larger than the continental phylogroups in four out of followed by periods of isolation. Although A. rutila is five measurements (Fig.  3). This is consistent with several allopatrically distributed in the Yucatán Peninsula and other avian species from Tres Marías (Grant 1965) and Central America, the studied individuals were placed is a clear example of the “Island Rule” (van Valen 1973; in a single group. This link between individuals close to Lomolino 1985; Clegg and Owens 2002; Boyer and Jetz the Yucatán Peninsula and Central America was also 2010, Benítez-López et  al. 2021). This shift in body size found in A. tzacatl (Miller et al. 2011), in which popula- has been attributed to changes in thermal ecology (Clegg tions in Tabasco (Mexico), Belize, Honduras, and Costa and Owens 2002), as the Tres Marías Islands are signifi - Rica shared several haplotypes despite the long distances cantly drier and warmer than the continental habitat of between them, and the most differentiated populations the species. However, dramatic divergence in body size were the southernmost populations in Panama. The can be the result of geographic isolation independent of periods of isolation mentioned above combined with environmental pressure (Seeholzer and Brumfield 2018). demographic decreases (see BSP results) could explain This divergence in body size with relatively similar beak the shallow substructure found in the West Nicaragua sizes suggest that beak morphology might be more tightly samples (pp = 0.99). However, our STRU CTU RE analysis constrained than body size due to this species’ feeding ecology. Vázquez‑López et al. Avian Res (2021) 12:61 Page 14 of 17 Implications for taxonomy name, but specimens from Huehuetán, Chiapas were We present a taxonomic proposal for the A. rutila com- described under A. cinnamomea saturata by Nelson plex based on the recovery of four independent clades in 1898, and this name was later considered a synonym with high support values corresponding to natural geo- of A. r. corallirostris (Ridgway 1901). As a consequence, graphic groups. the Chiapas lineage could be named A. saturata (Nelson Our results support the first differentiated clade, A. r. 1898). graysoni, distributed in the Tres Marías Islands (Nayarit, Mexico). This differentiation was also confirmed by our morphometric analyses, which showed that the Tres Conclusions Marías Islands population is the largest morphotype. A. The evolutionary hypothesis of diversification within this r. graysoni was first suggested as an independent species group is of geographic isolation (limited gene flow), dif - by Ridgway (1901), who described Amazilia graysoni ferences in current environmental conditions, and his- (Grayson’s Hummingbird) as a separate species similar torical habitat fragmentation promoted by past scenarios in coloration to A. rutila, but darker and much larger. (Pleistocene refuges). The A. rutila complex consists of Then following the alternative phylogenetic/evolution - four well-defined clades that, in our opinion, merit taxo - ary species taxonomy of the Mexican avifauna, Navarro- nomic reevaluation. Our data suggest that A. r. graysoni Sigüenza and Peterson (2004) proposed A. r. graysoni as (Tres Marías Islands) as well as A. r. rutila (Pacific coast) an independent evolutionary species, which is morpho- should be considered full species. The other two strongly logically distinct in size and coloration from the rest of supported clades are: (a) the Chiapas group (southern the subspecies within A. rutila. The proposal of elevating Mexico), and (b) the populations from Yucatán Peninsula this taxon to species level has been discussed recently and Central America. These clades belong to the A. rutila (de Silva et  al. 2020). The destruction of native habitat, corallirostris taxon as currently defined. We propose that the introduction of exotic species, and reduced popula- this group be split and renamed as two species: A. satu- tion sizes increase the vulnerability of birds on the Pacific rata for the Chiapas group and A. corallirostris for the Islands of Mexico and must be immediately addressed Yucatán Peninsula and Central America group. in conservation policies (Hahn et  al. 2012; de Silva et  al. 2020). Supplementary Information The second differentiated clade corresponds to the The online version contains supplementary material available at https:// doi. org/ 10. 1186/ s40657‑ 021‑ 00295‑0. group distributed along the Pacific coast. Within this second group, two subspecies have been described (A. r. Additional file 1: Table S1. Information of voucher numbers (ID), geo ‑ rutila and A. r. diluta), which, according to our results, graphic locations (longitude, latitude), and biological collections of tissue are not observable as separate groups. Therefore, we pro - samples used in this study. Table S2. Information of accession numbers pose that the Pacific coast phylogroup be represented by (AN) and geographic locations (locality, longitude, latitude) of GenBank’s sequences used in this study. A. r. rutila (the first named). The third phylogroup is distributed in Chiapas, Mex - Additional file 2: Figure S1. The Bayesian Inference and Maximum Likeli‑ hood of the mitochondrial DNA and nuclear DNA. ico (east of the Isthmus of Tehuantepec). This group Additional file 3: Figure S2. STRU CTU RE analysis of Amazilia rutila com‑ corresponds to A. r. corallirostris, which is distributed plex. Figure S3. Hierarchically structure analysis. from south and southeastern Mexico to Costa Rica. The Additional file 4: Figure S4. Statistics obtained for the selection of biocli‑ Chiapas phylogrouop is distributed in an important and matic variables (19 variables from WorldClim) that were used on Ecological well-studied main fauna region, considered an area of Niche Modelling projections, for the whole complex (Amazilia rutila), and each geographic group (Pacific, Chiapas, Yucatán_CA). endemism in western Mexico with one of the highest val- ues of avian species richness (Peterson and Navarro 2000; Additional file 5: Figure S5. Ecological Niche Modelling for Amazilia rutila under different past scenarios: LIG (Last Inter Glacial), LGM (Last Gla‑ García-Trejo and Navarro-Sigüenza 2004). cial Maxima), and MH (Mid Holocene) periods. Figure S6. Ecological Niche Finally, the fourth independent group corresponds to Modelling for each geographic group of Amazilia rutila under different the clade of the Yucatán Peninsula—Central America, the past scenarios: LIG (Last Inter Glacial), LGM (Last Glacial Maxima), and MH (Mid Holocene) periods. most genetically differentiated (7.52%), which also cor - Additional file 6: Figure S7. Bayesian skyline plots (BSP) for geographic responds to A. r. corallirostris. According to our results, groups (Pacific, Chiapas, Yucatán_CA), and for the A. rutila complex. A. r. corallirostris is a polyphyletic subspecies, since it is composed of two independent groups. A. r. corallirostris Acknowledgements was described by Bourcier and Mulsant (1846) with spec- We thank M. Robbins of the University of Kansas (Natural History Museum), imens from Escuintla, Guatemala. Therefore, the name A. and S. Birks of the University of Washington (Burke Museum of Natural History corallirostris should be applied to the Yucatán-CA clade and Culture) for providing tissue samples used in this study. We also thank Alejandro Gordillo Martínez, Fabiola Ramírez, Isabel Vargas‑Fernandez and (Ridgway 1901). This left the Chiapas group without a V ázquez‑López et al. Avian Res (2021) 12:61 Page 15 of 17 Bleiweiss R. Tempo and mode of hummingbird evolution. Biol J Linn Soc. Raúl Ivan Martínez Becerril for technical support. We acknowledge the General 1998;65:63–76. Direction of Wildlife (Instituto de Ecología, SEMARNAT, México) for providing Bleiweiss R, Kirsch JA, Matheus JC. DNA hybridization evidence for the collecting permits. L. Kiere reviewed the English. principal lineages of hummingbirds (Aves: Trochilidae). Mol Biol Evol. 1997;14:325–43. Authors’ contributions Booth TH, Nix HA, Busby JR, Hutchinson MF. BIOCLIM: the first species distribu‑ BEH‑B and MV ‑L designed the study. BEH‑B secured financial support. MV ‑L tion modelling package, its early applications and relevance to most carried out the laboratory work. MV‑L, BEH‑B, SMR‑B and LEZ ‑B analyzed the current MAXENT studies. Divers Distrib. 2014;20:1–9. data. BEH‑B, MV ‑L, NC‑R, AB‑H, LEZ ‑B and KR contributed to the writing and Bouckaert R, Vaughan TG, Barido‑Sottani J, Duchêne S, Fourment M, Gavry‑ improvement of the manuscript. All authors read and approved the final ushkina A, et al. BEAST 2.5: an advanced software platform for Bayesian manuscript. evolutionary analysis. PLoS Comput Biol. 2019;15:e1006650. Bourcier J, Mulsant E. Description de vingt espèces nouvelles d’oiseaux‑ Funding mouches. Ann Sci Phys Nat Agric. 1846;9:312–32. This research was supported by PAPIIT/DGAPA, Universidad Nacional Boyer AG, Jetz W. Biogeography of body size in Pacific island birds. Ecography. Autónoma de México (UNAM) through a grant to Blanca E. Hernández‑Baños 2010;33:369–79. (IN220620). LEZ‑B acknowledges the Postdoctoral scholarship provided by Braconnot P, Otto‑Bliesner B, Harrison S, Joussaume S, Peterchmitt JY, Abe ‑ DGAPA‑UNAM. Ouchi A, et al. Results of PMIP2 coupled simulations of the Mid‑Holo ‑ cene and Last Glacial Maximum‑Part 1: experiments and large ‑scale Availability of data and materials features. Clim Past. 2007;3:261–77. The sequences generated during the current study are available in GenBank Cavender‑Bares J, Gonzalez‑Rodriguez A, Pahlich A, Koehler K, Deacon N. (accession numbers: MZ998668‑MZ998740, OK624605‑ OK624657, OK614044‑ Phylogeography and climatic niche evolution in live oaks (Quercus OK614080, OK618334‑ OK618361). series Virentes) from the tropics to the temperate zone. J Biogeogr. 2011;38:962–81. Chaves AV, Clozato CL, Lacerda DR, Sari H, Santos FR. Molecular taxonomy Declarations of Brazilian tyrant‑flycatchers (Passeriformes, Tyrannidae). Mol Ecol. 2008;8:1169–77. Ethics approval and consent to participate Clegg SM, Owens PF. The “island rule” in birds: medium body size and its eco‑ This study was performed in line with the principles of the Declaration of logical explanation. Proc R Soc B Biol Sci. 2002;269:1359–65. Helsinki. Cobos ME, Osorio‑ Olvera L, Soberón J, Peterson AT. ellipsenm: ecological niche characterizations using ellipsoids. R package; 2020. https:// github. com/ Consent for publication marlo necob os/ ellip senm. Not applicable. 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: Competing interests low mitochondrial DNA differentiation reveals recent divergence. The authors declare that they have no competing interests. Condor. 2008;4:729–39. de Silva HG, Pérez Villafaña MG, Cruz‑Nieto J, Cruz‑Nieto MÁ. Are some of Author details 1 the birds endemic to the Tres Marías Islands (Mexico) species? Bull Br Museo de Zoología, Facultad de Ciencias, Universidad Nacional Autónoma de 2 Ornithol Club. 2020;140:7–37. México, Apartado Postal, 04510 Mexico City, Mexico. Department of Biol‑ 3 DeLattre PA. Oiseaux‑mouches nouveaux au peu connus, découverts au ogy, Ithaca College, Ithaca, NY 14850, USA. Museo de Zoología, Facultad de Guatimala. L’echo Du Monde Savant. 1843;7:1068–70. Estudios Superiores Zaragoza, Universidad Nacional Autónoma de México, 4 Drummond AJ, Rambaut A, Shapiro B, Pybus OG. Bayesian coalescent infer‑ Mexico City, Mexico. Department of Biology, Colorado State University, Fort ence of past population dynamics from molecular sequences. Mol Biol Collins, CO 80523, USA. Evol. 2005;22:1185–92. Earl DA, vonHoldt BM. STRU CTU RE HARVESTER: a website and program for Received: 19 February 2021 Accepted: 31 October 2021 visualizing STRU CTU RE output and implementing the Evanno method. Conserv Genet Res. 2012;4:359–61. Ellegren H. Molecular evolutionary genomics of birds. Cytogenet Genome Res. 2007;117:120–30. Evanno G, Regnaut S, Goudet J. Detecting the number of clusters of individu‑ References als using the software STRU CTU RE: a simulation study. Mol Ecol. Arango A, Villalobos F, Prieto‑ Torres DA, Guevara R. The phylogenetic diversity 2005;14:2611–20. and structure of the seasonally dry forest in the Neotropics. J Biogeogr. Excoffier LG, Schneider S. Arlequin v. 3.1: an integrated software package for 2021;48:176–86. population genetics data analysis. Evol Bioinform Online. 2006;1:47–50. Arbeláez‑ Cortés E, Navarro‑Sigüenza AG. Molecular evidence of the taxonomic Excoffier L, Smouse PE, Quattro JM. Analysis of molecular variance inferred status of western Mexican populations of Phaethornis longirostris (Aves: from metric distances among DNA haplotypes: application to human Trochilidae). Zootaxa. 2013;3716:81–97. mitochondrial restriction region. Genetics. 1992;131:479–91. Arbeláez‑ Cortés E, Roldán‑Piña D, Navarro ‑Sigüenza AG. Multilocus phylo ‑ Fu Y. Statistical tests of neutrality of mutations against population growth, geography and morphology give insights into the recent evolution hitchhiking and background selection. Genetics. 1997;147:915–25. of a Mexican endemic songbird: Vireo hypochryseus. J Avian Biol. García‑Deras GM, Cortés‑Rodríguez N, Honey M, Navarro ‑Sigüenza AG, García‑ 2014;45:253–63. Moreno J, Hernández‑Baños BE. Phylogenetic relationships within the Bandelt HJ, Foster P, Röhl A. Median‑joining networks for inferring intraspecific genus Cynanthus (Aves: Trochilidae), with emphasis on C. doubledayi. phylogenies. Mol Biol Evol. 1999;16:37–48. Zootaxa. 2008;1742:61–8. Barber BR, Klicka J. Two pulses of diversification across the Isthmus of García‑ Trejo EA, Navarro‑Sigüenza AG. Patrones biogeográficos de la riqueza Tehuantepec in a montane Mexican bird fauna. Proc R Soc B Biol Sci. de especies y el endemismo de la avifauna en el oeste de México. Acta 2010;277:2675–81. Zool Mex. 2004;20:167–85. Benítez‑López A, Santini L, Gallego ‑Zamorano J, Milá B, Walkden P, Huijbregts Gill F, Donsker D, Rasmussen P. IOC World Bird List (v.11.1). 2021. https:// www. MA, et al. The island rule explains consistent patterns of body size evo‑ world birdn ames. org/. lution in terrestrial vertebrates. Nat Ecol Evol. 2021;5:768–86. González C, Ornelas JF, Gutiérrez‑Rodríguez C. Selection and geographic Billerman SM, Keeney BK, Rodewald PG, Schulenberg TS. Birds of the world. isolation influence hummingbird speciation: genetic, acoustic and Ithaca: Cornell Laboratory of Ornithology; 2020. https:// birds ofthe world. org/ bow/ home. Vázquez‑López et al. Avian Res (2021) 12:61 Page 16 of 17 morphological divergence in the wedge‑tailed sabrewing (Campylop - Miller MA, Pfeiffer W, Schwartz T. Creating the CIPRES Science Gateway for terus curvipennins). BMC Evol Biol. 2011;11:38. inference of large phylogenetic trees. Gateway computing environ‑ Gordon CE, Ornelas JF. Comparing endemism and habitat restriction in Mesoa‑ ments workshop. New Orleans: IEEE; 2010. merican tropical deciduous forest birds: implications for biodiversity Montaño‑Rendón M, Sánchez‑ González LA, Hernández‑Alonso G, Navarro ‑ conservation planning. Bird Conserv Int. 2000;10:289–303. Sigüenza AG. Genetic differentiation in the Mexican endemic Rufous‑ Grant PR. The adaptive significance of some size trends in island birds. Evolu‑ backed Robin, Turdus rufopalliatus (Passeriformes: Turdidae). Zootaxa. tion. 1965;19:355–67. 2015;4034:495–514. Hahn IJ, Hogeback S, Römer U, Vergara PM. Biodiversity and biogeography Moore WS. Inferring phylogenies from mtDNA variation: mitochondrial‑ gene of birds in Pacific Mexico along an isolation gradient from mainland trees versus nuclear‑ gene trees. Evolution. 1995;49:718–26. Chamela via coastal Marías to oceanic Revillagigedo Islands. Vertebr Morrone JJ. Biogeographical regionalisation of the Neotropical region. Zool. 2012;62:123–44. Zootaxa. 2014;3782:1–110. Hall TA. BioEdit: a user‑friendly biological sequence alignment editor and Navarro‑Sigüenza AG, Peterson AT. An alternative species taxonomy of the analysis program for Windows 95/98/NT. Nucleic Acids Symp Ser. birds of Mexico. Biota Neotrop. 2004;4:1–32. 1999;41:95–8. Nelson EW. Descriptions of new birds from the Tres Marias Islands, Western Hernández‑Baños BE, Zamudio ‑Beltrán LE, Milá B. Phylogenetic relationships Mexico. Proc Biol Soc Washington. 1898;12:63‑4. and systematics of a subclade of Mesoamerican emerald humming‑ Ornelas JF, González C, de los Monteros AE, Rodríguez‑ Gómez F, García‑Feria birds (Aves: Trochilidae: Trochilini). Zootaxa. 2020;4748:581–91. LM. In and out of Mesoamerica: temporal divergence of Amazilia hum‑ Hijmans RJ, Cameron SE, Parra JL, Jones PG, Jarvis A. Very high resolution mingbirds pre‑ dates the orthodox account of the completion of the interpolated climate surfaces for global land areas. Int J Climatol. Isthmus of Panama. J Biogeogr. 2014;41:168–81. 2005;25:1965–78. Ortiz‑Ramírez MF, Sánchez‑ González LA, Castellanos‑Morales G, Ornelas JF, Hijmans RJ. raster: geographic data analysis and modeling. R package. version Navarro‑Sigüenza AG. Concerted Pleistocene dispersal and genetic dif‑ 3.0‑12. 2020. https:// CRAN.R‑ proje ct. org/ packa ge= raster . ferentiation in passerine birds from the Tres Marías Archipelago, Mexico. Howell SN, Webb S. A guide to the birds of Mexico and northern Central Auk. 2018;135:716–32. America. Oxford: Oxford University Press; 1995. Osorio‑ Olvera L, Barve V, Barve N, Soberón J. nichetoolbox: from getting Hubisz M, Falush D, Stephens M, Pritchard J. Inferring weak population struc‑ biodiversity data to evaluating species distribution models in a friendly ture with the assistance of sample group information. Mol Ecol Resour. GUI environment. R package. version 0.2.0.0. 2016. http:// shiny. conab io. 2009;9:1322–32.gob. mx: 3838/ niche toolb2/. Huelsenbeck JP, Ronquist F. MrBayes: a program for the Bayesian inference of Pacheco MA, Battistuzzi FU, Lentino M, Aguilar RF, Kumar S, Escalante AA. Evo‑ phylogeny. Bioinformatics. 2002;17:754–5. lution of modern birds revealed by mitogenomics: timing the radiation Janes JK, Miller JM, Dupuis JR, Malenfant RM, Gorrell JC, Cullingham CI, et al. and origin of major orders. Mol Biol Evol. 2011;28:1927–42. The K = 2 conundrum. Mol Ecol. 2017;26:3594–602. Ridgway R. Part V: Family Trochilidae in: The birds of North and Middle Jiménez A, Ornelas JF. Historical and current introgression in a Mesoamerican America: a descriptive catalogue. Washington: Bulletin of the United hummingbird species complex: a biogeographic perspective. PeerJ. States National Museum No. 50; 1901. 2016;4:e1556. Paynter RA. The ornithogeography of the Yucatán Peninsula. New Haven: Yale Koleff P, Urquiza‑Haas T, Contreras B. Prioridades de conservación de los University; 1955. bosques tropicales en México: reflexiones sobre su estado de conser ‑ Peterson AT, Navarro AG. Western Mexico: a significant center of avian end‑ vación y manejo. Ecosistemas. 2012;21:6–20. emism and challenge for conservation action. Cotinga. 2000;14:42–6. Lanfear R, Calcott B, Ho SY, Guindon S. PartitionFinder: combined selection Peterson AT, Navarro‑Sigüenza AG, Li X. Joint effects of marine intrusion of partitioning schemes and substitution models for phylogenetic and climate change on the Mexican avifauna. Ann Assoc Am Geogr. analyses. Mol Biol Evol. 2012;29:1695–701. 2010;100:908–16. Lanfear R, Frandsen PB, Wright AM, Senfeld T, Calcott B. PartitionFinder 2: new Phillips SJ, Dudík M, Schapire RE. Maxent software for modeling species niches methods for selecting partitioned models of evolution for molecular and distributions. version 3.4.1. 2021. and morphological phylogenetic analyses. Mol Biol Evol. 2017;34:772–3. Pritchard JK, Stephens M, Donnelly P. Inference of population structure using Lawrence GN. Annals of the Lycaeum of natural history of New York. New York: multilocus genotype data. Genetics. 2000;155:945–59. Lyceum of Natural History; 1866. R Core Team. R: a language and environment for statistical computing. Vienna: Lerdau M, Whitbeck J, Holbrook NM. Tropical deciduous forest: death of a R Foundation for Statistical Computing; 2020. https:// www.R‑ proje ct. biome. Trends Ecol Evol. 1991;6:201–2. org/. Lerner HRL, Meyer M, James HF, Hofreiter M, Fleischer RC. Multilocus resolution Rambaut A, Drummond AJ, Xie D, Baele G, Suchard MA. Posterior summariza‑ of phylogeny and timescale in the extant adaptive radiation of Hawai‑ tion in Bayesian phylogenetics using Tracer 1.7. Syst Biol. 2018;67:901–4. ian honeycreepers. Curr Biol. 2011;21:1838–44. Rambaut A, Drummond AJ. TreeAnnotator v.2.6.2:MCMC Output analysis. Librado P, Rozas J. DnaSP v5: a software for comprehensive analysis of DNA 2013; http:// beast. bio. ed. ac. uk. polymorphism data. Bioinformatics. 2009;25:1451–2. Rambaut A. FigTree v1.4.2: Tree figure drawing tool. 2014; http:// tree. bio. ed. Lomolino MV. Body size of mammals on islands: the island rule reexamined. ac. uk. Am Nat. 1985;125:310–6. Ramírez‑Barrera SM, Hernández‑Baños BE, Jaramillo ‑ Correa JP, Klicka J. Deep McCormack JE, Peterson AT, Bonaccorso E, Smith TB. Speciation in the high‑ divergence of Red‑ crowned Ant Tanager (Habia rubica: Cardinalidae), a lands of Mexico: genetic and phenotypic divergence in the Mexican jay multilocus phylogenetic analysis with emphasis in Mesoamerica. PeerJ. (Aphelocoma ultramarina). Mol Ecol. 2008;17:2505–21. 2018;6:e5496. McGuire JA, Witt CC, Altshuler DL, Remsen JV Jr. Phylogenetic systematics and Rodríguez‑ Gómez F, Ornelas JF. Genetic structuring and secondary contact biogeography of hummingbirds: Bayesian and maximum likelihood in the white‑ chested Amazilia hummingbird species complex. J Avian analyses of partitioned data and selection of an appropriate partition‑ Biol. 2018;49:e01536. ing strategy. Syst Biol. 2007;56:837–56. Ryan RM. The biotic provinces of Central America. Acta Zool Mex. 1963;6:1–54. McGuire JA, Witt CC, Remsen JV Jr, Corl A, Rabosky DL, Altshuler DL, et al. Rzedowski J. Vegetación de México. México: Limusa; 1978. Molecular phylogenetics and the diversification of hummingbirds. Curr Savage JM. The origins and history of the Central America herpetofauna. Biol. 2014;24:910–6. Copeia. 1966;4:719–66. Meave JA, Romero‑Romero MA, Salas‑Morales SH, Pérez‑ García EA, Gallardo‑ Seeholzer GF, Brumfield RT. Isolation by distance, not incipient ecological spe ‑ Cruz JA. Diversidad, amenazas y oportunidades para la conservación ciation, explains genetic differentiation in an Andean songbird (Aves: del bosque tropical caducifolio en el estado de Oaxaca, México. Ecosis‑ Furnariidae: Cranioleuca antisiensis, Line‑ cheeked Spinetail) despite near temas. 2012;21:85–100. threefold body size change across an environmental gradient. Mol Ecol. Miller MJ, Lelevier MJ, Bermingham E, Klicka JT, Escalante P, Winker K. Phyloge‑ 2018;27:279–96. ography of the Rufous‑tailed Hummingbird (Amazilia tzacatl). Condor. 2011;113:806–16. V ázquez‑López et al. Avian Res (2021) 12:61 Page 17 of 17 Smith H. Las provincias bióticas de México, según la distribución geográfica de Weller AA. Cinnamon Hummingbird. In: del Hoyo JJ, Elliott A, Sargatal J, edi‑ las lagartijas del género Sceloporus. Anales De La Escuela Nacional De tors. Handbook of the birds of the world. Barn‑ owls to hummingbirds. Ciencias Biológicas. 1941;2:103–10. Barcelona: Lynx Editions; 1999. p. 596‒597. Smith BT, Escalante P, Hernández‑Baños BE, Navarro ‑Sigüenza AG, Rohwer S, Werneck FP, Costa GC, Colli GR, Prado DE, Sites JW. Revisiting the historical Klicka J. The role of historical and contemporary processes on phylo‑ distribution of Seasonally Dry Tropical Forests: new insights based on geographic structure and genetic diversity in the Northern Cardinal, palaeodistribution modelling and palynological evidence. Global Ecol Cardinalis Cardinalis. BMC Evol Biol. 2011;11:136. Biogeogr. 2011;20:272–88. Sorenson MD, Ast JC, Dimcheff DE, Yuri T, Mindell DP. Primers for a PCR‑based West RC. The natural regions of Middle America. In: West RC, editor. Handbook approach to mitochondrial genome sequencing in birds and other of Middle American Indians. Austin: University of Texas Press; 1964. p. vertebrates. Mol Phylogenet Evol. 1999;12:105–14. 363–83. Stamakis A. RaxML version 8: a tool for phylogenetic analysis and post‑analysis Wickham H, Francois R, Henry L, Müller K. dplyr: A Grammar of Data Manipula‑ of large phylogenies. Bioinformatics. 2014;30:1312–3. tion. R package. version 0.8.5. 2020. https:// CRAN.R‑ proje ct. org/ packa Stecher G, Tamura K, Kumar S. Molecular evolutionary genetics analysis ge= dplyr. (MEGA) for macOS. Mol Biol Evol. 2020;37:1237–9. Wilcox TP, Zwickl DJ, Heath TA, Hillis DM. Phylogenetic relationships of the Stephens M, Scheet P. Accounting for decay of linkage disequilibrium in dwarf boas and a comparison of Bayesian and bootstrap measures of haplotype inference and missing‑ data imputation. Am J Human Genet. phylogenetic support. Mol Phylogenet Evol. 2002;25:361–71. 2005;76:449–62. Zamudio‑Beltrán LE, Hernández‑Baños BE. Genetic and morphometric diver ‑ Stephens M, Smith N, Donelly P. A new statistical method for haplotype recon‑ gence in the Garnet‑ Throated Hummingbird Lamprolaima rhami (Aves: struction from population data. Am J Human Genet. 2001;68:978–89. Trochilidae). PeerJ. 2018;6:e5733. Tajima F. Statistical method for testing the neutral mutation hypothesis by Zamudio‑Beltrán LE, Ornelas JF, Malpica A, Hernández‑Baños BE. Genetic and DNA polymorphism. Genetics. 1989;123:585–95. morphological differentiation among populations of the Rivoli’s Hum‑ Thompson JD, Higgins DG, Gibson TJ. CLUSTAL W: improving the sensitivity of mingbird (Eugenes fulgens) species complex (Aves: Trochilidae). Auk. progressive multiple sequence alignment through sequence weight‑ 2020;137:ukaa032. ing, position‑specific gap penalties and weight matrix choice. Nucleic Zizka A, Silvestro D, Andermann T, Azevedo J, Duarte‑Ritter C, Edler D, et al. Acids Res. 1994;22:4673–80. CoordinateCleaner: standardized cleaning of occurrence records from van Rossem AJ. A northwest race of the Cinnamon Hummingbird. Condor. biological collection databases. Method Ecol Evol. 2019;10:744–51. 1938;40:226–7. van Valen LM. Pattern and the balance of nature. Evol Theor. 1973;1:31–49. 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. 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. Weir JT. Divergent timing and patterns of species accumulation in lowland and highland neotropical birds. Evolution. 2006;60:842–55. Re Read ady y to to submit y submit your our re researc search h ? Choose BMC and benefit fr ? Choose BMC and benefit from om: : fast, convenient online submission thorough peer review by experienced researchers in your field rapid publication on acceptance support for research data, including large and complex data types • gold Open Access which fosters wider collaboration and increased citations maximum visibility for your research: over 100M website views per year At BMC, research is always in progress. Learn more biomedcentral.com/submissions

Journal

Avian ResearchSpringer Journals

Published: Nov 11, 2021

Keywords: Amazilia rutila; Genetic structure; Genetic variation; Phylogeography; Tres Marías Islands; Trochilidae; Tropical deciduous forest

There are no references for this article.