Genetic and morphometric divergence in the Garnet-Throated Hummingbird Lamprolaima rhami (Aves: Trochilidae)

Taxon sampling and sequencing

We obtained tissues from 54 individuals of L. rhami from 14 localities across most of its geographic range (Table S1, Fig. 1, using field collecting permit from Instituto Nacional de Ecología, SEMARNAT: FAUT-0169). We defined five groups a priori to evaluate genetic and morphological variation among the five allopatric regions where the species occurs: (1) the Sierra Madre Oriental (SMO), (2) the highlands of Guerrero (GRO), (3) the Sierra of Miahuatlan in Oaxaca (MIA), (4) the highlands of Chiapas and Guatemala (CHIS), and (5) the region in Central America comprising the highlands of Honduras and El Salvador (CA). Tissue samples were obtained for four geographic groups (excluding the CA group) from the following collections: “Museo de Zoología Alfonso L. Herrera” (Universidad Nacional Autónoma de México), Museum of Natural Science (Lousiana State University), and Museum of Vertebrate Zoology (University of California, Berkeley).

DNA was extracted using the DNAeasy™ kit (Qiagen Inc., Valencia, CA, USA), following the manufacturer’s protocols. To evaluate genetic variation in the complex, two mitochondrial markers were obtained from the 54 samples (Control Region, CR; and subunits 6 and 8 from ATPase gene, ATPase 6 & 8). To evaluate phylogenetic relationships and for species delimitation analyses between groups, two additional mitochondrial markers and four nuclear regions were surveyed in a subsample of 31 individuals (NADH dehydrogenase subunit 2, ND2; NADH dehydrogenase subunit 4, ND4; the 7th intron of the beta fibrinogen gene, BFib, the regions between exons 4 and 5 of the Muscle Skeletal Receptor Tyrosine Kinase gene, MUSK; a segment comprising the end of exon 6 and the beginning of exon 8 of the Ornithine Decarboxylase gene, ODC, and intron 5 of adenylate kinase gene, AK1). We included sequences from the same molecular markers for Eugenes fulgens, E. spectabilis, and Tilmatura dupontii, to serve as outgroups (McGuire et al., 2007; Zamudio-Beltrán & Hernández-Baños, 2015).

We amplified these molecular markers via polymerase chain reaction (PCR) using specific primers and protocols (Table S3). Reactions contained 10×  buffer (1.25 µL), 10 mM dNTP (0.19 µL), 50 mM MgCl₂ (0.38 µL), 10 µM of each primer (0.25 µL), 0.1 µL of Taq (INVITROGEN), and 0.5 µL of genomic DNA (12.5 µL total volume). PCR products were visualized on a 1% agarose gel, and DNA sequencing was performed by the High-Throughput Genomics Unit Service of the University of Washington. We edited and aligned chromatograms with Sequencher v4.8 (GeneCodes Corporation, Ann Arbor, MI, USA). The multilocus alignment can be found in FigShare (https://figshare.com/projects/Lamprolaima_rhami_project_/34487).

Population structure

To evaluate the number of haplotypes and their relationships, a statistical parsimony haplotype network was constructed for the concatenated dataset of mitochondrial markers (CR and ATPase 6 & 8), using the program TCS v1.21 (Clement, Posada & Crandall, 2000).

To analyze genetic diversity and genetic structure, we obtained values of haplotype diversity, nucleotide diversity, mean number of pairwise differences, and population F_ST values. These analyses were performed with 1,000 replicates, using the program Arlequin v3.11 (Excoffier, Laval & Schneider, 2005). Using the same program, we conducted an analysis of molecular variance (AMOVA; Excoffier, Smouse & Quattro, 1992) to detect structure between populations based on comparisons between geographically defined groups. According to our results, regions on both sides of the Isthmus of Tehuantepec were also evaluated (IT, east and west).

To evaluate the isolation by distance among geographic regions, we performed a Mantel Test with 1,000 iterations, comparing matrices of genetic and geographic distances, using the program zt v1.1 (Bonnet & De Peer, 2002). Statistical analyses were not performed in the MIA group because of limited number of samples (n = 2), but were considered in haplotype networks and phylogenetic analyses.

Demographic analyses

To evaluate demography and population stability, we obtained Tajima’s D and Fu’s Fs values, in Arlequin v2.11 (Excoffier, Laval & Schneider, 2005), with 1,000 replicates (mtDNA database). Using the same program, parameters, and database, we further evaluated the historical demography of each group under an expansion model with a MISMATCH distribution test and estimated its significance with the raggedness index (Slatkin & Hudson, 1991; Rogers & Harpending, 1992; Harpending, 1994). To analyze variation in effective population size over time, we used Bayesian skyline plots (BSP; Drummond et al., 2005) performed in BEAST v1.6.0 (Drummond & Rambaut, 2007), with 10 million steps for mtDNA, using a mean rate of 0.023 substitutions per site per lineage per million years (s/s/l/My), under Control Region and ATPase estimates (Lerner et al., 2011).

Evolutionary models and phylogenetic analyses

For each molecular marker (mtDNA and nuclear DNA), we calculated the evolutionary model that best fit the data based on the Akaike Information Criterion AIC (Akaike, 1987) using jModelTest (Posada, 2008). We performed a phylogenetic reconstruction using the Bayesian Inference (BI) approach in Mr. Bayes v3.0 (Huelsenbeck & Ronquist, 2002). We assigned different evolutionary models to each gene partition. We ran four simultaneous chains for each Monte Carlo Markov Chain analysis for 50 million generations, sampling every 1000 generations. The results were visualized to ensure ESS (effective sample sizes) values higher than 200, and the burn-in value was determined using Tracer v1.6.0 (Rambaut, Suchard & Drummond, 2013). The initial 20% of generations were eliminated. The remaining trees were used to construct a majority rule consensus tree with posterior probability distributions, which was visualized using the program FigTree v1.2.3 (http://tree.bio.ed.ac.uk/software/figtree/).

Morphological variation

To examine morphological variation between groups of L. rhami, we took five measurements from 208 voucher specimens (Table S2) corresponding to four of the five geographic groups defined a priori (SMO, GRO, CHIS, CA). These specimens were available from the following collections: Museo de Zoología “Alfonso L. Herrera” (MZFC, UNAM), Museum of Comparative Zoology (MCZ), American Museum of Natural History (AMNH), Donald R. Dickey Bird and Mammal Collection (BMC), and the Moore Lab of Zoology (MLZ). Measurements for bill length (from the base to the tip of the upper mandible), bill width (width at the nostrils), bill depth (from the upper mandible to the base of the bill at the nostrils), and wing chord (distance from the carpal joint to the tip of the longest primary) were taken with a dial calliper with a precision of 0.1 mm, while tail length (distance from the uropigial gland to the tip of the longest rectrix) was determined with a milimetric ruler. A single observer took all measurements. Subsets of individuals were measured twice to confirm consistency between measurements using correlation coefficients obtained in STATISTICA v7 (StatSoft, Inc., Tulsa, OK, USA). When variation among measurements was low or null, all voucher specimens were measured once and these values were used in further analyses. To test the normality of our data, we performed a Lilliefors (Kolmogorov–Smirnov) test using the R package Nortest v1.0-4 (Gross & Ligges, 2015). This test was conducted with raw and log-transformated data. Since the data were not normally distributed, a Wilcoxon/Mann–Whitney test (Bauer, 1972) was performed to evaluate differences between males and females. To evaluate differences among groups, and after confirming differences among sexes, we conducted two sets of Kruskall-Wallis tests (Hollander & Wolfe, 1973) to compare: (1) geographic groups (SMO, GRO, CHIS, CA), and (2) groups separated by the Isthmus of Tehuantepec (east, west). All tests were conducted for each variable, treating males and females separately. Statistical analyses were performed using RStudio v.1.1.447 (RStudio Team, 2016).

Species delimitation and divergence times

According to the results, we assess the limits between different groups based on: (1) disjunct populations (groups by geographic region: SMO, GRO, MIA, CHIS), (2) phylogroups (SMO/MIA, GRO, CHIS), and (3) groups separated by the Isthmus of Tehuantepec (east: SMO/MIA/GRO, and west: CHIS). We used the command line of coalescent approach implemented in Bayesian Phylogenetics and Phylogeography software (BP&P v3.4, (Rannala & Yang, 2003; Yang & Rannala, 2010). This method uses the multispecies coalescent model (MSC) to compare different models of species delimitation (Yang & Rannala, 2010; Rannala & Yang, 2013) and species phylogeny (Yang & Rannala, 2014; Rannala & Yang, 2017) in a Bayesian framework, accounting for incomplete lineage sorting due to ancestral polymorphism and gene tree-species tree discordance. We used the concatenated data set of eight molecular markers (mt DNA and nuclear DNA), but mitochondrial markers were treated as one locus, so the total number of this parameter was set to nloci = 5. To confirm consistency between runs, we performed multiple analyses using algorithms 0 and 1 (0: species tree given as fixed, 1: species tree given treated as the guide tree), selecting different seed number between runs and changing finetune parameters (ε, algorithm prior), as suggested by Yang & Rannala (2010). After confirming consistency between runs, the subsequent analyses were performed according to species delimitation using a user-specified guide tree (speciesdelimitation = 1, speciestree = 0), with values of ε = 5. The analyses were conducted using parameter finetune = 1, which allows the program to make automatic adjustments to prior parameters. Because different values of θ (ancestral population size, the product of effective population size N and mutation rate µ per site) can result in different posterior probabilities for the same guide tree (Leaché & Fujita, 2010; Yang, 2015), we used three different values: (1) low θ priors (0.0001, IG: 3, 0.0002), (2) medium θ priors (0.001, IG: 3, 0.002), and (3) high θ priors (0.01, IG: 3, 0.02). The inverse gamma prior to τ (species divergence times) was set to tauprior = 3, 0.03, 1.5% of sequence divergence. Each analysis was run with the reversible-jump Markov chain Monte Carlo algorithm (rjMCMC) for 100 thousand generations, sampling every five, and discarding 30 thousand generations as burn-in.

To estimate divergence times, we decided to evaluate the hypothesis of three independent groups based on full evidence (phylogroups: SMO/MIA, n = 12; GRO, n = 9; CHIS, n = 10). We used the concatenated data set (mt DNA and nuclear DNA) that included data from Eugenes fulgens, E. spectabilis, and Tilmatura dupontii as outgroups. For each partition, we assigned the previous selected evolutionary model. This analysis was performed using StarBeast (*Beast; Heled & Drummond, 2010). We employed an uncorrelated lognormal relaxed clock and a Yule process speciation model to model the tree prior. We assigned a calibration node based on a secondary calibration obtained for the split between the “Mountain Gems” clade (L. rhami, E. fulgens, and E. spectabilis) and “Bees” clade (T. dupontii; 12.5 Mya; (McGuire et al., 2014). We incorporated mean substitution rates reported previously (ATPase 6 and 8, ND2, ND4: Pacheco et al., 2011); CR: (Lerner et al., 2011); AK1, BFib, MUSK, ODC: (McGuire et al., 2014). Three independent analyses were run for 30 million generations, sampling every 1,000. The log and tree files from each analysis were combined in LogCombiner, and visualized in Tracer to confirm convergence and ensure acceptable ESS values (ESS > 200). We discarded the first 15% of trees as burn-in. We used TreeAnnotator v1.8.2 (Rambaut & Drummond, 2007) to summarize trees as a maximum clade credibility tree, and to obtain mean divergence times with 95% highest posterior density intervals. The resulting tree was visualized in FigTree v.1.2.3.

Genetic diversity and population structure

We obtained a concatenated dataset of 1,402 bp for 54 individuals (527 bp of CR and 875 bp of ATPase 6 & 8). The complementary dataset of five molecular markers for 31 individuals included 875 bp of ATPase 6 and 8, 527 bp of CR, 918 bp of ND2, 521 bp of ND4, 758 bp of BFib, 559 bp of MUSK, 495 bp of ODC, and 416 bp of AK1. The initial dataset included 33 haplotypes (24 found with CR and 15 with ATPase 6 & 8). Estimates of haplotype and nucleotide diversity are presented in Table 1. Overall, high haplotype diversity and low nucleotide diversity were observed within groups (SMO, GRO, CHIS).

The mtDNA network revealed significant population structure within the L. rhami complex (Fig. 1). There was a clear separation between populations on either side of the IT, which were separated by twelve mutational steps, while individuals of GRO were closely linked to those of the SMO. In general, the most frequent haplotype was present in populations from the Sierra Madre Oriental group (SMO). Haplotypes from GRO tended to separate from the main haplotype. In Table 2, F_ST values confirm high levels of geographic structure between regions. This further translates into a significant correlation between the genetic distance and the geographic distance matrices, according to the Mantel test, thus suggesting isolation by distance between groups (r = 0.87, p < 0.005).

AMOVA results indicated that the highest genetic variation was observed among rather than within populations, with similar percentages when grouping populations according to geographic groups or on either side of the IT: 76.41% and 78.74% respectively (P <  0.0001, Table 3).

Demographic analyses

Differing conclusions among the methods used to evaluate demographic history led to ambiguous results. The occurrence of historical population expansion was supported by negative and significant values of neutrality tests (Tajima’s D and Fu’s Fs), except for Tajima’s D statistic in SMO and GRO groups (Table 1). A mismatch distribution unimodal curve was recovered for CHIS population, but no significant values of the raggedness index indicated possible demographic expansion in all populations, as curves under the expansion model did not deviate from a unimodal distribution. BSP estimates revealed that effective population size was flat across time for GRO. This pattern was also found for CHIS, however, the higher posterior density low interval suggested a growing demographic tendency, and a subtle demographic expansion was recovered in SMO population (Fig. 2).

Evolutionary models and Phylogenetic analyses

We obtained a concatenated dataset of 5,069 bp. The best-fit models for each molecular marker were as follows: HKY (MUSK), HKY+I (ATPase 6 and 8, ND4), HKY+G (AK1), HKY+I+G (CR), TNR+G (ND2), TPM3uf (ODC), and TPM2uf +I (BFib). Phylogenetic relationships using the multilocus dataset resulted in one main monophyletic group corresponding to individuals from west of the IT (PP >  0.95, Fig. 3: Bayesian Inference). Most individuals from east of the IT were grouped into two well-supported separate clades, but no resolution was recovered for three individuals from this region. Moreover, one well-supported clade included most individuals from GRO group from west of the IT, with the rest of individuals merged in a polytomy with individuals from SMO region.

Morphometric variation

Normality of our data was rejected and there was dimorphism between males and females in all variables (Table S4). General comparisons between geographic groups resulted in significant differences in all variables for both males and females, while comparisons between groups separated by the Isthmus of Tehuantepec showed significant differences in bill width (χ² = 6.3669, p = 0.01163) and wing chord (χ² = 8.0642, p = 0.004515) for males; in females all traits were statistically different except for bill length (Table 4). Paired differences between geographic groups (comparing each group against the others), revealed that CA and GRO present significant differences in several traits (Table S4).

Species delimitation and divergence times

Three different species hypotheses were assessed: (A) species delimited by disjunct populations (groups by geographic region: SMO, GRO, MIA and CHIS), B) phylogroups (SMO/MIA, GRO and CHIS), and (C) groups separated by the Isthmus of Tehuantepec (east: SMO/MIA/GRO, and west: CHIS), using different values of θ priors (population size parameters) to test the sensitivity of the species delimitation results. BP&P analyses testing the first hypothesis (that allopatric populations are different species) resulted in low statistical support for the split of groups SMO and MIA, and suggested the split of GRO group as a valid species only when using low θ priors. The high probabilities for an ancestral node suggest splitting into multiple species, as well as the validity of CHIS as an independent group (Fig. 4A). Analysis of the three lineages corresponding to phylogroups resulted in high speciation probabilities in all cases, except when using high θ priors (0.01, IG: 3, 0.02) for the splitting of GRO and SMO/MIA groups (Fig. 4B). Analysis of the two-lineage hypothesis (i.e., two groups separated by the Isthmus of Tehuantepec) suggested that these groups are separate species (Fig. 4C). Varying the θ priors affected the resulting speciation probabilities; higher θ priors resulted in lower probabilities than lower θ priors. Our divergence time estimates (Fig. 4D) showed that the split between L. rhami complex and its sister group (genus Eugenes) was around 10.55 Mya (8.28–13.14 Mya). The estimate for the first split within the L. rhami complex was dated ca. 0.207 Mya (0.091–0.317 Mya), corresponding to the divergence between populations on either side of the IT. The split between the groups from west of the IT (SMO/MIA and GRO) was dated at ∼0.087 Mya (0.035–0.146 Mya). The independence of groups was supported by high values of posterior probability.