# High dispersal ability versus migratory traditions: Fine-scale population structure and post-glacial colonisation in bar-tailed godwits

Jesse R Conklin, Yvonne I Verkuil, Margaux J M Lefebvre, Phil F Battley, Roeland A Bom, Robert E Gill, Chris J Hassell, Job ten Horn, Daniel R Ruthrauff, T. Lee Tibbitts, et al.

<figure><img src="image_1.png" /></figure>

To cite this version:

Jesse R Conklin, Yvonne I Verkuil, Margaux J M Lefebvre, Phil F Battley, Roeland A Bom, et al. "High dispersal ability versus migratory traditions: Fine-scale population structure and post-glacial colonisation in bar-tailed godwits." Molecular Ecology, 2024, 33 (15), pp.e17452. 10.1111/mec.17452. hal-04779401

HAL Id: hal-04779401

https://hal.science/hal-04779401v1

Submitted on 13 Nov 2024

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_2.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
# High dispersal ability versus migratory traditions: Fine-scale population structure and post-glacial colonisation in bar-tailed godwits

Jesse R. Conklin<fnref n="1" />,<fnref n="2" /> ID

Yvonne I. Verkuil<fnref n="1" />,<fnref n="2" /> ID

Margaux J. M. Lefebvre<fnref n="3" /> ID

Phil F. Battley<fnref n="4" /> ID

Roeland A. Bom<fnref n="2" />,<fnref n="5" /> ID

Robert E. Gill Jr<fnref n="6" /> ID

Chris J. Hassell<fnref n="7" /> ID

Job ten Horn<fnref n="5" /> ID

Daniel R. Ruthrauff<fnref n="6" /> ID

T. Lee Tibbitts<fnref n="6" /> ID

Pavel S. Tomkovich<fnref n="8" /> ID

Nils Warnock<fnref n="9" /> ID

Theunis Piersma<fnref n="1" />,<fnref n="2" />,<fnref n="5" /> ID

Michaël C. Fontaine<fnref n="3" />,<fnref n="10" /> ID

**Correspondence**

Jesse R. Conklin, Conservation
Ecology Group, Groningen Institute for
Evolutionary Life Sciences (GELIFES),
University of Groningen, P.O. Box 11103,
9700 CC Groningen, The Netherlands.

Email: conklin.jesse@gmail.com

**Funding information**

NIOZ; Research Council (TRC) of the
Sultanate of Oman, Grant/Award Number:
ORG/EBR/12/002; Dutch Research
Council (NWO), Grant/Award Number:
824.01.001

**Handling Editor:** Michael M. Hansen

**Abstract**

In migratory animals, high mobility may reduce population structure through increased dispersal and enable adaptive responses to environmental change, whereas rigid migratory routines predict low dispersal, increased structure, and limited flexibility to respond to change. We explore the global population structure and phylo-geographic history of the bar-tailed godwit, *Limosa lapponica*, a migratory shorebird known for making the longest non-stop flights of any landbird. Using nextRAD sequencing of 14,318 single-nucleotide polymorphisms and scenario-testing in an Approximate Bayesian Computation framework, we infer that bar-tailed godwits existed in two main lineages at the last glacial maximum, when much of their present-day breeding range persisted in a vast, unglaciated Siberian-Beringian refugium, followed by admixture of these lineages in the eastern Palearctic. Subsequently, population structure developed at both longitudinal extremes: in the east, a genetic cline exists across latitude in the Alaska breeding range of subspecies *L. i. baueri*; in the west, one lineage diversified into three extant subspecies *L. i. lapponica*, *taymyrensis*, and *yamalensis*, the former two of which migrate through previously glaciated western

This is an open access article under the terms of the [Creative Commons Attribution License](https://creativecommons.org/licenses/by/4.0/), which permits use, distribution and reproduction in any medium, provided the original work is properly cited.

© 2024 The Author(s). *Molecular Ecology* published by John Wiley & Sons Ltd.

<footnote n="1">¹Conservation Ecology Group, Groningen Institute for Evolutionary Life Sciences (GELIFES), University of Groningen, Groningen, The Netherlands</footnote>

<footnote n="2">²BirdEyes, Centre for Global Ecological Change at the Faculties of Science & Engineering and Campus Fryslân, University of Groningen, Leeuwarden, The Netherlands</footnote>

<footnote n="3">³MiVEGEC, CNRS, IRD, University of Montpellier, Montpellier, France</footnote>

<footnote n="4">⁴Zoology and Ecology Group, School of Food Technology and Natural Sciences, Massey University, Palmerston North, New Zealand</footnote>

<footnote n="5">⁵Department of Coastal Systems, NIOZ Royal Netherlands Institute for Sea Research, Texel, The Netherlands</footnote>

<footnote n="6">⁶U.S. Geological Survey, Alaska Science Center, Anchorage, Alaska, USA</footnote>

<footnote n="7">⁷Global Flyway Network, Broome, Western Australia, Australia</footnote>

<footnote n="8">⁸Zoological Museum, Moscow MV Lomonosov State University, Moscow, Russia</footnote>

<footnote n="9">⁹Audubon Canyon Ranch, Cypress Grove Research Center, Marshall, California, USA</footnote>

<footnote n="10">¹⁰Groningen Institute for Evolutionary Life Sciences (GELIFES), University of Groningen, Groningen, The Netherlands</footnote>
Europe. In the global range of this long-distance migrant, we found evidence of both (1) fidelity to rigid behavioural routines promoting fine-scale geographic population structure (in the east) and (2) flexibility to colonise recently available migratory flyways and non-breeding areas (in the west). Our results suggest that cultural traditions in highly mobile vertebrates can override the expected effects of high dispersal ability on population structure, and provide insights for the evolution and flexibility of some of the world's longest migrations.

### KEYWORDS

bird migration, climate change, genetic differentiation, genotyping-by-sequencing, glacial refugia, *Limosa lapponica*, phylogeography, population genomics

# 1 | INTRODUCTION

Highly mobile species are expected to have a corresponding dispersal propensity, which should decrease the strength of physical or ecological barriers to gene flow, resulting in shallow genetic structure, or even panmixia, across large geographical areas (Palumbi, 1994; Theisen et al., 2008; Verkuil et al., 2012). Accordingly, a negative relationship between dispersal ability and population structure has been demonstrated in various taxa (Bohonak, 1999; Medina et al., 2018). Conversely, highly mobile vertebrates often display structured populations even in the apparent absence of barriers, for example in widely roaming marine species (Karl et al., 2011; Knutsen et al., 2003). This suggests that ecology and cultural or behavioural traditions can override the expected effects of mobility on dispersal and population structure (Rosenbaum et al., 2009; Sellas et al., 2005).

Seasonal migration may further reduce structure in mobile taxa by offering increased opportunities for dispersal among breeding populations (Moussy et al., 2013). This is particularly expected in birds, whose migrations may transcend the earth's most formidable ecological barriers, such as the Pacific Ocean (Piersma et al., 2022), the Sahara Desert (Schmaljohann et al., 2007) and the Himalayan mountain range (Hawkes et al., 2011). Avian migratory tendency and physiological flight capacity per se have been linked with increased dispersal and decreased population structure and speciation (Arguedas & Parker, 2000; Claramunt et al., 2012). Paradoxically, migration has also been linked with increased rates of speciation and diversification in birds (Rolland et al., 2014; Winger et al., 2014). Avian long-distance migration is often associated with high site-fidelity and flyway-specific behavioural and physiological adaptations that enable timely exploitation of predictable seasonal resources along the migratory route (Åkesson et al., 2017; Conklin et al., 2013; Winger et al., 2019). These traits predict low effective dispersal among migratory populations, through simple spatial separation, temporal isolation (Hendry & Day, 2005) or unsuccessful recruitment by immigrants (Nosil et al., 2005).

Additionally, because populations may have experienced very different climatic and geological conditions through time, present-day population structure is best understood with a phylogeographic

perspective (Avise, 2009; Knowles & Maddison, 2002). In birds, migration involves a highly integrated, multi-trait phenotype (including flight capacity, fuelling, phenology, navigation, moult, etc.; Åkesson et al., 2017; Piersma et al., 2005) that may imply limited flexibility to adapt to conditions outside those in which the migration evolved. Conversely, it is clear that some migratory bird populations have persisted through historical large-scale environmental perturbations, such as repeated Pleistocene glaciations (Batchelor et al., 2019; Hewitt, 2000); apparent responses included population sub-division and size changes, loss and gain of migration, and altered migration routes and distances (Avise & Walker, 1998; Conklin et al., 2022; Milá et al., 2006; Tan et al., 2023; Zink & Gardner, 2017). More empirical studies are needed to understand how dispersal and migration have interacted to generate and maintain population structure over time.

The bar-tailed godwit, *Limosa lapponica*, is a migratory shorebird (Order Charadriiformes, Family Scolopacidae) known for performing some of the longest migrations in the animal world. In particular, the migration of the Alaska-breeding subspecies (*L. i. baueri*) includes three separate trans-oceanic flights of >5000 km, including a southward flight of c. 12,000 km across the Pacific Ocean from Alaska to New Zealand (Figure 1), the longest non-stop flight recorded in any landbird (Conklin et al., 2017; Gill et al., 2009). There are six recognised subspecies of bar-tailed godwits (Bom et al., 2022; Engelmoer & Roselaar, 1998; Tomkovich, 2010), which vary greatly in migratory distances (Figure 1) and expected phylo-geographic histories. For example, the present-day range of bar-tailed godwits in the western Palearctic (Figure 1) was strongly affected by the glaciation of Europe during the last glacial maximum (LGM, c. 20,000–25,000 years before present, ybp; Hewitt, 2000). In particular, the entire present-day breeding range of the *lapponica* subspecies, which makes the shortest migration in the species (one-way distance of c. 2500–4500 km), was unavailable at the LGM. By contrast, the eastern half of the bar-tailed godwit breeding range falls within Beringia, an area of the Arctic that stretches from the Lena River in northeast Russia to the Mackenzie River in northwest Canada (Figure 1) and was largely ice-free at the LGM (Ehlers & Gibbard, 2007; Pielou, 1991) and previous glacial cycles (Batchelor et al., 2019). Unravelling the phylogeographic history of bar-tailed godwits may help understand the evolution of *baueri's*
<figure><img src="image_3.png" /><figcaption>**FIGURE 1** Global distribution and sampling of bar-tailed godwits. For each of six recognised subspecies (indicated by colour), arrows indicate general migration routes between northern breeding areas and southern non-breeding (i.e. boreal winter) areas (coloured areas). Numbers indicate total individuals sampled in each area (black circles); three shades of blue indicate sampling for *baueri* North, Central and South. Letters refer to sampling sites (in black; see Table S1 for details) and other locations mentioned in the text (in grey): A, Murmansk, Russia; B, Wadden Sea, The Netherlands; C, Banc d'Arguin, Mauritania; D, Taimyr Peninsula, Russia; E, Barr al Hikman, Oman; F, Roebuck Bay, Australia; G, SE Chukotka, Russia; H, Colville River, Alaska; I, Seward Peninsula, Alaska; J, Yukon-Kuskokwim Delta, Alaska; K, North Island, New Zealand; L, South Island, New Zealand; M, West Siberian Lake; N, West Siberian Plain; O, Lena River Delta; P, Bering Land Bridge; Q, Mackenzie River Delta. Dashed box indicates approximate extent of the Beringia region.</figcaption></figure>

trans-Pacific migration – for example, whether the physiological capabilities it requires arose in situ, or were products of eastward expansion from the Palearctic (Hedenström, 2010; Piersma et al., 2022).

In this study, we describe the history of divergences and degree of neutral genetic differentiation among global flyway populations of bar-tailed godwits. We use nextera-tagmented, reductively amplified DNA (nextRAD) sequencing (Russello et al., 2015) for de novo discovery of genome-wide single-nucleotide polymorphisms (SNPs) for population genetic analyses, and compare hypothesised evolutionary scenarios in an Approximate Bayesian Computation (ABC) framework using DIYABC-RF (Collin et al., 2021). By reconstructing the recent evolutionary history of bar-tailed godwits, we (1) describe the historical context for present-day population structure, (2) provide insights for the evolution and persistence of extreme non-stop flights, and (3) elucidate the origins and maintenance of population structure in long-distance migratory birds.

## 2 | MATERIALS AND METHODS

### 2.1 | Sampling and DNA extraction

We assembled DNA samples representing all recognised and hypothesised breeding populations within the global range of bar-tailed godwits (Figure 1, Table S1). Where possible, we used samples collected in known breeding areas; however, because bar-tailed

godwits breed in low densities in remote regions, for some populations there were few or no breeding samples available. In these cases, we used samples collected from non-breeding sites when sampled individuals could be confidently assigned to breeding areas either because they were remotely tracked to breeding areas using satellite-telemetry or light-level geolocation, or because long-term research programmes had established strong links between breeding and non-breeding areas (e.g. through mark-recapture/resight programmes). There were no verifiable samples available from the poorly studied *anadyrensis* population; however, we included eight samples collected from coastal southeast Chukotka, within 200km of the known *anadyrensis* breeding range (Figure 1; Tomkovich, 2010), in order to detect genetic variation that might indicate the distinctiveness of *anadyrensis*. All tissue samples were acquired from museum collections or collected by the authors and colleagues in the field under all requisite permits appropriate to their respective countries and institutions.

We extracted genomic DNA from samples using two methods. For blood or organ tissue samples preserved in 95% ethanol, we used the DNeasy Blood and Tissue Kit (Qiagen) following the manufacturer's instructions for tissue. For blood samples preserved in Queen's lysis buffer, we used the NucleoSpin Blood QuickPure Kit (Macherey-Nagel). Extract quality was first assessed on a 1.5% agarose gel to exclude extractions with insufficient yield or excessively degraded DNA. We then quantified DNA concentrations using a Qubit 3.0 fluorometer (Life Technologies), diluted extracts to achieve relatively even concentrations, and dried down samples
in a SpeedVac concentrator (ThermoFisher Scientific). We delivered 80–162 ng of DNA of 209 individual godwits (Table S1) for SNP discovery and genotyping.

## 2.2 | SNP genotyping using nextRAD sequencing

Genomic DNA was converted into nextRAD genomic fragment libraries (SNPsaurus, LLC, USA) following the method described by Russello et al. (2015). Genomic DNA was first fragmented with Nextera reagent (Illumina, Inc.), which also ligates short adapter sequences to the ends of the fragments. The Nextera reaction was scaled for fragmenting 10 ng of genomic DNA, although 20–30 ng of genomic DNA was used for input to compensate for degraded DNA in the samples and to increase fragment sizes. Fragmented DNA was then amplified, with one of the primers matching the adapter and extending 10 nucleotides into the genomic DNA with the selective sequence 'GTGTAGAGCC'. Thus, only fragments starting with a sequence that can be hybridised by the selective sequence of the primer were efficiently amplified. PCR amplification was done at 74°C for 27 cycles. The nextRAD libraries were sequenced on an Illumina HiSeq-4000 at the Genomics Core Facility, University of Oregon, USA.

For SNP-calling, reads were first trimmed using custom scripts (SNPsaurus, LLC) in bbduk (BBMap tools; Bushnell, 2016). Next, a de novo reference was created from abundant reads (after removal of low-quality (phred-scale quality <20) and very high-abundance reads) and reads that aligned to these. All 205,854,244 reads were mapped to the reference with an alignment identity threshold of 95% using bbmap (BBMap tools). Genotype calling was performed using SAMtools and BCFtools (Danecek et al., 2021; samtools mpileup -gu -Q 10 -t DP, DPR -f ref.fasta -b samples.txt | bcftools call -cv - > genotypes.vcf), applying a minimum read depth filter of 7×. The genotype table was then filtered using VCFtools v.01.14 (Danecek et al., 2011) to remove SNPs called in <80% of samples and putative alleles with a population frequency <3%, to exclude artefactual variants. The resulting VCF file included 6424 unique loci (150-bp sequences) containing 21,031 SNPs.

Additional cleaning was performed using VCFtools to remove indels (n=478 SNPs), low-quality SNPs (phred-scale quality [QUAL] <999; n=15 SNPs), and samples that failed to sequence (n=1 menzbieri individual) or were of unknown breeding population (n=9 individuals from New Zealand that could not be confidently assigned to a breeding region within Alaska based on tracking data). Next, SNPs deviating from Hardy-Weinberg (HW) proportions in all seven hypothesised populations (excluding *anadyrensis*) were identified and removed using VCFtools (p <.05; n=48 SNPs); the small sample size in *anadyrensis* made it statistically insensitive to HW deviations. To ensure independence of loci, we used PLINK v.1.9 (Chang et al., 2015) to identify and remove all SNPs in linkage disequilibrium (LD; r² > .20; n=6188 SNPs). The LD-pruned data set included 199 individuals and 14,318 unlinked SNPs on 6140 loci.

We then used VCFtools to calculate proportion of missing SNP calls per individual (range 1%–77%) and removed 15 individuals with >30% missing data. Because inclusion of related individuals can bias population genetic analyses (Anderson & Dunham, 2008; Rodríguez-Ramilo & Wang, 2012), we estimated individual pairwise relatedness in PLINK using the Identity-by-Descent estimator PI-HAT; we found no cases of relatedness involving half-siblings or closer (PI-HAT >0.20). The final data set included 184 individuals and 14,318 SNPs. We used VCFtools, PLINK and PGDspider v.2.1.0.3 (Lischer & Excoffier, 2012) to convert data to different formats required for further analysis.

## 2.3 | Inference of population structure and diversity

To assess major axes of genetic variation and structure, we performed a principal component analysis (PCA; Patterson et al., 2006) using the R (R Core Team, 2021) packages *gdsfmt* v.1.14.1 and *SNPRelate* v.1.12.2 (Zheng et al., 2012). We determined the number of ancestral populations or clusters (K) and ancestry proportions of these for each individual using ADMIXTURE v.1.3.0 (Alexander et al., 2009). We conducted 20 replicate runs (with random seeds) for each putative number of clusters (K) ranging from 1 to 8. We used cross-validation (CV) (Alexander & Lange, 2011) to assess how the CV error rate varied with increasing K values and which value provided the lowest CV error. Then, we used the CLUMPAK (Cluster Markov Packager Across K; Kopelman et al., 2015) web server ([http://clumpak.tau.ac.il/](http://clumpak.tau.ac.il/)) with default settings to summarise estimates of individual ancestry proportions to each cluster across replicate runs, inspect potential distinct solutions identified and visualise the most likely ancestry proportions at each value of K.

Globally and for each putative population, we characterized genetic diversity by calculating nucleotide diversity (π), heterozygosity and inbreeding coefficient (FIS) using VCFtools. We evaluated evidence for deviations from mutation-drift equilibrium, indicative of departures from a null hypothesis of constant population size, by estimating the per-locus Tajima's D values for each population using VCFtools. We estimated the degree of genetic differentiation among populations by calculating pairwise FST (Weir & Cockerham, 1984) using the *diffCalc* function in the R package *diversity* v.1.9.90 (Keenan et al., 2013), with 95% confidence intervals derived from 500 bootstraps. We calculated p-values for FST estimates using the *pairwiseTest* function (1000 permutations) in the R package *strataG* v.2.4.905 (Archer et al., 2017).

## 2.4 | Population evolutionary relationships and demographic history

To visualise evolutionary relationships among populations, we first constructed a midpoint-rooted neighbour-joining (NJ) tree based on Nei's genetic distance (Nei, 1972) using the R packages *poppr* v.2.8.5
(Kamvar et al., 2014) and *ape* v.5.3 (Paradis & Schliep, 2019), with missing data replaced by mean allele counts, and node support calculated from 1000 bootstraps.

To identify likely potential sources of admixture in the tree, we first conducted the three-population test for admixture (Reich et al., 2009) in R using ADMIXTOOLS 2.0.0 (Maier et al., 2023) for all possible population triads. Briefly, the *f*<sub>3</sub> statistic (Patterson et al., 2012) assesses shared genetic drift to test the null hypothesis that a target population (A) is related to two other populations (B, C) in a simple tree-like manner. We tested significance of the *f*<sub>3</sub> test using a blocked-jackknife resampling of SNPs, considering blocks of 100 SNPs. A significantly negative value of *f*<sub>3</sub> (A;(B,C)) indicates that the target population must have arisen from an admixture of the other two populations, or populations closely related to them.

To further investigate the genetic relationships among god-
wis populations, we estimated population networks using the ap-
proach implemented in ADMIXTOOLS 2.0.0 (Maier et al., 2023),
which uses the shared and private genetic ancestry components of
genome-wide allele frequency data to infer population branching
while accounting for potential historical migration and admixture
events among populations. To find the best network topology, we
used the function *find graphs* for 0–5 admixture events (m), each
with 100 replicates and a maximum of 300 generations. For each
replicate, we calculated the out-of-sample score to account for the
difference in likelihood score caused by an increasing number of
admixture events, which also increase the degrees of freedom. We
used the *baueri* South population as a root to fold the networks;
even if that population is not an actual outgroup, it allowed visualis-
ing population relationships in a consistent way. The actual popula-
tion branching order was later assessed statistically using the ABC
(Beaumont et al., 2002) Random Forest (RF) statistical framework
(Pudlo et al., 2016; Raynal et al., 2019). For the best-supported rep-
licates (i.e. those displaying the likelihood score closest to zero), we
assessed the goodness-of-fit with the R package *admixture-graph*
v.1.0.2 (Leppälä et al., 2017). This approach allows comparing the
observed and expected values of *f*<sub>4</sub> statistics among the different
alternatives and identifying the graph that best fits the data (i.e. the
one properly predicting the observed *f*<sub>4</sub> values for all possible four-
population combinations).

Finally, we evaluated the support for different possible sce-
narios of population divergence and admixture using the ABC-RF
statistical framework (Pudlo et al., 2016; Raynal et al., 2019) in
the R package DIYABC-RF v.1.2.1 (Collin et al., 2021). ABC-RF can
estimate posterior probabilities of historical scenarios, based on
massive coalescent simulations of genetic data. Simulations are
compared to observed SNP data using summary statistics to iden-
tify the best-fitting model by calculating the number of RF classifi-
cation votes and to approximate the posterior probability for the
best model (Pudlo et al., 2016). The best-fitting posterior param-
eter distribution values for the best model can be estimated using
a RF procedure applied in a regression setting (Raynal et al., 2019).
Estimated parameters include the effective size (N<sub>e</sub>) for each

population, admixture events and their rates (r), and split and ad-mixture times (t). We evaluated potential evolutionary scenarios in a stepwise manner, as follows.

As our primary question arising from the above analyses concerned the nature and relative timing of the admixed origin of *menzbieri* (see Section 3), we sought to first reduce the number of potential scenarios by fixing the topology of the relationships among the three *baueri* populations. So, in Step 1, we tested whether the apparently clinal variation among these populations was best represented by one of two possible simple branching topologies or with *baueri* Central as an admixture of *baueri* North and South (Figure S4, Table S2).

In Step 2, we created four possible scenarios (Figure S5c-f) for an admixed origin of *menzbieri* from neighbouring populations to the west and east, while allowing populations on the two main branches to diverge at different times relative to each other. Essentially, the admixed origin of *menzbieri* could occur: (1) before either branch began to diverge (Figure S5c); (2) after populations on one branch diverged but before the other branch began to diverge (Figure S5d,e) or (3) after populations on both branches had diverged (Figure S5f). To compete with these scenarios, we created two ‘null-hypothesis’ scenarios in which *menzbieri* was not admixed, but simply diverged from one of the two main branches (Figure S5a,b).

The scenario parameters were considered as random vari-
ables drawn from prior distributions (Tables S2 and S3). We used
DIYABC-RF to simulate 20,000 genetic data sets per scenario, with
the same properties as the observed data set (number of loci and
proportion of missing data). Simulated and observed data sets were
summarised using the whole set of summary statistics proposed
by DIYABC-RF for SNP markers (Table S4), describing population
genetic variation (e.g. heterozygosity and proportion of monomor-
phic loci), differentiation (e.g. FST and Nei's genetic distances) or
admixture (e.g. admixture coefficient, f3 and f4 statistics; Patterson
et al., 2012); see the full list and details in Table S4. Linear discrimi-
nant analysis components were also used as additional summary sta-
tistics (Estoup et al., 2012). The total number of summary statistics
was 50 and 1708 for Step 1 and 2 respectively.

We used the RF classification procedure (Pudlo et al., 2016) to
compare the likelihood of the competing scenarios at each step with
DIYABC-RF. RF is a machine-learning algorithm that uses hundreds
of bootstrapped decision trees to perform classification, using the
summary statistics as a set of predictor variables. Some simulations
are not used in decision tree building at each bootstrap (i.e. the out-
of-bag simulations), and are used to compute the ‘misclassification
error rate’, also known as the ‘prior error rate’; this metric provides
a direct method for estimating the CV error rate (Pudlo et al., 2016).
The classification performance, and thus the discrimination power
among the tested scenarios under the present simulation procedure,
can be directly visualized using the confusion matrix that is produced
at this step. The confusion matrix provides a global account of how
well simulated data sets generated under each scenario are properly
classified (or not) by the RF procedure to the scenario from which
they were generated (Schrider & Kern, 2018). The confusion matrix
can thus be used to approximate Type-I and Type-II error rates associated with the selected scenario.

At each step, we built a training set of 20,000 simulated data sets per scenario, with the same number of loci and individuals as the observed data set, and then grew classification forests of increasing size (100, 500 and 2000 trees) to test convergence. The RF computation provides a classification vote for each scenario (i.e. the number of times a model is selected from the decision trees). We selected the scenario with the highest classification vote as the most likely scenario, and we estimated its posterior probability following the recommendation of Pudlo et al. (2016). We assessed the global performance of the classification procedure and the confidence in the ABC-RF scenario choice procedure, by calculating the prior error rate based on the available out-of-bag simulations, together with the confusion matrix and the associated Type-I and Type-II error rates. We repeated the RF analysis (3 times in Step 1 and 10 times in Step 2) to ensure that the results converged.

For the best-supported model in Step 2, we estimated posterior distribution values of all parameters using a regression by RF methodology (Raynal et al., 2019) in DIYABC-RF. We used a classification forest of 500 decision trees, based on a training set of 100,000 simulations, to conduct the parameter estimations. We converted divergence time estimates to years assuming a generation time of 8 years (delayed maturity with adult annual survival c. 0.86; Méndez et al., 2018) and the genome-wide mutation rate calculated by Zhang et al. (2014) for Order Charadriiformes: $1.5 \times 10^{-9}$ substitutions per site per year. All of the steps involved in the ABC-RF analysis (simulations, computation of summary statistics, model checking, scenario

comparisons and estimations of parameter posterior distributions) were performed with DIYABC-RF.

# 3 | RESULTS

## 3.1 | Summary of nextRAD SNP data set

The final data set comprised 184 unrelated bar-tailed godwits genotyped at 14,318 unlinked high-quality SNPs; this included 15–31 individuals in each of seven hypothesised populations, but only three purported *anadyrensis* individuals (Table S1). On average, individuals were genotyped at 93.7% (range: 72%–99%) of SNPs and with a mean read depth of 50.4 (range: 15–116). Each SNP was genotyped in an average of 172.5 (range: 149–184 or 81%–100%) individuals. Globally, nucleotide diversity ($\pi$) was 0.245 and heterozygosity was 20.6%.

## 3.2 | Population structure and diversity

The first three axes (PC1–3) of the PCA explained 4.0% of the total genetic variation (Figure 2a). On PC1, individuals were essentially distributed on an axis of longitude in three clusters: a western Palearctic group containing all of *lapponica*, *taymyrensis* and *yamalensis*; an Alaska group containing all of *baueri*; and an intermediate group of *menzbieri*, suggesting an admixed origin (Figure 2a). PC2 separated 12 *lapponica* individuals from the main western Palearctic

<figure><img src="image_4.png" /><figcaption>FIGURE 2 (a) Population structure estimated by principal component analysis. The individual scores ($n=184$ individuals) for the first three principal components (PC1–3) are shown. The scree plot (inset) indicates the proportion of explained genetic variance by each PC, with PC1–3 explaining 4.0% of total variation (see Figure S1 for additional PC axes). (b) Individual genetic ancestries assigned to major clusters for $K=2-5$ estimated using ADMIXTURE. At each value of $K$, the ancestry proportions for the 184 individuals for the dominant solution were determined by CLUMPAK summary of 20 replicate runs. Numbers in parenthesis indicate proportion of replicate runs contributing to the dominant solution (see Figure S2 for minor clusters for $K=3-5$). Names below plot indicate sampling locations (see Figure 1 and Table S1); colours indicate a priori purported breeding populations.</figcaption></figure>
cluster. PC3 stretched *baueri* individuals along a latitudinal cline, partly separating the North, South and Central populations. The three purported *anadyrensis* individuals did not form a coherent cluster, one clustering with *menzbieri* and two with *baueri*. Higher-order axes of variation (PC4–6) explained an additional 2.2% of variation but revealed negligible further sub-structure (Figure S1).

Genetic ancestry analyses with ADMIXTURE for K=2 ancestral populations provided the lowest CV error criterion (Figure S2) and 100% support for the major clusters across 20 replicate runs. At K=2, populations were clearly divided between a western Palearctic cluster and an Alaska cluster, with *menzbieri* reflecting admixture between the two (Figure 2b). Additional structure identified in the PCA was apparent at K=3 and 4; that is, the clear structure within purported *lapponica* and the clinal variation within *baueri* respectively. At K=5, the identification of a predominantly *menzbieri* cluster may suggest significant genetic drift specific to that population, as opposed to an origin of simple and recent admixture. Higher values of K revealed no further important geographical structure (Figure S2). Again, the purported *anadyrensis* individuals did not form a coherent cluster, but appear to include one *menzbieri* and two *baueri* individuals (Figure 2b).

In both PCA and ADMIXTURE, only six of 27 purported *lapponica* individuals sampled on wintering grounds in The Netherlands clearly grouped with the four individuals sampled at known *lapponica* breeding areas in western Russia (see Figure 1, Table S1). Of the remaining 21, 19 clearly grouped with *taymyrensis* (sampled at breeding areas in central Russia and a non-breeding area in Mauritania), and two appear to reflect admixture between *lapponica* and *taymyrensis* (Figure 2a,b). Therefore, for further analysis, we considered a new population of ‘*taymyrensis* Europe’ (n=19), and a reduced *lapponica* population (n=10 ‘true’ *lapponica* plus the 2 admixed individuals). For clarity, we hereafter refer to the original *taymyrensis* samples as ‘*taymyrensis* Africa’ (n=29). As PCA and ADMIXTURE results do not

support an identifiable *anadyrensis* population, we excluded those samples (n=3) from further analysis. Therefore, we proceeded with 181 individuals in eight populations.

Globally, nucleotide diversity was 0.245, and this was similar for all populations (Table S5). Mean heterozygosity was also generally uniform (range of means: 0.198–0.215), with no evidence of inbreeding in any population (all mean $F_{IS}$<0.20; Table S5). Values of Tajima's D were low and similar among populations (range of means: 0.214–0.339; Table S5), suggesting no evidence of dramatic population-size changes or additional hidden structure.

Pairwise $F_{ST}$ values among the eight hypothesised populations ranged 0.001–0.048 (Table 1), and were greatest between the westernmost and easternmost breeding populations (*lapponica* vs. *baueri*; $F_{ST}$ values ranging from 0.046 to 0.048 for the three *baueri* subpopulations). The two groups that winter in Europe, all previously considered *lapponica* (*lapponica* and *taymyrensis* Europe), were moderately differentiated ($F_{ST}$=0.016). Differentiation among the three sub-populations of *baueri* ($F_{ST}$=0.006–0.009) was lower than for most recognised subspecies pairs, with 95% CIs that included zero. The lowest values were among the three groups breeding in central Russia (*taymyrensis* Europe, *taymyrensis* Africa and *yamalensis*; $F_{ST}$=0.001–0.003).

## 3.3 | Population evolutionary relationships and demographic history

The NJ tree based on Nei's distance (Figure 3a) again identified two main branches (western Palearctic and Alaska groups), with low node support within each group, and *menzbieri* occupying an intermediate position between them.

The three-population test (based on the Patterson's $f_3$ statistic of the form (A;(B,C)); Table S6) identified *menzbieri* as a likely product

<table><caption>TABLE 1 Mean population pairwise $F_{ST}$ and 95% confidence intervals of estimates.</caption><thead><tr><th rowspan="2"></th><th colspan="2">taymyrensis</th><th colspan="2">taymyrensis</th><th colspan="2">baueri</th></tr><tr><th>lapponica</th><th>Europe</th><th>Africa</th><th>yamalensis</th><th>menzbieri</th><th>North</th><th>Central</th></tr></thead><tbody><tr><td>taymyrensis Europe</td><td>0.016**<br>[0-0.042]</td><td></td><td></td><td></td><td></td><td></td><td></td></tr><tr><td>taymyrensis Africa</td><td>0.014**<br>[0-0.040]</td><td>0.001*<br>[0-0.014]</td><td></td><td></td><td></td><td></td><td></td></tr><tr><td>yamalensis</td><td>0.015**<br>[0.0003-0.040]</td><td>0.003**<br>[0-0.016]</td><td>0.002**<br>[0-0.010]</td><td></td><td></td><td></td><td></td></tr><tr><td>menzbieri</td><td>0.026**<br>[0.001-0.048]</td><td>0.011**<br>[0.002-0.024]</td><td>0.010**<br>[0.004-0.020]</td><td>0.013**<br>[0.006-0.022]</td><td></td><td></td><td></td></tr><tr><td>baueri North</td><td>0.048**<br>[0.031-0.073]</td><td>0.032**<br>[0.021-0.045]</td><td>0.031**<br>[0.023-0.042]</td><td>0.034**<br>[0.026-0.045]</td><td>0.016**<br>[0.007-0.026]</td><td></td><td></td></tr><tr><td>baueri Central</td><td>0.047**<br>[0.028-0.071]</td><td>0.032**<br>[0.018-0.053]</td><td>0.031**<br>[0.019-0.045]</td><td>0.033**<br>[0.022-0.048]</td><td>0.015**<br>[0.003-0.031]</td><td>0.007**<br>[0-0.027]</td><td></td></tr><tr><td>baueri South</td><td>0.046**<br>[0.030-0.070]</td><td>0.031**<br>[0.021-0.044]</td><td>0.030**<br>[0.024-0.040]</td><td>0.033**<br>[0.026-0.040]</td><td>0.014**<br>[0.007-0.024]</td><td>0.009**<br>[0-0.020]</td><td>0.006**<br>[0-0.023]</td></tr></tbody></table>

Note: All comparisons are significantly different from zero (1000 permutations; *p*=.038, **p**<.001).
<figure><img src="image_5.png" /><figcaption>FIGURE 3 (a) Unrooted neighbour-joining tree, based on Nei's minimum distance. Numbers indicate bootstrapped node confidence (%). (b) AdmixtureGraph showing best population graph with two migration reticulations (m = 2), arbitrarily rooted with *baueri* South population (see Section 2). Numbers indicate relative genetic drift (unlabelled branches = 0 distance). Dotted lines represent admixtures, with percentages indicating admixture proportions and arrows indicating direction of admixture.</figcaption></figure>

of admixture between the two main branches identified by the NJ tree (see Figure 3a). Of 168 possible population triads, 12 produced a significantly negative $f_3$ statistic; all of these included *menzbieri* as the target population (A) and one member of each of the two main branches (*lapponica/taymyrensis/yamalensis* or the *baueri* group respectively) as the source populations (B, C) (Table S6).

For admixture graphs estimated by ADMIXTOOLS 2, the greatest increase in likelihood occurred at $m=2$ admixture edges in the population graph (Figure S3a), but out-of-sample scores indicated the highest likelihood at $m=3$ (Figure S3b). Overall, the five best-supported replicates (i.e. those with likelihood scores closest to zero) included two distinct network topologies for $m=3$ and one for $m=2$. The best-supported topologies at $m=2$ and 3 (Figure S3c,d) each provided a perfect fit with the data as estimated by Admixture-graph (i.e. no outlying observed $f_4$ statistic values compared to those

theoretically expected from these graphs; Figure S3e,f). Therefore, according to parsimony we present the model at $m=2$ (Figure 3b). In this graph, *menzbieri* arose from the admixture of a population ancestral to present-day western Palearctic populations (86% contribution) and a population closely related to *baueri* North (14%). Additionally, *baueri* North has received 4% ancestry from a population closely related to *taymyrensis* Europe.

We further investigated the population evolutionary history, the branching order, admixture events and times of these events in the ABC-RF analyses. In Step 1, we tested whether the apparently clinal variation among the three *baueri* populations was best represented by a simple bifurcating topology or by historical admixture (Figure S4). The best-fitting scenario, representing *baueri* Central as an admixture of *baueri* North and South, received the greatest proportion of RF votes (44.5%), with a posterior probability of 51.0% (Figure S4, Table S7). This result was consistently recovered over nine independent RF analyses with 500, 1000 and 2000 RF decision trees (Table S7). This ABC-RF analysis had good power to discriminate among the three scenarios, as suggested by the low global *misclassification error rate* (also known as *prior error rate*) of 11.4% (Figure S4, Table S7). The confusion matrix (Table S8) further showed that our simulation design delivered good power to discriminate among the three alternative scenarios, producing an overall average of 11.4% misclassified simulations across scenarios. Following these Step 1 results, we fixed the topology of the main Alaskan branch accordingly in the global scenarios evaluated in Step 2.

In Step 2 (Figure 4 and Figure S5), we compared six possible scenarios for the origin of *menzbieri*: two 'null-hypothesis' scenarios in which *menzbieri* simply diverged from one of the two main branches (a,b), and four in which an admixed *menzbieri* arose at different time points relative to the diversification within each branch (c-f). The best-supported scenario was c, in which the historical admixture occurred prior to the main diversifications within both the western Palearctic and Alaska groups (Figure 4c); this scenario received (mean ± SD) 26.4 ± 0.9% of votes, with a posterior probability of 46.1 ± 2.6% (Figure 4 and Figure S5, Table S9). The ABC-RF analysis showed a very high level of convergence and consistency, selecting scenario c as the most probable scenario in nine of 10 replicates (Table S9). However, this scenario choice was closely challenged by scenario e, which received 25.3 ± 1.0% of the votes (Figure 4 and Figure S5, Table S9). Scenario e was indeed most similar to c, but with the admixture occurring before diversification among *baueri* populations and after *lapponica* diverged from the other western Palearctic populations (Figure 4e). The global discrimination power among the six scenarios based on our RF model choice procedure was good, but showed some difficulties in discriminating between scenarios c and e. While the overall *misclassification error rate* (or *prior error rate*) was 22.0 ± 0.0%, the confusion matrix obtained from the simulations under the priors of each scenario showed that 41% of the simulations generated under scenario c were correctly identified as originating from that scenario, but that 26% were wrongly assigned to scenario e (Table S10). This indicates a limited sensitivity to discriminate
<figure><img src="image_6.png" /><figcaption>FIGURE 4 Scenarios tested in Step 2 of DIYABC-RF analysis: Six possible scenarios for the origin of *menzbieri*, by either simple branching (a,b) or by admixture of the two main branches at relatively ancient (c), intermediate (d,e) or recent (f) time periods. Extant (sampled) populations are indicated by colours; inferred historical populations are shown in shades of grey. The best-supported scenario was (c), followed by (e), as inferred by the proportion of Random Forest classification votes; see Figure S5 and Tables S9 and S10 for details.</figcaption></figure>

between these two scenarios. Nevertheless, the specificity was good, as suggested by the low proportion (≤11.8%; Table S10) of simulations generated by other scenarios that were incorrectly classified as scenario c.

Each of the 20 demographic parameters of the best-supported scenario (c) was estimated within the ABC-RF framework (Table S11). Timing parameter estimates support the existence of three main lineages at the LGM, with the divergence of the western Palearctic and Alaskan branches occurring c. 38.2K ybp (modal estimate; 5th–95th quantiles: 22.0–93.1K), followed by an admixture giving rise to *menzbieri* c. 22.0K ybp (13.4–52.5K; Figure 5, Table S11). The post-LGM warming period featured *lapponica* diverging from the western Palearctic lineage (c. 16.4K ybp; 9.1–28.9K) and the diversification of *baueri* populations (c. 11.6–12.3K ybp); although we modelled the *baueri* diversification as a divergence followed by an admixture (see Step 1), the two time estimates were essentially indistinguishable, indicating a relatively contemporaneous radiation of three populations, rather than contact after a prolonged isolation (Figure 5). Most recently, the remaining western Palearctic lineage diversified into *taymyrensis* (Europe and Africa) and *yamalensis* populations (c. 4.5–5.2K ybp); these two divergence times were also indistinguishable, rather than clearly sequential (Figure 5).

## 4 | DISCUSSION

With their broad geographical distribution among relatively discrete migratory flyways (Figure 1) and their propensity for extended non-stop flights across formidable ecological barriers (Conklin et al., 2017; Piersma et al., 2022), bar-tailed godwits represent an elegant test of the consequences of high mobility for population structure. Within this species, we found evidence for both ends of a spectrum: in the west, negligible structure between two populations with adjacent breeding areas but separate ranges in the non-breeding season (*taymyrensis* and *yamalensis*), and in the east, clinal structure within a single subspecies, despite individuals freely mixing throughout the non-breeding season (*baueri*). As we discuss below, reconciling these seemingly incongruous patterns requires understanding the contrasting phylogeographic histories of these two regions.

### 4.1 | Population structure

Our results reveal a more complicated population structure in bar-tailed godwits than was previously known. Before the subspecific recognition of *taymyrensis* (Engelmoer & Roselaar, 1998) and
<figure><img src="image_7.png" /></figure>

**FIGURE 5** Best-supported Scenario in Step 2 of DIYABC-RF analysis (see Figure 4), scaled to relative time-parameter estimates (converted to years assuming a generation time of 8 years) for five divergence events (branches) and two admixture events (horizontal bars). On left, posterior distributions for each time parameter estimate, indicating mode, mean, median and interval of 5th–95th quantiles (vertical line). For reference, estimated timing of Last Glacial Maximum, Holocene Climatic Optimum, and most recent presence of Bering Land Bridge and glaciation of Europe are shown.

*tamalensis* (Bom et al., 2022), all bar-tailed godwits breeding in the western Palearctic (Figure 1) were considered differentially migrating populations of the nominate subspecies *L. l. lapponica*. Previous genetic work showed these three subspecies were indistinguishable at the mtDNA control region (Bom et al., 2022), consistent with a recent post-glacial expansion into three breeding areas and migration routes. However, in genome-wide SNPs, we found that *lapponica* was strongly differentiated from *taymyrensis* and *tamalensis*, which were weakly differentiated from each other (Table 1). This suggests a post-LGM colonisation by *lapponica* of recently deglaciated Europe, followed later by emerging structure along two flyways in the source population breeding in the central Russian Arctic.

We found unexpected structure among samples collected in the Dutch Wadden Sea during mid-winter, purported to belong to *lapponica* (Duijns et al., 2012): of these 27 samples, only six grouped with the four individuals sampled in the breeding area of *lapponica*, whereas 19 clearly clustered with *taymyrensis*, and another two appeared to represent admixture between the two subspecies (Figure 2). *Lapponica* and *taymyrensis* are generally considered ‘leap-frog’ migrants, with the entirely of *taymyrensis* wintering in west Africa, south of European-wintering *lapponica* (Figure 1; Duijns et al., 2012). However, our results strongly indicate that both *taymyrensis* and *lapponica* are present in western Europe in the winter. Currently, population estimates for *lapponica* (150,000–180,000 and increasing) and *taymyrensis* (380,000–420,000 and declining) are based entirely on winter census data (van Roomen et al., 2022), assuming no mixing of subspecies. Our

results from the Wadden Sea (indicating c. 70% *taymyrensis*) should be confirmed through genetic assignment of a geographically broader sample of the European wintering population, followed by a re-evaluation of the population sizes and conservation status of both *lapponica* and *taymyrensis*.

Within *baueri*, we confirm that behavioural and morphological differences observed along a latitudinal cline in the Alaska breeding range (Conklin et al., 2010, 2011; Conklin & Battley, 2011) are mirrored by a similar genetic cline (Figure 2). Bar-tailed godwits from the entire Alaska breeding range regularly mix at non-breeding sites in New Zealand (Battley et al., 2020; Conklin et al., 2011) and Australia (Wilson et al., 2007) and at migratory staging sites in the Yellow Sea and south-west Alaska (Battley et al., 2012). On northward migration, *baueri* also mixes freely with *menzbieri* in the Yellow Sea (see Figure 1; Battley et al., 2012; Choi et al., 2015). The clear genetic structure maintained among (sub-) populations that regularly meet in the non-breeding season implies that behavioural rigidity, in terms of strict migration timing and fidelity to a breeding region, is sufficient to preclude panmixia among and within subspecies. Within *baueri*, previous research found no clear link between migration timing, which generally follows the latitudinal genetic cline we detected, and variation at the gene Clock, thought to be associated annual-cycle scheduling in birds (Parody-Merino et al., 2019). Therefore, it remains unclear whether the phenotypic and genetic cline within *baueri* is primarily driven by neutral genetic drift or involves selection, and whether it is emerging or stably maintained over a long time period.
Due to the small number and high sequencing failure rate of purported *anadyrensis* samples in our study, we could not resolve the status or history of that subspecies, which has a distinct migration route through eastern Russia but an incompletely described non-breeding range (Chan et al., 2022; Tomkovich, 2010). The three samples in our final data set did not cluster together, grouping with either *menzbieri* (*n* = 1) or *baueri* (*n* = 2; Figure 2), suggesting that these individuals, sampled in southeast Chukotka, were most likely passage migrants from the neighbouring subspecies, rather than from the closest breeding population (Figure 1). Therefore, genetic evaluation of the status of *anadyrensis* awaits analysis of verifiable members of the subspecies.

## 4.2 | Phylogeographic history

Our demographic reconstruction in DIYABC-RF provides the first estimates for divergence times in bar-tailed godwits, and new insights regarding the historical origins of present-day flyway populations and migratory routes. The divergence time estimate of c. 38K ybp for the two main lineages suggests the ancestors of the western Palearctic and Alaskan populations existed in Siberia and Beringia well before the LGM (Figure 5). This corresponds to a period (c. 40K ybp) when the Bering Sea separated Asia and North America (Farmer et al., 2023) and much of the Chukotka Peninsula was glaciated (Batchelor et al., 2019); either or both of these factors could have caused or perpetuated isolation of the two lineages. We inferred a subsequent admixture of the two lineages c. 22K ybp, giving rise to present-day *menzbieri* (Figure 5); the wide confidence interval around this estimate (13.4–52.5K ybp), and the relatively similar support for two tested scenarios, means that this admixture could have occurred just before, during, or after the LGM. During this entire period, Beringia was largely unglaciated (Batchelor et al., 2019) and lower sea levels resulted in an intact Bering Land Bridge (c. 11–36K ybp; Farmer et al., 2023), a situation that could have created new breeding habitat available for colonisation from both the west and the east, or conceivably re-connected breeding populations in a continuous swath of tundra habitat stretching from Chukotka to Alaska. Indeed, reconstructed vegetation histories indicate that at least part of the Bering Land Bridge hosted graminoid and herbaceous tundra, the preferred breeding habitat of bar-tailed godwits, during the LGM (Wang et al., 2017). Presumably, subsequent re-flooding of the Bering Strait (by c. 11K ybp) and incremental reduction and isolation of tundra habitat associated with the warming period leading to the mid-Holocene Climatic Optimum (c. 8K ybp; Kraaijeveld & Nieboer, 2000; Stewart & Dalén, 2008; Wauchope et al., 2017) served to isolate these three extant lineages into their present-day breeding ranges (Figure 1). Similar structure across Beringia is found in numerous avian species (Winker et al., 2023), and the deeper biogeographical history of the region, which includes the periodic appearance and disappearance of the Bering Land Bridge, is reflected in the complex histories of diverse taxa, including plants and humans (Hoffecker et al., 2023; Wen et al., 2016).

Because significant portions of Alaska have remained unglaciated for hundreds of thousands of years (Batchelor et al., 2019; Winker et al., 2023) and our results suggest no temporal limit to *baueri*'s occupation of its present-day breeding range, we can infer no theoretical maximum age of that population's trans-Pacific migration (Gill et al., 2009). Therefore, we cannot help discern between scenarios of in situ evolution of these extremely long non-stop flights in the Pacific versus incremental up-scaling of flight distances through eastward expansion of godwits from Asia (Hedenström, 2010; Piersma et al., 2022). We can only infer that any such eastward expansion did not appear to occur recently (e.g. post-LGM), as *baueri* has apparently existed in Alaska for sufficient time for genetic and morphological structure to develop along a latitudinal cline (Conklin et al., 2011). Paleoclimatic modelling, to assess change in both atmospheric circulation and habitats in the Pacific Basin across time, may reveal more about how long *baueri* could theoretically have been migrating across the Pacific in the current manner (Piersma et al., 2022).

If bar-tailed godwits existed in a Siberian-Beringian refugium (stretching from the West Siberian Plain to Alaska) at the LGM and then expanded westward with the de-glaciation of Europe, then we would expect *lapponica* to have diverged quite recently from other western Palearctic lineages, particularly given that their present-day breeding range in Fennoscandia was not available until c. 9–12K ybp (Stroeven et al., 2016). Our somewhat older divergence time estimate for *lapponica* (c. 16K ybp) allows for the possibility that bar-tailed godwits were isolated in tundra breeding areas that were available south or southeast of glaciated Europe during or soon after the LGM. Such refugia, including in central/southern Europe, have been inferred for tundra-breeding shorebirds (Arcones et al., 2021) and other taxa (Hewitt, 2004; Sommer & Nadachowski, 2006; Taberlet et al., 1998). If *lapponica* then drifted northward towards Fennoscandia from central Europe as glaciers retreated, then its present-day migration route would re-trace that northward expansion, as has been hypothesised for numerous migratory birds that either originally evolved from sedentary populations (Rappole & Jones, 2002; Salewski & Bruderer, 2007) or oscillated between short- and long-distance migration as northern areas were impacted by cyclic glaciations (Bruderer & Salewski, 2008; Zink & Gardner, 2017).

At the LGM, the lineage that recently gave rise to *taymyrensis* and *yamalensis* could have existed east or south of the glaciers that extended from Europe to the Taimyr Peninsula of central Russia (Arcones et al., 2021; Batchelor et al., 2019). We may then expect that the present-day breeding range and migration of *yamalensis* (Figure 1) would most closely resemble the state of the western Palearctic lineage at the LGM. Consistent with this, there is evidence of tundra habitats existing south of the ice-dammed West Siberian Lake at the LGM (Binney et al., 2017; Mangerud et al., 2004), close to the present-day breeding area of *yamalensis* in the West Siberian Plain. Currently, *taymyrensis* and *yamalensis* have some degree of isolation in breeding areas along a latitudinal and longitudinal gradient (Figure 1; Bom et al., 2022). This may
suggest that a northeastward expansion of breeding range towards the Taimyr Peninsula was associated with discovery of the migration route of *taymyrensis* through western Europe (Figure 1) after glaciers retreated.

## 4.3 | Population structure in long-distance migrants

The very fine-scale structure we found within the Alaska breeding range of *baueri*, involving birds making the longest non-stop migratory flights on earth, demonstrates the weak association between dispersal *ability* (in terms of flight capacity and apparent strength of barriers) and actual dispersal. In fact, bar-tailed godwits are known for extremely rigid site-fidelity at non-breeding and migratory staging sites (Conklin et al., 2013; Rakhimberdiev et al., 2018), and are presumed to show similar fidelity to breeding sites as adults. Within-baueri structure further suggests that natal dispersal is also quite restricted (i.e. relative to distances of c. 400–800 km between our Alaska sampling sites; Figure 1). As we have no reason to expect weaker site-fidelity in other bar-tailed godwit populations, we presume the shallow differentiation between *taymyrensis* and *yamalensis*, which may meet at pre- and post-breeding areas and breed <500 km from each other (Bom et al., 2022), is a result of their recent divergence (c. 6000 ybp) rather than high levels of current gene flow. Our results support previous research showing that cultural traditions can be a stronger influence on population structure in social vertebrates than is the physiological capacity for movement (Carroll et al., 2015).

While not reliably associated with dispersal and gene flow, mobility of a species may be more indicative of its propensity for exploration and colonisation of new geographical regions. As illustrated in other long-distance migratory shorebirds, such as dunlin, *Calidris alpina* (Wenink et al., 1996) and red knot, *C. canutus* (Conklin et al., 2022), such species appear to respond relatively quickly to major climatic changes, such as retreat of glacial ice sheets, with establishment of new flyway populations (as we have demonstrated with the *lapponica* subspecies). Here, the flight capacity of a species may be directly related to the likelihood of discovering (and then exploiting) such opportunities for colonization. As migration of adult birds tends to be more conservative, such “innovations” in a population may result largely from inexperienced juveniles conducting their first migration (Cresswell, 2014; Verhoeven et al., 2022). Intriguingly, bar-tailed godwits also seem flexible (at least in evolutionary timescales) with regard to migration distance. For example, *lapponica* and *yamalensis*, among the youngest of bar-tailed godwit subspecies, migrate as little as 2500–5000 km between breeding and non-breeding grounds, whereas the other four subspecies travel one-way distances of 10,000–18,000 km (Figure 1). This shows that behavioural and physiological ‘upscaling’ to accommodate the world’s longest flights (Hedenström, 2010; Piersma et al., 2022) is not a one-way journey; it may be reversed,

with regression towards short-distance migration, if the ecological opportunity arises.

## AUTHOR CONTRIBUTIONS

J.R.C., Y.I.V., M.C.F. and T.P. conceived and designed the study. P.F.B., R.A.B., R.E.G., C.J.H., J.t.H., D.R.R., T.L.T., N.W. and P.S.T. organised and performed field sampling and/or maintained individual resight-history databases. J.R.C. and Y.I.V. curated samples and conducted the labwork. J.R.C., M.J.M.L. and M.C.F. analysed and interpreted the data, with assistance from Y.I.V. J.R.C. and M.C.F. wrote the manuscript, with major contributions from Y.I.V. and T.P. All authors provided feedback and helped edit the manuscript.

## ACKNOWLEDGEMENTS

We thank Per Palsbøll and Ritsert Jansen for guidance and help to develop the funding proposal, and to members of the Palsbøll and Fontaine laboratories at the University of Groningen and Montpellier for helpful discussions during project conception and analysis respectively. We thank Anneke Bol at NIOZ Royal Netherlands Institute for Sea Research for DNA extractions of Mauritania samples; Marco van der Velde for laboratory assistance in Groningen; Eric Johnson and Paul Etter (SNPsaurus) for sequencing, bioinformatics, and technical expertise; and the Centre for Information Technology of the University of Groningen (particularly Bob Dröge and Cristian Marocico) for technical support and access to the Peregrine high-performance computing cluster. For providing samples, we thank Massey University, Moscow Lomonosov State University Zoological Museum, NIOZ, Royal Ontario Museum (Oliver Haddrath and Mark Peck), U.S. Geological Survey Alaska Science Center, University of Groningen, and University of Washington Burke Museum. For godwit capture and field sampling, we thank Luke DeCicco, Jimmy de Fouw, Maksim Dementyev, Nick Hajdukovich, Jim Johnson, David Melville, Dan Mulcahy, Murray Potter, Adrian Riegen, Rob Schuckard, Craig Steed, and Cleland Wallace, plus numerous local volunteers. In Australia, we thank Broome Bird Observatory and the Australasian Wader Studies Group for logistical and field work support, and we acknowledge the Yawuru People via the offices of Nyamba Buru Yawuru Limited for permission to capture godwits on the shores of Roebuck Bay, traditional lands of the Yawuru people. Sampling was supported by the David and Lucile Packard Foundation (Alaska, Australia, and New Zealand), BirdLife Netherlands and WWF-Netherlands (Australia), and a Royal Society of New Zealand Marsden Fund grant to P.F.B. and Andrew Fidler (New Zealand). The Prins Bernhard Cultuurfonds grant to T.P. started up the long-term efforts in Mauritania, supported now by operational grants from NIOZ; we thank the authorities of the Parc National du Banc d'Arguin for access and support. Work in the Wadden Sea was supported by Waddenfonds (project Metawad, WF209925) and operational funding from NIOZ. Work in Oman was financially supported by the Research Council (TRC) of the Sultanate of Oman (ORG/EBR/12/002). The core of this project was supported by a Dutch Research Council (NWO) grant to T.P. (ALW-Open Programme grant,
‘Ecological drivers of global flyway evolution’ 824.01.001). Any use of trade, product or firm names is for descriptive purposes only and does not imply endorsement by the US Government.

## CONFLICT OF INTEREST STATEMENT

The authors declare no conflict of interest.

## DATA AVAILABILITY STATEMENT

Demultiplexed nextRAD short-read data with sample metadata were deposited in NCBI's SRA archives under BioProject ID PRJNA1118901 (accession numbers SAMN41630116-SAMN41630324). Associated files (VCF and metadata) were deposited in DRYAD (doi:10.5061/dryad.73n5tb352).

## BENEFIT-SHARING STATEMENT

All genetic samples were collected with collaboration and/or consultation with local communities and under applicable permits, and all principal collaborators who provided samples are included as co-authors. The research addresses an international priority concern, which is the conservation of the study organism. Benefits from this research accrue from the sharing of our data and results with the scientific and stakeholder communities, both directly and on public databases as described above.

## ORCID

Jesse R. Conklin https://orcid.org/0000-0002-1414-0587

Yvonne I. Verkuil https://orcid.org/0000-0001-7080-3026

Margaux J. M. Lefebvre https://orcid.org/0000-0002-9339-8657

Phil F. Battley https://orcid.org/0000-0002-8590-8098

Roeland A. Bom https://orcid.org/0000-0001-8180-1958

Robert E. Gill Jr https://orcid.org/0000-0002-6385-4500

Job ten Horn https://orcid.org/0000-0002-8767-5929

Daniel R. Ruthrauff https://orcid.org/0000-0003-1355-9156

T. Lee Tibbitts https://orcid.org/0000-0002-0290-7592

Pavel S. Tomkovich https://orcid.org/0000-0002-1563-2196

Nils Warnock https://orcid.org/0000-0002-8042-318X

Theunis Piersma https://orcid.org/0000-0001-9668-466X

Michaël C. Fontaine https://orcid.org/0000-0003-1156-4154

## REFERENCES

Åkesson, S., Ilieva, M., Karagicheva, J., Rakhimberdiev, E., Tomotani, B., & Helm, B. (2017). Timing avian long-distance migration: From internal clock mechanisms to global flights. *Philosophical Transactions of the Royal Society, B: Biological Sciences*, 372(1734), 20160252. https://doi.org/10.1098/rstb.2016.0252

Alexander, D. H., & Lange, K. (2011). Enhancements to the ADMIXTURE algorithm for individual ancestry estimation. *BMC Bioinformatics*, 12, 246. https://doi.org/10.1186/1471-2105-12-246

Alexander, D. H., Novembre, J., & Lange, K. (2009). Fast model-based estimation of ancestry in unrelated individuals. *Genome Research*, 19, 1655-1664. https://doi.org/10.1101/gr.094052.109

Anderson, E. C., & Dunham, K. K. (2008). The influence of family groups on inferences made with the program Structure. *Molecular Ecology Resources*, 8(6), 1219-1229. https://doi.org/10.1111/j.1755-0998.2008.02355.x

Archer, F. I., Adams, P. E., & Schneiders, B. B. (2017). *strataG: An R package for manipulating, summarizing and analysing population genetic data*. *Molecular Ecology Resources*, 17(1), 5–11. https://doi.org/10.1111/1755-0998.12559

Arcones, A., Ponti, R., Ferrer, X., & Vieites, D. R. (2021). Pleistocene glacial cycles as drivers of allopatric differentiation in Arctic shorebirds. *Journal of Biogeography*, 48(4), 747–759. https://doi.org/10.1111/jbi.14023

Arguedas, N., & Parker, P. (2000). Seasonal migration and genetic population structure in house wrens. *Condor*, 102(3), 517–528. https://doi.org/10.1093/condor/102.3.517

Avise, J. C. (2009). Phylogeography: Retrospect and prospect. *Journal of Biogeography*, 36(1), 3–15. https://doi.org/10.1111/j.1365-2699.2008.02032.x

Avise, J. C., & Walker, D. (1998). Pleistocene phylogeographic effects on avian populations and the speciation process. *Proceedings of the Royal Society B: Biological Sciences*, 265, 457–463.

Batchelor, C. L., Margold, M., Krapp, M., Murton, D. K., Dalton, A. S., Gibbard, P. L., Stokes, C. R., Murton, J. B., & Manica, A. (2019). The configuration of Northern Hemisphere ice sheets through the quaternary. *Nature Communications*, 10(1), 3713. https://doi.org/10.1038/s41467-019-11601-2

Battley, P. F., Conklin, J. R., Parody-Merino, Á. M., Langlands, P. A., Southey, I., Burns, T., Melville, D. S., Schuckard, R., Riegen, A. C., & Potter, M. A. (2020). Interacting roles of breeding geography and early-life settlement in godwit migration timing. *Frontiers in Ecology and Evolution*, 8, 52. https://doi.org/10.3389/fevo.2020.00052

Battley, P. F., Warnock, N., Tibbitts, T. L., Gill, R. E., Jr., Piersma, T., Hassell, C. J., Douglas, D. C., Mulcahy, D. M., Gartrell, B. D., Schuckard, R., Melville, D. S., & Riegen, A. C. (2012). Contrasting extreme long-distance migration patterns in bar-tailed godwits *Limosa lapponica*. *Journal of Avian Biology*, 43(1), 21–32. https://doi.org/10.1111/j.1600-048X.2011.05473.x

Beaumont, M. A., Zhang, W., & Balding, D. J. (2002). Approximate Bayesian computation in population genetics. *Genetics*, 162(4), 2025–2035. https://doi.org/10.1093/genetics/162.4.2025

Binney, H., Edwards, M., Macias-Fauria, M., Lozhkin, A., Anderson, P., Kaplan, J. O., Andreev, A., Bezrukova, E., Blyakharchuk, T., Jankovska, V., Khazina, I., Krivonogov, S., Kremenetski, K., Nield, J., Novenko, E., Ryabogina, N., Solovieva, N., Willis, K., & Zernitskaya, V. (2017). Vegetation of Eurasia from the last glacial maximum to present: Key biogeographic patterns. *Quaternary Science Reviews*, 157, 80–97. https://doi.org/10.1016/j.quascirev.2016.11.022

Bohonak, A. J. (1999). Dispersal, gene flow, and population structure. *Quarterly Review of Biology*, 74(1), 21–45. https://doi.org/10.1086/392950

Bom, R. A., Conklin, J. R., Verkuil, Y. I., Alves, J. A., de Fouw, J., Dekinga, A., Hassell, C. J., Klaassen, R. H. G., Kwarteng, A. Y., Rakhimberdiev, E., Rocha, A., ten Horn, J., Tibbitts, T. L., Tomkovich, P. S., Victor, R., & Piersma, T. (2022). Central-west Siberian-breeding bar-tailed godwits (*Limosa lapponica*) segregate in two morphologically distinct flyway populations. *Ibis*, 164(2), 468–485. https://doi.org/10.1111/ibi.13024

Bruderer, B., & Salewski, V. (2008). Evolution of bird migration in a biogeographical context. *Journal of Biogeography*, 35(11), 1951–1959. https://doi.org/10.1111/j.1365-2699.2008.01992.x

Bushnell, B. (2016). *BBMap*. University of California. http://sourceforge.net/projects/bbmap

Carroll, E. L., Baker, C. S., Watson, M., Alderman, R., Bannister, J., Gaggiotti, O. E., Gröcke, D. R., Patenaude, N., & Harcourt, R. (2015). Cultural traditions across a migratory network shape the genetic structure of southern right whales around Australia and New Zealand. *Scientific Reports*, 5, 16182. https://doi.org/10.1038/srep16182

Chan, Y. C., Tibbitts, T. L., Dorofeev, D., Hassell, C. J., & Piersma, T. (2022). Hidden in plain sight: Migration routes of the elusive
Anadry bar-tailed godwit revealed by satellite tracking. *Journal of Avian Biology*, 2022(8), e02988. https://doi.org/10.1111/jav.02988

Chang, C. C., Chow, C. C., Tellier, L. C. A. M., Vattikuti, S., Purcell, S. M., & Lee, J. J. (2015). Second-generation PLINK: Rising to the challenge of larger and richer datasets. *GigaScience*, 4(1), 1–16. https://doi.org/10.1186/s13742-015-0047-8

Choi, C. Y., Battley, P. F., Potter, M. A., Rogers, K. G., & Ma, Z. (2015). The importance of Yalujiang coastal wetland in the north Yellow Sea to bar-tailed godwits *Limosa lapponica* and great knots *Calidris tenuirostris* during northward migration. *Bird Conservation International*, 25(1), 53–70. https://doi.org/10.1017/S0959270914000124

Claramunt, S., Derryberry, E. P., Remsen, J. V., & Brumfield, R. T. (2012). High dispersal ability inhibits speciation in a continental radiation of passerine birds. *Proceedings of the Royal Society B: Biological Sciences*, 279(1733), 1567–1574. https://doi.org/10.1098/rspb.2011.1922

Collin, F. D., Durif, G., Raynal, L., Lombaert, E., Gautier, M., Vitalis, R., Marin, J. M., & Estoup, A. (2021). Extending approximate Bayesian computation with supervised machine learning to infer demographic history from genetic polymorphisms using DIYABC random Forest. *Molecular Ecology Resources*, 21(8), 2598–2613. https://doi.org/10.1111/1755-0998.13413

Conklin, J. R., & Battley, P. F. (2011). Contour-feather moult of bar-tailed godwits (*Limosa lapponica baueri*) in New Zealand and the northern hemisphere reveals multiple strategies by sex and breeding region. *Emu*, 111(4), 330–340. https://doi.org/10.1071/MU11011

Conklin, J. R., Battley, P. F., & Potter, M. A. (2013). Absolute consistency: Individual versus population variation in annual-cycle schedules of a long-distance migrant bird. *PLoS One*, 8(1), e54535. https://doi.org/10.1371/journal.pone.0054535

Conklin, J. R., Battley, P. F., Potter, M. A., & Fox, J. W. (2010). Breeding latitude drives individual schedules in a trans-hemispheric migrant bird. *Nature. Communications*, 1(5), ncomms1072. https://doi.org/10.1038/ncomms1072

Conklin, J. R., Battley, P. F., Potter, M. A., & Ruthrauff, D. R. (2011). Geographic variation in morphology of Alaska-breeding bar-tailed godwits (*Limosa lapponica*) is not maintained on their nonbreeding grounds in New Zealand. *Auk*, 128(2), 363–373. https://doi.org/10.1525/auk.2011.10231

Conklin, J. R., Senner, N. R., Battley, P. F., & Piersma, T. (2017). Extreme migration and the individual quality spectrum. *Journal of Avian Biology*, 48(1), 19–36. https://doi.org/10.1111/jav.01316

Conklin, J. R., Verkuil, Y. I., Battley, P. F., Hassell, C. J., ten Horn, J., Johnson, J. A., Tomkovich, P. S., Baker, A. J., Piersma, T., & Fontaine, M. C. (2022). Global flyway evolution in red knots *Calidris canutus* and genetic evidence for a Nearctic refugium. *Molecular Ecology*, 31(7), 2124–2139. https://doi.org/10.1111/mec.16379

Cresswell, W. (2014). Migratory connectivity of Palaearctic-African migratory birds and their responses to environmental change: The serial residency hypothesis. *Ibis*, 156(3), 493–510. https://doi.org/10.1111/ibi.12168

Danecek, P., Auton, A., Abecasis, G., Albers, C. A., Banks, E., DePristo, M. A., Handsaker, R. E., Lunter, G., Marth, G. T., Sherry, S. T., McVean, G., & Durbin, R. (2011). The variant call format and VCFtools. *Bioinformatics*, 27(15), 2156–2158. https://doi.org/10.1093/bioinformatics/btr330

Danecek, P., Bonfield, J. K., Liddle, J., Marshall, J., Ohan, V., Pollard, M. O., Whitwham, A., Keane, T., McCarthy, S. A., & Davies, R. M. (2021). Twelve years of SAMtools and BCFtools. *GigaScience*, 10(2), 1–4. https://doi.org/10.1093/gigascience/giab008

Duijns, S., Jukema, J., Spaans, B., van Horssen, P., & Piersma, T. (2012). Revisiting the proposed leap-frog migration of bar-tailed godwits along the East-Atlantic flyway. *Ardea*, 100(1), 37–43. https://doi.org/10.5253/078.100.0107

Ehlers, J., & Gibbard, P. L. (2007). The extent and chronology of Cenozoic global glaciation. *Quaternary International*, 164–165, 6–20. https://doi.org/10.1016/j.quaint.2006.10.008.

Engelmoer, M., & Roselaar, C. S. (1998). *Geographical variation in waders*. Kluwer Academic Publishers.

Estoup, A., Lombaert, E., Marin, J.-M., Guillemaud, T., Pudlo, P., Robert, C. P., & Cornuet, J.-M. (2012). Estimation of demo-genetic model probabilities with approximate Bayesian computation using linear discriminant analysis on summary statistics. *Molecular Ecology Resources*, 12(5), 846–855. https://doi.org/10.1111/j.1755-0998.2012.03153.x

Farmer, J. R., Pico, T., Underwood, O. M., Cleveland Stout, R., Granger, J., Cronin, T. M., Fripiat, F., Martínez-García, A., Haug, G. H., & Sigman, D. M. (2023). The Bering Strait was flooded 10,000 years before the last glacial maximum. *Proceedings of the National Academy of Sciences*, 120(1), e2206742119. https://doi.org/10.1073/pnas

Gill, R. E., Jr., Tibbitts, T. L., Douglas, D. C., Handel, C. M., Mulcahy, D. M., Gottschalck, J. C., Warnock, N., McCaffery, B. J., Battley, P. F., & Piersma, T. (2009). Extreme endurance flights by landbirds crossing the Pacific Ocean: Ecological corridor rather than barrier? *Proceedings of the Royal Society B: Biological Sciences*, 276(1656), 447–457. https://doi.org/10.1098/rspb.2008.1142

Hawkes, L. A., Balachandran, S., Batbayar, N., Butler, P. J., Frappell, P. B., Milsom, W. K., Tseveenmyadag, N., Newman, S. H., Scott, G. R., Sathiyaselvam, P., Takekawa, J. Y., Wikelski, M., & Bishop, C. M. (2011). The trans-Himalayan flights of bar-headed geese (*Anser indicus*). *Proceedings of the National Academy of Sciences*, 108(23), 9516–9519. https://doi.org/10.1073/pnas.1017295108

Hedenström, A. (2010). Extreme endurance migration: What is the limit to non-stop flight? *PLoS Biology*, 8(5), 1–6. https://doi.org/10.1371/journal.pbio.1000362

Hendry, A. P., & Day, T. (2005). Population structure attributable to reproductive time: Isolation by time and adaptation by time. *Molecular Ecology*, 14(4), 901–916. https://doi.org/10.1111/j.1365-294X.2005.02480.x

Hewitt, G. M. (2000). The genetic legacy of the Quaternary ice ages. *Nature*, 405, 907–913. https://doi.org/10.1038/35016000

Hewitt, G. M. (2004). Genetic consequences of climatic oscillations in the Quaternary. *Philosophical Transactions of the Royal Society, B: Biological Sciences*, 359(1442), 183–195. https://doi.org/10.1098/rstb.2003.1388

Hoffecker, J. F., Pitulko, V. V., & Pavlova, E. Y. (2023). Beringia and the settlement of the Western Hemisphere. *Proceedings of the Royal Society of London, Series B: Biological Sciences*, 67(3), 882–909. https://doi.org/10.21638/SPBU02.2022.313

Kamvar, Z. N., Tabima, J. F., & Grünwald, N. J. (2014). Poppr: An R package for genetic analysis of populations with clonal, partially clonal, and/or sexual reproduction. *PeerJ*, 2, e281. https://doi.org/10.7717/peerj.281

Karl, S. A., Castro, A. L. F., Lopez, J. A., Charvet, P., & Burgess, G. H. (2011). Phylogeography and conservation of the bull shark (*Carcharhinus leucas*) inferred from mitochondrial and microsatellite DNA. *Conservation Genetics*, 12(2), 371–382. https://doi.org/10.1007/s10592-010-0145-1

Keenan, K., McGinnity, P., Cross, T. F., Crozier, W. W., & Prodöhl, P. A. (2013). DiveRsity: An R package for the estimation and exploration of population genetics parameters and their associated errors. *Methods in Ecology and Evolution*, 4(8), 782–788. https://doi.org/10.1111/2041-210X.12067

Knowles, L. L., & Maddison, W. P. (2002). Statistical phylogeography. *Molecular Ecology*, 11(12), 2623–2635. https://doi.org/10.1046/j.1365-294X.2002.01637.x

Knutsen, H., Jorde, P. E., André, C., & Stenseth, N. C. (2003). Fine-scaled geographical population structuring in a highly mobile marine species: The Atlantic cod. *Molecular Ecology*, 12(2), 385–394. https://doi.org/10.1046/j.1365-294X.2003.01750.x

Kopelman, N. M., Mayzel, J., Jakobsson, M., Rosenberg, N. A., & Mayrose, I. (2015). Clumpak: A program for identifying clustering modes and packaging population structure inferences across *K. Molecular*
Ecology Resources, 15(5), 1179-1191. https://doi.org/10.1111/
1755-0998.12387

Kraaijeveld, K., & Nieboer, E. N. (2000). Late Quaternary paleogeography and evolution of arctic breeding waders. *Ardea, 88*(2), 193-205.

Leppälä, K., Nielsen, S. V., & Mailund, T. (2017). Admixturegraph: An R package for admixture graph manipulation and fitting. *Bioinformatics*, 33(11), 1738-1740. https://doi.org/10.1093/bioin
formatics/btx048

Lischer, H. E. L., & Excoffier, L. (2012). PGDSpider: An automated data conversion tool for connecting population genetics and genomics programs. *Bioinformatics*, 28(2), 298-299. https://doi.org/10.1093/
bioinformatics/btr642

Maier, R., Flegontov, P., Flegontova, O., Işıldak, U., & Changmai, P. (2023). On the limits of fitting complex models of population history to f-statistics. *eLife*, 12, e85492.

Mangerud, J., Jakobsson, M., Alexanderson, H., Astakhov, V., Clarke, G. K. C., Henriksen, M., Hjort, C., Krinner, G., Lunkka, J. P., Möller, P., Murray, A., Nikolskaya, O., Saarnisto, M., & Svendsen, J. I. (2004). Ice-dammed lakes and rerouting of the drainage of northern Eurasia during the last glaciation. *Quaternary Science Reviews*, 23(11-13), 1313-1332. https://doi.org/10.1016/j.quascirev.2003.12.009

Medina, I., Cooke, G. M., & Ord, T. J. (2018). Walk, swim or fly? Locomotor mode predicts genetic differentiation in vertebrates. *Ecology Letters*, 21(5), 638–645. https://doi.org/10.1111/ele.12930

Méndez, V., Alves, J. A., Gill, J. A., & Gunnarsson, T. G. (2018). Patterns and processes in shorebird survival rates: A global review. *Ibis*, 160(4), 723–741. https://doi.org/10.1111/ibi.12586

Milá, B., Smith, T. B., & Wayne, R. K. (2006). Postglacial population expansion drives the evolution of long-distance migration in a song-bird. *Evolution*, 60(11), 2403–2409. https://doi.org/10.1111/j.0014-3820.2006.tb01875.x

Moussy, C., Hosken, D. J., Mathews, F., Smith, G. C., Aegerter, J. N., & Bearhop, S. (2013). Migration and dispersal patterns of bats and their influence on genetic structure. *Mammal Review*, 43(3), 183–195. https://doi.org/10.1111/j.1365-2907.2012.00218.x

Nei, M. (1972). Genetic distance between populations. *American Naturalist*, 106(949), 283-292. https://doi.org/10.1086/282771

Nosil, P., Vines, T. H., & Funk, D. J. (2005). Perspective: Reproductive isolation caused by natural selection against immigrants from divergent habitats. *Evolution*, 59(4), 705–719. https://doi.org/10.1054/04-428

Palumbi, S. R. (1994). Genetic divergence, reproductive isolation, and marine speciation. *Annual Review of Ecology and Systematics*, 25(1), 547-572. https://doi.org/10.1146/annurev.es.25.110194.002555

Paradis, E., & Schliep, K. (2019). Ape 5.0: An environment for modern phylogenetics and evolutionary analyses in R. *Bioinformatics*, 35(3), 526-528. https://doi.org/10.1093/bioinformatics/bty633

Parody-Merino, Á. M., Battley, P. F., Conklin, J. R., & Fidler, A. E. (2019). No evidence for an association between *clock* gene allelic variation and migration timing in a long-distance migratory shorebird (*Limosa lapponica baueri*). *Oecologia*, 191(4), 843–859. https://doi.org/10.1007/s00442-019-04524-8

Patterson, N., Moorjani, P., Luo, Y., Mallick, S., Rohland, N., Zhan, Y., Genschoreck, T., Webster, T., & Reich, D. (2012). Ancient admixture in human history. *Genetics*, 192(3), 1065–1093. https://doi.org/10.1534/genetics.112.145037

Patterson, N., Price, A. L., & Reich, D. (2006). Population structure and eigenanalysis. *PLoS Genetics*, 2(12), 2074–2093. https://doi.org/10.1371/journal.pgen.0020190

Piersma, T., Gill, R. E., Jr, Ruthrauff, D. R., Guglielmo, C. G., Conklin, J. R., & Handel, C. M. (2022). The Pacific as the world's greatest theater of bird migration: Extreme flights spark questions about physiological capabilities, behavior, and the evolution of migratory pathways. *Ornithology*, 139(2), ukab086. https://doi.org/10.1093/ornithology/ukab086

Pielou, E. C. (1991). After the Ice Age: The return of life to glaciated North America. University of Chicago Press.

Piersma, T., Pérez-Tris, J., Mouritsen, H., Bauchinger, U., & Bairlein, F. (2005). Is there a "migratory syndrome" common to all migrant birds? *Annals of the New York Academy of Sciences*, 1046, 282-293. https://doi.org/10.1196/annals.1343.026

Pudlo, P., Marin, J. M., Estoup, A., Cornuet, J. M., Gautier, M., & Robert, C. P. (2016). Reliable ABC model choice via random forests. *Bioinformatics*, 32(6), 859-866. https://doi.org/10.1093/bioinformatics/btv684

R Core Team. (2021). *R: A language and environment for statistical computing*. R Foundation for Statistical Computing. https://www.r-project.org/

Rakhimberdiev, E., Duijns, S., Karagicheva, J., Camphuysen, C. J., V. R. S. Castricum., Dekinga, A., Dekker, R., Gavrilov, A., ten Horn, J., Jukema, J., Saveliev, A., Soloviev, M., Tibbitts, T. L., van Gils, J. A., & Piersma, T. (2018). Fuelling conditions at staging sites can mitigate Arctic warming effects in a migratory bird. *Nature Communications*, 9(1), 4263. https://doi.org/10.1038/s41467-018-06673-5

Rappole, J. H., & Jones, P. (2002). Evolution of old and new world migration systems. *Ardea*, 90(3), 525-537.

Raynal, L., Marin, J. M., Pudlo, P., Ribatet, M., Robert, C. P., & Estoup, A. (2019). ABC random forests for Bayesian parameter inference. *Bioinformatics*, 35(10), 1720-1728. https://doi.org/10.1093/bioin
formatics/bty867

Reich, D., Thangaraj, K., Patterson, N., Price, A. L., & Singh, L. (2009). Reconstructing Indian population history. *Nature*, 461(7263), 489-494. https://doi.org/10.1038/nature08365

Rodríguez-Ramilo, S. T., & Wang, J. (2012). The effect of close relatives on unsupervised Bayesian clustering algorithms in population genetic structure analysis. *Molecular Ecology Resources*, 12(5), 873-884. https://doi.org/10.1111/j.1755-0998.2012.03156.x

Rolland, J., Jiguet, F., Jönsson, K. A., Condamine, F. L., & Morlon, H. (2014). Settling down of seasonal migrants promotes bird diversification. *Proceedings of the Royal Society B: Biological Sciences*, 281(1784), 201404073. https://doi.org/10.1098/rspb.2014.0473

Rosenbaum, H. C., Pomilla, C., Mendez, M., Leslie, M. S., Best, P. B., Findlay, K. P., Minton, G., Ersts, P. J., Collins, T., Engel, M. H., Bonatto, S. L., Kotze, D. P. G. H., Meyer, M., Barendse, J., Thornton, M., Razafindrakoto, Y., Ngouessono, S., Vely, M., & Kiszka, J. (2009). Population structure of humpback whales from their breeding grounds in the South Atlantic and Indian oceans. *PLoS One*, 4(10), e7318. https://doi.org/10.1371/journal.pone.0007318

Russello, M. A., Waterhouse, M. D., Etter, P. D., & Johnson, E. A. (2015). From promise to practice: Pairing non-invasive sampling with genetics in conservation. *PeerJ*, 3, e1106. https://doi.org/10.7717/peerj.1106

Salewski, V., & Bruderer, B. (2007). The evolution of bird migration - A synthesis. *Naturwissenschaften*, 94(4), 268-279. https://doi.org/10.1007/s00114-006-0186-y

Schmaljohann, H., Liechi, F., & Bruderer, B. (2007). Songbird migration across the Sahara: The non-stop hypothesis rejected! *Proceedings of the Royal Society B: Biological Sciences*, 274(1610), 735-739. https://doi.org/10.1098/rspb.2006.0011

Schrider, D. R., & Kern, A. D. (2018). Supervised machine learning for population genetics: A new paradigm. *Trends in Genetics*, 34(4), 301-312. https://doi.org/10.1016/j.tig.2017.12.005

Sellas, A. B., Wells, R. S., & Rosel, P. E. (2005). Mitochondrial and nuclear DNA analyses reveal fine scale geographic structure in bottlenose dolphins (*Tursiops truncatus*) in the Gulf of Mexico. *Conservation Genetics*, 6(5), 715-728. https://doi.org/10.1007/s10592-005-9031-7

Sommer, R. S., & Nadachowski, A. (2006). Glacial refugia of mammals in Europe: Evidence from fossil records. *Mammal Review*, 36(4), 251-265. https://doi.org/10.1111/j.1365-2907.2006.00093.x

Stewart, J. R., & Dalén, L. (2008). Is the glacial refugium concept relevant for northern species? A comment on Pruett and Winker 2005. *Climatic Change*, 86, 19-22. https://doi.org/10.1007/s10584-007-9366-9
Stroeven, A. P., Hättestrand, K., Kleman, J., Heyman, J., Fabel, D., Fredin, O., Goodfellow, B. W., Harbor, J. M., Jansen, J. D., Olsen, L., Caffee, M. W., Fink, D., Lundqvist, J., Rosqvist, G. C., Strömberg, B., & Jansson, K. N. (2016). Deglaciation of Fennoscandia. *Quaternary Science Reviews*, 147, 91–121. https://doi.org/10.1016/j.quascirev.2015.09.016

Taberlet, P., Fumagalli, L., Wust-Saucy, A. G., & Cosson, J. F. (1998). Comparative phylogeography and postglacial colonization routes in Europe. *Molecular Ecology*, 7(4), 453–464. https://doi.org/10.1046/j.1365-294x.1998.00289.x

Tan, H. Z., Jansen, J. J. F. J., Allport, G. A., Garg, K. M., Chattopadhyay, B., Irestedt, M., Pang, S. E. H., Chilton, G., Gwee, C. Y., & Rheindt, F. E. (2023). Megafaunal extinctions, not climate change, may explain Holocene genetic diversity declines in *Numenius* shorebirds. *eLife*, 12, e85422.

Theisen, T. C., Bowen, B. W., Lanier, W., & Baldwin, J. D. (2008). High connectivity on a global scale in the pelagic wahoo, *Acanthocybium solandri* (tuna family Scombridae). *Molecular Ecology*, 17(19), 4233–4247. https://doi.org/10.1111/j.1365-294X.2008.03913.x

Tomkovich, P. S. (2010). Assessment of the Anadyr lowland subspecies of bar-tailed godwit *Limosa lapponica anadyrensis*. *Bulletin of the British Ornithologists' Club*, 130(2), 88–95.

van Roomen, M., Citegetse, G., Crowe, O., Dodman, T., Hagermeijer, W., Meise, K., & Schekkerman, H. (2022). *East Atlantic Flyway assessment 2020*. https://stats.sovon.nl/pub/publicatie/20593

Verhoeven, M. A., Loonstra, A. H. J., McBride, A. D., Kaspersma, W., Hooijmeijer, J. C. E. W., Both, C., Senner, N. R., & Piersma, T. (2022). Age-dependent timing and routes demonstrate developmental plasticity in a long-distance migratory bird. *Journal of Animal Ecology*, 91(3), 566–579. https://doi.org/10.1111/1365-2656.13641

Verkuil, Y. I., Piersma, T., Jukema, J., Hooijmeijer, J. C. E. W., Zwarts, L., & Baker, A. J. (2012). The interplay between habitat availability and population differentiation: A case study on genetic and morphological structure in an inland wader (Charadriiformes). *Biological Journal of the Linnean Society*, 106(3), 641–656. https://doi.org/10.1111/j.1095-8312.2012.01878.x

Wang, Y., Heintzman, P. D., Newsom, L., Bigelow, N. H., Wooller, M. J., Shapiro, B., & Williams, J. W. (2017). The southern coastal Beringian land bridge: Cryptic refugium or pseudorefugium for woody plants during the Last Glacial Maximum? *Journal of Biogeography*, 44(7), 1559–1571. https://doi.org/10.1111/jbi.13010

Wauchope, H. S., Shaw, J. D., Varpe, Ø., Lappo, E. G., Boertmann, D., Lanctot, R. B., & Fuller, R. A. (2017). Rapid climate-driven loss of breeding habitat for Arctic migratory birds. *Global Change Biology*, 23(3), 1085–1094. https://doi.org/10.1111/gcb.13404

Weir, B. S., & Cockerham, C. C. (1984). Estimating F-statistics for the analysis of population structure. *Evolution*, 38(6), 1358–1370. https://doi.org/10.2307/2409936

Wen, J., Nie, Z. L., & Ickert-Bond, S. M. (2016). Intercontinental disjunctions between eastern Asia and western North America in vascular plants highlight the biogeographic importance of the Bering land

bridge from late Cretaceous to Neogene. *Journal of Systematics and Evolution*, 54(5), 469–490. https://doi.org/10.1111/jse.12222

Wenink, P. W., Baker, A. J., Rösner, H. U., & Tilanus, M. G. (1996). Global mitochondrial DNA phylogeography of holarctic breeding dunlins (*Calidris alpina*). *Evolution*, 50(1), 318–330. https://doi.org/10.1111/j.1558-5646.1996.tb04495.x

Wilson, J. R., Nebel, S., & Minton, C. D. T. (2007). Migration ecology and morphometrics of two bar-tailed godwit populations in Australia. *Emu*, 107(4), 262–274. https://doi.org/10.1071/MU07026

Winger, B. M., Auteri, G. G., Pegan, T. M., & Weeks, B. C. (2019). A long winter for the Red Queen: Rethinking the evolution of seasonal migration. *Biological Reviews*, 94(3), 737–752. https://doi.org/10.1111/brv.12476

Winger, B. M., Barker, F. K., & Ree, R. H. (2014). Temperate origins of long-distance seasonal migration in New World songbirds. *Proceedings of the National Academy of Sciences*, 111(33), 12115–12120. https://doi.org/10.1073/pnas.1405000111

Winker, K., Withrow, J. J., Gibson, D. D., & Pruett, C. L. (2023). Beringia as a high-latitude engine of avian speciation. *Biological Reviews*, 98, 1081–1099. https://doi.org/10.1111/brv.12945

Zhang, G., Li, C., Li, Q., Li, B., Larkin, D. M., Lee, C., Storz, J. F., Antunes, A., Greenwold, M. J., Meredith, R. W., Zeng, Y., Xiong, Z., Liu, S., Zhou, L., Huang, Z., An, N., Wang, J., Zheng, Q., Xiong, Y., ... Wang, J. (2014). Comparative genomics reveals insights into avian genome evolution and adaptation. *Science*, 346(6215), 1311–1321.

Zheng, X., Levine, D., Shen, J., Gogarten, S. M., Laurie, C., & Weir, B. S. (2012). A high-performance computing toolset for relatedness and principal component analysis of SNP data. *Bioinformatics*, 28(24), 3326–3328. https://doi.org/10.1093/bioinformatics/bts606

Zink, R. M., & Gardner, A. S. (2017). Glaciation as a migratory switch. *Science Advances*, 3(9), e1603133. https://doi.org/10.1126/sciadv.1603133

## SUPPORTING INFORMATION

Additional supporting information can be found online in the Supporting Information section at the end of this article.

**How to cite this article:** Conklin, J. R., Verkuil, Y. I., Lefebvre, M. J. M., Battley, P. F., Bom, R. A., Gill, R. E. Jr, Hassell, C. J., ten Horn, J., Ruthrauff, D. R., Tibbitts, T. L., Tomkovich, P. S., Warnock, N., Piersma, T., & Fontaine, M. C. (2024). High dispersal ability versus migratory traditions: Fine-scale population structure and post-glacial colonisation in bar-tailed godwits. *Molecular Ecology*, 33, e17452. https://doi.org/10.1111/mec.17452