As the Goose Flies: Migration Routes and Timing Influence Patterns of Genetic Diversity in a Circumpolar Migratory Herbivore

Abstract

Migration schedules and the timing of other annual events (e.g., pair formation and molt) can affect the distribution of genetic diversity as much as where these events occur. The greater white-fronted goose (Anser albifrons) is a circumpolar goose species, exhibiting temporal and spatial variation of events among populations during the annual cycle. Previous range-wide genetic assessments of the nuclear genome based on eight microsatellite loci suggest a single, largely panmictic population despite up to five subspecies currently recognized based on phenotypic differences. We used double digest restriction-site associated DNA (ddRAD-seq) and mitochondrial DNA (mtDNA) sequence data to re-evaluate estimates of spatial genomic structure and to characterize how past and present processes have shaped the patterns of genetic diversity and connectivity across the Arctic and subarctic. We uncovered previously undetected inter-population differentiation with genetic clusters corresponding to sampling locales associated with current management groups. We further observed subtle genetic clustering within each management unit that can be at least partially explained by the timing and directionality of migration events along with other behaviors during the annual cycle. The Tule Goose (A. a. elgasi) and Greenland subspecies (A. a. flavirostris) showed the highest level of divergence among all sampling locales investigated. The recovery of previously undetected broad and fine-scale spatial structure suggests that the strong cultural transmission of migratory behavior restricts gene flow across portions of the species’ range. Our data further highlight the importance of re-evaluating previous assessments conducted based on a small number of highly variable genetic markers in phenotypically diverse species.

1. Introduction

Species typically exist as groups of populations or metapopulations, in which the degree of inter-population genetic connectivity or divergence varies [1] based on topographic features (e.g., mountain ranges or fragmented landscapes; [2]), species/community interactions [3], and behaviors (e.g., mobility and social organization; [4]). For migratory species, metapopulation dynamics are also influenced by site fidelity (in one or both sexes), the seasonal timing and location of pair formation [5], the latter is affected by the seasonal timing and direction of migration [6,7]. Thus, migration may reinforce population divergence when individuals breeding in adjacent regions use different migratory routes to occupy disparate non-breeding locales. So-called “migratory divides” can give rise to genetically differentiated populations, even where breeding environments appear homogeneous and/or lack obvious physical barriers [8,9,10]. Moreover, how migration patterns influence population genetic structure and potential dispersal and gene flow among breeding areas and management units is particularly important for modeling population persistence in the face of ongoing environmental changes [11,12,13,14]. Therefore, understanding how migration patterns have shaped the distribution of genetic diversity across the landscape and among populations is fundamental for deciphering the evolutionary history of a species, as well as for developing effective management and conservation strategies [15].

For migratory species, management practices are typically developed and implemented across large geographic scales (e.g., flyways), each of which may harbor multiple genetically divergent populations that could warrant separate management strategies. Thus, finer-scale assessment of population structure is essential for defining appropriate management units. Increasingly, genetic data are used as a critical complement to movement data for the delineation of populations. Movement studies that rely on recapture/resighting data alone may be insufficient or biased [16]. For example, banding of migratory birds often occurs along migration routes or on wintering grounds, or focuses on adults at breeding sites; in all these cases, the natal origin of sampled individuals is unknown. Furthermore, it may be difficult to collect resighting data at the appropriate time scale for long-lived species with delayed reproduction. Studies that use telemetry or geolocation can generate more detailed data on the movements of individuals, but periods of observation are often of short duration. Moreover, the attachment of transmitters or other devices may affect movement in waterfowl [17], limiting the efficacy of this approach. Finally, movement data typically provide uncertain information on gene flow, as there may be no opportunity to determine if dispersal is followed by reproduction. Thus, relying on movement data alone constrains our ability to link contemporary movement patterns with effective dispersal across broad temporal and spatial scales [18,19].

In contrast to methods based on direct observation, genetic data collected at appropriate temporal and geographic scales [20] can reveal both historical and contemporary patterns of population structure [21], estimate longer-term genetic connectivity beyond the dispersal capabilities of a single individual [22,23], and provide evidence of connectivity between the breeding and wintering sites of migratory species [24,25]. Population genomic data, therefore, provide an opportunity to address the shortcomings of individual movement data and are increasingly incorporated into management and conservation decisions for harvested species [26,27].

The greater white-fronted goose (Anser albifrons) is a long-lived, long-distance migratory herbivore with a circumpolar distribution [28]. Nesting in Arctic and subarctic habitats and wintering in temperate zones, the species exhibits complex patterns of migration [29,30,31,32,33], genetic and morphological differentiation [34,35], and variation in ecological (e.g., foraging) strategies (e.g., [36,37]). Greater white-fronted geese exhibit life history traits such as extended parent-offspring associations [38,39,40], strong long-term pair bonds (1.6% lifetime maximum divorce rate, [41]), solitary breeding [28], spring or summer pairing [5], and social learning of movement behavior ([42]; see [43] for summary of migratory culture), all of which suggest the potential for genetic differentiation at both regional and finer spatial scales. Variation in migratory strategies and differences in phenology of geographic areas used across the annual cycle may further promote the formation of genomic structure between different nesting areas, as migratory routes rather than potential mixing in wintering areas have been shown to explain population connectivity in other migratory birds [44].

Local demographic trends vary widely due to heterogeneity in stressors faced by greater white-fronted geese across their range, including habitat loss/degradation (e.g., on wintering grounds in China, [45]), climate change stimulating earlier departures from winter areas [46], warmer sea surface temperatures causing increased precipitation and spring snow cover on breeding areas and in turn reduced breeding success in west Greenland [47], and illegal hunting, resulting in overharvest relative to prevailing demographic parameters [48]. Nevertheless, most North American greater white-fronted goose populations have experienced substantial growth over the last 30 years [49]. Populations elsewhere, including the small, isolated populations in Cook Inlet, Alaska (~15,000 A. a. elgasi; [50]) and Greenland (<20,000 A. a. flavirostris; [48,51,52]), as well as geese nesting in the Kava River Valley west of Magadan in Russia (<1000 individuals; [53]), have remained stable or declined. Although the global population size exceeds three million [54], the maintenance of geographic variation in phenotypes [55] and genotypes [35] necessitates regional-level management, the effective implementation of which requires detailed information on population structure relative to anthropogenic and other environmental threats.

Previous genetic studies have documented significant differentiation in mitochondrial DNA (mtDNA) for both regional and local scales in greater white-fronted geese [34,35], but these studies also reported little or no differentiation at nuclear microsatellite loci, even when comparing phenotypically differentiated subspecies and populations. Limited genetic structure at nuclear loci might indicate that nesting areas are interconnected either directly or indirectly by male-mediated dispersal [56], as observed in other waterfowl species based on telemetry [57] or banding data. For some species, smaller “genetic” data sets (e.g., few to dozens of loci) and larger “genomic” data sets (e.g., hundreds or thousands of loci) have produced similar results (e.g., [58,59]. In other instances, genomic data have provided greater resolution, revealing cryptic biodiversity, the description of which can inform more effective conservation strategies [60,61]. Notably, genomic data revealed clear evidence of regional population structure in another wide-ranging circumpolar species, white-chinned petrel Procellaria aequinoctialis, that exhibits delayed reproduction and was previously thought to be panmictic [62]. Further, although gene flow can be estimated from standard measures of genetic differentiation (e.g., F_ST), the underlying assumption that populations are at genetic equilibrium is often not met, resulting in poor estimates of the current level of connectivity (effective dispersal) among populations [1]. Alternatively, genomic data sets combined with other analytical (individual-based) methods offer the power to detect more recently established patterns of spatial genetic structure despite the inherent time lag between a cessation in gene flow and subsequent increase in population divergence (F_ST; [63]). Consequently, genomic data have the potential to provide important insights, particularly for species such as greater white-fronted geese that exhibit within-species diversification in morphological and behavioral traits and for which smaller genetic datasets have failed to uncover differentiation.

Towards this end, we generated double digest restriction-site associated DNA (ddRAD-seq) and with previously published mitochondrial sequence data [35] used multiple analytical approaches to: (1) re-assess the geographic distribution of nuclear genomic diversity in greater white-fronted geese, and (2) investigate processes shaping the distribution of genomic variation across the Arctic and subarctic to achieve a deeper understanding of genetic connectivity among populations.

2. Materials and Methods

2.1. Sample Areas and Management Units

Genomic data were obtained from 234 greater white-fronted goose blood samples collected from 20 sites representing seven management units across their circumpolar distribution (Figure 1,

Table 1. General information on sampling areas and sample sizes (N) for the greater white-fronted geese included in this study including primary nesting (including habitat type: taiga † or tundra *) and wintering sites, management units and notes on conservation concerns. Numbers (#) in first column correspond to map locations in Figure 1.

Figure 1 Map #	Continent	Flyway	Region	Nesting Area † Taiga, * Tundra	Wintering Area	Management Unit Name and Population Size (year)source	Conservation Concerns	Subspecies	N
1	North America	Pacific	-	USA (Alaska): Cook Inlet †	USA: California	“Tule” ~15,000 (2019) 1	Small stable population, poor monitoring, hunting	elgasi	25
2	North America	Pacific	-	USA (Alaska): Bristol Bay a,†	Western Mexico	“Pacific” 730,000 (2018) 2	Increasing, low concern, over- abundance increasing causing conflict with “Tule”	frontalis/sponsabc	15
3	North America	Pacific	-	USA (Alaska): Yukon-Kuskokwim Delta *	USA: California			frontalis/sponsabc	21
4	North America	Central	Interior	USA (Alaska): Koyukuk †	USA: southern interior and Mexico	“Midcontinent” 2,000,000–3,000,000 (2016) 3	Increasing, low concern	frontalis	7
5	North America	Central	Interior	USA (Alaska): Selawik †	USA: southern interior and Mexico			frontalis	6
6	North America	Central	Interior	USA (Alaska): Kanuti †	USA: southern interior and Mexico			frontalis	7
7	North America	Central	Interior	Canada (Yukon): Old Crow †	USA: southern interior and Mexico			frontalis	20
8	North America	Central	Western Arctic	USA (Alaska): North Slope *	USA: southern interior and Mexico			frontalis	20
9	North America	Central	Western Arctic	Canada (Northwest Territories): MacKenzie River Delta *	USA: southern interior and Mexico			frontalis	10
10	North America	Central	Western Arctic	Canada (Northwest Territories): Anderson River *	USA: southern interior and Mexico			frontalis	10
11	North America	Central/Mississippi	Eastern Arctic	Canada (Nunavut): Kiillinnguyaq *	USA: southern interior and Mexico			frontalis	5
12	North America	Central/Mississippi	Eastern Arctic	Canada (Nunavut): Victoria Island *	USA: southern interior and Mexico			frontalis	5
13	North America	Central/Mississippi	Eastern Arctic	Canada (Nunavut): Queen Maud Gulf *	USA: southern interior and Mexico			frontalis	5
14	North America	Central/Mississippi	Eastern Arctic	Canada (Nunavut): Rasmussen Basin *	USA: southern interior and Mexico			frontalis	5
15	Asia	West Pacific Palearctic	-	Russia: Anadyr Lowlands *	Korea and Japan	“West Pacific” 402,000–424,000 (2020) 4	Increasing, low concern	albifrons	10
16	Asia	West Pacific Palearctic	-	Russia: Kolyma River Delta *	Korea, Japan, China			albifrons	10
17	Asia	West Pacific Palearctic	-	Russia: Magadan †	?	?	Data deficient	albifrons	4
18	Asia	East Asian Continental Palearctic	-	Russia: Lena River *	China	“East Asian Continental” 48,000 (2020) 4	Decreasing thought due to wintering habitat loss	albifrons	10
19	Asia/ Europe	Western Palearctic	-	Russia: Kanin to Taimyr Peninsula *	Europe and SW Asia	“European” 1,400,000 5	Over-abundance increasing causing conflict	albifrons	20
20	North America	Western Palearctic	Greenland	Greenland *	Ireland/United Kingdom	“Greenland” 20,200 6	Depressed reproductive success (climate related) failing to balance natural mortality	flavirostris	19

^a Nesting habitat in Bristol Bay is intermediate between wet tundra habitat and boreal forest. Typically consists of shrub covered tundra meadows with interspersed stunted spruce and birch trees. ^b Geese nesting areas in southwest Alaska have been proposed as a separate subspecies, A. a. sponsa, due to smaller body size [65]. ^c Subspecies frontalis is sometimes referred to as gambelli. ¹ [50], ² [66], ³ [67], ⁴ [45], ⁵ [68], ⁶ [69].

). Whenever possible, adults with reproductive status confirmed as “breeding” were preferentially included. Greenland samples were collected on the wintering grounds; many of these individuals were resighted in subsequent years on wintering and staging sites in Iceland used only by the Greenland population [64]. Since we did not observe these individuals on the breeding grounds in Greenland, we could not confirm that all individuals attempted to breed during their lifetime. Juveniles that were resighted as reproductive-age adults (including in some cases with their own offspring), and therefore with a known breeding location, were included to bolster the sample sizes for some populations.

Nearctic populations of greater white-fronted geese migrate along the Pacific, Central, and Mississippi Flyways, which are currently managed as three discrete populations in their respective geographical flyway (Figure 1). However, geese in the Pacific Flyway nest in three distinct breeding areas and represent two additional management units that correspond to different subspecies: (1) Cook Inlet (Tule goose, subspecies elgasi), and (2) “Pacific” populations nesting in the Yukon-Kuskokwim Delta and Bristol Bay region (subspecies frontalis/gambelli or sponsa; see [65]). The Central and Mississippi Flyways (also representing subspecies frontalis/gambelli) are comprised of both tundra and interior taiga nesting populations and are typically managed as a single unit referred to as the “Midcontinent” population [70]. Although there is overlap in wintering site use, geese from the western Arctic and interior Alaska are more likely to use the Central Flyway, whereas eastern Arctic tundra nesting birds more often use the Mississippi Flyway [30]. All management units in the Nearctic have population sizes that are considered either stable or increasing, with low conservation concern, except for the Tule Goose, which is considered at risk due to its low census size and a potential increase in competition due to increased spatial overlap with other Pacific geese during winter [50].

Greater white-fronted geese nesting in the Palearctic (subspecies albifrons) can generally be grouped into three main management units: (1) West Pacific (wintering in Japan and Korea), (2) East Asian Continental (wintering in China), and (3) European (wintering from the Caspian Sea, Iraq and Turkey up through Western Europe to the North Sea coasts) (Figure 1). Geese in the West Pacific and East Asian Continental management units display different habitat preferences in their respective wintering areas, which is thought to contribute to differing population demographic trajectories. West Pacific birds wintering in Korea and Japan prefer coastal areas [45], have a strong preference for feeding in agricultural land, and are less reliant on wetlands, all of which is thought to be correlated with an increasing population size, as observed in most North American goose species [71,72]. In contrast, the habitat use of East Asian Continental birds wintering in China is largely restricted to natural wetlands, a pattern that may be enforced by greater human activity on adjacent farmland, including illegal hunting and gleaning of spilled grain by domestic waterfowl [73]. In addition, wetland areas in China are decreasing and becoming more degraded, which may “trap” wintering geese in natural wetlands of decreasing size and quality, despite the availability of nearby farmland [74].

In Greenland, greater white-fronted geese nest in the thin strip of land on the western side of the island that becomes ice-free in summer (100–160 km wide between 66° and 72°N) and stage in Iceland during spring and autumn en route to and from wintering areas in Ireland, Scotland and Wales [52]. This small population represents a morphologically distinct subspecies (flavirostris) that rarely mixes with other populations in wintering areas; it is managed separately from other populations in both North America and Europe. The Greenland population has recently declined because of poor breeding success (likely due to a series of years with increased spring snow cover), which has failed to balance annual mortality, despite the closing of hunting throughout the range [75].

2.2. Library Preparation and Bioinformatics

Genomic DNA was extracted using a DNeasy Blood and Tissue following the manufacturer’s protocols (Qiagen, Valencia, CA, USA). Extractions were quantified using a Broad Range Quant-iT dsDNA Assay Kit (Thermo Fisher Scientific, Inc., Waltham, MA, USA). Double digest restriction-site associated DNA (ddRAD) library preparation and bioinformatic pipelines followed DaCosta and Sorenson ([76]; Python scripts available at http://github.com/BU-RAD-seq/ddRAD-seq-Pipeline, accessed on 18 July 2018). Single-end sequencing (150 bp) was completed on an Illumina HiSeq 2500 at the Tufts University Core Genomics Facility. Libraries were indexed with dual 6 base pair (bp) indices and demultiplexed using bcl2fastq-1.8.4 software (Illumina Inc., San Diego, CA, USA). After clustering of similar sequences into putative loci, genomic positions relative to the Gallus gallus reference genome (GenBank assembly GCA_000002315.2) were determined using blastn v.2 [77], allowing clusters representing the same locus to be combined, thereby improving the detection of polymorphic insertions and deletions (indels) [76]. Generally, individual genotypes at each locus were scored as heterozygous when two distinct haplotypes (i.e., alleles) each accounted for more than 29% of sequence reads for a given sample, and as homozygous when >93% of reads were consistent with a single haplotype. Putative heterozygotes with only 20–29% of reads representing a second allele were also allowed if that allele was present in other individuals. If not, the genotype was flagged as ambiguous. When only 7–20% of reads represented a second allele or there was evidence of more than two alleles in a given sample, genotypes were also flagged as ambiguous. Loci with a median per sample sequencing depth ≥10, <10% missing genotypes, and <10% flagged genotypes were retained for downstream analyses.

2.3. Population Divergence and Nucleotide Diversity

We calculated per population nucleotide diversity and composite pairwise estimates of relative population divergence (Φ_ST) using a Python script (out2phistA.py; available at http://github.com/BU-RAD-seq/Out-Conversions, accessed on 18 July 2018). These comparisons considered seven populations based on defined management units as outlined in Table 1, including: (1) the primary management units in North America (Tule Goose, non-Tule Pacific, Midcontinent), (2) Greenland, and (3) Palearctic flyways (West Pacific, East Asian Continental, and Western). For each pairwise comparison, we tested whether the Φ_ST value was significantly different from zero by comparing it to a null distribution generated by randomly reassigning individuals to populations (n = 101 replicates for each pair of population). We also calculated matrilineal genetic divergence (Φ_ST) among the same groups and mitochondrial nucleotide diversity using previously published mitochondrial (mtDNA) control region sequence data (376 base pairs; [35]) in ARLEQUIN v.3.5.1.2 [78].

2.4. Population Structure

Population structure was analyzed using four methods based on the ddRAD-seq data using individual-based approaches: (1) principal components analysis (PCA) to characterize overall patterns of genetic variation; (2) maximum likelihood clustering analysis using ADMIXTURE [79,80] to test for the presence of multiple genetic clusters; (3) fineRADstructure [81] to assess genetic relationships based on recent shared co-ancestry; and (4) Estimated Effective Migration Surface (EEMS, [82]) to identify geographic regions that deviate from a null model of isolation-by-distance (IBD).

First, we visualized genetic structure with PCA following the approach of Novembre and Stephens [83]. We coded individual genotypes as 0, 1, or 2 (i.e., the number of copies of the alternate allele for each biallelic polymorphism) and analyzed the data in R v. 4.0.2 [84]. Missing or ambiguous alleles (comprising ~2% of the data matrix) were assigned a score equal to the alternate allele frequency in the overall data set for the respective single nucleotide polymorphism (SNP) or indel. We excluded rare SNPs with a global frequency <1% (i.e., 4 or fewer copies among the 217 individuals analyzed = 434 alleles). Finally, because the genetic similarity of close familial relatives can generate a stronger signal in PCA analyses than population-level similarities and differences, we tested for the presence of close relatives within each of the sampled populations using COLONY v.2.0.6.8 [85], scoring each unique haplotype at each locus as an allele. Given the large sample of loci in the ddRAD-seq data, COLONY runs using a range of priors for allelic dropout and genotyping errors produced identical results. For the PCA analysis, we retained one individual from each set of close familial relatives, resulting in the exclusion of 17 individuals (see Supplemental Table S1 for summary of COLONY results). For other analyses, the exclusion of close relatives did not have an appreciable effect on results, so all samples were included in the remaining analyses.

Second, we used SNP data to obtain maximum likelihood estimates of population assignments for each individual with ADMIXTURE v.1.3. Biallelic SNPs were first formatted using plink [86]; SNPs with a rare allele observed only in a single individual were excluded from the analysis. ADMIXTURE analyses were run without a priori assignment of individuals to populations with a 10-fold cross validation (CV). A quasi-Newton algorithm was employed to accelerate convergence [87] using 100 iterations for each number of populations (K) from 1 to 15; the optimum K was estimated based on the average of CV-errors per K. Results for other values of K were visualized to explore alternative summaries of population structure consistent with the biology of this species [88]. Final admixture proportions for each value of K and per sample assignment probabilities (Q estimates; the log likelihood of group assignment) were based on CLUMPP v.1.1 [89] summaries of all 100 replicates per K value. The R program PopHelper [90] was used to convert ADMIXTURE outputs into CLUMPP input files for each K value.

Third, we used the fineRADstructure program to infer population structure based on recent shared ancestry. By emphasizing the most recent coalescent events, the co-ancestry coefficient calculated by fineRADstructure is particularly sensitive to patterns of recent co-ancestry, and thus to subtle patterns of population structure that may have been more recently established. As the analysis is based on the sharing of identical or nearest-neighbor haplotypes at each locus, all SNPs within each locus are concatenated to define alleles. Missingness (the proportion of missing alleles per individual, which included both missing and ambiguous genotypes) ranged from 0.007 to 0.168 with a median value of 0.021 and the median values per population ranged from 0.013 to 0.047. Samples were assigned to populations using 5,000,000 iterations sampled every 1000 steps with a burn-in of 500,000. We used 1,000,000 iterations of the tree-building algorithm to define and assess the relationships among clusters. Finally, the output was visualized using the R scripts, fineradstructureplot.r and finestructurelibrary.r (available at http://github.com/millanek/fineRADstructure, accessed on 10 November 2019).

Fourthly, we implemented the spatial method EEMS to estimate effective gene flow (m) and genetic diversity (q) relative to geographic distance in order to reveal patterns across the landscape. EEMS uses a stepping stone model to identify geographic areas where genetic dissimilarity decays more quickly or more slowly than expected under a model of isolation by distance (IBD) based on the estimated migration rate (m). A gene flow surface that correlates genetic variation with geographic distance is interpolated to visualize potential barriers or corridors to movement. In addition, an effective diversity parameter (q) estimates the expected within-deme coalescent time and is proportional to average heterozygosity. A coordinate file was constructed using Google Maps API v3 (http://www.birdtheme.org/useful/v3tool.html, accessed on 19 November 2019) to define an outer boundary that included the entire Arctic nesting distribution of greater white-fronted geese. As the distribution encompasses the international date line, we added 360 to negative longitude coordinates in order to have a closed polygon to fulfill the requirements of the program and allow for movement across the Bering Strait. We used the same SNP data set as in the ADMIXTURE analysis. Parameters were adjusted during preliminary runs until the accepted proportion of steps in the Markov chain Monte Carlo (MCMC) chain was at least 10%. We ran three independent analyses using 5,000,000 burn-in steps followed by 25,000,000 MCMC iterations sampled every 2000 steps for each of three deme sizes (100, 250, 500). We checked for convergence and visualized effective migration and genetic diversity surfaces using the R package rEEMSplots [82].

3. Results

3.1. Bioinformatics

We obtained a median of 682,280 raw sequence reads per sampled individual (range: 344,091 to 1,643,653). A total of 4155 clusters (i.e., putative single-copy loci) met the median depth per sample threshold of ≥10, and among these, 3888 loci passed thresholds for genotype and alignment quality, yielding 36,522 SNPs or unique indels (only 5 loci were invariant). Median sequencing depth of 79 reads per locus per individual for this set of loci resulted in generally robust inference of genotypes, with <1% missing data and ~1% ambiguous genotypes in the overall data matrix.

We detected pairs or trios of close familial relatives, including full- and half-sibs as well as parent-offspring relationships, in several of the sampled populations (Supplemental Table S1). The Greenland samples included a trio and pair of full sibs as well as three parent-offspring pairs, all confirmed by post-marking behavioral observations. The Cook Inlet sample set also included a trio of full sibs, along with their inferred father; one of the siblings was captured in Oregon during spring migration the year after the other samples were collected and was tracked to nesting locations in Cook Inlet in two subsequent years.

3.2. Population Divergence and Molecular Diversity

Based on ddRAD-seq loci, we observed a low but statistically significant level of genetic differentiation across the management units (Φ_ST = 0.021 across all populations and loci; with pairwise comparisons ranging from 0.005 to 0.078; Supplemental Figure S1). We recovered relatively higher differentiation in pairwise comparisons, including either the Cook Inlet (Tule; all Φ_ST ≥ 0.025) or Greenland nesting areas (all Φ_ST ≥ 0.046). Lower levels of genetic differentiation were found within North America excluding Cook Inlet (Φ_ST = 0.005) and the Palearctic (Φ_ST = 0.007 across all Palearctic populations; pairwise Φ_ST = 0.007–0.009; Supplemental Figure S1). Similarly, we observed a low level of intercontinental differentiation, and particularly between Eurasian and North American continental populations (

Table 2. Pairwise Φ_ST values calculated from 3888 ddRAD-seq loci (below diagonal) with highest value for an individual locus in parentheses, and Φ_ST values for the mtDNA control region (above diagonal) for greater white-fronted goose management units. Values within boxes indicate estimates within the same continent. Significant values for mtDNA (α = 0.05) and ddRAD-seq loci (α = 0.01) are indicated in bold text. See Supplemental Figure S1 divergence estimates among individual nesting areas.

	Western Palearctic	East Asian Continental Palearctic	West Pacific Palearctic	Greenland	Midcontinent	Pacific	Cook Inlet Basin
Western Palearctic	-	0.041	0.184	0.103	0.076	0.160	0.474
East Asian Continental Palearctic	0.007 (0.251)	-	0.248	0.204	0.084	0.249	0.599
West Pacific Palearctic	0.007 (0.351)	0.009 (0.515)	-	0.237	0.111	0.181	0.159
Greenland	0.053 (0.610)	0.057 (0.593)	0.057 (0.711)	-	0.060	0.135	0.524 a
Midcontinent	0.007 (0.381)	0.009 (0.689)	0.008 (0.214)	0.046 (0.663)	-	0.134	0.248
Pacific	0.010 (0.485)	0.011 (0.476)	0.010 (0.535)	0.052 (0.574)	0.005 (0.153)	-	0.303
Cook Inlet Basin (Tule goose)	0.036 (0.487)	0.038 (0.523)	0.035 (0.539)	0.078 (0.698)	0.025 (0.432)	0.027 (0.302)	-

^a Value from Wilson et al. [35].

). Likewise, on a locus-by-locus basis, there were few loci showing evidence of elevated divergence (defined here as Φ_ST > 0.1) among management units (Supplemental Figure S2). The Greenland sample had the largest fraction of loci showing relatively greater divergence (11–23% of loci with Φ_ST > 0.1 in comparisons with other populations), whereas the Tule population showed elevated divergence (Φ_ST > 0.1) from non-Greenland samples at 3.7–6.4% of loci. In all other pairwise comparisons, less than 1% of loci had Φ_ST values exceeding 0.1. Nucleotide diversity for ddRAD-seq loci was similar for all management units and nesting areas/regions (range: 0.0058–0.0066;

Table 3. Nucleotide diversity for the mtDNA control region (this study and [35]) and 3888 ddRAD-seq loci for management units and subpopulations within the Midcontinent and Pacific Units. Standard deviations are in parentheses.

Continent	Management Unit (Nesting Region/Area)	mtDNA	ddRAD-Seq
Palearctic
	West Pacific Palearctic	0.0198 (0.0104)	0.0065 (0.0067)
	East Asian Continental Palearctic	0.0047 (0.0034) a	0.0063 (0.0066)
	Western Palearctic	0.0154 (0.0084) a	0.0066 (0.0066)
Greenland
	Greenland	0.0084 (0.0050) a	0.0058 (0.0067)
North America	Midcontinent	0.0188 (0.0098)	0.0066 (0.0066)
	Interior Alaska	0.0206 (0.0107)	0.0064 (0.0066)
	Interior Yukon/Old Crow	0.0186 (0.0102) a	0.0065 (0.0066)
	Western Arctic	0.0156 (0.0084)	0.0065 (0.0067)
	Eastern Arctic	0.0143 (0.0077)	0.0064 (0.0067)
	Pacific	0.0186 (0.0099)	0.0066 (0.0067)
	Bristol Bay	0.0109 (0.0064) a	0.0064 (0.0068)
	Yukon-Kuskokwim Delta	0.0217 (0.0116) a	0.0065 (0.0068)
	Cook Inlet Basin (Tule goose)	0.0094 (0.0054) a	0.0063 (0.0068)

^a Values from Wilson et al. [35].

); the highest percentages of loci lacking any within population variation were observed in Greenland (Greenland management unit; 17.9%), Lena River (East Asian Continental management unit; 13.7%), and Cook Inlet (Tule management unit; 11.8%; Supplemental Figure S2).

The level of mtDNA divergence among populations covered a much broader range (Φ_ST = 0.041–0.599, Table 2), with the highest values observed in comparisons involving Cook Inlet. For both mtDNA and the ddRAD-seq data, low to moderate levels of divergence were observed within management units, with neighboring nesting areas showing the lowest values for both data types (Supplemental Figure S1). In contrast to the ddRAD-seq data, for which nucleotide diversity was similar across management units, we found somewhat lower mtDNA diversity at Lena River and in Greenland, whereas the highest levels of mtDNA diversity were in the Interior Alaska and YK-Delta nesting areas (Table 3).

3.3. Nuclear Population Structure

The first two components of the PCA captured the distinctiveness of the Greenland and Cook Inlet populations from each other and all of the other sampled populations (Figure 2A). Additional axes captured divergence between North America and Eurasia (PC3) and among populations within each region (PC4 and PC5, respectively), but divergence is easier to visualize in a separate PCA excluding the Greenland and Cook Inlet samples (Figure 2B). In the latter analysis, PC1 separates Palearctic and North American populations, whereas PC2 suggests east–west IBD among the North American populations. Similarly, PC3 for this analysis (results not shown) captured a similar east–west gradient in the Palearctic. Generally, individuals clustered with other members of their respective sampling location and/or management unit with one notable exception; one Tule Goose (Tule B_40) was included within the Palearctic cluster, a result also observed in other types of clustering analyses (see below).

As with PCA, ADMIXTURE generally grouped individuals by continent and management unit (Figure 1 and Supplemental Figure S3). The CV-error was lowest for K = 2 (CV-error = 0.2550) which grouped Greenland and Cook Inlet as a single cluster with all remaining individuals in the other grouping. However, exploration of higher K-values that had CV-error similar to K = 2 (within 0.0129) revealed additional biologically and geographically relevant partitioning. At K = 3 (CV-error = 0.2564), individuals sampled in Greenland and Cook Inlet formed nesting-area specific clusters (Supplemental Figure S3), with the exception of one Tule goose sampled on the nesting grounds (the same individual noted above as having a Palearctic genotype noted in the PCA). At K = 4 (CV-error = 0.2613), individuals representing the Palearctic (subspecies albifrons) formed a distinct cluster from all other non-Tule Holarctic samples. At K = 5 (CV-error = 0.2679), evidence of genetic structure across the Pacific and Midcontinent management units emerged, with the majority of individuals from the Interior region (see Table 1) having an admixed signature between the two management units (Figure 1).

Results from fineRADstructure (Figure 3) were broadly consistent with our other analyses, but also highlighted additional patterns of population structure not evident in the PCA or ADMIXTURE results. Individuals generally clustered by nesting areas, which in turn clustered by geographic proximity and migration routes (see map in Figure 1). Pairs of samples with high co-ancestry values (e.g., >30) are close familial relatives confirmed by the COLONY analysis; clusters of close relatives generally fall within a larger cluster of samples from the same population except in the case of three related individuals from the YK-Delta, which form a distinct cluster. Despite being placed outside the main YK-Delta cluster, there is low support (0.18) for the relevant node in the tree and these three YK-Delta individuals share greater co-ancestry with samples from the YK-Delta than with any other populations. Three clusters included samples from multiple areas; a cluster comprising samples from Interior Alaska also included one sample each from the Mackenzie River Delta and Bristol Bay, along with seven samples from Old Crow. Another cluster comprised samples from the Lena River (East Asia Continental Palearctic) along with three samples each from the West Pacific (Kolyma River Delta) and Western Palearctic (Taimyr Peninsula) flyways along with the one Cook Inlet (B_40) individual noted above. Finally, one North Slope Alaska goose clustered with eastern Canadian Arctic birds. Recent co-ancestry was particularly high within the Tule and Greenland samples, respectively, reflecting their somewhat lower level of genetic diversity and presumably relatively greater isolation. More significantly, the fineRADstructure analysis showed that the Greenland population shared relatively greater co-ancestry with samples from the Canadian Eastern Arctic than with samples from other regions, which is also reflected in mtDNA-based Φ_ST calculations (Φ_ST = 0.049 for this comparison versus 0.086–0.257 for other comparisons, see Supplemental Figure S1). Similarly, the Tule samples have elevated co-ancestry with samples from the Pacific management unit and Alaska Interior nesting areas. These results likely reflect some level of historical connectivity between these respective populations.

The EEMS analysis highlighted regions with lower gene flow than expected under a model of IBD (Figure 4). Gene flow surfaces were similar across assumed deme sizes with one exception noted below. Regions indicated as having significantly reduced gene flow (posterior probability > 95%) included: (1) the West Pacific Palearctic (Anadyr Lowlands and Magadan); (2) south-central Alaska extending into interior boreal Alaska, a region that includes both Pacific and Central Flyway populations along with the Tule Goose; (3) Kiillinnguyaq/Victoria Island in eastern Canada; and (4) Greenland. In addition, the Mackenzie River Delta showed evidence of reduced gene flow but only with an assumed deme size of 500; Supplemental Figure S5). These results roughly correspond to the population clusters identified by ADMIXTURE and fineRADstructure. High connectivity (i.e., genetic decay was slower than expected under IBD; blue areas in Figure 4) was characteristic of most of the Arctic extending from the western Palearctic to the Kolyma River Delta in Eurasia (>95% posterior probability in Kolyma River Delta and Lena River), and across most of the North American Arctic (west and east of the Mackenzie River Delta with >95% posterior probability in North Slope Alaska). The EEMS analysis also identified lower levels of genetic diversity (q) for Cook Inlet, Magadan, and Greenland.

4. Discussion

Greater white-fronted geese exhibit more complex patterns of population genetic structure than previously understood. The distribution of genetic diversity across the landscape coincides, at least partially, with differences in the timing and direction of migration combined with the timing of behavioral events such as pair formation. By surveying a large number of nuclear loci, we recovered subtle patterns of genetic structure within and across flyways that were not detected in previous studies that assayed fewer loci (microsatellite loci; [34,35]). While the results presented above are concordant with results based on mtDNA (Supplemental Figure S1; [35]), mtDNA data does not offer the same level of resolution. Analyses leveraging insights from coalescent theory (i.e., fineRADstructure, EEMS) reveal fine-scale genomic structure and provide insight into both historical and contemporary processes that have likely shaped patterns of connectivity and divergence among greater white-fronted goose populations.

Greater white-fronted geese are genetically clustered into five main groups (Figure 1): (1) Greenland, (2) Cook Inlet, Alaska (Tule Goose), and (3) the Palearctic. The remaining North American populations across the Pacific and Central/Midcontinent Flyways show clear evidence of further fine-scale structuring that distinguishes the (4) Pacific and (5) Midcontinent populations in North America (Figure 3). Nonetheless, the overall level of genetic differentiation was generally low across the range of the greater white-fronted goose, with the exception of Cook Inlet and Greenland, which harbor distinct and more highly differentiated populations. The more subtle structure uncovered within flyways in the Nearctic and Palearctic is also biologically relevant. Notably, across- and within-flyway structuring as well as signatures of admixture across portions of the North American and Palearctic ranges, respectively, is consistent with previously described biogeographic provinces defined by historical (e.g., glaciation events) and contemporary (e.g., migration strategy and spatio-temporal habitat use) processes as well as with morphology (see below).

4.1. Historical Biogeography

The population genetic structure of Arctic species is often influenced by vicariance during the Pleistocene, as ice sheets fragmented the landscape and displaced species to lower latitudes and/or high-latitude refugia [91]. Glacial refugia for greater white-fronted geese have been proposed in Alaska and eastern Russia (Beringia), the High Canadian Arctic, Newfoundland (Atlantic Shelf), Svalbard, and Greenland/Iceland, whereas remaining nesting areas were glaciated [92]. Presence of an ice shield along the Brooks Range and nunataks in northern and eastern Beringia contributed to additional genetic structure within species occupying this region during the Last Glacial Maximum [93,94]. Indeed, the heterogeneous landscape within Beringia along with smaller refugia in southeastern Alaska [95] may have promoted genetic differentiation among nesting areas in this region. Geese nesting north of the Brooks Range are genetically differentiated from birds breeding in other parts of Alaska, and morphological variation between geese nesting in western and south-central Alaska further suggests isolation in multiple refugia. Morphological similarity between greater white-fronted geese in western Alaska and the Anadyr Lowlands [55], suggests that the geese in these areas may have expanded from the same refugium in Beringia. Moreover, lower genetic diversity in the Western Palearctic and East Asian Continental Palearctic than in the West Pacific Palearctic, observed in both mtDNA [35] and ddRAD-seq loci, supports the inference of east to west expansion across the Palearctic from a single refugium in Beringia (see also collared lemmings Dicrostonyx torquatus [96]).

The origin of the Greenland subspecies of greater white-fronted goose cannot be ascertained based on fossil evidence [97]. However, subfossil evidence suggests the greater white-fronted goose was the most widespread goose species in low-arctic west Greenland from historical times (at least c. 1200–1400 A.D.) until the last few decades, when populations of other goose species (snow goose Anser caerulescens, Canada goose Branta canadensis) began to expand [98]. Subfossil evidence also shows that, in historic times, a small Branta species, possibly brant (Branta bernicla) or a small-bodied white-cheeked goose (B. hutchinsii), coexisted with greater white-fronted geese [98]. Prolonged isolation of greater white-fronted geese nesting in Greenland is suggested by both their genetic and morphological distinctiveness (i.e., body and bare parts color, head and bill size and proportions, body size and mass; [55]). Ploeger [92] hypothesized that the Greenland geese originated from western European stocks and were isolated in a North Sea refugium during the Last Glacial Maximum. A North American origin was dismissed due to: (1) the lack of regular wintering grounds in North America; (2) the considerable morphological differences between the Greenland subspecies and other North American birds; and (3) the lack of continuous nesting areas across eastern North America and Greenland. Our genetic analyses, however, indicate that the Greenland geese share more recent co-ancestry with North American than Eurasian populations (Figure 3; see also [35]). The disjunct, present-day nesting distribution of greater white-fronted geese (see Figure 1) may be attributable to the more recent retreat of the Canadian Shield, remnants of which remained over portions of eastern Arctic Canada until at least 5.7 cal ka BP, well after western areas of North America were already free of ice [99]. Indeed, greater white-fronted geese are not known to nest in northeastern Canada. One potential scenario that explains the spatial patterns in allelic and haplotypic diversity is that following climate amelioration, greater white-fronted geese occupying High Canadian Arctic refugia colonized habitats made available by retreating ice sheets in Greenland, which commenced deglaciation by 10.1 cal ka BP [99]. As numerous ice divides and domes were still present in the High Canadian Arctic Archipelago, dispersing over the Greenland Icecap into Iceland, which is currently a key stopover site in the present-day migration route, may have been more successful. While hypotheses for the origin of Greenland migratory flyways typically evoke expansion from Europe via Iceland, the establishment of Greenland to Europe migration originating from an ancestral population in North America has also been proposed for the red knot Calidris canutus, which nests in both the High Canadian Arctic Archipelago and Northern Greenland and shows a pattern of mtDNA structure similar to that of greater white-fronted geese [100].

4.2. Effects of Timing and Location of Annual Cycle Events and Location on Genetic Exchange

Movement patterns throughout the annual cycle influence population demographic parameters within and between populations [101]. Although long recognized that a strong spatial barrier between allopatric populations can limit genetic intermixing, temporal isolation among populations can also restrict gene flow, such that “when” individuals or populations occupy particular sites can influence genetic differentiation just as much as “where” individuals are distributed [101,102,103,104]. Staggered use of northern staging areas in autumn and spring by greater white-fronted geese occupying different nesting areas has been documented in North America [29,30], resulting in population-specific chronology of key events such as the timing of spring migration and nesting. Marked variation in when individuals from different nesting areas occupy non-breeding sites likely contributes to the subtle yet detectable genetic structure observed within both the Pacific and Central Flyways in North America as pairing likely occurs on these sites during the spring [105,106]. Further, the greater white-fronted goose is a highly social species that maintains family bonds longer than most other goose species (up to 13 years in the case of the Greenland population; [40]). Social connectivity throughout the annual cycle and across years coupled with cultural transmission of behaviors such as the timing of migration and breeding across generations may reduce the propensity for gene flow, while also reducing the fitness of any dispersing individuals (e.g., barnacle geese Branta leucopsis [107]; greylag geese Anser anser [108]).

The spatial scale of population genetic structure varied among flyways. Within the Pacific Flyway, genetic differences were evident among the three major nesting areas in Alaska: Bristol Bay, Yukon-Kuskokwim Delta, and most distinctively, Cook Inlet (Tule). Spanning a much larger geographic area, birds from nesting areas within the Midcontinental management unit, which use the Central and Mississippi Flyways, formed four genetic clusters (Figure 3): (1) western Arctic—Mackenzie River Delta area and west; (2) eastern Arctic—east of Mackenzie River Delta area; (3) a mixed cluster comprising birds from Interior Alaska and Old Crow Flats; and (4) a smaller cluster comprising birds exclusively from Old Crow Flats. As indicated by the EEMS analysis (Figure 4), gene flow among Alaska nesting areas is limited, resulting in greater differentiation among populations on a relatively small spatial scale, likely reflecting a higher level of effective isolation among nesting areas within this region. Specifically, nesting areas within the Pacific Flyway exhibit both spatial and temporal segregation, with little temporal overlap among geese from different nesting areas across most of the annual cycle [29]. Conversely, birds from nesting areas within the Midcontinent management unit have a higher degree of spatial and temporal overlap in their use of stopover sites during migration [30], which may explain, in part, the lower level of genetic structure found within this unit.

Although populations are largely segregated (spatially and/or temporally) throughout the year, notably on nesting and wintering areas, synchronous occupancy of molting or staging grounds may be important in maintaining genomic connectivity between some flyways [33]. Geese nesting within interior and southern Alaska, respectively, utilize shared molting sites in Holy Cross and Innoko, Alaska [35], where geese that migrate along the Pacific and Central Flyways, respectively, intermix. This may provide opportunities for the formation of either temporary or longer-term social bonds that influence migration behavior and in turn lead to the exchange of individuals between these populations. Indeed, geese nesting in Interior Alaska (Central Flyway) have greater shared co-ancestry with Pacific Flyway birds, including those that nest in similar habitat (i.e., Bristol Bay), than with geese that nest in other North American regions. At Old Crow Flats, molting geese from elsewhere and potentially from Interior Alaska intermix with locally breeding geese [109]; this likely explains why birds sampled at Old Crow Flats fall into two distinct clusters in the fineRADstructure analysis, one that includes geese originating from Interior Alaska and one comprising birds from Old Crow Flats only. While we hypothesize that this result is attributable to the temporary mixing of geese from different nesting areas, the similar cultural behaviors (e.g., timing and location of staging sites along migration routes) and nesting preferences in these areas may help maintain their greater shared co-ancestry. Additionally, different nesting habitat preferences may be a barrier to dispersal, as no non-origin assignments of individuals were observed between tundra (e.g., Yukon-Kuskokwim Delta) and taiga (e.g., Interior Alaska, Bristol Bay) nesting geese that share molting sites. However, it should be noted that one tundra nesting goose sampled on the North Slope, Alaska clustered with the Interior Alaska geese. These two regions do not share primary molting sites but there is some temporal overlap in other areas along migratory routes [30], which may provide opportunities for intermixing. Restricted dispersal between tundra and taiga nesting locations suggests that even in complete sympatry during parts of the annual cycle, conspecific cues such as body size may be used to discriminate among potential mates as either tundra and taiga nesters or alternatively, high level of philopatry to nesting grounds may help to maintain genetic distinctiveness [110].

Less information is available for Palearctic populations but telemetry data suggest that some Kolyma River geese may move between populations and intermix with both Lena River and Anadyr geese [45]. In our analyses, a few birds sampled at both Kolyma and the Taimyr Peninsula cluster with the Lena River samples (Figure 3), consistent with dispersal of Lena River birds in both directions. A similar pattern was observed for mtDNA, with haplotypes shared between Lena River and sample sites to both the west and east, but with no haplotypes shared between Taimyr and Anadyr [35]. Our results for the Palearctic, including the EEMS analysis (Figure 4), suggest few barriers to gene flow between Kolyma and populations to the east suggesting that isolation by distance is sufficient to explain a longitudinal gradient in genomic variation that corresponds to a proposed morphological cline across the Palearctic [55].

5. Future Directions

The warming of the Arctic and habitat changes along migratory routes may disrupt potential isolating mechanisms, including temporal segregation among greater white-fronted goose populations. Indeed, the increase in agriculture and concomitant loss of wetlands in the past 150 years has undoubtedly eroded ecological segregation of greater white-fronted geese in wintering areas. This is most evident in the Pacific Flyway of North America, where loss of wetlands and the introduction of grain crops, particularly rice, has led to increased mixing during winter of birds from allopatric nesting areas (i.e., Cook Inlet and Yukon-Kuskokwim Delta geese; see [34]). An increased level of interaction between Cook Inlet and Yukon-Kuskokwim Delta geese could lead to higher levels of genetic exchange (introgression; [111]), which might bolster genetic diversity within the Tule goose and ultimately its adaptive potential in response to future environmental change. Our analyses, however, suggest that gene flow has not yet increased to a level that would modify the genetic integrity of this distinctive population. Similar changes in the spatial distribution of greater white-fronted geese during winter have occurred elsewhere, as multiple phenotypes (and genotypes) co-occur in wintering areas in the western Palearctic (i.e., geese from both Greenland and Russia), and possibly in the eastern Palearctic on wintering grounds in China. Opportunities to bolster genomic diversity are especially important for species that nest at high latitudes, which are expected to experience the most rapid increases in temperatures and thus may be most vulnerable to impending climate change [112].

Despite the potential threats of climate change, many of the larger continental populations of greater white-fronted geese in both Eurasia and North America have experienced substantial increases in population size over the last few decades. In contrast, smaller and more isolated populations, including the distinctive Tule and Greenland subspecies as well as populations wintering in China and the Caspian Sea, have remained static or declined [49]. A potential mismatch with the environment may be developing in the Greenland population, where geese have advanced their departure from wintering areas by 15 days over 43 years [46], whereas their arrival on the nesting grounds has not changed since the 1880s. Thus, despite being capable of some plasticity in response to environmental change during certain stages of the annual cycle, there may be limitations at other stages, including constraints on changes in movement patterns and physiological processes related to breeding. If so, stronger cultural transmission and extended familial bonds in greater white-fronted geese as compared to other North American goose species [5] and the low percentage of successfully reproducing adults [51] may impair the capability of the Greenland population to adapt to new environmental stressors more than other forms. Other goose species have already exhibited substantial behavioral and distributional changes, including the loss or reduction of migratory behavior, in response to environmental change (e.g., brant [113]; barnacle goose [114]). In the barnacle goose, reduced transmission of behavioral traditions associated with a decrease in the duration of familial bonds may have led to increased exploratory behavior, enabling this species to respond to a changing environment by establishing new traditions via developmental plasticity [114]. Although no intercontinental movement based on banding data has been reported [35], our data set included an individual goose with a Palearctic genotype sampled in the Cook Inlet, Alaska. This is the first evidence of potential movement across the Bering Strait for greater white-fronted geese. Nevertheless, because the Palearctic subspecies albifrons and Nearctic frontalis are morphologically similar [55], it is difficult to determine whether this sample represents an isolated vagrant or is indicative of a more regular low-frequency occurrence. Given its widespread distribution encompassing different breeding and wintering habitat types, along with potentially variable responses to environmental changes across populations (see [115]), the greater white-fronted goose may be a valuable model for comparative studies of species’ responses to rapid climate change in the Holarctic. Finally, this study provides an important baseline for future analyses in exploring how population structure and distinctiveness of isolated populations (e.g., the Tule Goose and Greenland subspecies) is either augmented or degraded as individuals and populations respond through behavioral modifications and shifts in distribution to environmental and landscape alteration over time.