Hybridization among Arctic white-headed gulls (Larus spp.) obscures the genetic legacy of the Pleistocene

We studied the influence of glacial oscillations on the genetic structure of seven species of white-headed gull that breed at high latitudes (Larus argentatus, L. canus, L. glaucescens, L. glaucoides, L. hyperboreus, L. schistisagus, and L. thayeri). We evaluated localities hypothesized as ice-free areas or glacial refugia in other Arctic vertebrates using molecular data from 11 microsatellite loci, mitochondrial DNA (mtDNA) control region, and six nuclear introns for 32 populations across the Holarctic. Moderate levels of genetic structure were observed for microsatellites (FST= 0.129), introns (ΦST= 0.185), and mtDNA control region (ΦST= 0.461), with among-group variation maximized when populations were grouped based on subspecific classification. Two haplotype and at least two allele groups were observed across all loci. However, no haplotype/allele group was composed solely of individuals of a single species, a pattern consistent with recent divergence. Furthermore, northernmost populations were not well differentiated and among-group variation was maximized when L. argentatus and L. hyberboreus populations were grouped by locality rather than species, indicating recent hybridization. Four populations are located in putative Pleistocene glacial refugia and had larger τ estimates than the other 28 populations. However, we were unable to substantiate these putative refugia using coalescent theory, as all populations had genetic signatures of stability based on mtDNA. The extent of haplotype and allele sharing among Arctic white-headed gull species is noteworthy. Studies of other Arctic taxa have generally revealed species-specific clusters as well as genetic structure within species, usually correlated with geography. Aspects of white-headed gull behavioral biology, such as colonization ability and propensity to hybridize, as well as their recent evolutionary history, have likely played a large role in the limited genetic structure observed.

Population and species divergence in Arctic breeding species often reflect the climatic oscillations of the Pleistocene (Hewitt 2004). Glacial activity had profound effects on the distribution of Arctic taxa through these major climatic shifts. During glacial maxima, ice sheets subdivided ancestral populations into temperate or high-latitude ice-free areas, often resulting in the formation of phenotypically similar species with shallow genetic differentiation (Schmitt 2007). Warming during interglacial periods allowed species to expand their distributions into newly available habitat, resulting in clinal variation in genetic diversity (Hewitt 2004). This pattern of expansion and contraction of species distributions occurred many times in the Pleistocene; more than 20 glacial cycles have been recorded (Williams et al. 1998). Following the glacial maxima, many species that were isolated in southern and northern refugia expanded, came into secondary contact, and hybridized to varying extents (Hewitt 2001). This likely resulted in cycles of isolation and hybridization throughout the Pleistocene. Concordance of secondary contact zones (suture zones) has been observed across several Arctic species, which suggests a commonality in the location and persistence of glacial refugia during the last glacial maximum (Hewitt 2000).
Over the past decade, molecular markers have aided in the identification of cryptic glacial refugia and substantiated previously hypothesized refugia (Waltari and Cook 2005;Schafer et al. 2010;Sonsthagen et al. 2011). Previously, identification of the location of Pleistocene glacial refugia required knowledge of species distributions prior to glaciations or were inferred from paleoecological data. However, population contractions and expansions as a result of glacial cycling left predictable genetic signatures (Avise 2000). Populations arising via postglacial colonization of a region through successive founder events are expected to show reduced genetic diversity relative to populations residing in nonglaciated areas and to exhibit a genetic signature of population expansion from low-diversity founder populations (Galbreath and Cook 2004). However, current or past hybridization among closely related taxa may make it difficult to assess genetic relationships among Arctic populations. Introgression would likely have maintained or increased genetic diversity when a recently deglaciated area was colonized, and, therefore, would not be expected to produce a genetic signature of population expansion.
White-headed gulls (Larus spp.) are a geographically widespread clade of 18 species (Liebers et al. 2004;Olsen and Larsson 2004;Pons et al. 2005). Included in this clade is a subclade of 13 very closely related species of large white-headed gulls (Pons et al. 2005), which present particularly vexing problems to biologists. Some species within the white-headed gull complex have a circumpolar distribution (L. argentatus, L. canus, and L. hyperboreus; Fig. 1), while others are restricted to more circumscribed areas at high latitude (L. glaucoides, L. schistisagus, and L. thayeri;Fig. 1;Olsen and Larsson 2004). Previous assessments of relationships among white-headed gull populations revealed low to moderate genetic structure based on mitochondrial DNA (mtDNA), microsatellite, and amplified fragment length polymorphism (AFLP) loci (Liebers et al. 2001;Liebers and Helbig 2002;Crochet et al. 2003;Liebers et al. 2004;Pons et al. 2005;Gay et al. 2007;Vigfúsdóttir et al. 2008;Sternkopf et al. 2010). However, these studies focused largely on populations in Europe where climatic oscillations of the late Quaternary were significantly different than those in North America (Hewitt 1996), notably in the extensive glacial advances in North America (Velichko et al. 1997) and the likely absence of long-term high-latitude glacial refugia in Europe (Schmitt 2007). The presence of several white-headed gull species restricted to northern latitudes, coupled with the relatively low genetic differentiation observed among taxa, suggests that glacial oscillations associated with the late Pleistocene may have played a large role in the diversification of this group.
We studied the influence of glacial oscillations on the genetic structure of seven species of white-headed gull that breed at high latitudes (L. argentatus, L. canus, L. glaucescens, L. glaucoides, L. hyperboreus, L. schistisagus, and L. thayeri) using microsatellite genotypes from 11 autosomal loci, intron sequences from six autosomal nuclear genes, and mtDNA sequences from the control region. We evaluated Holarctic localities that have been hypothesized as ice-free areas or glacial refugia in other Arctic vertebrates, including the southern edge of the Bering Land Bridge, northern Beringia, Haida Gwaii, Newfoundland Bank, Spitsbergen Bank, and northwest Norway. Specifically, we employed traditional frequency-based and coalescent-based analyses to test if populations residing at high latitudes have the genetic signature of refugia. Populations formed through postglacial colonization are characterized by lower levels of nucleotide and haplotype diversity (Avise 2000), later times of expansion relative to other sampled populations, and genetic signatures of population growth based on the coalescent (Lessa et al. 2004). Inclusion of multiple taxa that occupy the Arctic allow us to examine whether geographically concordant contact regions suggestive of secondary contact are observed among populations expanding out of different Pleistocene refugia (i.e., suture zones) (Avise 2000).

Sampling and DNA extraction
Tissue samples of breeding-season adults, representing seven species and 32 populations of white-headed gulls, were collected or obtained through tissue loans ( Fig. 1; Appendix S1): L. argentatus argentatus (Norway), L. a. argenteus (France, Iceland, and United Kingdom), L. a. smithsonianus (Canada and United States), L. canus brachyrhynchus (Canada and United States), L. c. canus (Sweden and United Kingdom), L. c. kamtschatschensis (Russia), L. glaucescens (Canada and United States), L. glaucoides kumlieni (Canada), L. hyperboreus barrovianus (United States), L. h. hyperboreus (Canada, Greenland, Iceland, and Norway), L. h. pallidissimus (Russia), L. schistisagus (Russia and United States), and L. thayeri (Canada and United States). Because of the limited number of breeding individuals of L. thayeri in tissue collections, nonbreeding adults of this species were included in this study. Care was taken to ensure that plumage characteristics were consistent with pure species, given the tendency for hybridization in this group (Pierotti 1987;Olsen and Larsson 2004 and citations therein). Species classifications follow the American Ornithologists' Union Checklist of North American Birds (Banks et al. 2008 (Tha). Extent of the most recent glacial ice sheets is illustrated in white, and unglaciated regions are illustrated in gray (Hewitt 2004). Sample sizes are in parentheses. See Appendix S1 for physical descriptions of the localities. 1280 c 2012 The Authors. Published by Blackwell Publishing Ltd. (Olsen and Larsson 2004). Total genomic DNA was extracted from each sample using an AutoGen animal tissue extraction kit (AutoGen, Holliston, Maine). Genomic DNA concentrations were quantified using spectrophotometry and diluted to 50 ng μL -1 working solutions.

Microsatellite genotyping
Twenty-four individuals were screened at 30 microsatellite loci known to be variable for gull species (Laridae). Of these, 11 polymorphic loci containing dinucleotide repeat motifs were selected for further analyses of all tissue samples: Hg16, Hg18, Hg25 (Crochet et al. 2003), K16 (Tirard et al. 2002), LarZAP12, LarZAP19, LarSNX24, LarZAP26 (Gregory andQuinn 2006), Rbg13, Rbg18, andRbg29 (Given et al. 2002). Polymerase chain reaction (PCR) amplifications followed Sonsthagen et al. (2007) with two modifications. The forward primer was end-labeled with one of two fluorescent phosphoramidite dyes (FAM or HEX). Fluorescently labeled PCR products were electrophoresed on an automated DNA sequencer (ABI 3130XL; Applied Biosystems, Foster City, CA) and sized using GENEMAPPER R version 4.0 (Applied Biosystems) with a universal ROX-labeled size standard (DeWoody et al. 2004). Ten percent of the samples were amplified and sized in duplicate for quality control purposes.

MtDNA and nuclear intron sequencing
We followed Liebers et al. (2001) and amplified a 2500 base pair (bp) fragment of the mtDNA genome to avoid amplifying nuclear pseudogenes observed in this group. From this, we directly sequenced 430 bp of domain I of the control region. Twenty nuclear autosomal introns were screened for variability and six selected for further analysis: α-enolase intron 8, ghrelin (ghrel) intron 3, ornithine carboxylase (od) intron 7, clathrin heavy-chain (chc) intron 5, myelin proteolipid protein (mpp) intron 4, and glyceraldehyde-3-phosphate dehydrogenase (gapdh) intron 11. Ghrel had two insert/deletions. To obtain sequence information from the entire fragment for individuals that are heterozygous for both insert/deletions, two internal sequencing primers were developed. See Appendix S2 for primer information. PCR amplifications, cycle sequencing, and postsequencing protocols followed Sonsthagen et al. (2007). ExoSAP-IT R (USB Corporation, Cleveland, OH) was used to remove excess primers and dNTPs in PCR products. Sequences are accessioned in GenBank (JQ708216-JQ710335).

Estimation of genetic diversity
Allelic phases of nuclear introns were inferred from diploid sequence data using PHASE 2.0 (Stephens et al. 2001). PHASE uses a Bayesian approach to reconstruct haplotypes from population genotypic data and allows for recombination and the decay of linkage disequilibrium with distance. The accuracy of haplotypes reconstructed by PHASE has been validated and shown to be greater than that of cloning with large datasets (Harrigan et al. 2008). The PHASE analysis (1000 iterations with a 1000 iteration burn-in period) was repeated three times and was consistent across runs. MtDNA and nuclear intron sequences were analyzed in NETWORK 4.5.0.0 (Fluxus Technology Ltd. 2008) using the median joining method (Bandelt et al. 1995), to illustrate possible reticulations in the gene trees because of homoplasy or recombination.
We calculated allelic frequencies, inbreeding coefficient (F IS ), and expected and observed heterozygosities for each microsatellite locus, mtDNA, and the six nuclear introns in FSTAT 2.9.3 (Goudet 1995). Hardy-Weinberg equilibrium (HWE) and linkage disequilibrium were tested in FSTAT 2.9.3 for microsatellite and nuclear intron loci, adjusting for multiple comparisons using Bonferroni corrections (α = 0.05). We verified the selective neutrality for mtDNA control region sequence data using Tajima's D (Tajima 1989), implemented in ARLEQUIN 3.11 (Schneider et al. 2000).

Detecting spatial genetic structure
Levels of population structure among sampled sites were assessed with pairwise F ST , R ST , ST , overall F-statistics, and R-statistics calculated in ARLEQUIN, adjusting for multiple comparisons using Bonferroni correction (α = 0.05). We used Hedrick's (2005) method, implemented in Recode-Data version 1.0 (Meirmans 2006) to calculate the maximum value of F ST obtainable for our suite of microsatellite loci. Interallelic and interhaplotypic sequence divergences were used to calculate pairwise ST (Excoffier et al. 1992), and nuclear intron alleles were paired by individual. MODEL-TEST 3.06 (Posada and Crandall 1998) was used to determine the minimum parameter nucleotide substitution model that best fit the nuclear intron and mtDNA sequence data under Akaike's information criterion (Appendix S2; Akaike 1974).
Genotypic nuclear data (microsatellite and intron) were analyzed in STRUCTURE 2.2.3 (Pritchard et al. 2000) to detect the occurrence of population structure without a priori knowledge of putative populations. A series of analyses were performed (1) among large white-headed gull individuals (excludes L. canus individuals; Pons et al. 2005), and (2) within species represented by multiple populations. Data were analyzed using an admixture model assuming correlated frequencies to probabilistically assign individuals to putative populations (parameters: burn-in 10,000 iterations; 500,000 Markov chain Monte Carlo iterations) with the possible populations (K ) ranging from 1 to 15. Analyses were repeated five times and were consistent across runs. We used the method c 2012 The Authors. Published by Blackwell Publishing Ltd.
of Evanno et al. (2005) to determine the most likely number of clusters.
Hierarchical analyses of molecular variance (AMOVA) were conducted in ARLEQUIN to test for significance of geographic and taxonomic (subspecies and species) partitioning of a priori hypothesized genetic units using microsatellite, nuclear intron, and mtDNA loci. Populations were grouped to test (1) specific designations, (2) subspecific designations, (3) geographic proximity irrespective of species status, and (4) geographic proximity and species status. We assumed that groupings that maximized the among-group variance ( CT ) and were significantly different from random distributions constituted the most probable subdivision (Sonsthagen et al. 2011). L. glaucoides, L. schistisagus, and L. thayeri were not included in AMOVA comparisons because these species were represented by a single population.

Estimation of population demography
Evidence for historical fluctuations in population size was evaluated for 11 microsatellite loci using BOTTLENECK 1.2.02 (Cornuet and Luikart 1996) and for sequence data using LAMARC 2.1.2b (Kuhner 2006;Kuhner and Smith 2007). Fluctuations in population size inferred from microsatellite data were assessed using a Wilcoxon sign-rank test in BOTTLENECK. The probability distribution was established using 1000 permutations under two models: stepwise mutation model (SMM) and two-phase model of mutation (TPM; parameters: 79% SMM, variance 9; Rousset 1996). Heterozygote deficiency relative to the number of alleles indicates recent population growth, whereas heterozygote excess indicates a recent population bottleneck (Cornuet and Luikart 1996). It is important to note that BOTTLENECK compares heterozygote deficiency and excess relative to genetic diversity, not to HWE expectation (Cornuet and Luikart 1996). LAMARC was run using Bayesian search parameters: 10 short chains (1000 trees used out of 20,000 sampled) and three long chains (10,000 trees used out of 2,000,000 sampled). Data were analyzed three times and parameters converged across runs. Finally, mismatch distributions of mtDNA haplotype data were calculated in ARLEQUIN to gain further insight into historical population demography (Rogers and Harpending 1992).

Genetic diversity
Multilocus microsatellite genotypes were collected from 343 individuals representing seven species. The number of alleles per locus ranged from 4 to 18. Allelic richness ranged from 1.62 to 2.86 with a mean of 2.62 across all populations (Table 1). Observed heterozygosities ranged from 22.7% to 70.7%; the mean across all populations was 52.0% (Table 1).
In general, L. hyperboreus and populations of L. argentatus from Europe had lower levels of heterozygosity than did other species. Two populations (AeFrc and CanNWT) exhibited heterozygote deficiency and did not conform to HWE (Table 1). The population of L. canus from the Northwest Territories (CanNWT) was in linkage disequilibrium at nine loci pairs (Lar24xK16, Lar24xHg18, Lar12xLar19, Lar12xK16, Lar12xHg18, Lar12xRbg29, Lar26xHg18, Lar19xHg18, K16xHg18), but the overall comparison was not significant. The remaining populations and loci were in linkage equilibrium and HWE and all loci were retained for subsequent analyses. The inbreeding coefficient (F IS ) ranged from -0.098 to 0.286 across populations; mean value was 0.112. The inbreeding coefficient for CanNWT was significantly larger than expected (α > 0.05) ( Table 1).
Nuclear introns were 323-665 bp in length and contained 15-46 variable sites (Appendix S2). PHASE reconstructed 24-117 alleles for the individual introns ( Fig. 2A-F; Appendix S2). Probabilities of reconstructed haplotypes ranged primarily from 0.80 to 1.00, although a minority of individuals had probabilities ranging from 0.50% to 0.78 (5% of individuals for chc, 2% for enolase, 12% for gapdh, 2% for ghrel, 17% for mpp, and 2% for od7). We attribute the lower probabilities of reconstructed haplotypes for gapdh and mpp to the high occurrence of autapomorphies (single novel polymorphisms occurring on one allele in one individual) in these loci; 30% and 18% of individuals had novel polymorphisms for gapdh and mpp, respectively. Private alleles were observed for most species at most loci; however, private alleles were only observed in two (enolase and gapdh) or three (chc, enolase, and gapdh) nuclear introns, respectively, for L. glaucoides and L. thayeri ( Fig. 2A-F). Observed heterozygosities ranged from 47.7% to 77.3%, with a mean of 61.3% across all populations (Table 1). Populations and loci were in linkage equilibrium and HWE. The inbreeding coefficient (F IS ) ranged from -0.171 to 0.311, mean value was 0.088. Inbreeding coefficients for CanAnc and CanRus were significantly larger than expected (α > 0.05) ( Table 1). Haplotype (h) and nucleotide (π) diversity ranged from 0.957 to 1.000 and 0.000 to 0.009, respectively (Table 1).
We assayed a 392 bp fragment of the mtDNA control region characterized by 54 variable sites among 134 unique haplotypes (Appendix S2; Fig. 2G). Number of haplotypes per population ranged from 2 to 14 (mean = 6.48; Table 1). Private haplotypes were observed for all species (Fig. 2G). Haplotype (h) and nucleotide (π) diversity were high for most populations, with values ranging from 0.286 to 1.000 and 0.001 to 0.024, respectively (Table 1)   Spatial genetic structure

Analyses among species
Overall estimates of population subdivision were significant across all marker types ( Table 2). RecodeData calculated an upper F ST limit of 0.415 for the microsatellite data. Therefore, the overall F ST of 0.129 accounts for 31.1% of the maximum possible level of genetic structure. Moderate levels of population structure were observed among species at the 11 microsatellite loci and six nuclear loci (Table 2; Appendix S3). Most of the significant interpopulation comparisons among species were observed between L. glaucescens populations and all other populations and between L. canus populations and all other populations; values were typically higher when a mutation model was applied to the dataset (R ST and ST ; Appendix S3). In contrast to the nuclear data, high levels of variance in mtDNA haplotypic frequency were observed between most population pairs. As with the nuclear data, most interpopulation comparisons between L. canus and all other species were significant (Appendix S3). Our STRUCTURE analyses indicated that K was maximized among all large white-headed gull individuals when K equaled 2 ( K = 327.5) and 4 ( K = 16.4) for the microsatellite and nuclear intron loci, respectively (Fig. 3). Patterns of individual assignment differed slightly between marker types: L. glaucescens and L. schistisagus individuals were assigned predominantly to one cluster and L. argentatus and L. hyperboreus individuals to the other cluster based on microsatellite data (Fig. 3A). Gulls representing L. thayeri, L. glaucoides, and northern populations of L. argentatus (AaTrm, SmiAK, SmiMN, VegCh, VegYa) and L. hyperboreus (HypYKD) were assigned in approximately equal proportions to both clusters based on microsatellite data (Fig. 3A), a signal consistent with hybridization at the northern locales. Similarly, L. glaucescens individuals were predominantly assigned to a single cluster (blue) based on the nuclear intron data, with L. hyperboreus and some individuals representing L. argentatus (AeFrc, AeIc, AaTrm) predominantly assigned to another cluster (red) (Fig. 3B). The remaining individuals were assigned to the four clusters in approximately equal portions (Fig. 3B).
Partitions in the nuclear and mitochondrial genomes appear to conform to subspecific classifications (Fig. 4). Among-group variation (F CT ) based on the microsatellite data was maximized when populations were grouped by subspecies or a combination of subspecies and geographic proximity for F ST and R ST estimates, respectively (Fig. 4). Similarly, among-group variation based on nuclear intron and mtDNA data (F ST and ST ) was maximized by subspecies. Because the high level of structure between L. glaucescens and L. canus populations and all other populations is likely driving variance estimates, we ran additional AMOVAs including only populations of L. argentatus and L. hyperboreus, the other two species represented by multiple populations, to gain additional insight into the partitioning of genetic variation in these species. Among-group variation, as estimated from microsatellite (F ST and R ST ) and nuclear intron frequency data (F ST ; Fig. 4), was maximized when L. argentatus and L. hyperboreus populations were grouped by geographic proximity regardless of specific or subspecific classification. This is suggestive of contemporary hybridization. In contrast, among-group variance estimates calculated from the mtDNA (F ST and ST ) and nuclear intron ( ST ) sequence data were maximized when populations were grouped by subspecies (Fig. 4).

Analyses within species
Variance in microsatellite and nuclear intron allelic frequencies (F ST , R ST , and ST ) ranged from 0.030 to 0.114 (Table 2). In contrast to interspecies comparisons, F ST values were typically greater than R ST values for the microsatellite data, with a few exceptions: between northern and southern populations of L. glaucescens, between the L. argentatus population from Tromsø and other L. argentatus populations, and between populations of L. canus from North America and from Russia populations (Appendix S3). However, we observed few significant interpopulation comparisons based on the nuclear intron data (Appendix S3). Variance in mtDNA haplotypic frequencies (F ST ) calculated over all populations ranged from 0.080 to 0.537, with larger values observed for most species when the nucleotide substitution model was applied to the data set ( ST = 0.076-0.406; Table 2). Across all within taxa analyses, the likelihood generated for the nuclear intron genotypic data was maximized when K equaled 1. In contrast, genetic partitioning was observed based on microsatellite data. Within L. argentatus, K was maximized when K equaled 3 ( K = 223.7) (Fig. 5A). Larus argentatus individuals from Alaska and Russia clustered together (blue), individuals breeding on the eastern coast of North America grouped together (yellow), and individuals from Iceland, Tromsø, and France clustered together (red) (Fig. 5A). Larus argentatus individuals from the Northwest Territories were assigned to the blue and yellow clusters in approximately equal frequencies (Fig. 5A). Within L. canus, K was maximized when K equaled 2 ( K = 368.4); these two groups corresponded to subspecific classifications, with individuals from North America assigned to one cluster and those from Sweden and Russia to the other (Fig. 5B). K also was maximized when K equaled 2 ( K = 452.9) within L. glaucescens, although the assignment of individuals did not correspond to sample locality (Fig. 5C). For L. hyperboreus, K was maximized when K equaled 4 ( K = 96.9) (Fig. 5D). The assignment of individuals appeared to correspond loosely to sample locality, although the signal was not strong; individuals from the Yukon-Kuskokwim Delta were assigned predominately to the green cluster, those from Iceland to the yellow cluster, those from northern Alaska and Baffin Island to the blue cluster, and individuals from Greenland and Svalbard to the red cluster (Fig. 5D).

Population demography
Evidence for significant fluctuations in historical population demography was detected based on microsatellite genotypes. Under the SMM (Table 3), population growth (heterozygote deficiency) was observed for populations of L. schistisagus from Kamchatka Peninsula; L. glaucoides; L. argentatus from Tromsø, France, and Iceland; L. hyperboreus from Greenland, Iceland, Chukotka Peninsula, and Yukon-Kuskokwim Delta; L. glaucescens from Aleutian Islands, Homer, and Middleton Island; and L. canus from south-central Alaska and Northwest Territories. Similar results were observed under the TPM, though fewer populations had signatures of population growth (Table 3).
Significant population growth based on nuclear intron sequences was detected, using LAMARC, for all populations, with theta (4N e μ) ranging from 0.002 to 0.014 (Table 3). In contrast, all populations had a signal of population stability when g was estimated from mtDNA sequence data, consistent with a pattern of populations located in glacial refugia (Lessa et al. 2004). Theta (2N f μ)   ( Table 3). Mismatch distributions estimated from mtDNA sequence data failed to reject the sudden expansion model, based on Harpending's raggedness index (Harpending 1994), for any population. Parameter estimates for time of expansion (τ ) ranged from 0.0 to 16.5 (Table 3). The smallest values were observed for populations of L. hyperboreus from Iceland, L. argentatus from Minnesota, Tromsø, and France, and L. glaucoides; the largest estimates for populations of L. glaucescens from Middleton Island, L. hyperboreus from Greenland, L. argentatus from Maryland and Iceland, and L. canus from south-central Alaska (Table 3).

Multilocus genetic structure within and among species
Despite extensive allele and haplotype sharing among whiteheaded gull species, genetic substructure was observed within and among species across all marker types. Species and populations at high latitude exhibited lower genetic differentiation then their southern counterparts. Furthermore, individuals breeding at northern latitudes clustered together regardless of species designation, consistent with contemporary hybridization. At least two haplotype/allele groups were observed at each locus; however no haplotype/allele group was represented by a single species at any of the loci. Private alleles and haplotypes were observed for most species at most loci; however, private alleles were only observed in two or three nuclear introns, respectively, for L. glaucoides and L. thayeri (Fig. 2).
Genetic evidence for contemporary hybridization among northern populations of Arctic white-headed gulls is corroborated by field reports (e.g., L. glaucescens × L. hyperboreus; L. argentatus × L. hyperboreus; L. argentatus × L. glaucescens; L. glaucescens × L. schistisagus; L. argentatus × L. schistisagus; L. glaucoides × L. thayeri; Olsen and Larsson 2004 and  Table 3. Results of demographic analyses for 11 microsatellite loci under the stepwise mutation model (SMM), and two-phased model of mutation (TPM), and for sequence data from six introns and the mitochondrial control region calculated from assayed Arctic gull populations. Parameter estimates θ (4N e μ for nuclear DNA, 2N f μ for mtDNA), exponential growth rate (g), and time of expansion (τ ) calculated from mismatch distributions with standard deviations (SD) are provided for each Arctic gull population.  citations therein). Hybridization would be expected to homogenize allelic frequencies by locality, as neutral loci will remain similar because of introgression and recombination (Mallet 2005). Species appear to have been isolated long enough to have accumulated unique mutations, as indicated by the partitions in the nuclear and mtDNA genomes. Therefore, we contend that hybridization has occurred only recently in Arctic white-headed gull evolutionary history, likely from secondary contact following contemporary range expansion. Introgression of species-specific alleles may be maintained through local adaptation to intermediate habitat types where species coexist, as hybrids have been reported to display adaptive traits of both parental species (L. glaucescens × L. occidentalis; Good et al. 2000). Recent speciation and contemporary hybridization likely both play a role in the magnitude of allele and haplotype sharing observed among white-headed gulls. Of particular interest is the extent of introgression/hybridization occurring at the northern limits of species' ranges and among white-headed gulls that breed exclusively at high latitudes. Long-term stable hybrid zones have been reported in temperate areas for several white-headed gull taxa (L. occidentalis × L. glaucescens, Bell 1996Bell , 1997Good et al. 2000;L. argentatus × L. marinus, Crochet et al. 2003) and may be maintained by hybrid superiority at the hybrid zone (Moore 1977). In contrast, white-headed gull species appear to hybridize pervasively throughout northern latitudes (L. argentatus × L. hyperboreus, Vigfúsdóttir et al. 2008;L. argentatus × L. glaucescens, Williamson and Peyton 1963;L. argenatus × L. schistisagus, Olsen and Larsson 2004;L. glaucoides × L. thayeri, Weir et al. 2000;L. glaucescens × L. schistisagus, Olsen and Larsson 2004) and discrete hybrid zones appear to be absent. Differences in the degree of hybridization may be attributable to the stability of the habitat where these areas of secondary contact occur. Stable secondary contact zones for gulls are observed at temperate latitudes, where presumably habitat has remained relatively stable throughout the last glacial maximum, allowing species to diverge in allopatry without coming into secondary contact during interglacial periods. Conversely, Arctic species reside in more stochastic environments, where suitable habitat repeatedly contracted and expanded during the Pleistocene glacial cycles (Hewitt 2004). These highly variable climatic conditions likely resulted in a cycle of isolation during glacial periods and secondary contact during interglacial periods, potentially limiting species divergence and development of preand postzygotic isolating mechanisms.
Comparatively higher estimates of population structure observed for mtDNA than nuclear DNA markers are consistent with Haldane's rule. Haldane's rule states that hybrids of the heterogametic sex will experience reduced fitness (i.e., greater inviability or sterility) relative to those of the homogametic sex (Coyne and Orr 2004). In birds, females are heterogametic; therefore, if hybrid females were experiencing a strong disadvantage relative to hybrid males, observed genetic differentiation would be greater in mtDNA than nuclear markers. In a study of mainly European white-headed gull populations, researchers proposed that the large discrepancy in interspecific comparisons (mtDNA estimates were 3.3-14.5 times greater than estimates using microsatellites) between marker types could be attributed to a strong disadvantage for hybrid females (Crochet et al. 2003). We did observe high levels of interpopulation comparisons among species for mtDNA ( ST = 0.130-0.821; Appendix S3) with no significant comparisons at nuclear markers. However, Haldane's rule would presumably have the greatest influence at secondary contact zones. In these areas, 23% (5/22; Appendix S3) of the comparisons had significant mtDNA estimates. Of those five, only two comparisons (HypYKD × VegCh, HypIc × AeIc) had proportionally greater degree of divergence than can be explained by differences in the effective population size between genomes (Zink and Barrowclough 2008) and maximum possible F ST value (Meirmans 2006), although differences were slight (HypYKD × VegCh, F ST Exp. = 0.039, F ST Obs. = not significant; HypIc × AeIc, F ST Exp. = 0.051, F ST Obs. = not significant). Therefore, Haldane's rule does not appear to apply to the white-headed gull species studied here.
The lower F ST values observed at nuclear fragment assays (i.e., microsatellites) relative to mtDNA sequence data among species may be attributable to fragment length homoplasy, through not identifying unique alleles because fragments of the same length may have different sequences or may have mutated back to the ancestral state. Both types of homoplasy could pose problems when assessing population structure with fragment analysis based on detecting allelic frequency differences among populations. Most interpopulation comparisons among species have higher R ST than F ST estimates, indicating that the mutation process, and therefore homoplasy, is having an effect on estimators of population subdivision. Caution should be taken when interpreting pairwise population comparisons of allelic variance among species. Rousset (1996) showed that there are no simple effects of homoplasy on estimators of population differentiation (F ST and R ST ) for loci evolving under the SMM and island model of migration, making it difficult to assess potential biases in estimates. However, we observed similar signals of population structure based on microsatellite and nuclear intron data. Therefore, the differences in the degree of population structure observed between the genomes studied here may be more attributable to introgression reducing the rate of lineage sorting in the nuclear genome and, to a lesser extent, fragment length homoplasy, although experimental evidence is needed to test this hypothesis.

Pleistocene refugia and comparisons with other taxa
Two distinct mtDNA haplotype groups were observed within L. argentatus, L. canus, and L. hyperboreus, the three sampled taxa with circumpolar distributions, consistent with other studies on white-headed gulls ( Fig. 2G; Liebers et al. 2004;Sternkopf et al. 2010). A pattern of at least two allele groups was also observed for the nuclear intron loci ( Fig. 2A-F). Concordance in haplotype and allele groups suggests that white-headed gulls were subdivided into at least two refugia that persisted for extended periods of time during the Pleistocene. Furthermore, substructuring observed within mtDNA corresponds to locality. The small central haplotype group is represented by L. argentatus individuals from Iceland and Tromsø and L. hyperboreus individuals from Iceland, Greenland, Svalbard, and a single individual from northern Alaska (Fig. 2G, population data not shown). All populations of L. argentatus are represented in the large haplotype group; however, only North American and Greenland L. hyperboreus individuals (and a single individual from Iceland) are observed within this group (Fig. 2G, population data not shown). These findings differ from those of Sternkopf et al. (2010), who identified genetic similarity between L. hyperboreus and L. a. smithsonianus (the North American subspecies); we did not observe a haplotype group restricted to North America in L. argentatus. The presence of a primarily Scandinavian/Greenland/Iceland haplotype group indicates the restriction of at least L. argentatus and L. hyperboreus into a high-latitude refugium in the North Atlantic/Arctic Ocean, possibly Spitsbergen Bank or northwest Norway. Furthermore, given the restricted geographical distribution of L. hyperboreus haplotypes within the central clade, the Scandinavian/Greenland/Iceland refugium was likely isolated from other L. hyperboreus populations and did not substantially contribute to the postglacial colonization of North America and Europe.
Despite the presence of species with distributions restricted to northern latitudes, suggestive of restriction to Pleistocene refugia, we were unable to identify glacial refugia for the white-headed gulls studied here based on the coalescent. White-headed gulls are characterized by strong dispersal ability and a propensity to hybridize in areas of secondary contact, as is reflected in contemporary accounts of long-range colonization and subsequent hybridization (Pierotti 1987;Olsen and Larsson 2004;Vigfúsdóttir et al. 2008). The tendency for hybridization at areas of secondary contact is very strong: 16 of the 18 species (89%) are reported to hybridize in nature (Pierotti 1987;Olsen and Larsson 2004 and citations therein) and appear to be free of postzygotic barriers to hybridization (Shields 1987;Snell 1991). Hybridization and subsequent introgression may have erased the genetic legacy of the Pleistocene for white-headed gulls; reductions in effective population sizes associated with the restriction to glacial refugia were not observed, even for populations currently located in glaciated areas. Alternatively, historic Arctic white-headed gull populations that were restricted south and north of the ice sheets likely followed habitat made available by the retreating glacial ice sheets to present day locations. Short movements from refugia would have allowed these historical populations to retain genetic diversity because effective population sizes would not be reduced (Hewitt 1996), especially if colonization occurred over a long period. However, it is unlikely that all populations studied here colonized slowly subsequent to glacial retreat, given the dispersal and colonization ability of Arctic white-headed gulls. Therefore, a more parsimonious explanation is that the strong tendency for hybridization in this group erased the genetic signature of Pleistocene refugia. Holarctic species typically exhibit shallow but clear phylogeographic signal, due to fragmentation of the distributions of these species into high-latitude glacial refugia, resulting in genetic and morphological differentiation across ranges (Hewitt 2004). Although Arctic white-headed gulls have morphologically distinct forms across their distribution (i.e., subspecies), limited population genetic subdivision among northern Arctic white-headed gull populations was observed. The combination of limited sorting among species (Liebers et al. 2004;Pons et al. 2005) and morphological diversity of taxa, as is evident in gulls, is suggestive of recent allopatric fragmentation and restriction to multiple glacial refugia. Indeed, four populations of gulls (L. hyperboreus from Greenland, L. glaucescens from Middleton Island, L. argentatus from Iceland, and L. canus from south-central Alaska) have larger time-of-divergence estimates than other sampled populations and coincide with proposed high-latitude glacial refugia for other taxa (Ploeger 1968). Historical and contemporary hybridization among regions, coupled with small effective population sizes, could have overwhelmed the genetic uniqueness accumulated by northern populations during the last glacial maximum, producing a signal of low genetic structure.
Arctic white-headed gulls are less genetically differentiated at the northern limits of their distribution; this pattern was observed both among species that breed exclusively at high latitudes (L. schistisagus, L. glaucoides, L. thayeri, and L. hyperboreus) and within species (e.g., L. argentatus and L. glaucescens). Liebers and Helbig (2002) observed this pattern between L. fuscus and L. cachinnans; L. fuscus exhibited less genetic structure than its southern counterpart L. cachinnans. The authors proposed that differences in the degree of population subdivision indicated that northern gulls are phylogenetically younger than their southern counterparts and were more affected by glacial cycles. In contrast, their southern relatives would have been able to maintain larger long-term stable population sizes during glacial periods. Although Liebers and Helbig's (2002) hypothesis is consistent with observations within Europe, it may not apply to the North American species studied here. In contrast to Europe, high-latitude glacial refugia, notably Beringia, were present in North America during the last glacial maximum and likely promoted genetic diversification in Arctic taxa (Hewitt 2004). As glacial ice sheets retreated, the ranges of temperate species likely expanded northward and came into contact with northern "refugial" populations. Hybridization between "newly arriving" temperate species and northern "refugial" gull species provides an alternate hypothesis for the genetic pattern of decreasing genetic differentiation with increasing latitude.