Molecular Data.—

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 4.6.1.0 (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⁻⁸ 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⁻⁹ substitutions per site per year for autosomal loci (Ellegren 2007), and for the Z-linked loci a rate of 1.45 × 10⁻⁹ 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.

Morphometrics.—

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.

Molecular Data.—

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 F_ST 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.

FIG. 2. 

Allele networks of the six nuclear loci sequenced, where black circles represent the north group and grey circles represent the south group. Note that all loci show allele sharing between geographic groups.

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).

FIG. 3. 

Posterior probabilities for (a) effective population size, (b) gene flow; 2NM, and (c) time since divergence using Isolation with Migration Analytical (IMa).

Morphological Data.—

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).

TABLE 3.

Untransformed means and standard deviations of morphological characters measured (in millimeters) from I. chrysater.

TABLE 4.

Principal component analysis of the morphometric characters. Factor loadings for PC1 and PC2 explain 64% of the variance. Also, it is shown that wing, tarsus, and tail length are the higher values.

FIG. 4. 

Plotting of the Principal Components of the six variables for the four subspecies. Note there is a gradual cline from the subspecies from north (+, ▪), to south (▵, ▴), but there is no observable diagnosable character difference between them.

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 F_ST 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.