Geographic barriers between populations of a species can result in divergence of genes, morphology, or behaviors that can lead to speciation. The Yellow-backed Oriole (Icterus chrysater) is distributed from southern Mexico to Colombia, but with a major range disjunction of 600 km in Costa Rica. We examined molecular and morphological data for differences between northern and southern populations. We sequenced the mitochondrial control region and six nuclear introns. Genetic data show strong north–south population structure with evidence of gene flow. The evidence of gene flow between populations is surprising because of the large geographic break between populations. We also measured six morphological characters from specimens collected along the species’ distribution and found shallow north–south divergence.
Studies of variation in organisms across their distributional ranges are important for understanding the first steps along the path to speciation (Rising 1970, Price 2008). Color, plumage patterns, vocalizations, morphology, and DNA sequences can all vary geographically in natural populations. Divergence between populations can build up over time and may manifest as major morphological, genetic, or behavioral differences between them (de Queiroz 2005, Leaché et al. 2009). It is crucial to use multiple criteria, including reproductive isolation, ecological or niche divergence, phenotypic distinctiveness, genetic clustering, and coalescence of alleles to fully understand the patterns and processes of early evolutionary divergence (Hart 2011). In particular, use of multiple criteria in species delimita-tion increases the likelihood of detecting recent divergence. However, using only one criterion at a time can yield conflicting views (de Queiroz 2005, Leaché et al. 2009).
In recent years, mitochondrial DNA (mtDNA) has been the most commonly used marker in studies of early divergence. Even though mtDNA can be an excellent marker, studies of early divergence and speciation based on it alone are not ideal for several reasons (Joseph and Omland 2009). First, mtDNA is a single genetic locus and thus represents a limited sample of the whole genome; second, it is maternally inherited so some processes may bias mtDNA evolution (e.g., sex-biased dispersal, hybridization and/or incomplete lineage sorting); third, different parts of the genome have different evolutionary rates. Thus, relying on one gene can generate biased inferences about the population history of a species (Peters et al. 2008, Joseph and Omland 2009). The use of multiple independent loci will result in a more complete understanding of population histories since mutation, genetic drift, and natural selection occur independently in unlinked loci (Peters et al. 2008, Peterson and Nyári 2008).
One of the factors that may drive lineage separation is habitat change because of climatic oscillations (e.g., in temperature, precipitation, and sea levels). Interactions between climate change and the complex topography of Central/South America affected the ranges of many species during the Pleistocene (Hewitt 2000, Zarza et al. 2008). Tropical forests were fragmented by drier habitats like savannas and deserts that expanded during cooler times (Hewitt 2000), leading to diversification of numerous North American birds, mammals and squamates (Hewitt 2004). In birds, high levels of intraspecific morphological variation can result from sexual selection on characteristics such as coloration (Hoekstra and Price 2004) and/or viability selection (perhaps related to foraging), on size and shape (Burns et al. 2002), or other ecological adaptations (James 1983). Within the genus Icterus, extensive color and pattern variation occurs, with no two species identical in plumage (Omland and Lanyon 2000).
Icterus chrysater (Yellow-backed Oriole) is a tropical species found from southern Mexico to South America (Fig. 1). It occurs from sea level up to 2,900 m and inhabits a wide variety of habitats, from scrubby woodlands to pine-oak forest and evergreens (Howell and Webb 1995, Jaramillo and Burke 1999). The distribution of this species is remarkable in that it is separated by a 600 km gap into two disjunct ranges: (1) north, from the Isthmus of Tehuantepec south to Nicaragua, and (2) south, from Panama to northern Colombia (Fig. 1). Differences in the overall size of individuals among different populations lead to recognition of four subspecies—two in the north (I. c. chrysater and I. c. mayensis), and two in the south (I. c. giraudii, and I. c. hondae). Furthermore, some authors have proposed that Icterus chrysater comprises multiple species (Navarro-Sigüenza and Peterson 2004), and some subspecies were once classified as species (Jaramillo and Burke 1999).
No study to date has used molecular phylogeography to inform taxonomy or knowledge of species limits within this complex. The main goal of this study was to use a combination of multi-locus genetic approaches and morphometric characters to determine whether this species has diverged across its major range disjunction. We investigated (1) whether degree of genetic variation corresponds to morphological variation, (2) using coalescent methods, the approximate time of splitting between populations, and (3) evidence of gene flow between the two groups. This study furthers understanding of phylogeographic patterns and speciation of Neotropical bird species and of the overall biogeographic history of this complex region of Central America and northern South America.
We sampled 32 individuals of the Yellow-backed Oriole complex (19 fresh tissue specimens, and toepads from 13 museum skins) to cover as much of the species’ distribution as possible (Table 1). For comparison and tree rooting, we used Icterus graduacauda, the sister species of I. chrysater (Omland et al. 1999, Jacobsen et al. 2010). DNA was extracted using the DNeasy extraction Tissue Kit (QIAGEN N.V., Hilden, Germany). To avoid contamination of the toepads, we extracted and amplified these samples in a laboratory that does not deal with modern bird DNA. We modified the DNeasy protocol when dealing with old toepads by adding 25 µl of 100 mg/ml solution of DTT (dithiothreitol) to the lysis buffer and the 20 µl of Proteinase K.
We amplified a 345-bp fragment of the mitochondrial control region (CR-Domain I) and six nuclear introns (α-enolase, GADPH, Rhodopsin 2, TGFB, SLC9 and MUSK; Cortés-Rodríguez et al. 2013) via the polymerase chain reaction (PCR). To obtain enough intron PCR product for sequencing, we performed a “re-amplification” of the toepad PCR product, using 2–4 µl as a template. We conducted PCRs on the Gene Amp PCR System 9700 (Applied Biosystems Inc., Foster City, CA, USA) and verified products on a 1% agarose gel with ethidium bromide. We cleaned samples by using Exosap – it (Affymetrix Inc., Santa Clara, CA, USA) and sequenced by using ABI Prism Dye Terminator Cycle Sequencing ready reaction with AmpliTaq DNA Polymerase FS. We removed excess dye terminators by the EDTA/ethanol precipitation protocol recommended by ABI, and we ran sequences on an ABI 3100 automated sequencer, located in the Department of Biological Sciences at the University of Maryland, Baltimore County.
Population Genetic and Coalescent Analyses
We edited and aligned sequences using Sequencher 4.1 (Gene Codes Corp., Ann Arbor, MI, USA). We constructed haplotype networks using Network 22.214.171.124 (Fluxus Technology Ltd., Suffolk, UK; Bandelt et al. 1999). The populations were first analyzed according to the four subspecies, and then we compared the two main groups (north and south of Costa Rica). To estimate population divergence between the groups at the mitochondrial level, we performed an analysis of molecular variance (AMOVA; Excoffier et al. 1992), which estimates variance components like Φ-statistics, analogs of Wright's F-statistics as implemented in Arlequin v3.01 (Bowie et al. 2004, Excoffier et al. 2005). Significance values for AMOVA were tested with 1,000 permutations. In addition, we generated mitochondrial Bayesian skyline plots, which estimate historic patterns of population size using their genealogies without a prior demographic model, in the program BEAST v. 1.6.2 (Drummond and Rambaut 2007, Ho and Shapiro 2011). Each run was performed for 10 million generations with a relaxed clock model and a coalescent Bayesian skyline tree. Trees were sampled every 1,000 generations with a burn-in of 10,000 generations. The time since divergence for the control region was calculated using the geometric mean of mutation rates which was 0.025 s/s/l/m (substitution/site/lineage/million years; Kondo et al. 2004). We used TRACER 1.5 (Rambaut and Drummond 2009) to visualize the results and to ensure convergence and ESS values. In addition, we used STRUCTURE 2.3.X (Pritchard et al. 2000) to look for signs of population differentiation between north and south ranges; we include all six introns and the control region in this analysis and we set K to 2 and 4 populations.
We used coalescent analysis of multiple independent loci to get a more complete picture of evolutionary history and population genetics (Peters et al. 2008, Peterson and Nyári 2008). We used the Bayesian model implemented in Isolation with Migration Analytical (IMa) to capture the population dynamics of the species in early stages of differentiation (Hey 2007). To ensure mixing of the Markov chains (to facilitate convergence), we ran multiple heated chains, and we checked the autocorrelation and estimates of ESS (effective sample sizes; Hey 2007). We calculated the effective number of migrants per population per generation (2NM) directly in IMa 2.0 (Hey and Nielsen 2007).
To convert the estimated parameters into biologically informative values, we used previously published information for passerine birds. We assumed an average generation time of 1.7 years (Sæther et al. 2005) and a substitution rate of 5.0 × 10−8 substitutions per site per year (s/s/y) for the mitochondrial control region (Milá et al. 2007). We used a mutation rate of 1.35 × 10−9 substitutions per site per year for autosomal loci (Ellegren 2007), and for the Z-linked loci a rate of 1.45 × 10−9 substitutions per site per year (Axelsson et al. 2004, Ellegren 2007). We used three replicate runs with different seed numbers to check for convergence of the posterior probabilities using the same priors (-q1 20 -q2 10 -qA 5 -m1 15 -m2 10 - t 1.5 -b 1000000 -L 4.0 -N 30 -s125699 -p245 -fg -g1 0.95 -g2 0.9 -j1). Because we were running multiple loci and multiple chains, we used the “heated mode” to ensure better swapping of chains.
We measured 140 adult male specimens of Icterus chrysater: 55 from its northern distribution and 85 from its southern range (see Appendix 1). We measured six morphological characters: wing length (WL), bill length (BL), bill width (BW), bill depth (BD), tarsus length (TAR), and tail length (TAIL). We took tail and wing measurements by using a ruler to the nearest 1 mm; all other measurements were taken with digital calipers (Mitutoyo CD-6” CX, Mitutoyo Corp., Kawasaki, Japan) with 0.01mm precision. All measurements were taken three times and the means used for further analyses (Yezerinac et al. 1992).
We first compared data among the four named subspecies, and then we grouped them into south and north groups to make comparisons across the range disjunction in Costa Rica. All statistical analyses were performed using SPSS v. 21.0 (IBM Corp., Armonk, NY, USA). Before conducting the analyses, we tested the variables for homoscedasticity, both untransformed and after applying log-transformations; however, all of the variables failed the equal variances test, so it was necessary to use non-parametric analyses. We performed Principal Components Analysis (PCA) to control for correlation among the variables. Principal components with eigenvalues >1.0 were then analyzed using a two-sample Mann-Whitney U test as a non-parametric alternative to ANOVA.
Control region sequences for Icterus chrysater had a nucleotide diversity (π) of 0.00249, while haplotype diversity (h) was 0.496 (six haplotypes). The haplotype network shows no sharing of haplotypes between north and south (Fig. 1); however, populations in the north and south regions are not reciprocally monophyletic. The average nucleotide difference for the control region between north and south populations was 1.2%. STRUCTURE did not reveal any population structure, showing only one panmictic population, and the historical demographic test (Bayesian skyline-plots) showed no evidence of a population expansion. However, AMOVA showed that 61% of the observed variation is between north and south groups, with an FST of 0.71 (Table 2).
The six nuclear introns sequenced totaled 2,789 base pairs, of which 53 (2%) were polymorphic sites. The allele networks show similar star-like patterns (Fig. 2) with some sharing of common internal alleles, while other alleles are not shared. Also, the networks show that populations from the north have a higher allelic diversity, whereas populations from the south generally show very low allelic diversity. We conducted AMOVA for each locus (except for SLC9 because we only had one population from the south and the analysis could not be performed). The average of those AMOVAs (analyses for separate loci not shown) resulted in 9.5% of the variation explained by north versus south groups, 24.3% because of variation among subspecies within groups, and 66.0% because of variation within subspecies.
Isolation with Migration Analytical (IMa) allowed us to have a better understanding of the demographic history of the populations. Inheritance scalars were defined depending on the mode of inheritance for each locus (autosomal = 1 and Z-linked = 0.75), and the infinite-sites mutation model was used. IMa estimated the Ne for the north group as 500,000 individuals (90% HPD = 170,000–200,000 individuals) and the Ne for the south group as 170,000 individuals (90% HPD = 40,000–750,000 individuals), making the north Ne almost three times larger (Fig. 3). The estimated migration between the two populations was higher than one migrant per generation: 2NM was 1.87 from the south to the north (95% HPD = 0.32–24 migrants per population per generation) while from north to south it was 1.07 (95% HPD = 0.12–16), to suggest the presence of some gene flow between the two populations across Costa Rica (Fig. 3b). Finally, as can be observed in Fig. 3c, the approximate estimated time of divergence between these two groups is about 150,000 years ago (90% HDP = 67,000–650,000 years; mid to late Pleistocene).
Table 3 shows the untransformed means and standard deviations for each of the six morphological characters measured in I. chrysater. PCA scores indicate significant differentiation between north and south groups; PC I and PC II together explained 64% of the variation distributed respectively as 48% and 16% (Table 4). The characters that contributed the most to PC I were wing length, tarsus length, and tail length (“overall size”); PC II factor loadings were stronger for bill length and bill depth (“bill”). Also, figure 4 shows some morphometric differentiation between the north and south subspecies. Because these birds are found in a wide variety of habitats and altitudes, we separated the data into lowlands (>1,000 m) and highlands (<1,000 m) to see if there would be a pattern between these two groupings; however, we did not find any cluster that separates lowland and highland populations. Because our data could not be normalized, we performed a Mann-Whitney test of PC1, and it showed that the north vs. south groups differ significantly (U = 718, P < 0.001).
Early Stages of Genetic Differentiation and Gene Flow across a Wide Geographic Break
Many organisms show strong geographic variation within named species in morphology, behavior, ecology, or genetics (Zink and Dittmann 1993, Greenberg et al. 1998, Milá et al. 2000, Baker et al. 2003). Our results show low genetic variability of the control region within Icterus chrysater (π = 0.0024), which has also been observed in other members of the genus Icterus (Baker et al. 2003, Cortés-Rodríguez et al. 2008, Kondo et al. 2004, Sturge et al. 2013) to demonstrate deep divergences between Caribbean Island populations. The AMOVA shows that a large percentage of the genetic variation is explained by variation among groups (61.8%) and with a very high FST value of 0.71 (Table 2). However, even though the mtDNA network shows no haplotypes shared between north and south, there are no shared base pairs that distinguish all north and south populations which implies shallow genetic divergence (Joseph and Omland 2009). The evidence of gene flow between populations is surprising given the 600 km geographic break separating them and our lack of known seasonal migration. IMa analyses indicate roughly two migrants per generation in each direction (2NM).
This finding of gene flow is in striking contrast to what we have found in previous comparisons within the same genus (Icterus). For example, in a companion study, we compared Icterus chrysater to its sister species, Icterus graduacauda (2.6% mtDNA divergent); I. graduacauda is separated from I. chrysater by only the narrow Isthmus of Tehuantepec, a lowland valley that separates these two highland species and with a minimum distance between the species of ~50 km (Cortés-Rodríguez et al. 2013). Despite this very narrow barrier between them (and their very close relationship), IMa analysis showed little or no evidence of gene flow between the two species in either direction (95% HPD 0–0.002 and 0–0.004; Cortés-Rodríguez et al. 2013).
Shallow Morphometric Divergence between Populations of Icterus chrysater.—We found a statistically significant mean difference in overall morphology between north and south populations, but the range of variation within those groups was widely overlapping. Thus, we found apparent lack of clear morphological divergence, even though there are clear habitat differences among populations (Cortés-Rodríguez et al. 2013).
In conclusion, our study highlights the importance of using multiple lines of evidence, such as mtDNA, nuDNA and morphometrics, when studying population divergence within species. The case of Icterus chrysater is especially interesting, because this species has two disjunct populations separated by a 600 km gap. The presence of gene flow across the wide geographic gap is interesting, because it is enough to slow or prevent divergence between the populations. This result contrasts with a companion study comparing Icterus chrysater to its sister species, Icterus graduacauda. The two species are only separated narrowly by the Isthmus of Tehuantepec (~50 km), yet our coalescent analysis showed no gene flow between them (Cortes-Rodriguez et al. 2013). Thus, we suggest that one of the factors that may be crucial in facilitating early genetic divergence between populations is plumage differentiation. In the absence of plumage differences, ecological differences and geographic distance are likely not sufficient for preventing gene flow during the early stages of divergence. Clearly, the north and south populations of I. chrysater, which show no known plumage divergence, represent the earliest stages of intraspecific divergence.
We would like to thank the curators, collectors, and staff of the following institutions for tissue loans: Marjorie Barrick Museum of Natural History, University of Nevada; Burke Museum of Natural History and Culture, University of Washington, Seattle; Peabody Museum of Natural History, Yale University; Louisiana State University Museum of Natural Sciences; Delaware Museum of Natural History; Academy of Natural Sciences, Philadelphia; Museo de Zoología Facultad de Ciencias, UNAM, Mexico City; Colección Nacional de Aves Instituto de Biología, Universidad, Nacional Autónoma de México; and University of Kansas Natural History Museum. We extend special thanks to members of the Omland lab and the Leips lab at UMBC. Blanca E. Hernandez-Baños, Frode Jacobsen, Jeffrey Peters, Flor Rodríguez-Gómez and César A. Ríos Muñoz provided assistance during data analyses. Also, Samuel López de Aquino (Museo de las Aves, Saltillo, México) and the staff at Museo de Zoología Facultad de Ciencias (UNAM) assisted during fieldwork, especially M. en C. Alejandro Gordillo Martínez and Iván Anuar López López. This research was supported by a U.S. National Science Foundation grant to KEO (DEB-1119506), and from Mexico PAPIIT-DGAPA (IN225611) and Consejo Nacional de Ciencia y Tecnología (CONACyT).