Flyway structure in the circumpolar greater white‐fronted goose

Abstract Dispersal and migratory behavior are influential factors in determining how genetic diversity is distributed across the landscape. In migratory species, genetic structure can be promoted via several mechanisms including fidelity to distinct migratory routes. Particularly within North America, waterfowl management units have been delineated according to distinct longitudinal migratory flyways supported by banding data and other direct evidence. The greater white‐fronted goose (Anser albifrons) is a migratory waterfowl species with a largely circumpolar distribution consisting of up to six subspecies roughly corresponding to phenotypic variation. We examined the rangewide population genetic structure of greater white‐fronted geese using mtDNA control region sequence data and microsatellite loci from 23 locales across North America and Eurasia. We found significant differentiation in mtDNA between sampling locales with flyway delineation explaining a significant portion of the observed genetic variation (~12%). This is concordant with band recovery data which shows little interflyway or intercontinental movements. However, microsatellite loci revealed little genetic structure suggesting a panmictic population across most of the Arctic. As with many high‐latitude species, Beringia appears to have played a role in the diversification of this species. A common Beringian origin of North America and Asian populations and a recent divergence could at least partly explain the general lack of structure at nuclear markers. Further, our results do not provide strong support for the various taxonomic proposals for this species except for supporting the distinctness of two isolated breeding populations within Cook Inlet, Alaska (A. a. elgasi) and Greenland (A. a. flavirostris), consistent with their subspecies status.

routes; however, individual migratory strategies can be flexible and, in some instances, can alter rapidly (Jonker et al., 2013;Pulido, 2007;Rolshausen et al., 2013;Sutherland, 1998). How this flexibility influences the genetic composition of populations will largely depend not only on the relative frequency of flyway switching (abmigration) but ultimately if these observed migratory irregularities lead to homogenization of previously separated populations (Rockwell & Barrowclough, 1987).
In northern high latitudes, most waterfowl species are highly mobile and migrate seasonally from nesting areas at higher latitudes during the summer months to areas at lower latitudes during winter months. Banding, telemetry, bird counts throughout the year, and morphological data have led to the identification of major migratory flyways, which are an integral part of management strategies, particularly in North America. However, the boundaries between these flyways are not always discrete and fidelity to these migratory flyways varies within and across taxonomic groups (Baldassarre, 2014;Ely & Scribner, 1994;Guillemain, Sadoul, & Simon, 2005;Lavretsky, Miller, Bahn, & Peters, 2014;Madsen, Tjørnløv, Frederiksen, Mitchell, & Sigfússon, 2014).
Although observational data, such as the distribution of band recoveries, frequently suggest high migratory connectivity, it is relatively unknown in many waterfowl species if fidelity to migration flyway reflects philopatry (natal-and breeding-site fidelity), which would promote genetic structure. Contrasting patterns in structure ascertained from genetic information and observational data have been uncovered for many avian species (e.g., Koenig, van Vuren, & Hooge, 1996;Kraus et al., 2014;Liu, Keller, & Heckel, 2012;Pearce et al., 2014), such that although observational data revealed little or no interchange, genetic data showed limited (or no) genetic signal of flyway structure. Lack of correspondence between genetic structure and observational data has been attributed to mainly male-biased dispersal in cases where annual pair formation occurs on the winter grounds, and this has provided a mechanism enabling genetic interchange among breeding locales (see, e.g., Peters & Omland, 2007;Wilson, Gust, Petersen, & Talbot, 2016). Waterfowl are harvested largely in wintering areas, yet population counts (upon which management decisions are based) are typically conducted on breeding areas. Thus, an understanding of the strength of the relationship between philopatry and fidelity to flyway is integral for species management and maintaining the future viability of populations (Baldassarre, 2014).
Geese occupying northern high latitudes exhibit life-history traits that may facilitate population structure and restrict interflyway genetic exchange, such as high philopatry in both sexes along with long-term pair bonds and familial associations and delayed reproduction (Ely & Scribner, 1994;Scribner et al., 2003;Ely, Wilson, & Talbot, 2017). In contrast to ducks, pair-bonding in some goose species is thought to occur primarily during the spring and summer when genetically similar individuals are segregated (Ely & Scribner, 1994;Leafloor, Moore, & Scribner, 2013;Weegman et al., 2015), and this would provide an additional mechanism to further limit gene flow among breeding areas. The greater white-fronted goose (Anser albifrons, Figure 1) is only one of two goose species with a nearly circumpolar distribution (the other being brant, Branta bernicla), and is comprised of populations that utilize five major flyways ( Figure 2). Across Eurasia and North America, there is considerable phenotypic variation (Ely et al., 2005), which has led to the naming of up to six morphological subspecies (Banks, 2011;Delacour, 1954;Mooij & Zöckler, 2000); Mooij and Zöckler (2000) also include an ecological component to their subspecies attributions. However, geographic distribution of some subspecies remains uncertain, and it is unclear whether there is correspondence between subspecies designations and genetic partitioning. The maintenance of subspecies boundaries and lack of interflyway banding recoveries suggest restricted gene flow among subspecies that loosely corresponds to flyway (see Ely & Dzubin, 1994 for summary). However, distinct phenotypic variants occur in sympatry at certain periods during the annual cycle. Within the Pacific Flyway, for example, three morphologically distinct populations, including the largest and smallest forms, overlap in migratory pathways and winter distribution (Ely & Takekawa, 1996;Ely et al., 2005;Orthmeyer, Takekawa, Ely, Wege, & Newton, 1995). Several mechanisms have been proposed in the maintenance of reproductive isolation of sympatric wintering populations, such as microgeographic and behavioral barriers (see Ely et al., 2017). Furthermore, greater white-fronted goose populations, at least in North America, tend to be spatially and temporally segregated during migration (Ely, Neiman, Alisauskas, Schmutz, & Hines, 2013;Ely & Takekawa, 1996). Further, diverse topography, as characterized by mountains and river drainages, promote temporal variation in timing of breeding and hence fall migration for breeding chronology. This is particularly evident within Alaskan breeding populations, which exhibit greater latitudinal variation in breeding chronology than other regions (Ely et al., 2005) that may promote structure both across and within flyways.
Previous population genetic studies on greater white-fronted geese have focused on either within flyway (Ely et al., 2017;Volkovsky, Fisenko, Gerasimov, & Zhuravlev, 2016;Volkovsky, Kulikova, Gerasimov, & Zhuravlev, 2013) or between Eurasia and Greenland populations (Eda, Shimada, Ushiyama, Mizota, & Koike, 2013). Here we present the most comprehensive population genetic assessment of the greater white-fronted goose to date across their entire breeding range, using genotypic data from eight autosomal microsatellite loci and sequence data from the mitochondrial DNA (mtDNA) control region, as well as contemporary observational data in the form of band recovery distribution through 2015. Specifically, we determined the level of genetic partitioning (continent, flyway, or within flyway) and used a Bayesian model of isolation with migration to test for the presence of intercontinental gene flow. As greater white-fronted geese exhibit high levels of migratory connectivity (>50% return rates to both winter and breeding grounds; Fox & Stroud, 1988;Wilson, Norriss, Walsh, Fox, & Stroud, 1991;Ely & Dzubin, 1994;Alisauskas & Lindberg, 2002) and have complex long-standing family structure, we predict there would be significant genetic structuring not only among flyways but also among breeding areas.

| Sample collection
Blood, feather, muscle tissue, or eggshell membranes were collected from greater white-fronted geese throughout most of their known breeding range, including North America (eight localities representing two major flyways and two to three subspecies in Alaska and seven localities representing midcontinent breeding areas in Canada), Asia (seven localities representing 1-2 subspecies) and Greenland (one subspecies; Figure 2).

| Study species
Within North America, the Pacific Flyway population of greater white-fronted geese is comprised of three distinct breeding regions representing two subspecies ( Figure 2): Cook Inlet (tule goose, A. a. elagsi; also known as A. a. gambeli or A. a. gambelli; see Banks, 2011 for summary of taxonomic history) and Yukon-Kuskokwim Delta and Bristol Bay region (A. a. frontalis). The two western Alaskan breeding locales have been recently proposed but not currently accepted as comprising a separate subspecies due to their small body size (A. a. sponsa;Banks, 2011) with the Bristol Bay population being of intermediate body size (Orthmeyer et al., 1995). The midcontinent breeding population is located across the taiga and arctic habitats of central and northern Alaska and Canada with migratory routes mainly within the Central Flyway but to a lesser extent in the Mississippi Flyway ( Figure 2). The midcontinent population is most often considered a single subspecies (A. a. frontalis).
Within Eurasia, the Khatanga River has been proposed as the geographic break between western and eastern Palearctic populations (Mooij & Zöckler, 2000); however, Delacour (1954) proposed this biogeographic break was further east, at the Kolyma River. F I G U R E 2 Sampling localities (A-W) of greater white-fronted goose used in this study (refer to Table 1 for location names). The species' breeding range is highlighted with shaded color corresponding to flyway designation and haplotypes presented in Figure 3. Suture zones at Lena River, Russia, and Mackenzie River Delta, Canada, are indicated by the letter S and J, respectively. These two locations also indicate the proposed western and eastern boundary of Beringia near the Khatanga River  Within western Beringia, the Lena River and Yana River are of particular interest in that these areas tend to have high biodiversity given they are situated in between Atlantic and Pacific flyways (Gilg et al., 2000). Although, this region is often grouped within the eastern Palearctic, results from a recent movement study (Li, Si, Ji, & Gong, 2017) showed these populations winter farther inland in China and possibly represent a more central Palearctic Flyway as Far East populations (e.g., Anadyr and Kolyma) tend to winter primarily in coastal regions of Japan and Korea. However, in agreement with Delacour (1954), observations of migrating geese suggested that the Lena River area might be partially composed of western Palearctic geese (Syroechkovskiy, 2006). Thus, the subspecific designation of geese in these regions, which include much of western Beringia, is still in debate. While some authorities consider that all Eurasian populations comprise a single subspecies (A. a. albifrons;Owen, 1980;Portenko, 1989;Banks, 2011), others (e.g., Delacour, 1954) place the eastern Palearctic populations within the North American subspecies (A. a. frontalis), based on similarity in body size. Mooji (2000) and Mooij and Zöckler (2000) further propose designating the eastern Palearctic population as a separate subspecies (A. a. albicans), based on wintering distribution, migratory routes, and slightly larger body size than western Palearctic areas. Geese nesting in Greenland and wintering primarily in Ireland and United Kingdom are designated as a separate subspecies (A. a. flavirostris) due to their large size, dark coloration and nonoverlapping breeding and wintering distribution with other potential subspecies.

| DNA isolation and sequencing
Genomic DNA was extracted from blood, muscle, feather, or eggshell membranes using a "salting out" procedure described by Medrano, Aasen, and Sharrow (1990), with modifications described in Sonsthagen, Talbot, and White (2004) for blood and muscle and in Talbot et al. (2011) for feathers and eggshell membranes. Genomic DNA concentrations were quantified using fluorometry and diluted to 50 ng mL -1 working solutions.

| Population subdivision
The degree of population subdivision among breeding areas was assessed by calculating pairwise F ST for mtDNA and microsatellite in ARLEQUIN, adjusting for multiple comparisons using Benjamini and Yekutieli-modified false discovery rate (α = 0.05; Benjamini & Yekutieli, 2001;Narum, 2006). Because the upper possible F ST value for a set of microsatellite loci is usually <1.0 (Hedrick & Goodnight, 2005), we used RECODEDATA, version 1.0 (Meirmans, 2006), to calculate the uppermost limit of F ST for a given data set.
We used two approaches to explore the genetic partitioning of genetic variation among and within breeding groups. We first used an analysis of molecular variance (AMOVA) in ARLEQUIN to test for significance of geographic partitioning of a priori hypothesized genetic units using mtDNA and microsatellite loci with statistical significance tested by 16,000 permutations. Populations were grouped to test (see Figure 1): (a) current and prior subspecific groupings, (b) geographic/flyway division, (c) nesting habitattundra vs. taiga and (d) proposed refugia (see Ploeger, 1968). As Greenland and Cook Inlet both only represented a single subspecies, we excluded these populations from subspecies groupings, as having a group represented by a single population will bias within-and between-group variance. As well, Greenland was excluded from major flyway groupings, as it is the only population utilizing its flyway. We assumed that groupings that maximize among-group variation (Φ CT or F CT ) and significantly different than random distribution of individuals were the most probable geographic divisions (p < 0.05).
Secondly, we used a Bayesian-clustering program, STRUCTURE 2.2.3 (Pritchard, Stephens, & Donnelly, 2000), to determine the level of population structure in the autosomal microsatellite data set without providing a priori information on the geographic origin of the individuals. If no structure was observed, the LOCPRIOR option was used as this model is able to detect population structure in datasets with a weak signal of structure not detectable under standard models (Hubisz, Falush, Stephens, & Pritchard, 2009 five times for each K and based on these preliminary results, 10 additional replicated were performed for K = 1-10 for a total of fifteen independent runs. We used the ∆K method of Evanno, Regnaut, and Goudet (2005) and evaluated the estimate of the posterior probability of the data given K, Ln P(D), to determine the most likely number of groups at the uppermost level of population structure.

| Demographic history and gene flow
Demographic histories of the greater white-fronted goose based on mtDNA sequence data were evaluated using two approaches: standard qualitative test statistics, Tajima's D and Fu's F s , and coalescent-based estimation implemented in IMa2 (Hey & Nielsen, 2004). To test for genetic signatures of recent effective population size changes, we calculated Fu's F s (Fu, 1997) and Tajima Table 1 for population groupings) as the breeding locales generally grouped by continent; albeit with low support (Supporting Information Appendix S1).
Preliminary results using different evolutionary divergence scenarios between major flyways yielded differing rates and directionality of gene flow, further suggesting that the simplified model of Old vs.
New World was most appropriate given the dataset.
Here, we estimated three different evolutionary parameters including ϴ (where ϴ = 4N e μ, and N e is the effective population size and μ is the mutation rate per site per generation) for each contemporary population and the ancestral populations, 2Nm (which is ϴm im /2, here m im is the migration rate relative to the mutation rate estimated in IMa2, N is the population size, and m is the proportion of the population consisting of immigrants per generation), and t (where t = T/μ, and T is the number of years since divergence).
We ran 20 Markov chains under a geometric heating scheme (option -hfg), with chain temperatures of 0.90 and 0.75. We ran preliminary analyses using large, flat priors for each parameter. From the results of these runs, we defined narrower upper bounds for each param- s is the expected adult survival rate (Saether et al., 2005). Greater white-fronted geese reach sexual maturity at age 3 (Baldassarre, 2014;Campbell & Goodwin, 1985), and adult survival rates averaged 0.749 of the Pacific flyway population (Schmutz & Ely, 1999). Using those values, we estimated G to be approximately 5.9 years, which corresponds to generation time estimates of other goose species (Dillingham, 2010).

| Compilation of band recovery data
We obtained banding and recovery data from the U. S. Geological

| Genetic diversity
One hundred and sixty-two unique mtDNA haplotypes were found     with an average of 8.38 (SD 4.14) alleles per locus (Table 1). Molecular diversity indices were similar across regions with allelic richness ranging from 2.83 (Greenland) to 3.73 (Anderson River, Canada; Table 1).
Observed heterozygosity ranged from 51.3% (Lena River) to 66.8% (Koyukuk, Alaska) with an overall value of 59.9% (SD 6.7%). All populations and loci were in Hardy-Weinberg and linkage equilibrium.
No significant correlation with allelic richness and longitude was detected within Old or New World regions or within each flyway.

| Low microsatellite differentiation but high mtDNA differentiation
We found significant overall differentiation for mtDNA (F ST = 0.131, Φ ST = 0.295, p < 0.05). However, not all sampling sites were significantly differentiated from each other (Table 2). Notably the few nonsignificant pairwise F ST comparisons were primarily restricted to either between Palearctic populations or within the North American midcontinent. In the case of the Magadan population, sample size was low (n = 5) so nonsignificant pairwise F ST values associated with Magadan may reflect Type II error.
An analysis of variance (AMOVA) revealed that the best partition of genetic variance was when all populations were grouped primarily by major flyway (western Palearctic, eastern Palearctic, Pacific, and midcontinent, Table 3) with Lena River and Yana River representing a potential central Palearctic Flyway (Φ CT = 0.117, p = 0.001). Alternatively, if Lena River and Yana River was placed within the eastern or western Palearctic, 9.3% and 10.8% of the genetic variation was explained with these groupings (Table 3).
Most geographic proximity groupings with the inclusion of Greenland also explained a significant portion of variance (5.5%-11.9%; Table 3 and Supporting Information Appendix S4). The best geographic proximity grouping (11.9% of genetic variation explained; Φ CT = 0.119, p = 0.001) was the same as best flyway grouping with the only difference was Greenland was included within North American midcontinent population. In general, the inclusion of Greenland within the western Palearctic for hypothesized geographic groupings generally lowered the amount of genetic variance explained compared to when Greenland was placed within North America (Supporting Information Appendix S4). Groupings based on nesting habitat did not explain a significant amount of generic variation while groupings based on proposed refugia (Ploeger, 1968) explained 5.2% of the variance (Table 3).
When testing the different subspecies classifications (excluding Greenland and tule goose), the best grouping was proposed by Delacour (1954) with North America populations (currently classified as A. a. frontalis) and eastern Palearctic (with break at Lena River) grouped together and the western Palearctic as a second group; although it was marginally significant (Φ CT = 0.135, p = 0.06, Table 3). If the break between Palearctic regions was placed at the Kolyma River Delta, 9.9% (p = 0.027) of the genetic variation was explained. In general, while flyway, geographic or subspecies groupings typically accounted for approximately F I G U R E 4 Plot of mtDNA control region nucleotide diversity (π) by longitude. Boundaries of Beringia as defined by Hultén (1937) are indicated by dashed line which also represent two wellknown suture zones at Lena River, Russia and Mackenzie River Delta, Canada In contrast to mtDNA sequence data, we found very little structure within the microsatellite dataset. Significant pairwise F ST estimates were restricted to comparisons with Greenland and Cook Inlet populations (Table 2). This general lack of genetic differentiation was also reflected in the STRUCTURE analysis where the most likely number of genetic clusters was K = 3 deemed by both Evanno's method (ΔK = 5.0) and assessing Ln P(K) (K = 3; −13,317.7 vs. K = 1; −13414.6) using the LOCIPRIOR (r < 1). When no location information prior was used, the most likely number of genetic clusters was one. The three genetic clusters corresponded with F ST estimates whereby Greenland and Cook Inlet individuals were grouped in separate clusters with the remaining populations from Asia and North America representing a single cluster or showing admixture ( Figure 5). Further, AMOVA groupings revealed small but significant genetic partitioning (<1% of variance) explained by multiple groupings (e.g., flyway, geographic proximity, Old vs. New World, and subspecies; Table 3 and Supporting Information Appendix S4); while >95% of total genetic variation was partitioned within populations.

| Demographic history and gene flow based on mtDNA
In general, the population size parameter (ϴ) was larger than the ancestral effective population size ( Figure 6). The population size parameter for New World populations was approximately five times larger than Old World (not including Greenland), suggesting a clear population expansion in North America. Although the population size was twice as large in the Old World relative to ancestral population size, suggestive of population expansion, the HPD95 ranges overlapped, potentially indicating population stasis after divergence. This lack of clear demographic expansion was also evident in the neutrality tests, Fu's F s and Tajima's D. Significant negative test statistics, indicating TA B L E 2 Pairwise F ST values for microsatellite data (above diagonal) and mtDNA control region (below diagonal). F ST values in bold text are significant after Benjamini and Yekutieli-modified false discovery rate correction. Values outlined in thick lines include within continent comparisons and thin line and shaded comparisons indicate are within flyway comparisons. Eastern and western Palearctic Flyway designations follow Mooij and Zöckler (2000). Sampling locations are coded as letters (A-W) and refer to codes given in Table 1  TA B L E 3 Hierarchical analysis of molecular variance of haplotypic and allelic frequencies to test hypotheses associated with (a) subspecies classification schemes; (b) flyway designation; (c) geographic proximity; (d) nesting habitat; (e) putative refugia for greater white-fronted goose populations. For complete list of hypothesized groupings, see Supporting Information Appendix S4. Some populations were not included in all groupings as it would be the sole representative for that group. For example, Cook Inlet and Greenland populations were excluded from subspecies grouping analysis as they are the only members of their respective subspecies. Significant fixation indices (p < 0.05) are indicated in bold. Please refer to Figure 1 and Table 1

| Band recovery data
Band recovery data showed that greater white-fronted geese had a high degree of fidelity to the flyway and the continent where they were initially captured (Supporting Information Appendix S2). We observed only four instances of intercontinental dispersal and less than 1% of flyway switching within North America or Eurasian flyways.

| D ISCUSS I ON
We found significant differentiation based on mtDNA control region sequence data, and more limited differentiation based on a Population codes given in parentheses are in reference to locality names and locations indicated in Figure 1 and Table 1 Cook Inlet (tule goose) being the most distinctive. The high mtDNA structure, and general lack of nuclear DNA structure is a common finding in waterfowl; however, this is one of the first studies in waterfowl to our knowledge to show that proposed flyway delineation corresponds to a significant partitioning in mtDNA genetic variance (Kraus et al., 2014;Liu et al., 2012;Peters & Omland, 2007).
Unlike for brant, the other circumpolar goose species, and some other Arctic bird species (see Henningson & Alerstam, 2008), migration routes used by greater white-fronted geese are maintained within the same continent, which is reflected by the lack of intercontinental band recoveries and low mtDNA haplotype sharing across continents. Within waterfowl, geese and swans are unique in their long-lasting parental associations, and some greater whitefronted geese have familial associations lasting for 8 years or more (Ely, 1993;Warren, Fox, Walsh, & O'Sullivan, 1993). Thus cultural transmission of migratory tendencies likely plays a role in the maintenance and integrity of traditional migratory routes as well as the distribution of genetic variation.

| Maintenance of flyway structure
Our analysis based on mtDNA showed significant structure which can at least be partially be explained by flyway designation (Table 3).
The general lack of haplotype sharing across major flyway desig- American waterfowl in eastern Eurasia, especially brant, lesser snow geese, Chen caerulescens, and northern pintail, Anas acuta (Flint et al., 2009;Henny, 1973). This suggests that broad-scale dispersal may be relatively lower in greater white-fronted geese compared to other  CAFF, 2018). Incomplete lineage sorting, therefore, may also play a role in the observed pattern, and hypothesis is supported by the ubiquity and high frequency of the most common allele(s) for each microsatellite locus across regions. Given more recent coalescent times associated with microsatellite loci and the potential connectivity via male dispersal, we suggest that a SNP-based approach may be needed to achieve the necessary statistical power to investigate within flyway structure (Jonker et al., 2012); however, is likely a combination of both gene flow (e.g., within flyway) and incomplete lineage sorting (e.g., across continents) is influencing the pattern of genetic diversity of greater white-fronted geese.
Greater white-fronted geese possess many behavioral and ecological characteristics that might restrict genetic interchange among populations, not only at the broad flyway (macrogeographic) scale but also on a more local (microgeographic) scale (see Ely et al., 2017). As seen in Pacific Flyway white-fronted geese and other goose species, the timing of pairing can have pronounced implications on genetic structuring in geese (Ely & Scribner, 1994;Ely et al., 2017).  (Ely & Takekawa, 1996;Ely et al., 2013). This temporal and spatial segregation is more pronounced in the Pacific Flyway with segregation of site use observed at the population level (Ely et al., 2017), while the midcontinent population tends to be more temporally segregated by general breeding areas (Ely et al., 2013). In the midcontinent, Ely et al. (2013) found that taiga-nesting birds (including interior Alaska and Old Crow Flats, Canada)

| Demographic history
The greater white-fronted goose is the only representative of the "gray" geese in North America and is thought to have a Eurasian origin (Ottenburghs et al., 2016). Our isolation-with-migration analysis suggests that North American and Eurasian populations began diverging at least 50,000 years ago, with post-divergence gene flow biased toward movement into North America. Directionality in gene F I G U R E 6 Estimates of the population size (Ɵ; top) and number of migrants per generation (population migration rate; bottom) for Old World (Eurasia) and New World (North America) greater white-fronted goose populations with lower and upper bound of the estimated 95% higher posterior interval (HPD) indicated. Please see Figure 2 and Table 1  flow and pronounced population expansion in North America further suggests an Old World origin for this species complex as observed in many other species where Beringia was used as an entrance into North America (e.g., Waltari, Hoberg, Lessa, & Cook, 2007). In North America, at the start of the interglacial period approximately 14,000 years ago, an ice-free corridor was established between the Cordilleran and Laurentide ice sheets (Adams, 1997;Pielou, 1991 although this remains to be tested. Ploeger (1968)  The Lena River is generally accepted to be the westernmost boundary of Beringia and is a well-known suture zone (Hewitt, 2004 The same genetic pattern is seen in North America at the Mackenzie River (the easternmost boundary of Beringia; Hultén, 1937) and Anderson River populations. These two locales also lie near a known suture zone (Hewitt, 2004)

| Taxonomic implications
The subspecies concept and its usefulness has been a subject of controversy in ornithology for decades (Mayr, 1942;Winker, 2010; see review in Ornithological Monographs No. 67). The taxonomy of the greater white-fronted goose is a prime example of this, as its taxonomic history has been highly debated and confusing, with respect to the number of subspecies and their nomenclature (see Banks, 2011). As subspecies rank can often be used to inform and be the basis of conservation policies (Winker, 2010), understanding the relationships between the distribution of genetic diversity and subspecies designations can be crucial. As mentioned by Zink and Dittman (1992) and reiterated by Cicero and Johnson (2006), analyses of geographic variation and potential subsequent delineation into subspecies must be done with known breeding origins. Much of the confusion surrounding the subspecific status of greater whitefronted geese likely centers from the use of wintering specimens for taxonomic determinations (Ely & Dzubin, 1994) and mean morphological differences in subspecies descriptions, where more rigorous testing is needed (Cicero & Johnson, 2006;Patten & Unitt, 2002).
The classification of Cook Inlet (tule; A. a. elagsi) and Greenland (A. a. flavirostris) populations as valid subspecies has been little challenged, given their large body size and darker coloration, with slight overlap with other populations (Ely et al., 2005), and especially their isolated breeding areas. These multiple lines of evidence suggest these two populations are on independent evolutionary trajectories, and therefore, it is not surprising that only geese from these two locales showed significant genetic differentiation at both marker types. Both subspecies have restricted ranges throughout their life cycle and reduced population sizes, relative to the more broadly distributed proposed subspecies with which they overlap in distribution to some degree during winter. As such, the genetic data support morphological, ecological and behavioral evidence that provided the foundation of the subspecies hypotheses associated with these two discrete breeding populations. We note that the tule goose is considered to be at risk by the International Waterfowl Research Bureau (Alaska Department of Fish and Game, 2015;Green, 1996), and the Greenland greater white-fronted goose is one of the first subspecies to be red listed in the United Kingdom by the Birds of Conservation Concern assessment (Eaton et al., 2009 (Figure 2; letters Q-T, W) show a close affiliation to North American midcontinent populations (Figure 2; letters D-O), to which they show morphologically similarities (Delacour, 1954). However, the amount of mtDNA genetic variation explained when eastern Palearctic are grouped by themselves (flyway AMOVA grouping) explained approximately the same amount of mtDNA genetic variation (11% vs. 13%) as that under the subspecies model (subspecies AMOVA grouping, see Table 3). In addition, we found no genetic support for the proposed subspecies (A. a. sponsa; Banks, 2011; but see Orthmeyer et al., 1995), based on their smaller average body size, encompassing the Yukon-Kuskokwim Delta and Bristol Bay regions. Due to the lack of clear distinction in subspecies attributions, based on mtDNA sequence data and the presence of multiple body sizes along with a clear phenotypic body size cline in Eurasia, a genomic approach may be needed to determine if the different phenotypes represent genetically discrete populations, and to test taxonomic designations used to make management prescriptions.

ACK N OWLED G M ENTS
We are thankful for the tremendous support of colleagues from around the world who collected samples for us at various locations.
USGS laboratory assistance was provided by J. Gust, A. Palmer, and G. K. Sage. We appreciate the very helpful reviews of the manuscript by J. Pearce, S. Sonsthagen, and two anonymous reviewers and to K.

CO N FLI C T O F I NTE R E S T
None declared.

AUTH O R CO NTR I B UTI O N S
REW, CRE, and SLT were involved in study design. REW collected and analyzed genetic data, and CRE collected data and analyzed banding data. REW lead the writing. All authors contributed to writing and interpretation of results. All authors approve the final version of the manuscript.

DATA ACCE SS I B I LIT Y
Microsatellite genotype data and sample information are available