Intercontinental genetic structure and gene flow in Dunlin (Calidris alpina), a potential vector of avian influenza

Waterfowl (Anseriformes) and shorebirds (Charadriiformes) are the most common wild vectors of influenza A viruses. Due to their migratory behavior, some may transmit disease over long distances. Migratory connectivity studies can link breeding and nonbreeding grounds while illustrating potential interactions among populations that may spread diseases. We investigated Dunlin (Calidris alpina), a shorebird with a subspecies (C. a. arcticola) that migrates from nonbreeding areas endemic to avian influenza in eastern Asia to breeding grounds in northern Alaska. Using microsatellites and mitochondrial DNA, we illustrate genetic structure among six subspecies: C. a. arcticola,C. a. pacifica,C. a. hudsonia,C. a. sakhalina,C. a. kistchinski, and C. a. actites. We demonstrate that mitochondrial DNA can help distinguish C. a. arcticola on the Asian nonbreeding grounds with >70% accuracy depending on their relative abundance, indicating that genetics can help determine whether C. a. arcticola occurs where they may be exposed to highly pathogenic avian influenza (HPAI) during outbreaks. Our data reveal asymmetric intercontinental gene flow, with some C. a. arcticola short-stopping migration to breed with C. a. pacifica in western Alaska. Because C. a. pacifica migrates along the Pacific Coast of North America, interactions between these subspecies and other taxa provide route for transmission of HPAI into other parts of North America.

migratory behavior of many virus hosts (Kilpatrick et al. 2006). Although most migratory movements occur within continents, intercontinental migration can also occur. For example, up to three million birds and thousands of infected individuals cross the Bering Strait from Asia into Alaska each year (Winker and Gibson 2010).
The likelihood that an individual bird species may contribute to the intercontinental spread of avian influenza depends in part on the details of its seasonal migratory patterns. Thus, migratory connectivity studies of birds can be used to define important migratory pathways and identify the population of origin of individuals at all stages of the annual cycle (Webster et al. 2002). Such studies take on new importance in the age of widespread disease transfer by birds (e.g., Rappole et al. 2000;Ishiguro et al. 2005;Morshed et al. 2005;Fergus et al. 2006;Gilbert et al. 2006;Dusek et al. 2014). If the identity and origin of avian disease carriers can be determined and if their migratory pathways are understood, it may be possible to predict the next occurrence of a virulent disease near human population centers, implement precautionary measures to limit human-bird contact, and adopt practices to try to minimize the potential for further spread of the disease to other geographic regions.
The Dunlin (Calidris alpina) is a circumpolar migratory shorebird that breeds throughout arctic and subarctic tundra regions and winters in the southern portion of the Northern Hemisphere (Del Hoyo et al. 1996). There are up to 11 described subspecies that show varying degrees of morphological variation (Greenwood 1986;Tomkovich 1986;Nechaev and Tomkovich 1987;Browning 1991;AOU 2013). These purported subspecies are believed to use separate breeding grounds, but their migratory flyways and nonbreeding areas may overlap (Warnock and Gill 1996;Lappo et al. 2012;Gill et al. 2013). Five subspecies of Dunlin breed in the East Asia and Alaska region known as Beringia ( Fig. 1): Calidris alpina actites, C. a. kistchinski, and C. a. sakhalina breed in the Russian Far East while Figure 1 Breeding distribution of six subspecies of Dunlin (Calidris alpina) sampled for genetic analysis. Sites codes are congruent with those listed in Table 1.

Dunlin genetic structure
Miller et al.
C. a. arcticola and C. a. pacifica breed in Alaska (Warnock and Gill 1996;AOU 2013;Fig. 1). A sixth North American subspecies exists (C. a. hudsonia), but breeds in central and eastern Canada and winters along the Atlantic Coast and Gulf of Mexico (Fern andez et al. 2008). Potential interactions among the Beringia subspecies are complex: C. a. pacifica breeds in western Alaska and migrates south along the Pacific Coast of North America to winter in the western United States and Mexico (Fern andez et al. 2008;Gill et al. 2013). Calidris alpina arcticola breeds in northern Alaska, but migrates across the Bering Strait to winter along the Pacific Coast of Asia where it potentially intermixes with the three East Asia subspecies (Fern andez et al. 2008;Lanctot et al. 2009;Gill et al. 2013). Dunlin were ranked the second-highest of 26 priority taxa to be routinely monitored for HPAI in Alaska when extensive sampling was initiated during the H5N1 HPAI outbreak in 2006 (U.S. Fish and Wildlife Service and U.S. Geological Survey 2007). The rankings were based on each taxon's distribution in Asia, proximity to locations where HPAI has been previously identified, general habitat use patterns, ease of sampling, and population size in Alaska (Alaska Interagency HPAI Bird Surveillance Working Group 2006; Ip et al. 2008). Dunlin ranked high primarily because they winter in areas where outbreaks of HPAI occur in Asia and because so many individuals (300 000-700 000 birds; Andres et al. 2012) migrate from Asia to Alaska each year. Dunlin are also highly susceptible to HPAI H5N1 (Hall et al. 2011). Mortality is likely common among infected juveniles (Hall et al. 2011), but infected adults may survive and transmit viruses. Surveys of wild-caught Dunlin in Alaska between 2006 and 2007 revealed that 0.22% were positive for AI based on RT-PCR analyses of cloacal swabs or fecal samples (Ip et al. 2008), indicating that active shedding of AI viruses was occurring at the time of sampling. This value likely underestimates the true infection rate, as Hall et al. (2011) found that RT-PCR detection of H5N1 in experimental challenges was longer lasting and more consistent from oropharyngeal samples as opposed to cloacal samples. Furthermore, Pearce et al. (2012) found that 2.6% of Dunlin sampled in Alaska during the late summer of 2010 demonstrated evidence for prior AI exposure based on serologic assays. While actual numbers are likely to vary substantially from year to year based on the dynamics of viral outbreaks in Asia, these studies nominally suggest that between 1540 and 18 200 (based on estimated population sizes) infected Dunlin could be in Alaska in any given year. Collectively, this information indicates that C. a. arcticola is an important subspecies to consider when evaluating potential routes and mechanisms by which Asian influenza strains can be transmitted to North America.
Although all Dunlin subspecies show some phenotypic variation (Tomkovich 1986;Nechaev and Tomkovich 1987;Browning 1991), it is difficult to separate them outside of the breeding grounds using commonly employed morphological characters such as plumage or culmen, head, wing, and tarsus measurements (Warnock and Gill 1996;Wennerberg et al. 1999; but see Gates et al. 2013). This is particularly true in eastern Asia where four subspecies are thought to intermix during the nonbreeding season (Lanctot et al. 2009;Gates et al. 2013). In these circumstances, hypervariable molecular markers and DNA sequences may be useful for illuminating patterns of population connectivity and movements of individuals throughout the annual cycle . We therefore used mitochondrial DNA sequences (mtDNA) from the cytochrome b gene and control region along with eight nuclear microsatellite loci to address multiple questions associated with the differentiation of Dunlin subspecies and the extent of gene flow and interactions among groups from Asia and North America.
(i) Do genetic data provide evidence for differentiation among Dunlin subspecies and breeding populations from the region? While prior work has examined phylogeographic patterns in the Northern Hemisphere, most studies were based on small sample sizes and had limited (or no) sampling within Beringia-associated subspecies (e.g., Wennerberg et al. 1999;Marthinsen et al. 2007). (ii) Does genetic differentiation among subspecies provide a basis for the probabilistic identification of subspecies where they co-occur (sensu Patten and Unitt 2002)? In particular, we were interested in determining whether genetic data can distinguish C. a. arcticola from the three other Dunlin subspecies that winter in East Asia (C. a. actites, C. a. kistchinski, and C. a. sakhalina). If distinguishable, then nonbreeding populations of C. a. arcticola could be more easily identified, leading to a better understanding of the likelihood of this subspecies becoming infected and transmitting HPAI into North America. (iii) Can genetic data characterize the extent of gene flow and interaction among the three proximate Beringian subspecies (C. a. sakhalina, C. a. arcticola, and C. a. pacifica)? Given the geographic locations of their breeding ranges (Fig. 1), opportunities for gene flow among subspecies may occur. Furthermore, a portion of the C. a. arcticola and C. a. pacifica populations intermix during postbreeding staging in western Alaska (Gill et al. 2013), but the extent of gene flow between these groups is not well known. If gene flow is extensive, then the data may point to greater-than-expected interactions between these two subspecies. Because C. a. pacifica winters along the Pacific Coast of North America, interactions with C. a. arcticola during the breeding or postbreeding season may increase the risk of transmission of Asian influenza strains from Alaska into other parts of North America.

Sample collection and molecular methods
We collected 370 Dunlin blood or tissue samples from 18 breeding areas during the 2003 to 2009 breeding seasons ( Fig. 1, Table 1). Samples included putative representatives from the five subspecies that inhabit eastern Asia and Alaska (C. a. actites, C. a. kistchinski, C. a. sakhalina, C. a. pacifica, and C. a. arcticola) and were our primary focus for this study. However, we also included samples from three C. a. hudsonia breeding populations in eastern North America to help provide greater genetic and spatial context to our analyses. Individual birds were captured with bownets at nest sites (most subspecies) or lethally collected (C. a. kistchinski samples) on breeding territories. Live-captured birds had up to 0.3 mL of blood collected into a heparinized tube via brachial puncture with a 26-to 27.5-gauge needle. Additional breeding season tissues were obtained from the University of Washington Burke Museum to augment the Russian populations (UWBM Accession Numbers 43910, 44120, 44121, 51684, 51687, 51693, 51694, 51695, and 69903). Blood or tissue samples were preserved in Longmire buffer (Longmire et al. 1997) until used for genetic analyses. DNA was extracted as described in Haig et al. (2004). We used polymerase chain reaction (PCR) to amplify partial sequences of the mitochondrial cytochrome b gene (cyt b) and control region (D-loop) in 234 samples (Table 1). Primer pairs, including L14996-H15646 (http://people. bu.edu/msoren/primers.html, accessed 15 January 2015) and TS96L-TS778H (Wenink et al. 1994), were used to amplify the mitochondrial cyt b and D-loop sequences, respectively. All primer sequences and annealing temperatures are shown in Appendix A. PCR amplifications were performed in 20 lL reactions containing 2.5 mM MgCl 2 , 1 lM of primers, 100 lM of each dNTP, 19 PCR buffer (Perkin Elmer, Waltham, MA, USA), and 1 U AmpliTaq Gold DNA polymerase (Perkin Elmer). Thermal-cycling parameters included initial denaturation at 94°C followed by 35 cycles of denaturing at 94°C (30 s), the annealing temperature listed in Appendix A (30 s), and extension at 72°C (60 s). PCR products were bidirectionally sequenced with BigDye â Terminator 3.1 Cycle Sequencing chemistry (Life Technologies, Grand Island, NY, USA) and resolved on an ABI 3730 automated DNA sequencer, with resulting chromatograms aligned, edited, and trimmed using the program SeqMan ver. 8.0.2 (DNAStar Inc., Madison, WI, USA). The final 1112-bp alignment contained concatenated sequences from each individual and included 633 bp of cyt b and 479 bp from the D-loop.
Nuclear microsatellite genotypes were obtained at eight loci for 370 individuals (Table 1; Appendix A). We obtained primers for loci CALP2 and 4A11 from Wennerberg (2001a), and for loci Cme2, Cme10, and Cme12 from van Treuren et al. (1999), whereas loci D25, D26, and D110 were characterized de novo for this specific investigation during an Illumina GAIIx Genome Analyzer paired-end 80 run (sensu Jennings et al. 2011). Library construction followed recommended Illumina protocols with the exception that index sequencing 'bar-coded' adapters (Craig et al. 2008;Cronn et al. 2008) were substituted for standard paired-end adapters. Primer sequences and annealing temperatures for all microsatellites are provided in Appendix A. PCRs were performed in a 10 lL reaction volume with the following reagent concentrations: 19 PCR buffer (Promega Inc., Madison, WI, USA), 0.5 lM of each primer, 2.5 mM MgCl 2 , 100 lM of each dNTP, and 1 U Taq DNA polymerase (Promega, Inc.). Thermal-cycling parameters included 2 min denaturation at 93°C, followed by 30 cycles of 30 s at 93°C, 30 s at the appropriate annealing temperature, and elongation at 72°C for 1 min. Amplification products were analyzed on an ABI 3100 capillary DNA automated sequencer. ABI GENESCAN software was used to size fragments based on internal lane standard GeneScan 500 [Rox]. ABI GENEMAPPER software was used to score alleles sizes.

Differentiation among subspecies
We characterized the mitochondrial and microsatellite data to provide heuristic indicators of differences among subspecies. For the mtDNA data, we used FaBox (Villesen 2007) to identify unique haplotypes in the data set and create tables reflecting haplotype frequencies and shared haplotypes among groups. ARLEQUIN version 3.1 (Excoffier et al. 2005) was used to quantify gene diversity (H) and nucleotide diversity (p) in mtDNA data within each subspecies. Tables documenting microsatellite allele frequency variation among subspecies were created using CONVERT (Glaubitz 2004). Likewise, program GDA version 1.1 (Lewis and Zaykin 2002) was used to calculate allelic richness and observed and expected heterozygosity (H O and H E , respectively). HP-Rare (Kalinowski 2005) was used to obtain rarefied estimates of allelic richness that accounted for differences in sample size.
We used phylogenetic analyses to examine differentiation of subspecies based on the mtDNA data. The program PhyML 3.0 (Guindon et al. 2010) was used to infer phylogenetic relationships among haplotypes using the maximum-likelihood (ML) criterion. The best-fit nucleotide substitution model was identified using jModeltest2 (Darriba et al. 2012). One thousand bootstrap replicates were used to evaluate clade support. Bayesian phylogenetic analyses were performed using MRBAYES version 3.1.2 (Huelsenbeck and Ronquist 2001), where four concurrent chains were run for 6 9 10 6 generations. Trees were sampled every 2000 generations and 'burn in' included the initial 25% of samples. jModeltest 2 was also used to identify nucleotide substitution models for Bayesian analyses, but was restricted to the subset of models supported by MRBA-YES when performing model selection. Resulting phylogenetic trees from both analyses were visualized and annotated using MEGA 5.2 (Tamura et al. 2011).
We used STRUCTURE version 2.2.3 (Pritchard et al. 2000) to analyze the microsatellite data to identify the number of genetic clusters and to probabilistically assign each analyzed individual to one of the identified clusters. Analyses assumed numbers of clusters (K) ranging from one through seven and were based on the uncorrelated allele frequency model and no admixture. Ten replicate analyses were performed for each value of K with each replicate using an initial 10 6 burn-in steps followed by 10 7 analysis replicates. We evaluated the outcome of analyses in two different ways: by identifying the value of K that produced the highest average likelihood score over replicates and through the use of the DK procedure of Evanno et al. (2005). In both cases, results were summarized over replicates using the program CLUMPP (Jakobsson and Rosenberg 2007). Prior to all microsatellite analyses, we used GDA version 1.1 (Lewis and Zaykin 2002) to identify deviations from Hardy-Weinberg genotypic proportions and test for linkage disequilibrium between pairs of loci within each subspecies. Composite test results for Hardy-Weinberg disequilibrium within each subspecies were obtained by combining P-values from locus-specific analyses using the Z-transform test (Whitlock 2005).
ARLEQUIN was used to perform an analysis of molecular variance (AMOVA; Excoffier et al. 1992) and quantify genetic structure among Dunlin subspecies. In this analysis, Φst (for mtDNA), F ST , and R ST (both for microsatellite data, the latter assuming a strict stepwise mutation; Slatkin 1995) were calculated to determine the overall and pairwise levels of differentiation among different subspecies. P-values associated with these statistics were obtained using 10 000 randomization replicates.
Distinguishing C. a. arcticola from other subspecies that winter in Asia Results from STRUCTURE analyses (described above) were further evaluated to determine whether the microsatellite data could be used to probabilistically distinguish among Dunlin subspecies that winter in Asia. If STRUC-TURE identified more than one cluster, then assignment values for individuals within each cluster may facilitate accurate subspecific diagnoses of individual birds from mixed groups on the nonbreeding grounds. We also used the individual assignment approach encapsulated in GeneClass2 (Piry et al. 2004), where we determined whether birds could be assigned to one of the predefined Dunlin subspecies with a high degree of confidence. Analyses used the Bayesian computation criterion of Rannala and Mountain (1997) and probability computations as described in Cornuet et al. (1999) using 10 000 simulated individuals. After analyses, we determined the proportion of individuals that were correctly reassigned to their respective subspecies and the average probability associated with correct assignments.
The diagnostic utility of the mtDNA data was also evaluated. Results from the phylogenetic analyses initially suggested that our mtDNA could be used to distinguish C. a. arcticola from other subspecies that winter in Asia (see Results and Discussion). Specifically, haplotypes from C. a. arcticola and C. a. pacifica (hereafter referred to as clade I haplotypes) formed a clade that was largely distinct from haplotypes detected in the Asian subspecies C. a. kistchinski, C. a. sakhalina, and C. a. actites (see Results and Fig. 2). The sole exception to this pattern was the detection of seven C. a. sakhalina individuals that possessed clade I haplotypes (Fig. 2). Therefore, to more formally quantify the diagnostic potential of the mtDNA, we applied a simple formulation of Bayes' theorem (Sokal and Rohlf 1995) to estimate P(arcticola|I): the probability that an individual sampled on the nonbreeding grounds with a haplotype from clade I is actually C. a. arcticola rather than C. a. sakhalina. The probability is calculated as and relies on the following quantities: the probability of detecting clade I haplotypes in C. a. arcticola: P(I|arcticola) = 60/60 = 1.0; the probability of detecting clade I haplotypes in C. a. sakhalina: P(I|sakhalina) = 7/54 = 0.13; the probability of selecting a bird that is C. a. arcticola: P(arcticola); and the probability of selecting a bird that is C. a. sakhalina: P(sakhalina) = 1 À P(arcticola). Calidris alpina arcticola and C. a. sakhalina are believed to use similar areas during the winter, primarily Japan, coastal mainland China, Taiwan, and South Korea (Lanctot et al. 2009;Clements et al. 2013;Gill et al. 2013). Because P(arcticola) and P(sakhalina) reflect the probability of randomly selecting an individual from each subspecies, these quantities therefore depend on the abundance of each subspecies on the wintering grounds. Based on population estimates, there are 100 000 to 1 000 000 C. a. sakhalina individuals (Bamford et al. 2008), whereas 300 000 to 700 000 C. a. arcticola winter in East Asia (Andres et al. 2012). We therefore calculated P(arcticola|I) using the upper bound, lower bound, and approximate midpoint of each population size estimate in calculations.

Quantifying gene flow among Beringian subspecies
We used MIGRATE-N version 3.5.1 (Beerli and Palczewski 2010) to obtain Bayesian estimates of mutation-scaled effective population sizes and asymmetric migration rates among the three proximate Beringian subspecies (C. a. sakhalina, C. a. arcticola, and C. a. pacifica) that were most likely to exhibit gene flow. Limiting our analyses to three subspecies substantially reduced the number of parameters that needed to be simultaneously estimated, thereby providing a more tractable computational problem with a greater likelihood of success relative to analysis of the full data set (analysis required estimation of three as opposed to six effective population size parameters and six rather than thirty gene flow parameters) (Beerli 2009). MIGRATE-N estimates long-term effective population sizes as h = xN e l, where l is the mutation rate and x is an inheritance scaling factor that takes on values of 1 for mtDNA and 4 for codominant nuclear markers such as microsatellites. Long-term migration patterns are estimated over the time scales reflected by the set of sampled gene genealogies using the mutation-scaled quantity M = m/l, where m is the proportion of immigrants. Note that the product of the parameter estimates divided by the scaling factor (hM/x) provides a basis for estimating N e m, the effective number of immigrants into a population per generation.
Analysis parameter values and settings for MIGRATE-N were selected after preliminary exploratory analyses and with input from the program's developer (P. Beerli, personal communication). mtDNA analyses used the basic DNA sequence model, and priors for h were specified as a uniform distribution with minimum and maximum values of 0 and 0.03, respectively. Uniform priors with minimum and maximum values of 0 and 10 000 were likewise specified for M. Two independent runs based on random starting trees were performed to ensure convergence and consistency of parameter estimates. Each run was based on 10 6 recorded steps with a recording interval of 50 steps. Four concurrent chains were implemented during each run, with each chain using a static heating scheme based on temperature values of 1.0, 1.5, 3.0, and 10 5 . Microsatellite analyses were performed using the Brownian motion model. Lower and upper bounds for the uniform prior on h were specified as 0 and 10.0, whereas uniform priors for M were bound by 0 and 500. Two completely independent runs using starting UPGMA trees were performed, with each run based on 20 concurrent chains with 1000 recording steps made at 100 step intervals. The same heating scheme used for the mtDNA was applied to the microsatellites.     Table 2). Most haplotypes were restricted to a single subspecies (84 of 94; Appendix B).

Differentiation among subspecies
jModeltest2 identified the TrN+I+G model as most appropriate for ML analyses. The unrooted ML tree grouped the 94 unique haplotypes into four clades that included (1) a C. a. actites group, (2) a C. a. hudsonia group, (3) a C. a. kistchinski/sakhalina group, and (4) a group comprised primarily of C. a. arcticola/pacifica specimens (Fig. 2). With the exception of the detection of four C. a. arcticola/pacifica haplotypes in seven C. a. sakhalina specimens, there was no additional evidence of haplotype sharing among groups ( Fig. 2 and Appendix B). jModel-test2 indicated that the HKY model was most appropriate of those supported by MRBAYES. Trees from Bayesian analyses showed clear signs of convergence across the four runs (scale reduction factor of estimated parameters ranged from 0.99 to 1.01; standard deviation of split frequencies = 0.0093) and were virtually indistinguishable from the ML tree. Consequently, only the ML tree is presented here (Fig. 2).
There was highly significant differentiation among subspecies based on the mitochondrial data (Φ ST = 0.773,   Table 3A). All pairwise comparisons among subspecies were highly significant (Table 3A). Consistent with the phylogenetic analysis (Fig. 2), the lowest Φ ST values were detected in the C. a. pacifica/arcticola contrast and the C. a. sakhalina/kistchinski contrast-the two subspecies pairs that were not phylogenetically distinct in our analyses. As with the mtDNA data, C. a. sakhalina demonstrated the highest microsatellite allelic richness and H E values. However, the lowest microsatellite diversity was detected in C. a. actites ( Table 2). The microsatellites demonstrated no evidence for significant deviations from Hardy-Weinberg genotypic proportions after sequential Bonferroni corrections. Likewise, the 168 linkage disequilibrium tests performed (28 locus-pair analyses per subspecies * 6 subspecies) revealed only five significant results at the 0.05 level. These significant tests were detected across several subspecies (C. a. actites, C. a. pacifica, and C. a. sakhalina) and could have been observed by chance alone given the large number of individual tests that were performed.
The microsatellite analyses provided varying insights regarding genetic differentiation patterns in Dunlin. STRUCTURE suggested no evidence of differentiation among subspecies. Although the greatest average likelihood score was observed for the K = 5 case and the DK procedure suggested that there were K = 2 clusters, individual assignment probabilities to individual clusters were low and nearly uniform across clusters (Appendix C). This outcome indicates that the analysis procedure overestimated the true number of clusters and that subspecies-level subdivisions cannot be resolved with this analytical approach. In contrast, the global estimate of F ST from the microsatellite data indicated that significant genetic structure existed (F ST = 0.032, P < 0.001) (Table 3B). However, in comparison with the mitochondrial analysis, the microsatellite differentiation was generally small and reflected subtle differences in allele frequencies among subspecies (Appendix D). Most pairwise subspecific measures of differentiation were significant, with the exception of the comparison of C. a. sakhalina and C. a. kistchinski (F ST = 0.004, P = 0.129) ( Table 3B). The pairwise R ST values and their associated P-values were similar to those of F ST estimate, with the added finding of nonsignificant differentiation between C. a. arcticola and C. a. kistchinski (Table 3C).

Distinguishing C. a. arcticola from other subspecies that winter in Asia
Our STRUCTURE analyses suggested that the microsatellites possessed little utility for diagnosing subspecies (Appendix C). The GeneClass2 assignment tests provided similar insights. In general, success of the assignment approach was poor, with only 128 (34.6%) of the 370 individuals successfully assigned to the correct subspecies and only 31 of the 144 C. a. arcticola specimens (21.5%) correctly assigned. The average assignment probability of a properly assigned C. a. arcticola was only 0.576, indicating that there was low confidence in the correct assignments that were observed.
By contrast, our application of Bayes' theorem indicated a greater potential for genetic identification of C. a. arcticola if mtDNA data were used. In this case, the probability of a correct identification depends in part on the relative population sizes of C. a. arcticola and C. a. sakhalina (eqn 1; Table 4): the two subspecies that winter in Asia and that also can possess a type I haplotype. Using the upper and lower bounds of population size estimates for each subspecies, our calculations suggest that, under the extreme case where the ratio of C. a. sakhalina to C. a. arcticola is 1 000 000:300 000, the probability that a bird possessing a clade I haplotype is a C. a. arcticola individual is 0.698 (Table 4). This probability increases to 0.885 when population sizes are assumed to be equal and is as high as 0.982 when the population size of C. a. arcticola is assumed to be the upper extent of its estimated range and C. a. sakhalina is assumed to be at the lower extent of its range (Table 4).

Gene flow among Beringian subspecies
Results of MIGRATE-N analyses were comparable between independent runs for each data set, indicating that Table 4. Outcomes of calculations to infer P(arcticola|I): the probability that a randomly selected nonbreeding bird in East Asia with a haplotype from the main C. a. arcticola/pacifica group (Fig. 2) is actually C. a. arcticola as opposed to C. a. sakhalina. Calculations depend on the relative abundance of C. a. arcticola and C. a. sakhalina and are described in the Materials and methods (eqn 1). This table presents outcomes that evaluated upper, lower, and approximate midpoint population size estimates given by Bamford et al. (2008) and Andres et al. (2012). was also detected, albeit at lower levels than the rate into C. a. pacifica (Table 5).

Discussion
Migratory birds may facilitate the spread of HPAI from Asia to North America (Winker and Gibson 2010). In this investigation, we used large sample sizes and two genetic data sources (mitochondrial DNA and microsatellites) to determine genetic structure patterns among six Dunlin subspecies that reside in and migrate through eastern Asia and North America. We specifically focused on determining whether the four subspecies of Dunlin that winter in Asia can be differentiated and if genetic evidence for gene flow among Beringian subspecies exists. We suggest that our results may be useful for documenting potential HPAI transmission routes and the pathways that may facilitate the spread of disease across continents. Birds have reduced genetic structure relative to many other organisms, likely due to their capacity for flight and long distance movement (Greenwood and Harvey 1976;Zink et al. 1997). Many Arctic avian species, particularly migratory species, show lower levels of population genetic structure as a result of these high dispersal tendencies (Crochet 1996). For example, most shorebirds migrate long distances between breeding and nonbreeding areas (Brown et al. 2001), which may result in high gene flow and reduced genetic differentiation (e.g., Baker et al. 1994;Wenink et al. 1994;Haig et al. 1997;Wennerberg 2001b;Draheim et al. 2010;Miller et al. 2012). In contrast to past genetic studies of Dunlin that included limited sampling of Beringia-associated subspecies (e.g., Wenink et al. 1993;Wennerberg et al. 1999Wennerberg et al. , 2008, genetic analyses from our investigation revealed marked genetic differentiation among some Dunlin subspecies based on mtDNA analyses. Phylogenetic analysis revealed four separate phylogroups with high levels of statistical support (Fig. 2). Two of these groups consisted of samples from only C. a. hudsonia or C. a. actites, which occur in the most eastern and western regions of our study area. The other two groups contained mixtures of birds from more than one subspecies. The latter groups largely corresponded to birds that breed in relatively close proximity to one another, either in Asia (C. a. sakhalina and C. a. kistchinski) or in Alaska (C. a. arcticola and C. a. pacifica), although a few C. a. sakhalina birds from sites O and Q (Fig. 1) possessed haplotypes from the C. a. arcticola/ C. a. pacifica group (Fig. 2). The lack of clear structure between the C. a. sakhalina/kistchinski and C. a. arcticola/ pacifica groups suggests, in part, that the taxonomic status of these subspecies may require revision, although we recognize that other factors are important for defining subspecies (e.g., morphology, behavior, etc.; Haig et al. 2006).
Differentiation among subspecies was less pronounced based on the microsatellites, but significant structure was nonetheless detected between most subspecies pairs (Table 3). Male-biased gene flow (Clark et al. 1997; Gibbs Table 5. Bayesian estimates of mutation-scaled effective population sizes (h) and asymmetric migration rates (M) among the Dunlin subspecies C. a. arcticola, C. a. pacifica, and C. a. sakhalina. 95% credibility intervals are reported for each parameter, as is the derived parameter N e m reflecting the effective number of migrants per generation. See text for more details. Posterior distributions of estimated parameters are illustrated in Appendix E.  (Brown 1983) are plausible hypotheses that may explain differences between data sets. However, adult male Dunlin usually exhibit higher breeding site fidelity relative to females (Soikkeli 1967(Soikkeli , 1970Jackson 1994;Tomkovich 1994;Hill 2012). Thus, the lower effective population size and greater strength of genetic drift associated with maternally inherited haploid genomes may be the most reasonable explanation for the greater differentiation identified in the mtDNA data. Regardless of data set, the genetic structure patterns that we detected are likely the result of some degree of breeding site fidelity (Warnock and Gill 1996;Hill 2012) and reasonably strong population-specific migratory connectivity exhibited by some subspecies (Fern andez et al. 2008;Gill et al. 2013;S. Yezerinac and R. Lanctot, unpublished data).
Assuming that our sample of individuals and subspecies is representative of Dunlin from East Asia, our analysis suggests that we can use our data to obtain rudimentary estimates of the probability of correctly distinguishing Asian-versus Alaskan-breeding birds with mtDNA when sampling takes place in the East Asian nonbreeding areas. With the exception of seven C. a. sakhalina individuals, our representative mtDNA sequences from C. a. sakhalina, C. a. kistchinski, and C. a. actites (n = 107 total) were phylogenetically distinct from the haplotypes identified in Alaskan breeders (C. a. arcticola: n = 60; C. a. pacifica: n = 51; Fig. 2). Thus, if an individual possessed a haplotype associated with the C. a. actites or C. a. kistchinski/ sakhalina groups, the probability that the individual also breeds in Asia approaches 100% because no Alaska breeders possessed haplotypes from those groups. By contrast, if a bird sampled on the East Asia nonbreeding grounds possesses a haplotype from the main C. a. arcticola/pacifica group, our results suggest that the individual may either be C. a. sakhalina or C. a. arcticola ( Fig. 2; C. a. pacifica can be excluded from consideration given that this subspecies is entirely restricted to western North America). In this case, our application of Bayes' theorem indicates that there is nominally a~70% chance that a randomly selected bird possessing a haplotype from group I is C. a. arcticola (Table 4). The probability of a correct inference becomes even larger as the population size ratio of C. a. arcticola to C. a. sakhalina increases (Table 4). These probabilities are higher than the 53-60% correct assignment rates found by Gates et al. (2013) when using morphology to differentiate subspecies. Future analyses that combine genetic and morphological data may increase the likelihood of identifying C. a. arcticola in the Asian nonbreeding areas.
An unexpected outcome of our analyses included the detection of asymmetric gene flow from C. a. arcticola into C. a. pacifica and to a lesser extent also into C. a. sakhalina (Table 5). After considering potential reasons for this pattern, we highlight the simple fact that C. a. arcticola performs the longest spring migration out of all of the subspecies examined and that its northbound migration pathway crosses over part of the C. a. sakhalina and C. a. pacifica breeding areas (Fig. 1). It is feasible that some C. a. arcticola individuals 'short-stop' their migration in eastern Russia before crossing the Bering Sea to breed with C. a. sakhalina, and even more stop in western Alaska rather than continuing on to northern Alaska. Most reported cases of migratory short-stopping are associated with fall migrations en route to nonbreeding grounds, with the increased availability of supplemental food from agricultural systems (Wilson 1999;Jefferies et al. 2003) or climate change (Austin and Rehfisch 2005;La Sorte and Thompson 2007;Visser et al. 2009;Charmantier and Gienapp 2013) commonly invoked as possible explanations. In our case, we suggest that the frequency of short-stopping during spring migration may instead be correlated with poor weather conditions, resource limitations encountered during migration, or with the overall health and condition of the short-stopping individuals themselves. Evidence for migratory short-stopping during northbound breeding migrations has also been identified in lesser snow geese (Chen caerulescens caerulescens;Shorey et al. 2011). Given the shallow mitochondrial differentiation of C. a. pacifica and C. a. arcticola (Fig. 2), we also cannot rule out the possibility that the signal of asymmetric gene flow is the result of recent divergence of the two subspecies. However, a recent divergence does not preclude the possibility of ongoing gene flow, especially considering the geographic proximity of the breeding ranges of the two subspecies, the long migration flight undertaken by C. a. arcticola, and the fact that the northbound migratory path leads directly over C. a. pacifica's breeding range. In contrast, the signal of asymmetric gene flow from C. a. arcticola into C. a. sakhalina is most likely not the result of a recent divergence. The mtDNA-based phylogenetic tree illustrates that the two subspecies are reasonably well differentiated (Fig. 2), thereby leaving gene flow as a more tenable explanation for the analysis outcome.
Our finding of asymmetric gene flow indicates that, in addition to C. a. arcticola's usual northern Alaska breeding grounds, the western Alaska breeding grounds for C. a. pacifica need to be considered as a possible secondary entry point for Dunlin to carry AI into North America. This may be especially relevant if the migratory shortstopping behavior is influenced by an individual's health status, particularly if ill due to a viral disease. Because western Alaska and northern Alaska do not possess the same avian assemblages (Gabrielson and Lincoln 1959;Johnson and Herter 1989), the introduction of AI into western Alaska could lead to outbreaks in an additional and different suite of species than would an outbreak centered in northern Alaska.
The inference of asymmetric gene flow also implies the occurrence of direct interactions between C. a. arcticola and C. a. pacifica that could facilitate virus transmission between subspecies. Prior studies indicated that the two subspecies intermix during the fall after the breeding season (Taylor et al. 2011;Gill et al. 2013). If this were the only period of interaction, then the likelihood of HPAI spreading between subspecies would be low because any C. a. arcticola individuals harboring the virus would have had to (i) be infected on the wintering grounds and then (ii) live for 3-4 months with an active infection prior to intermixing with C. a. pacifica in the fall. However, our results suggest that individuals of the two subspecies sexually reproduce and thus likely share incubation duties for about 20 days (Warnock and Gill 1996). The breeding period occurs not long after migration and may coincide with the time when active shedding of HPAI by infected individuals is occurring.
Although our new findings do not specifically identify strategies for preventing the transmission of HPAI into North America, they nonetheless reveal a mechanism by which Dunlin could facilitate the spread of HPAI into North America and Mexico. This is particularly pertinent given that Dunlin are highly susceptible to infection with the H5N1 HPAI, and that some individuals may live to spread the disease, possibly after undergoing a migration (Hall et al. 2011). Although only a few Dunlin sampled in western North America have been documented with actively shedding AI (Ip et al. 2008;Iverson et al. 2008;US-FWS and USGS 2011), the continued emergence of new HPAI strains (e.g., H5N8, H7N9) and the fact that most efforts to date have detected prior exposure (i.e., antibodies, see Pearce et al. 2012;Johnson et al. 2014) indicates that the evolution of new strains remains problematic and that Dunlin are a potential route for HPAI to reach and spread within North America.