Comparative population genomics reveals key barriers to dispersal in Southern Ocean penguins

The mechanisms that determine patterns of species dispersal are important factors in the production and maintenance of biodiversity. Understanding these mechanisms helps to forecast the responses of species to environmental change. Here, we used a comparative framework and genomewide data obtained through RAD‐Seq to compare the patterns of connectivity among breeding colonies for five penguin species with shared ancestry, overlapping distributions and differing ecological niches, allowing an examination of the intrinsic and extrinsic barriers governing dispersal patterns. Our findings show that at‐sea range and oceanography underlie patterns of dispersal in these penguins. The pelagic niche of emperor (Aptenodytes forsteri), king (A. patagonicus), Adélie (Pygoscelis adeliae) and chinstrap (P. antarctica) penguins facilitates gene flow over thousands of kilometres. In contrast, the coastal niche of gentoo penguins (P. papua) limits dispersal, resulting in population divergences. Oceanographic fronts also act as dispersal barriers to some extent. We recommend that forecasts of extinction risk incorporate dispersal and that management units are defined by at‐sea range and oceanography in species lacking genetic data.


| INTRODUC TI ON
During the last 2.58 million years, periods of global warming and cooling have shaped the evolutionary trajectories of species around the globe (Hewitt, 2004), with those inhabiting the polar regions having been subject to the most extreme climatic shifts. The Anthropocene (Lewis & Maslin, 2015) will be characterized by changes outside of this natural variability, contributing to the sixth global mass extinction event (Barnosky et al., 2011) and altering the evolutionary pressures acting on species. The movement of individuals among populations within species can distribute adaptive genetic variants, potentially aiding resilience to changing conditions, and preventing populations from diverging through genetic drift (Slatkin, 1987).
Barriers to dispersal can isolate populations, resulting in allopatric speciation if both populations are viable, or extirpation if they are not. For successful biodiversity conservation, understanding dispersal is key to distinguishing breeding populations and defining management units (Funk, McKay, Hohenlohe, & Allendorf, 2012), and to enable accurate forecasts of local or global extinction risk in response to habitat change.
Seabirds are highly threatened (Croxall et al., 2012), mobile and capable of dispersing large distances with few apparent abiotic barriers to movement, yet most are characterized by high levels of philopatry (Coulson, 2002). This appears to be the case for Antarctic and sub-Antarctic penguins, but they present a significant logistical challenge for studying dispersal (defined here as movement away from natal colonies to alternate breeding sites), as the vast majority of colonies are in remote locations. Banding studies initially suggested a high degree of philopatry in many species (Weimerskirch, Jouventin, Mougin, Stahl, & Van, 1985), and, until recently (Jenouvrier, Garnier, Patout, & Desvillettes, 2017), forecasts of extinction risk had not considered the potential buffering effect of dispersal (Cimino, Lynch, Saba, & Oliver, 2016;Jenouvrier et al., 2014). Genetic analyses (Clucas, Younger et al., 2016;Freer et al., 2015;Roeder et al., 2001;Younger, Clucas, et al., 2015, observations of colony movements (LaRue, Kooyman, Lynch, & Fretwell, 2015) and fluctuations in colony size (Kooyman & Ponganis, 2017) indicate that dispersal may be common. However, hydrographic features are thought to act as barriers to dispersal in a handful of sub-Antarctic and temperate penguin species (see Munro & Burg, 2017, for a review). A comprehensive study of dispersal barriers in Antarctic and sub-Antarctic penguins is therefore needed to clarify their potential evolutionary responses to increasing threats , improve the accuracy of estimates of extinction risk, identify effective management strategies for their conservation and shed light upon the mechanisms that prevent or promote dispersal in these charismatic seabirds.
Here, we examine the relative importance of different ecological and evolutionary factors in determining dispersal and population differentiation in the Aptenodytes and Pygoscelis genera using a comparative population genomic framework. By comparing range-wide patterns of genetic differentiation in ecologically divergent species with overlapping distributions and shared ancestry, we aim to tease apart the mechanisms which have led to the distributions of genetic diversity that we see today. Many factors are known to influence the patterns of dispersal in seabirds (Friesen, Burg, & McCoy, 2007), and our comparative framework is designed to examine the relative importance of several of these factors to dispersal: oceanographic fronts, ephemerality of breeding habitat, geographic continuity of breeding habitat, and at-sea range. We generated robust single nucleotide polymorphism (SNP) data sets for five species of penguin, covering the majority of each species' range (Figure 1), and compared the levels of dispersal within each species.

| Study species
Our study focused on Aptenodytes and Pygoscelis, which are sister genera within the penguin (Spheniscidae) family (Gavryushkina et al., 2017). We included all species within these genera: emperor (Aptenodytes forsteri), king (A. patagonicus), Adélie (Pygoscelis adeliae), chinstrap (P. antarctica), and gentoo penguins (P. papua). The breeding distributions of all species are shown in Figure 1. Emperor penguins are restricted to the Antarctic continent and breed primarily on sea ice, with a relatively continuous breeding distribution around the continent (Fretwell et al., 2012). Adélie penguins have a breeding distribution that encompasses ice-free areas of the Antarctic continent's coastline, along with several islands, all south of the Antarctic Polar Front (Schwaller, Southwell, & Emmerson, 2013 patagonicus), Adélie (Pygoscelis adeliae) and chinstrap (P. antarctica) penguins facilitates gene flow over thousands of kilometres. In contrast, the coastal niche of gentoo penguins (P. papua) limits dispersal, resulting in population divergences. Oceanographic fronts also act as dispersal barriers to some extent. We recommend that forecasts of extinction risk incorporate dispersal and that management units are defined by at-sea range and oceanography in species lacking genetic data.

K E Y W O R D S
Aptenodytes, genetic differentiation, Polar Front, population genomics, Pygoscelis, RAD-Seq penguins similarly breed on ice-free areas but are not widespread on the Antarctic continent-colonies are found only at the Antarctic Peninsula and various islands south of the Antarctic Polar Front (Borboroglu & Boersma, 2013). Gentoo and king penguins have more northerly distributions, and both of these species have colonies both north and south of the Antarctic Polar Front, a potential dispersal barrier (Clucas, Younger et al., 2016;Friesen, 2015;Munro & Burg, 2017). King penguins are found exclusively on islands, whereas gentoo penguins also breed on the Antarctic Peninsula (Borboroglu & Boersma, 2013). We aimed to encompass as much of these breeding ranges as possible in our study design ( Figure 1).

| Sampling and sequencing
Blood or tissue samples were collected from up to 16 individuals per colony across a large part of the range of each of the study species ( Figure 1). We sampled a single representative colony for those islands or archipelagos with multiple colonies. Colony names, collection dates and tissue types are provided in Supporting information Table S1. Further details of the tissues collected from Adélie penguins at Béchervaise Island, Welch Island, Blakeney Point and Pétrels Island can be found in Ref. (Younger, Emmerson, Southwell, Lelliott, & Miller, 2015). Details of the tissue samples collected from emperor penguins at Halley Bay, Fold Island, Auster, Amanda Bay and Pointe Géologie can be found in Ref. (Younger, van den Hoff, Wienecke, Hindell, & Miller, 2016;Younger, Clucas, et al., 2015). All other samples were blood samples. To take blood, penguins were held with the flippers restrained and the head placed under the arm of the handler, or they were wrapped in cushioned material covering the head and preventing movement, to minimize stress during handling (Le Maho et al., 1992). A second handler took up to 1 ml blood from the brachial, intertarsal or jugular vein using a 25-G or 23-G needle and 1-ml syringe, after cleaning the area with an alcohol swab. Total restraint time was generally 2 to 3 min. All field activities were conducted under appropriate permits and were subject to independent ethical review. Samples were either stored frozen for transport, stored frozen in ethanol or Queen's Lysis buffer (Seutin, White, & Boag, 1991) for transport, or stored in RNAlater (Life Technologies) and transported at ambient temperature. All samples were then stored frozen at −20°C in Australia or the UK. DNA was extracted from blood and tissue samples using QIAGEN DNeasy Blood and Tissue Kits. The digestion step was modified to include 40 μL proteinase K (at 20 mg/ml) and extended to 3 hr for blood samples. Details of the modifications made to the protocols for tissue samples are available in Younger, Emmerson, et al. (2015) and Younger, Clucas, et al. (2015). All samples were treated with 1 μL RiboShredder (Epicentre) to reduce RNA contamination, and DNA was visualized on a 1% agarose gel to confirm high molecular DNA was present. DNA concentration and purity were measured on a Qubit and NanoDrop (Thermo Fisher Scientific), respectively.
To identify genomewide SNPs for each species, we used standard restriction site-associated DNA sequencing (RAD-Seq) (Baird et al., 2008) with individual barcoding and the Sbf1 restriction enzyme.
The NERC Biomolecular Analysis Facility at Edinburgh Genomics (https://genomics.ed.ac.uk) performed the library preparation and sequencing as described by Gonen et al. (2014) following Etter, Bassham, Hohenlohe, Johnson, and Cresko (2011). In short, 250 ng of DNA per individual was digested with Sbf1-HF (NEB) and then ligated to barcoded P1 adapters. Individuals were multiplexed into 18 libraries consisting of 19-23 barcoded individuals and were sheared into fragments of <300-400 bp. The number of individuals to multiplex into each lane was estimated using the cutting frequency of the Sbf1 enzyme and a genome size of 1.2 Gb, as per the published Adélie and Emperor penguin reference genomes at the time of sequencing. Subsequent analyses have shown their genomes are likely to be 1.25 and 1.39 Gb, respectively (Li et al., 2014). Size selection was performed by gel electrophoresis. Libraries were blunt-ended (NEB Quick Blunting Kit) and A-tailed before P2 adapters (IDT) were ligated. Enrichment PCR and purification with AMPure beads were performed before libraries were checked for size and quantity using Qubit and qPCR assays. Each library was then sequenced in a lane of an Illumina HiSeq 2500 using 125 base paired-end reads in high output mode (v4 chemistry).

| Bioinformatics
Read quality was assessed with FastQC (http://www.bioinformatics.babraham.ac.uk/projects/fastqc). Demultiplexing, removal of reads with adapter contamination and trimming to 113 bp were performed with process_radtags from the stacks pipeline v1.35 (Catchen, Amores, Hohenlohe, Cresko, & Postlethwait, 2011;Catchen, Hohenlohe, Bassham, Amores, & Cresko, 2013). We also used process_radtags to remove any read pairs with uncalled bases, a low quality score and/or a barcode or cut site with more than one mismatch. King and emperor penguin reads were aligned to the published emperor penguin reference genome (http://gigadb.org/ dataset/100005; scaffold-level assembly) while Adélie, chinstrap and gentoo penguin reads were aligned to the published Adélie penguin reference genome (http://gigadb.org/dataset/100006; scaffold-level assembly) using bwa-mem (Li, 2013). Terminal alignments were prevented by enforcing a clipping penalty of 100, and reads with more than five mismatches, multiple alignments and/or more than two indels were removed using a custom python script F I G U R E 1 Breeding distributions (grey) and the locations of colonies sampled in this study (coloured triangles) for (a) king penguins, which were sampled both north and south of the Antarctic Polar Front (Borboroglu & Boersma, 2013); (b) emperor penguins (LaRue et al., 2015); (c) chinstrap penguins (Borboroglu & Boersma, 2013); (d) Adélie penguins (Schwaller et al., 2013); and (e) gentoo penguins, which were sampled both north and south of the Antarctic Polar Front (Borboroglu & Boersma, 2013). The number of individuals genotyped at each colony is indicated [Colour figure can be viewed at wileyonlinelibrary.com] (filter.py, available from https://doi.org/10.5061/dryad.7c0q8). PCR duplicates were removed with Picard Tools (http://broadinstitute. github.io/picard).
We called and filtered SNPs separately for each species using the Stacks pipeline, following many of the suggestions outlined in Benestan et al. (2016). All of the settings and filters were applied in the same way for each species. We ran the modules of the pipeline separately: pstacks -cstacks -sstacks -rxstacks -cstacks -sstacks -populations. Briefly, in pstacks we required a minimum stack depth (-m) of six reads mapping to the same location and used the bounded SNP model (significance level of α = 0.05, upper bound = 0.1, lower bound = 0.00041 corresponding to the highest sequencing error rate recorded by phiX spikes). We found that setting the stack depth to six gave sufficient numbers of polymorphic loci in all species after downstream steps had been performed, while maintaining sufficient (≥ 6X) depth at each locus to reliably call heterozygotes. In cstacks, we used all the individuals in each species to build the catalog of loci.
In rxstacks, we removed confounded loci with a conservative confidence limit of 0.25, removed excess haplotypes from individuals and also removed any loci with a mean log likelihood <−10. We filtered SNPs in the populations module to retain a single random SNP per RAD-tag, remove any SNPs with a minor allele frequency (MAF) < 0.01 or heterozygosity > 0.5, remove any loci that were not present in all colonies and remove any SNPs that were not genotyped in at least 80% of individuals per colony. In the adegenet package (Jombart, 2008;Jombart & Ahmed, 2011) in R (R Core Team, 2013), we calculated whether SNPs were in Hardy-Weinberg equilibrium (HWE) in each colony and used vcftools v0.1.13 (Danecek et al., 2011) to calculate mean coverage for each SNP. SNPs were removed if they were out of HWE in > 50% of the colonies or had a mean coverage greater than twice the standard deviation for the species. For the number of SNPs before and after filtering, please refer to Supporting information Table S3. pgdspider v2.0.8.2 (Lischer & Excoffier, 2012) was used to convert the vcf file into other formats for further analyses.

| Outlier loci detection
We identified SNPs that were potentially under selection in each species using the F ST outlier method in bayescan v2.1 (Foll & Gaggiotti, 2008). These loci were removed prior to coalescent-based or population genetic analyses that assume loci are evolving neutrally. With this aim, the high false-positive rate associated with BayeScan (Lotterhos & Whitlock, 2014) was not a concern, and its power to detect loci genuinely under selection under a range of demographic scenarios was advantageous (Lotterhos & Whitlock, 2014). We set a conservative prior on the odds of neutrality (for every five loci, our prior expectation is that one is under selection) to identify all loci that could potentially be under selection. We deemed q-values <0.1 to be significant, meaning that one in ten loci identified was expected to be a false-positive neutral locus (Lotterhos & Whitlock, 2014;Storey & Tibshirani, 2003). SNPs that were identified to be putatively under selection were removed using VCFtools. We refer to the remaining SNPs as "neutral SNP data sets," and those with the full complement TA B L E 1 Summary statistics of the data sets generated for each species before outliers were removed as the "total SNP data set," although these definitions are applicable only to gentoo penguins, as these were the only species in which a large number of outliers were identified (Table 1).

| Contemporary population structure and summary statistics
The number of private alleles in each colony was calculated with the populations module in Stacks, using the total SNP data sets for each species. We calculated the observed (H O ) and expected (H S ) heterozygosity for each colony with the neutral SNP data sets using genodive v2.0b27 (Meirmans & Van Tienderen, 2004). We also used GenoDive to calculate AMOVA-based F ST estimates for each species using 999 permutations to calculate significance, and the Weir and Cockerham unbiased weighted F ST estimator (Weir & Cockerham, 1984) between all pairs of colonies within each species. This measure has been shown to be robust to small sample sizes when F ST is low (Willing, Dreyer, & van Oosterhout, 2012).
Significance was calculated using 10,000 permutations of the data and corrected for multiple tests using the sequential goodness of fit method (SGoF+) (Carvajal-Rodriguez & de Uña-Alvarez, 2011). For gentoo colonies, we measured F ST using both the neutral SNP data set and the total SNP data set, whereas for the other four species, the neutral data sets were used, since the number of outlier loci was very small (Table 1).
We used three different clustering methods for each species using the neutral SNP data sets: the Bayesian clustering algorithm In each case, we first ran the model for 100,000 generations, discarding the first 50,000 as burn-in, setting K (the number of clusters) to one but allowing lambda to vary in order to estimate the species-specific value of lambda to use. For subsequent runs, the species-specific value of lambda was set and the number of clusters was allowed to vary from K = 1 to K = N, where N was the number of colonies sampled for that species. Each analysis was run for 150,000 generations, discarding the first 50,000 as burn-in, and repeated ten times from a different random seed. We used structure harvester web v0.6.94 (Earl, 2012) to compare replicates and prepare files for clumpp v1.1.2 (Jakobsson & Rosenberg, 2007), which aligns the results from replicate runs of structure to check for multimodality and calculates the average membership coefficients of each individual to each cluster, ready for visualization with distruct v1.1 (Rosenberg, 2004). Mostly, we did not use the results from the Evanno method for estimating the "true" number of clusters in the data (Evanno, Regnaut, & Goudet, 2005), as it is not defined for K = 1 and it was often hard to find biological meaning in the results. This is not unexpected, as the Evanno method has been shown to perform poorly for scenarios of moderate to low genetic differentiation (Waples & Gaggiotti, 2006). There is a large body of literature suggesting that it is unrealistic to expect that a "true" value of K exists for any given data set (Benestan et al., 2016;Gilbert et al., 2012;Janes et al., 2017) and that in order to gain insight into different levels of genetic structure, it is better to view and report multiple K-values. Therefore, we discuss our results for multiple values of K for each taxon, to fully understand the levels of structure in the data.
Secondly, we performed PCAs for each species using the neutral SNP data sets. Allele frequencies were scaled and centred, and missing values were replaced with the species' mean allele frequency using the scaleGen function from adegenet. PCA was computed with the dudi.pca function from the ade4 v1.7-11 package. The first three PCs were plotted against one another, but only the first two are shown, as population structure was not visible beyond PC2 in all species.
Finally, we used DAPC, which can be used to describe genetic clusters by creating synthetic variables (discriminant functions) that maximize the variance among clusters while minimizing the variance within them. When genetic differentiation is moderate to strong, individuals can be assigned to clusters using successive Kmeans clustering-the find.clusters function in adegenet (Jombart, 2008;Jombart & Ahmed, 2011)-before DAPC, thus negating the a priori assignment of individuals to groups determined by their sampling location. However, in all species other than gentoo penguins, the successive K-means clustering suggested K = 1 was most likely (the Bayesian Inference Criterion was at its minimum at K = 1) and so DAPC was performed when individuals were grouped by their colony of origin. Cross-validation with 1,000 replicates was used to determine the appropriate number of PCs to retain, and the posterior membership probability of each individual to each colony was plotted to determine how well individuals were assigned back to their colony of origin.

| Phylogenetics and species delimitation of gentoo penguins
To investigate the phylogeographic relationships among gentoo penguin colonies, we used the coalescent species tree approach implemented in the snapp package (Bryant, Bouckaert, Felsenstein, Rosenberg, & RoyChoudhury, 2012) in beast v.2.4.0 . SNAPP infers species trees from unlinked biallelic markers, such as SNPs. The method calculates species tree likelihoods directly from the data by estimating the probability of allele frequency change across nodes, thus avoiding the necessity of finding and combining individual gene trees. Nevertheless, the method is highly computationally demanding; therefore, we selected two random individuals (i.e., four haplotypes) per colony to include in the analysis and repeated the analysis twice with different individuals to ensure reproducibility. The neutral data set was used, and loci that were no longer polymorphic in the reduced set of individuals were removed, leaving 6,868 and 6,754 SNPs. The forward and backward mutation rates (u and v) were calculated from the data rather than estimated as part of the MCMC (note that in SNAPP analyses, the mutation rate (μ) is fixed at 1 and divergence times are not estimated). The MCMCs were run for three million generations with the first 10% discarded as burn-in. We monitored the traces for convergence using tracer v1.6 (Rambaut & Drummond, 2007), and when ESSs for all parameters were large (> 300) and the traces had reached stationarity, we concluded the analyses. densitree v2.0.1 was used to visualize the posterior distributions of topologies as cloudograms, hence allowing for a clear depiction of uncertainty in the topology.
To investigate whether described (Stonehouse, 1970) and puta- in the data set. We conducted 20 independent maximum-likelihood tree inferences and then drew bootstrap supports from 1,000 replicates onto the best scoring topology. All searches were conducted under the GTRGAMMA nucleotide substitution model.
To further investigate the validity of the proposed subspecies divisions within gentoo penguins (de Dinechin et al., 2012;Stonehouse, 1970), we compared species delimitation models using the Bayes factor delimitation method BFD* (Leaché, Fujita, Minin, & Bouckaert, 2014) as implemented within the snapp package (Bryant et al., 2012) in BEAST 2.4.3 . The BFD* method estimates the marginal likelihoods of competing species delimitation models using path sampling, such that the models can be ranked by their marginal likelihoods, with Bayes factors then used to assess support for the model rankings. A total of 16 individuals were included in the analysis, such that each putative subspecies unit within our species delimitation models had a minimum of four representatives.
Therefore, we included four individuals from each of the following: (d) all colonies grouped except for Kerguelen. We used the neutral SNP data set for the BFD* analysis, as required for coalescent-based methods. SNPs that were no longer polymorphic within the reduced data set of individuals were removed, leaving 8,116 SNPs. In SNAPP, we specified a gamma prior distribution for the theta parameter, with alpha = 1 and beta = 1000, for a prior mean of 0.001 on theta.
This value was chosen to reflect the 0.1% sequence divergence observed among gentoo penguin alleles known to belong to a single subspecies. The speciation rate, lambda, was fixed at 495, as calculated from the tree height and number of tips using the Python script yule.py (https://github.com/joaks1/pyule). Initial exploratory path sampling analyses were conducted to determine an appropriate number of steps to produce stability of the marginal likelihood, with 72 deemed more than sufficient. Path sampling analyses of 72 steps were then conducted for each of four species delimitation models with 100,000 MCMC generations following 10,000 pre-burn-in.
Finally, to estimate divergence times among the major gen- A maximum clade credibility tree with mean node heights was estimated from each posterior after removing the first 10% of samples as burn-in.

| Genotyping
We genotyped 376 individuals across five penguin species, representing 32 breeding colonies (Figure 1). To ease interspecific comparisons, we refer to island or archipelago names rather than colonies (Supporting information Tables S1, S2). RAD-Seq yielded an average of 11.6 million reads per individual, with 97.1% retained after quality control. Alignment to reference genomes, SNP calling and filtering generated high-coverage, neutral SNP data sets (Table 1).
Median sequencing depth per SNP ranged from 26.7 in chinstrap penguins to 32.1 in emperor penguins (Table 1, Supporting information Figure S1). Sequencing depth of individuals was similar among species (Supporting information Figure S2, Table S3), as was the distribution of minor allele frequencies (Supporting information Figure   S3). There was no evidence for lane effects in our data: Sequencing depth did not vary significantly among lanes, the locations of SNPs called within reads were even and similar across lanes, and the patterns of population structure that we found showed no relationship to the lanes on which individuals were sequenced.
Given the large number of outliers in the gentoo penguin data set, we investigated whether loci putatively under selection were influencing patterns of genetic differentiation by comparing pairwise F ST values from both the neutral and total SNP data sets (Supporting information   Tables S4, S5). The differentiation patterns were the same; therefore, we used the neutral SNP data set for all subsequent analyses.

| Levels of intraspecific variation among species
The mean minor allele frequency (MAF) did not differ among species (

| King penguins
Genetic differentiation among king penguin colonies at the Falkland Islands, South Georgia, the Crozet Islands and Macquarie Island ( Figure 1a) was subtle, despite being separated by thousands of kilometres of open ocean (Clucas, Younger et al., 2016). Most surprisingly, the Falkland Islands were subtly genetically differentiated from nearby South Georgia, ca. 1,400 km away but on the opposite side of the Polar Front, but indistinguishable from the Crozet Islands ca. 7,500 km away and on the same side of the Polar Front. Our analyses showed South Georgia to be the most divergent population, followed by Macquarie Island (Figure 4a, Supporting information Figure S5; Table S6; see also Clucas, Younger et al., 2016). Genetic differentiation across the range of the king penguin was subtle, as evidenced by small pairwise F ST values (Supporting information Table S6); invariant genetic diversity among colonies (Supporting information Table S2); the selection of K = 1 as the best clustering solution by both the successive K-means clustering algorithm and structure (Supporting information Figure S4e Figure S6).

| Emperor penguins
Among the eight sampled emperor penguin colonies (Figure 1b), there are four genetically differentiated metapopulations (Younger et al., 2017). The Ross Sea metapopulation appears the most divergent from the rest, consistent with findings from mitochondrial DNA (Younger, Clucas, et al., 2015).  Table S7) and posterior membership probabilities of individuals to their colony of origin were relatively high with DAPC (Supporting information Figure S11b).

| Chinstrap penguins
Population differentiation among chinstrap penguin colonies across their range was extremely low. Only three of ten pairwise F ST s were significant (Supporting information Table S8), successive K-means clustering could not detect any clusters, and structure had the highest mean posterior probability when K = 1 (Supporting information Figure S4c) and could not discern clusters without location priors  Figure S8). This differentiation was not as obvious with DAPC (Figure 3, Supporting information Figure S11c) or PCA (Supporting information Figure S5c), suggesting that they are only subtly differentiated. It should be noted that at the South Orkney Islands, we only sampled three individuals and conclusions regarding that colony should be considered preliminary.

| Adélie penguins
The Adélie penguin colonies sampled can be divided into the "western colonies" of the Antarctic Peninsula and Scotia Arc, and the "eastern colonies" along the Mawson Coast and East Antarctica ( Figure 1d). We uncovered subtle genetic differentiation coincident with this geographic division (Figure 3d, 4d). Despite separations in excess of 4,000 km, the genetic divergence was very slight, suggesting ongoing gene flow. Successive K-means clustering could not detect any genetic differentiation and the highest posterior mean log likelihood was achieved at K = 1, both with and without location priors in structure (Supporting information Figure S4d). Without  Table S9), and PCA showed only subtle differentiation between eastern and western colonies (Supporting information Figure S5).
Finally, a high level of admixture among colonies was visible when we plotted individual membership probabilities to each colony with DAPC (Supporting information Figure S11d).

| Gentoo penguins
For gentoo penguins, we designed our sampling scheme (Figure 1e) to include the two currently recognized subspecies: the northern gentoo (the nominate subspecies, Pygoscelis papua papua (Forster, 1781)), which is formally distributed north of 60°S; and the southern gentoo (Pygoscelis papua ellsworthii), formally distributed on the Antarctic Peninsula and maritime Antarctic islands south of 60°S (Clements et al., 2017;Murphy, 1947;Stonehouse, 1970); and the putative Indian Ocean subspecies (de Dinechin et al., 2012), which is still formally regarded as P. p. papua (Clements et al., 2017).
Our analyses showed four distinct groupings of gentoo pen-   Table S11). It should be noted that our estimates of divergence time differ markedly from those based on mitochondrial data (Clucas et al., 2014;Levy et al., 2016), most likely because mitochondrial data alone were unable to resolve the topology completely, resulting in more recent co- log likelihood occurred at K = 2; and individuals from South Georgia were fully assigned to a separate cluster (Supporting information Figure S10). For the results of all the hierarchical structure analyses for gentoo penguins, see Supporting information Figure S10.
Pairwise F ST values among the four clades ranged 0.127 to 0.298 and were all highly significant (p < 0.001, Supporting information Table S4). F ST values were two orders of magnitude greater than observed within the other penguin species studied, even though geographic distances among colonies were similar (Figure 2). The genetic diversity of the four clades was significantly different and greatest for Kerguelen (Figure 5c).
The southern gentoo colonies on the South Shetland Islands and western Antarctic Peninsula were clearly differentiated from one another in a structure analysis with and without location priors  Table S4; all other species range = −0.008 to 0.008, Supporting information Tables S6-S9). Given the geographic proximity of these colonies (50-400 km separations), this level of genetic differentiation is in stark contrast to the other species studied.

| Factors influencing patterns of genetic variation in penguins
To identify key factors that influence dispersal in penguins, we com-  (Weimerskirch et al., 1985;Woehler, 1989), and the recent formation of several new king penguin colonies (Delord, Barbraud, & Weimerskirch, 2004 ing either dispersal or a high incidence of skipped breeding. Overall, there is a growing body of evidence indicating that dispersal in many species of penguins is a regular occurrence. Many factors have been previously identified as potential drivers of dispersal patterns in seabirds (Friesen, 2015;Friesen et al., 2007), and we will discuss our genomic results with respect to the most relevant of these for Southern Ocean penguins.

| At-sea range
The most important factor in determining patterns of intraspecific genetic variation in penguins appears to be their at-sea range. The four species for which we found evidence of dispersal over large spatial scales are all considered pelagic, spending at least a portion of their life history in the open ocean far from their colonies.
Adult emperor, Adélie and chinstrap penguins travel up to 1,400 km , 2,200 km (Dunn, Silk, & Trathan, 2011) and 3,900 km (Hinke et al., 2015) away from their colonies during the nonbreeding period, respectively. Juvenile emperor penguins are even more mobile than adults, with recorded journeys covering in excess of 7,000 km in just eight months (Thiebot, Lescroël, Barbraud, & Bost, 2013), and individuals documented in the vicinity of other breeding colonies (Kooyman, Kooyman, Horning, & Kooyman, 1996;Wienecke, Raymond, & Robertson, 2010). In the winter, king penguins travel up to 1,800 km to forage in the marginal ice zone (Bost, Charrassin, Clerquin, Ropert-Coudert, & Le Maho, 2004;Charrassin & Bost, 2001) and juveniles have been observed at breeding colonies up to 5,600 km from their natal colonies (Weimerskirch et al., 1985). With the exception of the king penguin, the at-sea ranges of these pelagic penguins exceed the average distances between colonies. This wide-ranging behaviour is likely to facilitate dispersal, as evidenced by overall low genetic differentiation within all the pelagic species.
In contrast, gentoo penguins have a coastal lifestyle, rarely ranging beyond the continental shelf, and forage inshore on locally available prey (Lescroël & Bost, 2005;  have discrete breeding grounds (Clapham, 1996), but undertake long foraging journeys during the non-breeding season that bring them into contact with individuals from other breeding locations (Amaral et al., 2016). This appears to facilitate dispersal, with direct observations of individuals moving between breeding sites in different oceans (Stevick et al., 2011) and evidence of gene flow across ocean basins (Rosenbaum et al., 2009). Similarly, Atlantic Bluefin tuna (Thunnus thynnus) show site fidelity to two disparate spawning grounds on either side of the Atlantic Ocean, but individuals from both populations mix on their pelagic foraging grounds across the North Atlantic (Block et al., 2005). This migratory behaviour and intermixing of stocks likely facilitate gene flow, which may explain the low levels of genetic differentiation between the Gulf of Mexico and Mediterranean populations (Rooker et al., 2007).
However, the influence of wide-ranging behaviour on dispersal patterns cannot be generalized to all marine taxa. Atlantic salmon (Salmo salar) are highly migratory with individuals from both sides of the Atlantic mixing on foraging grounds off western Greenland, yet high levels of natal homing have led to genetic differences allowing individuals to be assigned back to their natal river with high success (King, Kalinowski, Schill, Spidle, & Lubinski, 2001 (Irwin, Irwin, & Smith, 2011;Ruegg et al., 2014).

| Natal philopatry
We observed genetic differentiation between gentoo penguin colonies separated by less than 50 km. This finding cannot be explained solely by its coastal lifestyle, because the species is known to visit other colonies within this range . Given the very small spatial scale over which population differentiation was observed, it is possible that natal philopatry also plays a role in limiting gene flow in gentoo penguins. On Possession Island (Crozet Archipelago), almost all the first breeding attempts occurred in a range of 2-5 km of the natal colony (C. A Bost, pers. comm.). Natal philopatry is thought to be common among seabirds (Coulson, 2002) and has been identified as a barrier to gene flow in other species, although it usually acts in combination with other isolating mechanisms (Friesen, 2015).
Paradoxically, the range of the gentoo penguin is expanding southwards coincident with sea ice decline (Lynch, Naveen, Trathan, & Fagan, 2012). The genetic differentiation found here would suggest that immigration rates at new colonies should be low. Instead, high rates of breeding success and recruitment may explain rapid colony growth after establishment.

| Oceanographic fronts
Both gentoo and king penguins have breeding distributions spanning the Antarctic Polar Front (Figure 1a, e), and our results indicate that the front may be a barrier to dispersal in both species.
The Antarctic Polar Front is the convergent boundary between cold Antarctic waters and warmer sub-Antarctic waters and constitutes an important feature for seabird communities (Bost et al., 2009).
King penguins from South Georgia, which is the only breeding population situated south of the front in our data set, were the most genetically divergent, albeit subtly. Furthermore, king penguins from the Falkland Islands were genetically indistinguishable from those at Crozet, ca. 7,500 km away but also north of the Polar Front,

| Breeding habitat quality, continuity and ephemerality
Emperor penguins have a relatively continuous distribution around Antarctica with most colonies being situated within the range of individuals foraging from adjacent colonies (Fretwell et al., 2012) ( Figure 1b). Their fast-ice breeding habitat is highly ephemeral, leading to changes in colony locations over years (Fretwell, Trathan, Wienecke, & Kooyman, 2014;LaRue et al., 2015;Trathan, Fretwell, & Stonehouse, 2011) and millennia (Younger, Clucas, et al., 2015;Younger et al., 2016). The low levels of genetic differentiation among emperor colonies likely reflect the need for flexibility in breeding location. The case of the Adélie penguin is similar, in that its breeding habitat is somewhat ephemeral, with access blocked periodically by sea ice or icebergs (Dugger et al., 2010), and with several large discontinuities in its circumpolar range (Figure 1d) where ice-free habitat does not exist. We found that Adélie penguins were subtly genetically differentiated across a gap of several thousand kilometres in their breeding distribution (Figure 1d)

| Cryptic diversity within gentoo penguins
The currently recognized taxonomy of gentoo penguins is for two subspecies, the northern (Pygoscelis papua papua) distributed north of 60°S, and the southern (Pygoscelis papua ellsworthii) distributed on the Antarctic Peninsula and maritime Antarctic islands south of 60°S (Clements et al., 2017;Forster, 1781;Murphy, 1947). Our data support the existing classification of a northern gentoo subspecies; however, contrary to current taxonomic limits (Clements et al., 2017), we found that South Georgian gentoos are more closely related to the southern subspecies than the northern, a conclusion that is supported by previous studies of morphology (de Dinechin et al., 2012), mitochondrial DNA (Clucas et al., 2014) and microsatellites (Levy et al., 2016). We recommend formal taxonomic revision of the boundary between northern and southern gentoo penguins to reflect this.
The degree of genetic divergence of gentoo penguins at Kerguelen points to a need for morphological and ecological study to determine whether these are a distinct species worthy of formal description. The case for revision has been based until now on mitochondrial DNA (Clucas et al., 2014;de Dinechin et al., 2012;Vianna et al., 2017) and microsatellites (Levy et al., 2016;Vianna et al., 2017), and we have now confirmed deep lineage divergences using genomewide data. In the light of these results, there is also an urgent need to characterize gentoo penguins breeding at other archipelagos using genomic data, particularly Crozet Archipelago and Macquarie Island, as there is likely to be more cryptic diversity. Accurate species boundaries and the recognition of cryptic species are crucial for the conservation of biodiversity, particularly in the light of the challenges  that will face Southern Ocean biota in the Anthropocene. The three lineages of gentoo penguins are on separate evolutionary trajectories. By conserving their full spectrum of genetic variation, the evolutionary and adaptive potential of gentoo penguins can be maximized.

| Predicting population structure
Understanding the mechanisms behind patterns of species dispersal has never been more important. Climate change is dramatically altering the marine environment, leading to changes in habitat availability, quality and ephemerality, as well as shifting oceanographic features.
Understanding current barriers to dispersal is essential for forecasting how species might respond to changes in their environment and for implementing ecologically meaningful conservation strategies.
Our findings show that at-sea range and oceanography are likely predictive of population structure in penguins. For species that journey into pelagic waters and range further than the average distance between colonies, we observed very little population differentiation. For colonies that are separated by oceanographic fronts, we observed greater genetic divergence than would be expected based on distance alone. We suggest that for colonies or species of penguins for which genetic data are unavailable, these predictive factors could be used to guide estimates of management units.

DATA AVA I L A B I L I T Y
The Illumina reads are available from the NCBI Sequence Read

CO M PE TI N G FI N A N CI A L I NTER E S TS
The authors do not have any competing financial interests.

E TH I C A L A PPROVA L
Ethical approval for all research involving the handling of birds was provided by the following: the University of Oxford, the University