Extensive introgressive hybridization within the northern oriole group (Genus Icterus) revealed by three-species isolation with migration analysis

Until recently, studies of divergence and gene flow among closely-related taxa were generally limited to pairs of sister taxa. However, organisms frequently exchange genes with other non-sister taxa. The “northern oriole” group within genus Icterus exemplifies this problem. This group involves the extensively studied hybrid zone between Baltimore oriole (Icterus galbula) and Bullock's oriole (I. bullockii), an alleged hybrid zone between I. bullockii and black-backed oriole (I. abeillei), and likely mtDNA introgression between I. galbula and I. abeillei. Here, we examine the divergence population genetics of the entire northern oriole group using a multipopulation Isolation-with-Migration (IM) model. In accordance with Haldane's rule, nuclear loci introgress extensively beyond the I. galbula–I. bullockii hybrid zone, while mtDNA does not. We found no evidence of introgression between I. bullockii and I. abeillei or between I. galbula and I. abeillei when all three species were analyzed together in a three-population model. However, traditional pairwise analysis suggested some nuclear introgression from I. abeillei into I. galbula, probably reflecting genetic contributions from I. bullockii unaccounted for in a two-population model. Thus, only by including all members of this group in the analysis was it possible to rigorously estimate the level of gene flow among these three closely related species.


Introduction
The study of phylogeography and population divergence is challenged by the difficulty of distinguishing between shared retained ancestral polymorphism and introgressed alleles (Nielsen and Wakeley 2001). The much larger effective population size of nuclear DNA relative to mtDNA slows the rate of lineage sorting and the speed at which reciprocal monophyly is achieved at the nuclear level (Avise 2004). Slow and stochastic lineage sorting thus causes a pattern of allele sharing among species that can seem similar to patterns caused by gene flow through introgressive hybridization, even between long-divergent lineages (Hudson and Turelli 2003). Multiple unlinked loci and coalescent-based "divergence population genet-ics" methods are needed to distinguish between these two processes (e.g., Kliman et al. 2000;Machado et al. 2002).
The Isolation-with-Migration (IM) model Nielsen 2004, 2007) provides a powerful tool for the study of population divergence and for testing different speciation models and demographic scenarios (Hey 2006;Peters et al. 2007;Kondo et al. 2008;Pinho and Hey 2010). Shared alleles allow multilocus coalescent approaches such as IM to reconstruct the divergence history of closely related species (Nielsen and Wakeley 2001;Nielsen 2004, 2007). The original IM model was designed for two-population studies between sister taxa that have not exchanged genes with other taxa since they split Nielsen 2004, 2007). However, organisms frequently exchange genes with other non-sister taxa (Arnold 2006;Nosil 2008;Petit and Excoffier 2009). Genetic contributions from other source populations not included in the analysis (i.e., "ghost" populations) could thus severely bias or invalidate inference about the history of closely related taxa (Beerli 2004;Slatkin 2005;Strasburg and Rieseberg 2010). Thus, until recently, IM analyses were limited to pairs of taxa that share more recent common ancestry and gene flow with each other. However, this limitation was removed with the release of IMa2, allowing the simultaneous consideration of up to ten populations granted a population tree (i.e., phylogeny) is available for the taxa involved (IMa2: Hey 2010b).
In this study, we used IMa2 to analyze the three species in the "northern oriole" group within the New World orioles (Jacobsen and Omland 2011;Fig. 1). This group includes the well-known hybrid zone between the eastern Baltimore oriole (Icterus galbula) and Bullock's oriole (I. bullockii) (Sutton 1938(Sutton , 1968Sibley and Short 1964;Rising 1969Rising , 1970Rising , 1983Rising , 1996Corbin et al. 1979). I. galbula breeds across the eastern United States south to Louisiana, and across most of Canada from Alberta in the west to Nova Scotia in the east (Fig. 2). The breeding range of I. bullockii extends across western North America from British Columbia, Canada in the north to Durango, Mexico in the south. Both I. galbula and I. bullockii are long-distance migrants and mainly overwinter in southern Mexico and Central America (Rising et al. 1998(Rising et al. , 1999. Their winter ranges overlap with the range of the third member of the group, I. abeillei, a Mexican endemic limited to the central Mexican plateau (Jaramillo and Burke 1999) and traditionally considered a close relative and subspecies within I. bullockii due to their overall phenotypic similarity (Jaramillo and Burke 1999).
The parapatric breeding ranges of I. galbula and I. bullockii overlap along a long and narrow transect (150 km in Kansas, Allen 2002) that stretches over 1800 km from southern Alberta to central Texas (Rising 1983). The two orioles interbreed extensively in this contact zone, which falls within the "Rocky Mountain-Great Plains" suture zone (Remington et al. 1968), a hotspot of avian hybrid zones between eastern and western species pairs (e.g., Carling and Brumfield 2008;Flockhart and Wiebe 2009;Mettler and Spellman 2009). Documentation of viable hybrid offspring within this contact zone (Sutton 1938;Sibley and Short 1964;Sutton 1968;Rising 1970) led to the lumping of I. galbula and I. bullockii (including I. abeillei) into a single species, the northern oriole (AOU 1973). Later studies found that the I. galbula-I. bullockii hybrid zone was stable (Rohwer and Manning 1990;Rohwer and Johnson 1992;Rising 1996;Allen 2002), and full species status was restored for all three members of the group (AOU 1995). Hybridization has never been firmly documented between I. abeillei and the other two mem-bers of the group, but there is some evidence suggesting hybridization with I. bullockii where their ranges abut in northern Durango, Mexico (Miller 1906; but see Rising 1973). We are unaware of any documented instances of hybridization between members of the northern oriole group and other oriole species.
Despite the extensive research on the northern oriole group, we have limited knowledge about the extent of nuclear gene flow within this species complex. Allen

2414
(2002) examined clinal variation of phenotypic and molecular markers across the I. galbula-I. bullockii hybrid zone in Kansas, and found that neither male plumage traits nor mtDNA haplotypes introgressed beyond the zone of contact. These findings were recently corroborated by Carling et al. (2011), who estimated the mitochondrial cline to be 328 km wide. In contrast, of 152 neutral AFLP markers examined by Allen (2002), there were no fixed differences between the two species, and only ten AFLP bands differed more than 40% in frequency across the hybrid zone. Allen (2002) interpreted this apparent lack of differentiation in AFLP markers between I. galbula and I. bullockii as evidence of unhampered nuclear introgression since the two species came into secondary contact. However, such extensive allele sharing among closely related species could also result from the slow sorting of retained ancestral polymorphism (Moore 1995;Hudson and Turelli 2003).
The goals of this study were to (1) estimate levels of introgression among all three species in the northern oriole group while accounting for their recent shared ancestry using a three-population IMa2 model and (2) examine the effects of genetic exchange with unsampled taxa in two-population analyses.

Taxonomic sampling
We sought a broad representation from throughout the ranges of all three species whilst avoiding known areas of sympatry (see Appendix1 for sampling locations). Our sampling scheme provided a total of 21 to 24 individuals of each of the three species, I. abeillei, I. bullockii and I. galbula. In addition, a single individual of western meadowlark (Sturnella neglecta) was included as an outgroup taxon.

Multilocus nuclear data
This study included sequence data from eight intronspanning loci that map to six autosomes and the Z sex chromosome in the zebra finch (Taeniopygia guttata) genome (see Table 1 for details). The two Z-linked loci are well separated on the zebra finch Z chromosome (>7 Mb apart) and were therefore treated as independent (Backström et al. 2006(Backström et al. , 2010. All eight loci were amplified and sequenced following procedures described in previous phylogenetic studies on Icterus (Jacobsen et al. 2010;Jacobsen and Omland 2011). Resulting DNA chromatograms were aligned and edited in SEQUENCHER 4.7 (Gene Codes Corp., Ann Arbor, MI). Edited sequences for all loci were imported into MEGA5 (Tamura et al. 2011) and aligned using MUSCLE (Edgar 2004) with final manual adjustments. The haplotypes of individuals heterozygous at more than one position at a given locus were determined using the software PHASE v.2.1.1 (Stephens et al. 2001;Stephens and Scheet 2005), with the following parameter settings: burn-in = 1000, number of iterations = 10,000, and thinning interval = 1. Each analysis was repeated 10 times using different random starting seed and the run that received the highest Hatched areas indicate range overlap between I. bullockii and I. galbula. Digital distribution maps for each species were downloaded from NatureServe (Ridgely et al. 2007).
log-likelihood was used to infer the haplotypes. To avoid systematic bias in downstream coalescent analyses, all haplotypes resolved at a probability of 0.5 or higher were included in the study (Garrick et al. 2010). An ambiguously aligned single-nucleotide repeat (poly-A) region within locus FGB4 was excluded from the alignment prior to analyses. Linkage analysis of the two Z-linked loci (ALDOB5 and BRM15) using HAPLOVIEW v.4.2 (Barrett et al. 2005) revealed no evidence of LD among polymorphic sites (results not shown). We tested for intralocus recombination using the difference in sums of squares (DSS) test in TOPALi v.2.5 (Milne et al. 2004), with 500 bootstrap replicates to determine statistical significance. A significant DSS peak was detected at locus b-ACT2, and only the largest non-recombinant sequence block (558 bp) was retained for downstream analysis. Accession numbers for sequences obtained from previous studies on this group are listed in Appendix 2. All new sequences collected in this study have been deposited in GenBank under accession numbers JX403068-JX403342 and JX412950-JX412955.
Tests for selection were conducted using a multilocus HKA test (Hudson et al. 1987) implemented in the HKA program available at http://lifesci.rutgers.edu/ $ heylab/ HeylabSoftware.htm. We chose western meadowlark (Sturnella neglecta) as outgroup for this test, the most distant relative within the family Icteridae available for comparison across the eight loci examined in this study. Kimura 2P distances between the three northern oriole species and S. neglecta were calculated using DNASP v.5.10.0.1 (Librado and Rozas 2009). Statistical significance of the HKA test was determined through 10,000 coalescent simulations based on the number of sampled alleles and observed levels of polymorphism and divergence across loci. These simulations also allowed us to conduct a multilocus test of Tajima's D statistic (Tajima 1989). Fu's (1997) Fs and Ramos-Onsins and Rozas (2002)R 2 neutrality statistics were also calculated using DnaSP. Confidence intervals and statistical significance of test values were calculated using 10,000 random permutations and coalescent simulations.

Multilocus coalescent analyses
We examined patterns of gene flow among all three members of the northern oriole group using the Bayesian IM model implemented in the program IMa2 (Hey 2010b), which assumes random mating, constant population sizes, selective neutrality, free recombination among loci and no recombination within loci Nielsen 2004, 2007;Hey 2010b).
By default, all model parameters estimated by IMa2 are scaled to the population mutation rate. The parameter estimates were converted into demographic units (i.e., effective population sizes in number of individuals and divergence times in years) using a neutral mutation rate of 1.35*10 À9 substitutions per site per year for autosomal loci (Ellegren 2007) and 1.45*10 À9 substitutions per site per year for Z-linked loci (Axelsson et al. 2004;Ellegren 2007). To reflect the different modes of inheritance, auto- Table 1. Locus information. Chromosome: intron location determined based on BLAST searches of zebra finch reference genome assembly, with start coordinates for Z-linked loci provided in parentheses; Length: size (in base pairs) of largest non-recombinant alignment segment; N: number of non-recombinant alleles included in the study. Full locus information is as follows: Beta-actin gene, intron 2 (b-ACT2); alpha-enolase gene, intron 8 (a-ENO8); fibrinogen B beta polypeptide gene, intron 4 (FGB4); glyceraldehyde-3-phosphate dehydrogenase gene, intron 11 (GAPDH11); rhodopsin gene, intron 2 (RDP2); and transforming growth factor beta-2, intron 5 (TGFb5); aldolase-B fructose-biphosphate intron 5 (ALDOB5); Brahma protein gene, intron 15 (BRM15).
somal (1.0) and Z-linked (0.75) loci were assigned respective inheritance scalars in the IMa2 input file. Population sizes were converted into effective number of individuals using a generation time of 1.7 years, the average age at reproduction documented for multiple passerine lineages (Saether et al. 2005). The actual generation time in orioles may well be slightly longer, as most oriole species exhibit delayed plumage maturation, which could delay the onset of reproduction (for males in particular) (e.g., Flood 1984;Rohwer and Manning 1990;Butcher 1991;Richardson and Burke 1999). Effective population sizes reported in this study are therefore likely conservative estimates of the actual effective population sizes. The HKY mutation model was specified for all loci.

Three-population IMa2 analysis
First, we estimated levels of gene flow among all three species within the northern oriole group simultaneously in a three-population IMa2 analysis. We conducted separate analyses with two alternative population trees: (1) the nuclear species tree inferred by our recent species tree analysis of Icterus (Jacobsen and Omland 2011): (I. galbula, (I. abeillei, I. bullockii)) and (2) the mtDNA gene tree inferred by Omland et al. (1999): (I. bullockii, (I. abeillei, I. galbula)). Short preliminary runs were conducted to determine appropriate upper bounds of the demographic parameters. Subsequently, two identical long runs of 3-5*10 6 MCMC steps with a burnin period of 1*10 6 steps were started with different starting seeds and run until minimum effective sample sizes (ESS) of 200 were achieved. We assessed convergence by monitoring the mixing properties of the MCMC during each run and by ensuring that similar parameter estimates were obtained from the two independent runs. Acceptable chain mixing (i.e., absence of long-term trends in plots of L[P] and t over the course of each run) and low autocorrelations (<0.1) were achieved by running 100 chains under a geometric heating scheme (with heating parameters set to a = 0.975 and b = 0.75). Population-specific parameter priors were provided in a separate file. We used exponential priors on migration rates (option-j7), and determined appropriate prior distribution means based on the peak posterior parameter estimates from preliminary runs.

Two-population IMa2 analyses
To assess the effects of violating the crucial assumption of no allelic contributions from other source (i.e., "ghost") populations, we also conducted three separate pairwise IMa2 analyses of the three northern oriole species. Run settings were kept unchanged, except for slight adjustment of parameter priors where posterior density distributions exceeded the upper bounds determined by preliminary runs. Each MCMC (M) mode analysis was replicated using different random seeds to ensure convergence onto similar posterior distributions, using the same criteria for assessing convergence as described above. Finally, we conducted nested model testing by running IMa2 in "Load genealogies" (L) mode to specifically compare the fit of reduced models that do not allow for post divergence gene flow to the full IMa2 model (i.e., including all six parameters Θ 1 , Θ 2 , Θ a , m 1 , m 2 , t). This was performed by sampling 2.5*10 5 genealogies evenly from the two MCMC runs performed on each population pair and compared the log-likelihood scores of all possible nested models provided in a separate file in the IMa2 distribution package.

Polymorphism and divergence
The three species within the northern oriole group were polymorphic at all loci used in this study, except for I. abeillei that was monomorphic at ALDOB5. Sample sizes and data for individual loci are summarized in  (Table 2), but not significantly so (single factor ANOVA: h: F 2,21 = 0.656, P = 0.529; Hd: F 2,21 = 0.585, P = 0.566; p: F 2,21 = 0.943, P = 0.405; respectively). The differences in mean genetic diversity between I. abeillei and the two other members of the group are mainly driven by the lack of sequence variation at the sex-linked locus ALDOB5 ( Table 2). The haplotype networks revealed extensive allele sharing and few fixed differences among the three northern oriole species (Fig. 3). Some loci were characterized by a distinct star-shaped topology (e.g., a-ENO8, ALDOB5 and BRM15), with a common central (and presumed ancestral) haplotype surrounded by several rare, derived haplotypes one or two mutational steps away from the ancestral haplotype (Fig. 3b,g,h).
The three species showed varying degrees of genetic differentiation across loci, but generally exhibited extensive haplotype sharing across loci (Table 3, also see Fig. 3). Strongest differentiation was observed between the highly disjunct I. abeillei and I. galbula, with highly significant

2418
F ST values at six of the eight loci examined. The two widely hybridizing I. bullockii and I. galbula were strongly differentiated at only half (4/8) of the loci examined. The two sister species, I. abeillei and I. bullockii, were only weakly differentiated and only two loci (FGB4 and ALDO5) remained significant after Bonferroni corrections ( Table 3). The multilocus HKA test indicated no significant departures from neutrality in our dataset. The multilocus Tajima's D test was significant only for I. bullockii, a result that is mainly driven by significantly more negative than expected D values at loci a-ENO8 and BRM15 (Table 4). Fs values were negative at the majority of loci, indicative of recent population growth, but this statistic deviated only significantly from neutral expectations at less than half the loci following Bonferroni corrections ( Table 4). The R 2 test was only significant at a single locus a-ENO8 (Table 4).

Three-population IMa2 analysis
The inclusion of all three members of this species complex in a single IMa2 analysis revealed some interesting patterns of gene introgression (Fig. 4). The two competing population trees (nDNA vs. mtDNA) yielded similar parameter estimates (with one important exception discussed below). We therefore only report the results based on the nuclear population tree (I. galbula, (I. abeillei, I. bullockii)).
First, we found clear evidence of substantial introgression across the hybrid zone between I. galbula and I. bullockii, (2Nm I. galbula > I. bullockii = 3.388, 95% HPD = 0.41-12.30, vs. 2Nm I. bullockii > I. galbula = 2.125, 95% HPD = 0.193-6.504, Fig. 4a). Interestingly, the gene flow appears to be slightly skewed in the direction of I. bullockii. In contrast, we did not detect any gene flow between the sister species I. bullockii (Fig. 4b). The gene flow estimates virtually peaked at zero (2Nm = 0.001 in both directions) and the 95% HPD intervals ranged from 0 to 1.148 and 0 to 0.495 in the direction of I. bullockii and I. abeillei, respectively. Similarly, the three-population analysis did not detect any gene flow between the highly allopatric I. galbula and I. abeillei (Fig. 4c), with gene flow estimates peaking near zero and 95% HPD intervals ranging from 0 to 0.631 and 0 to 0.719 in the direction of I. galbula and I. abeillei, respectively. There was no evidence of historical gene flow between the I. galbula lineage and the ancestral I. abeillei/ I. bullockii lineage.
IMa2 provided very concise estimates of the timing of divergence events within the group. The analysis based on the nuclear species tree suggested that the initial divergence between the I. galbula lineage and the ancestor of I. abeillei and I. bullockii appears to have occurred approximately 350,000 years ago (95% HPD = 215,838-568,125), whereas I. abeillei and I. bullockii split much more recently approximately 95,000 years ago (41,865-159,087, Fig. 5a). In contrast, the analysis based on the mtDNA population tree (I. bullockii, (I. abeillei, I. galbula)) resulted in perfectly overlapping estimates of the two splitting events around 215-250,000 years ago, with nearly perfectly overlapping HPD intervals (Fig. 5b). Effective population sizes (N e ) were approximately 398,000 blackbacked orioles, 557,000 Baltimore orioles, and 921,000 Bullock's orioles (Fig. 6).

Two-population IMa2 analyses
A summary of the three separate two-population analyses is presented in Table 3. Estimates of gene flow were mostly consistent between the pairwise two-population analyses and the three-population analysis. However, the two-population analysis of I. abeillei and I. galbula suggested a small but significant amount of introgression from I. abeillei into I. galbula (2Nm I. abeillei > I. galbula = 0.477 (Fig. 4f). Although the 95% HPD interval included zero (see Table 5), nested model testing conducted in IMa2 indicated that a reduced model with zero migration was a significantly worse fit to the data than the full model allowing for migration (results not shown). Furthermore, the two-population analysis of I. galbula and I. bullockii indicated lower but equal rates of introgression in both directions (2Nm I. bullockii > I. galbula = 1.331 and 2Nm I. galbula > I. bullocki = 1.275, see Fig. 4d). The divergence time estimates from the pairwise twopopulation analyses were comparable to those produced by the three-population IMa2 analysis, except that both divergence events (I. galbula-I. abeillei/I. bullockii split, and I. abeillei-I. bullockii split) were estimated to have occurred roughly 70,000 years earlier than that indicated by the three-population analysis (see previous section, Fig. 5). In agreement with the nuclear species tree of this group (Jacobsen and Omland 2011), the most recent divergence in this group between I. abeillei and I. bullockii was estimated at approximately 163,000 years ago, well within the 95% HPD interval of the three-population   estimate (see previous section). Similarly, the older split between the I. galbula lineage and the sister taxa I. abeillei and I. bullockii were estimated at 411,000 and 424,000 years ago, respectively. Effective population size (N e ) estimates from the pairwise two-population analyses were also roughly similar to those estimated by the three-population analysis. The peak estimates were slightly higher in the two-population analyses, but the 95% HPD intervals overlap greatly with the estimates from the three-population analysis. The most striking deviations in population size estimates arose when the non-sister species, I. galbula and I. abeillei, were analyzed separately, which resulted in overestimation of N I. galbula and a slight underestimation of N I. abeillei (Table 5). Similarly, N I. bullockii was substantially smaller when paired with the non-sister I. galbula than when correctly paired with its sister species I. abeillei (Table 5, Fig. 6).

Discussion
The recent development of a multispecies coalescent framework for the study of divergence and gene flow among three or more taxa simultaneously in the program IMa2 (Hey 2010b) has proven useful for understanding the population dynamics of complex taxonomic groups consisting of several recently diverged species (e.g., Hey 2010a; Illera et al. 2011;Schoville et al. 2011;So et al. 2011;Li et al. 2012). In this study, we were able to estimate levels of nuclear introgression between allopatric populations of the three species in the northern oriole group while also accounting for their recent shared ancestry.
The hybrid zone between the non-sister species, I. galbula and I. bullockii, represents perhaps one of the most blatant case examples for the need of a multipopulation framework to disentangle the demographic history of larger groups of taxa. The reproductive barriers that so effectively prevent mtDNA and phenotypes of both I. galbula and I. bullockii from introgressing beyond the   (HPD95Lo) and upper (HPD95Hi) bounds of 95%HPD intervals for the effective population size of population 1 (N 1 ), population 2 (N 2 ), and ancestral population (N a ), splitting time (t), and population migration rate from population 2 into population 1 (2Nm1 > 2) and from population 1 into population 2 (2Nm2 > 1) forward in time. Population sizes are given in 1000 individuals and time since divergence in 1000 years. stable hybrid zone appear to be "leaky", allowing nearneutral nuclear alleles to introgress and successfully spread in a heterospecific genetic background without erasing the diagnostic species differences. Such a substantial rate of introgression evident in populations several hundred miles away from the nearest known areas of sympatry has major implications for our understanding of the population dynamics of all three species within the group. Limited mtDNA introgression (Allen 2002;Carling et al. 2011) but extensive nDNA introgression (Allen 2002; this study) follows the predictions of Haldane's rule (Haldane 1922), wherein reduced female hybrid fitness prevents maternally inherited mtDNA, but not bi-parentally inherited nDNA from introgressing (Brookfield 1993;Turelli and Orr 1995;Orr 1997). Sex-biased (i.e., male driven) gene flow due to female heterogamety and the effects of Haldane's rule have increasingly been documented in a wide range of avian groups (e.g., Borge et al. 2005;Carling and Brumfield 2008;Carling et al. 2010;Storchová et al. 2010;Backström and Väli 2011). A question of interest to researchers studying divergence and gene flow is to determine the timing of migration events between diverging populations. If possible, knowing when gene flow occurred in the speciation process would be of great importance to distinguish between an allopatric speciation model with recent secondary contact following range expansion (gene flow late) and a parapatric speciation model where divergence occurred in the presence of gene flow with subsequent decrease, then cessation of gene flow once reproductive barriers evolved. However, it is now clear that IMa2 and similar methods are unable to distinguish between these different modes of speciation (Sousa et al. 2011;Strasburg and Rieseberg 2011). We therefore cannot infer with confidence whether the currently hybridizing I. galbula and I. bullockii diverged in allopatry and only recently came back into secondary contact following postglacial range expansion, or if they diverged in parapatry while still exchanging genes throughout the early stages of divergence. Nonetheless, the fact that these two orioles are nearly 5% divergent in mtDNA (Omland et al. 1999) indicates that they have been evolving independently for an extended period of time, consistent with allopatric speciation followed by one or more episodes of gene flow during the last glacial cycles of the Pleistocene epoch. Pliocene or Pleistocene divergence followed by postglacial expansion has been proposed to explain the current distribution of many parapatric species found in North America and Eurasia (Klicka and Zink 1999;Hewitt 2004;Johnson and Cicero 2004). A recent exchange of nuclear variation through introgressive hybridization during secondary contact in the late Pleistocene would explain the wide gap between the mtDNA and nDNA divergence time estimates for these two species, even though the TMRCA for a single locus (e.g., mtDNA) will necessarily always predate the divergence between two species (Edwards and Beerli 2000).
In contrast, the IMa2 analyses revealed no evidence of introgression between I. abeillei and I. bullockii. The region of parapatry in Durango, Mexico where hybridization between these two species is alleged to have taken place was unfortunately not represented in this study. It is therefore possible that some hybridization could occur within this region, but restricted enough to prevent widespread introgression into the core ranges of either species. Thus, the ubiquitous nuclear haplotype sharing observed between these sister species most likely represents ancestral polymorphism retained in both descendent lineages since their recent split in the late Pleistocene (Hudson and Turelli 2003) The low rate of unidirectional introgression from I. abeillei into I. galbula that was indicated by the two-population IMa2 analysis is somewhat puzzling, given that no evidence of gene flow between these two species was revealed in the three-population analysis. Previous studies on these two species using IM based on mtDNA and two nuclear introns also found no evidence of gene flow between these two mitochondrial sister taxa (Kondo et al. 2004(Kondo et al. , 2008. The fact that they are reciprocally monophyletic at the mtDNA level separated by a single fixed substitution also strongly suggests that there is no recent or ongoing mitochondrial gene flow between these highly allopatric species (Kondo et al. 2008).
Rather, we suspect that this signature of gene flow is a spurious outcome arising when violating the assumption of the IM model that no genetic exchange has occurred between the two focal populations and ghost populations (Beerli 2004;Slatkin 2005). As I. abeillei shares a more recent evolutionary history with I. bullockii than with I. galbula (Jacobsen et al. 2010;Jacobsen and Omland 2011), retained ancestral polymorphisms shared between I. abeillei and I. bullockii could introgress into I. galbula via I. bullockii over the extensive hybrid zone. This study revealed just how easily neutral nuclear alleles can cross this species boundary and freely introgress across the entire range of the other species. By excluding I. bullockii from the analysis, such shared ancestral alleles could mistakenly be inferred by IMa2 to have introgressed directly from I. abeillei into I. galbula (Fig. 7). Furthermore, the many alleles introgressing from I. bullockii into I. galbula have probably contributed to the high nucleotide diversity in I. galbula, causing IMa2 to overestimate the current effective population size of I. galbula in analyses that do not account for I. bullockii.
Thus, exclusion from the analysis of other taxa that have exchanged genes with either of the two focal taxa through shared ancestry or gene flow can bias estimation 2422 of population sizes and gene flow in IM, and the severity of the bias depends on how large the genetic contribution from the third taxon has been (Strasburg and Rieseberg 2010).

Conclusions
Our study reveals the complexity of evolutionary processes (recent divergence, slow sorting of ancestral variation, and introgressive hybridization) that have contributed to produce the observed pattern of rampant allele sharing among the three species within the northern oriole group. Furthermore, our findings highlight the importance of accounting for all potentially interbreeding taxa when interpreting the divergence history of a group of taxa, and the extent to which gene flow may retard the rate of divergence at neutral loci. Thus, this study provides an interesting perspective on the broad-scale impacts of introgressive hybridization, which has important implications for everything from understanding adaptive evolution to inferring species trees. The ease with which genome-scale datasets can now be obtained will soon enable researchers to study even more complex clades with four or more taxa that may interbreed.  The full three-population analysis inferred gene flow only between the hybridizing Icterus galbula and I. bullockii, whereas (b) two-population analyses also indicated a low level of introgression (highlighted in red) from I. abeillei into I. galbula. The introgression between these highly allopatric species is best ascribed to ancestral polymorphisms shared between I. abeillei and I. bullockii introduced to I. galbula across the I. galbula-I. bullockii hybrid zone. This study clearly illustrates how introgression from other taxa not included in the analysis (in this case I. bullockii) can bias gene flow estimation if not accounted for.