Estimation of population allele frequencies from next‐generation sequencing data: pool‐versus individual‐based genotyping

Molecular markers produced by next‐generation sequencing (NGS) technologies are revolutionizing genetic research. However, the costs of analysing large numbers of individual genomes remain prohibitive for most population genetics studies. Here, we present results based on mathematical derivations showing that, under many realistic experimental designs, NGS of DNA pools from diploid individuals allows to estimate the allele frequencies at single nucleotide polymorphisms (SNPs) with at least the same accuracy as individual‐based analyses, for considerably lower library construction and sequencing efforts. These findings remain true when taking into account the possibility of substantially unequal contributions of each individual to the final pool of sequence reads. We propose the intuitive notion of effective pool size to account for unequal pooling and derive a Bayesian hierarchical model to estimate this parameter directly from the data. We provide a user‐friendly application assessing the accuracy of allele frequency estimation from both pool‐ and individual‐based NGS population data under various sampling, sequencing depth and experimental error designs. We illustrate our findings with theoretical examples and real data sets corresponding to SNP loci obtained using restriction site–associated DNA (RAD) sequencing in pool‐ and individual‐based experiments carried out on the same population of the pine processionary moth (Thaumetopoea pityocampa). NGS of DNA pools might not be optimal for all types of studies but provides a cost‐effective approach for estimating allele frequencies for very large numbers of SNPs. It thus allows comparison of genome‐wide patterns of genetic variation for large numbers of individuals in multiple populations.     

Minimum sample sizes for population genomics: an empirical study from an Amazonian plant species

High‐throughput DNA sequencing facilitates the analysis of large portions of the genome in nonmodel organisms, ensuring high accuracy of population genetic parameters. However, empirical studies evaluating the appropriate sample size for these kinds of studies are still scarce. In this study, we use double‐digest restriction‐associated DNA sequencing (ddRADseq) to recover thousands of single nucleotide polymorphisms (SNPs) for two physically isolated populations of Amphirrhox longifolia (Violaceae), a nonmodel plant species for which no reference genome is available. We used resampling techniques to construct simulated populations with a random subset of individuals and SNPs to determine how many individuals and biallelic markers should be sampled for accurate estimates of intra‐ and interpopulation genetic diversity. We identified 3646 and 4900 polymorphic SNPs for the two populations of A. longifolia, respectively. Our simulations show that, overall, a sample size greater than eight individuals has little impact on estimates of genetic diversity within A. longifolia populations, when 1000 SNPs or higher are used. Our results also show that even at a very small sample size (i.e. two individuals), accurate estimates of F_ST can be obtained with a large number of SNPs (≥1500). These results highlight the potential of high‐throughput genomic sequencing approaches to address questions related to evolutionary biology in nonmodel organisms. Furthermore, our findings also provide insights into the optimization of sampling strategies in the era of population genomics.     

SNP development from RNA‐seq data in a nonmodel fish: how many individuals are needed for accurate allele frequency prediction?

Single nucleotide polymorphisms (SNPs) are rapidly becoming the marker of choice in population genetics due to a variety of advantages relative to other markers, including higher genomic density, data quality, reproducibility and genotyping efficiency, as well as ease of portability between laboratories. Advances in sequencing technology and methodologies to reduce genomic representation have made the isolation of SNPs feasible for nonmodel organisms. RNA‐seq is one such technique for the discovery of SNPs and development of markers for large‐scale genotyping. Here, we report the development of 192 validated SNP markers for parentage analysis in Tripterygion delaisi (the black‐faced blenny), a small rocky‐shore fish from the Mediterranean Sea. RNA‐seq data for 15 individual samples were used for SNP discovery by applying a series of selection criteria. Genotypes were then collected from 1599 individuals from the same population with the resulting loci. Differences in heterozygosity and allele frequencies were found between the two data sets. Heterozygosity was lower, on average, in the population sample, and the mean difference between the frequencies of particular alleles in the two data sets was 0.135 ± 0.100. We used bootstrap resampling of the sequence data to predict appropriate sample sizes for SNP discovery. As cDNA library production is time‐consuming and expensive, we suggest that using seven individuals for RNA sequencing reduces the probability of discarding highly informative SNP loci, due to lack of observed polymorphism, whereas use of more than 12 samples does not considerably improve prediction of true allele frequencies.     

Inferring the mode of origin of polyploid species from next‐generation sequence data

Many eukaryote organisms are polyploid. However, despite their importance, evolutionary inference of polyploid origins and modes of inheritance has been limited by a need for analyses of allele segregation at multiple loci using crosses. The increasing availability of sequence data for nonmodel species now allows the application of established approaches for the analysis of genomic data in polyploids. Here, we ask whether approximate Bayesian computation (ABC), applied to realistic traditional and next‐generation sequence data, allows correct inference of the evolutionary and demographic history of polyploids. Using simulations, we evaluate the robustness of evolutionary inference by ABC for tetraploid species as a function of the number of individuals and loci sampled, and the presence or absence of an outgroup. We find that ABC adequately retrieves the recent evolutionary history of polyploid species on the basis of both old and new sequencing technologies. The application of ABC to sequence data from diploid and polyploid species of the plant genus Capsella confirms its utility. Our analysis strongly supports an allopolyploid origin of C. bursa‐pastoris about 80 000 years ago. This conclusion runs contrary to previous findings based on the same data set but using an alternative approach and is in agreement with recent findings based on whole‐genome sequencing. Our results indicate that ABC is a promising and powerful method for revealing the evolution of polyploid species, without the need to attribute alleles to a homeologous chromosome pair. The approach can readily be extended to more complex scenarios involving higher ploidy levels.     

Restriction site‐associated DNA sequencing,genotyping error estimation and de novo assembly optimization for population genetic inference

Restriction site‐associated DNA sequencing (RADseq) provides researchers with the ability to record genetic polymorphism across thousands of loci for nonmodel organisms, potentially revolutionizing the field of molecular ecology. However, as with other genotyping methods, RADseq is prone to a number of sources of error that may have consequential effects for population genetic inferences, and these have received only limited attention in terms of the estimation and reporting of genotyping error rates. Here we use individual sample replicates, under the expectation of identical genotypes, to quantify genotyping error in the absence of a reference genome. We then use sample replicates to (i) optimize de novo assembly parameters within the program Stacks, by minimizing error and maximizing the retrieval of informative loci; and (ii) quantify error rates for loci, alleles and single‐nucleotide polymorphisms. As an empirical example, we use a double‐digest RAD data set of a nonmodel plant species, Berberis alpina, collected from high‐altitude mountains in Mexico.     

Exploring Approximate Bayesian Computation for inferring recent demographic history with genomic markers in nonmodel species

Approximate Bayesian computation (ABC) is widely used to infer demographic history of populations and species using DNA markers. Genomic markers can now be developed for nonmodel species using reduced representation library (RRL) sequencing methods that select a fraction of the genome using targeted sequence capture or restriction enzymes (genotyping‐by‐sequencing, GBS). We explored the influence of marker number and length, knowledge of gametic phase, and tradeoffs between sample size and sequencing depth on the quality of demographic inferences performed with ABC. We focused on two‐population models of recent spatial expansion with varying numbers of unknown parameters. Performing ABC on simulated data sets with known parameter values, we found that the timing of a recent spatial expansion event could be precisely estimated in a three‐parameter model. Taking into account uncertainty in parameters such as initial population size and migration rate collectively decreased the precision of inferences dramatically. Phasing haplotypes did not improve results, regardless of sequence length. Numerous short sequences were as valuable as fewer, longer sequences, and performed best when a large sample size was sequenced at low individual depth, even when sequencing errors were added. ABC results were similar to results obtained with an alternative method based on the site frequency spectrum (SFS) when performed with unphased GBS‐type markers. We conclude that unphased GBS‐type data sets can be sufficient to precisely infer simple demographic models, and discuss possible improvements for the use of ABC with genomic data.     

Sampling related individuals within ponds biases estimates of population structure in a pond‐breeding amphibian

Effective conservation and management of pond‐breeding amphibians depends on the accurate estimation of population structure, demographic parameters, and the influence of landscape features on breeding‐site connectivity. Population‐level studies of pond‐breeding amphibians typically sample larval life stages because they are easily captured and can be sampled nondestructively. These studies often identify high levels of relatedness between individuals from the same pond, which can be exacerbated by sampling the larval stage. Yet, the effect of these related individuals on population genetic studies using genomic data is not yet fully understood. Here, we assess the effect of within‐pond relatedness on population and landscape genetic analyses by focusing on the barred tiger salamanders (Ambystoma mavortium) from the Nebraska Sandhills. Utilizing genome‐wide SNPs generated using a double‐digest RADseq approach, we conducted standard population and landscape genetic analyses using datasets with and without siblings. We found that reduced sample sizes influenced parameter estimates more than the inclusion of siblings, but that within‐pond relatedness led to the inference of spurious population structure when analyses depended on allele frequencies. Our landscape genetic analyses also supported different models across datasets depending on the spatial resolution analyzed. We recommend that future studies not only test for relatedness among larval samples but also remove siblings before conducting population or landscape genetic analyses. We also recommend alternative sampling strategies to reduce sampling siblings before sequencing takes place. Biases introduced by unknowingly including siblings can have significant implications for population and landscape genetic analyses, and in turn, for species conservation strategies and outcomes.     

WFABC: a Wright–Fisher ABC‐based approach for inferring effective population sizes and selection coefficients from time‐sampled data

With novel developments in sequencing technologies, time‐sampled data are becoming more available and accessible. Naturally, there have been efforts in parallel to infer population genetic parameters from these data sets. Here, we compare and analyse four recent approaches based on the Wright–Fisher model for inferring selection coefficients (s) given effective population size (N_e), with simulated temporal data sets. Furthermore, we demonstrate the advantage of a recently proposed approximate Bayesian computation (ABC)‐based method that is able to correctly infer genomewide average N_e from time‐serial data, which is then set as a prior for inferring per‐site selection coefficients accurately and precisely. We implement this ABC method in a new software and apply it to a classical time‐serial data set of the medionigra genotype in the moth Panaxia dominula. We show that a recessive lethal model is the best explanation for the observed variation in allele frequency by implementing an estimator of the dominance ratio (h).     

The influence of locus number and information content on species delimitation: an empirical test case in an endangered Mexican salamander

Perhaps the most important recent advance in species delimitation has been the development of model‐based approaches to objectively diagnose species diversity from genetic data. Additionally, the growing accessibility of next‐generation sequence data sets provides powerful insights into genome‐wide patterns of divergence during speciation. However, applying complex models to large data sets is time‐consuming and computationally costly, requiring careful consideration of the influence of both individual and population sampling, as well as the number and informativeness of loci on species delimitation conclusions. Here, we investigated how locus number and information content affect species delimitation results for an endangered Mexican salamander species, Ambystoma ordinarium. We compared results for an eight‐locus, 137‐individual data set and an 89‐locus, seven‐individual data set. For both data sets, we used species discovery methods to define delimitation models and species validation methods to rigorously test these hypotheses. We also used integrated demographic model selection tools to choose among delimitation models, while accounting for gene flow. Our results indicate that while cryptic lineages may be delimited with relatively few loci, sampling larger numbers of loci may be required to ensure that enough informative loci are available to accurately identify and validate shallow‐scale divergences. These analyses highlight the importance of striking a balance between dense sampling of loci and individuals, particularly in shallowly diverged lineages. They also suggest the presence of a currently unrecognized, endangered species in the western part of A. ordinarium's range.     

Linkage disequilibrium network analysis (LDna) gives a global view of chromosomal inversions,local adaptation and geographic structure

Recent advances in sequencing allow population‐genomic data to be generated for virtually any species. However, approaches to analyse such data lag behind the ability to generate it, particularly in nonmodel species. Linkage disequilibrium (LD, the nonrandom association of alleles from different loci) is a highly sensitive indicator of many evolutionary phenomena including chromosomal inversions, local adaptation and geographical structure. Here, we present linkage disequilibrium network analysis (LDna), which accesses information on LD shared between multiple loci genomewide. In LD networks, vertices represent loci, and connections between vertices represent the LD between them. We analysed such networks in two test cases: a new restriction‐site‐associated DNA sequence (RAD‐seq) data set for Anopheles baimaii, a Southeast Asian malaria vector; and a well‐characterized single nucleotide polymorphism (SNP) data set from 21 three‐spined stickleback individuals. In each case, we readily identified five distinct LD network clusters (single‐outlier clusters, SOCs), each comprising many loci connected by high LD. In A. baimaii, further population‐genetic analyses supported the inference that each SOC corresponds to a large inversion, consistent with previous cytological studies. For sticklebacks, we inferred that each SOC was associated with a distinct evolutionary phenomenon: two chromosomal inversions, local adaptation, population‐demographic history and geographic structure. LDna is thus a useful exploratory tool, able to give a global overview of LD associated with diverse evolutionary phenomena and identify loci potentially involved. LDna does not require a linkage map or reference genome, so it is applicable to any population‐genomic data set, making it especially valuable for nonmodel species.     

Genomic distribution and estimation of nucleotide diversity in natural populations: perspectives from the collared flycatcher (Ficedula albicollis) genome

Properly estimating genetic diversity in populations of nonmodel species requires a basic understanding of how diversity is distributed across the genome and among individuals. To this end, we analysed whole‐genome resequencing data from 20 collared flycatchers (genome size ≈1.1 Gb; 10.13 million single nucleotide polymorphisms detected). Genomewide nucleotide diversity was almost identical among individuals (mean = 0.00394, range = 0.00384–0.00401), but diversity levels varied extensively across the genome (95% confidence interval for 200‐kb windows = 0.0013–0.0053). Diversity was related to selective constraint such that in comparison with intergenic DNA, diversity at fourfold degenerate sites was reduced to 85%, 3′ UTRs to 82%, 5′ UTRs to 70% and nondegenerate sites to 12%. There was a strong positive correlation between diversity and chromosome size, probably driven by a higher density of targets for selection on smaller chromosomes increasing the diversity‐reducing effect of linked selection. Simulations exploring the ability of sequence data from a small number of genetic markers to capture the observed diversity clearly demonstrated that diversity estimation from finite sampling of such data is bound to be associated with large confidence intervals. Nevertheless, we show that precision in diversity estimation in large outbred population benefits from increasing the number of loci rather than the number of individuals. Simulations mimicking RAD sequencing showed that this approach gives accurate estimates of genomewide diversity. Based on the patterns of observed diversity and the performed simulations, we provide broad recommendations for how genetic diversity should be estimated in natural populations.     

Why to account for finite sites in population genetic studies and how to do this with Jaatha 2.0

With the advent of next‐generation sequencing technologies, large data sets of several thousand loci from multiple conspecific individuals are available. Such data sets should make it possible to obtain accurate estimates of population genetic parameters, even for complex models of population history. In the analyses of large data sets, it is difficult to consider finite‐sites mutation models (FSMs). Here, we use extensive simulations to demonstrate that the inclusion of FSMs is necessary to avoid severe biases in the estimation of the population mutation rate θ, population divergence times, and migration rates. We present a new version of Jaatha, an efficient composite‐likelihood method for estimating demographic parameters from population genetic data and evaluate the usefulness of Jaatha in two biological examples. For the first application, we infer the speciation process of two wild tomato species, Solanum chilense and Solanum peruvianum. In our second application example, we demonstrate that Jaatha is readily applicable to NGS data by analyzing genome‐wide data from two southern European populations of Arabidopsis thaliana. Jaatha is now freely available as an R package from the Comprehensive R Archive Network (CRAN).     

A fast and reliable computational method for estimating population genetic parameters

The estimation of ancestral and current effective population sizes in expanding populations is a fundamental problem in population genetics. Recently it has become possible to scan entire genomes of several individuals within a population. These genomic data sets can be used to estimate basic population parameters such as the effective population size and population growth rate. Full-data-likelihood methods potentially offer a powerful statistical framework for inferring population genetic parameters. However, for large data sets, computationally intensive methods based upon full-likelihood estimates may encounter difficulties. First, the computational method may be prohibitively slow or difficult to implement for large data. Second, estimation bias may markedly affect the accuracy and reliability of parameter estimates, as suggested from past work on coalescent methods. To address these problems, a fast and computationally efficient least-squares method for estimating population parameters from genomic data is presented here. Instead of modeling genomic data using a full likelihood, this new approach uses an analogous function, in which the full data are replaced with a vector of summary statistics. Furthermore, these least-squares estimators may show significantly less estimation bias for growth rate and genetic diversity than a corresponding maximum-likelihood estimator for the same coalescent process. The least-squares statistics also scale up to genome-sized data sets with many nucleotides and loci. These results demonstrate that least-squares statistics will likely prove useful for nonlinear parameter estimation when the underlying population genomic processes have complex evolutionary dynamics involving interactions between mutation, selection, demography, and recombination.     

Unraveling hierarchical genetic structure in a marine metapopulation: A comparison of three high‐throughput genotyping approaches

Marine metapopulations often exhibit subtle population structure that can be difficult to detect. Given recent advances in high‐throughput sequencing, an emerging question is whether various genotyping approaches, in concert with improved sampling designs, will substantially improve our understanding of genetic structure in the sea. To address this question, we explored hierarchical patterns of structure in the coral reef fish Elacatinus lori using a high‐resolution approach with respect to both genetic and geographic sampling. Previously, we identified three putative E. lori populations within Belize using traditional genetic markers and sparse geographic sampling: barrier reef and Turneffe Atoll; Glover's Atoll; and Lighthouse Atoll. Here, we systematically sampled individuals at ~10 km intervals throughout these reefs (1,129 individuals from 35 sites) and sequenced all individuals at three sets of markers: 2,418 SNPs; 89 microsatellites; and 57 nonrepetitive nuclear loci. At broad spatial scales, the markers were consistent with each other and with previous findings. At finer spatial scales, there was new evidence of genetic substructure, but our three marker sets differed slightly in their ability to detect these patterns. Specifically, we found subtle structure between the barrier reef and Turneffe Atoll, with SNPs resolving this pattern most effectively. We also documented isolation by distance within the barrier reef. Sensitivity analyses revealed that the number of loci (and alleles) had a strong effect on the detection of structure for all three marker sets, particularly at small spatial scales. Taken together, these results illustrate empirically that high‐throughput genotyping data can elucidate subtle genetic structure at previously‐undetected scales in a dispersive marine fish.     

Comparison of 454 pyrosequencing methods for characterizing the major histocompatibility complex of nonmodel species and the advantages of ultra deep coverage

Characterization and population genetic analysis of multilocus genes, such as those found in the major histocompatibility complex (MHC) is challenging in nonmodel vertebrates. The traditional method of extensive cloning and Sanger sequencing is costly and time‐intensive and indirect methods of assessment often underestimate total variation. Here, we explored the suitability of 454 pyrosequencing for characterizing multilocus genes for use in population genetic studies. We compared two sample tagging protocols and two bioinformatic procedures for 454 sequencing through characterization of a 185‐bp fragment of MHC DRB exon 2 in wolverines (Gulo gulo) and further compared the results with those from cloning and Sanger sequencing. We found 10 putative DRB alleles in the 88 individuals screened with between two and four alleles per individual, suggesting amplification of a duplicated DRB gene. In addition to the putative alleles, all individuals possessed an easily identifiable pseudogene. In our system, sequence variants with a frequency below 6% in an individual sample were usually artefacts. However, we found that sample preparation and data processing procedures can greatly affect variant frequencies in addition to the complexity of the multilocus system. Therefore, we recommend determining a per‐amplicon‐variant frequency threshold for each unique system. The extremely deep coverage obtained in our study (approximately 5000×) coupled with the semi‐quantitative nature of pyrosequencing enabled us to assign all putative alleles to the two DRB loci, which is generally not possible using traditional methods. Our method of obtaining locus‐specific MHC genotypes will enhance population genetic analyses and studies on disease susceptibility in nonmodel wildlife species.     

The aggregate site frequency spectrum for comparative population genomic inference

Understanding how assemblages of species responded to past climate change is a central goal of comparative phylogeography and comparative population genomics, an endeavour that has increasing potential to integrate with community ecology. New sequencing technology now provides the potential to perform complex demographic inference at unprecedented resolution across assemblages of nonmodel species. To this end, we introduce the aggregate site frequency spectrum (aSFS), an expansion of the site frequency spectrum to use single nucleotide polymorphism (SNP) data sets collected from multiple, co‐distributed species for assemblage‐level demographic inference. We describe how the aSFS is constructed over an arbitrary number of independent population samples and then demonstrate how the aSFS can differentiate various multispecies demographic histories under a wide range of sampling configurations while allowing effective population sizes and expansion magnitudes to vary independently. We subsequently couple the aSFS with a hierarchical approximate Bayesian computation (hABC) framework to estimate degree of temporal synchronicity in expansion times across taxa, including an empirical demonstration with a data set consisting of five populations of the threespine stickleback (Gasterosteus aculeatus). Corroborating what is generally understood about the recent postglacial origins of these populations, the joint aSFS/hABC analysis strongly suggests that the stickleback data are most consistent with synchronous expansion after the Last Glacial Maximum (posterior probability = 0.99). The aSFS will have general application for multilevel statistical frameworks to test models involving assemblages and/or communities, and as large‐scale SNP data from nonmodel species become routine, the aSFS expands the potential for powerful next‐generation comparative population genomic inference.     

Population structure of Atlantic mackerel inferred from RAD‐seq‐derived SNP markers: effects of sequence clustering parameters and hierarchical SNP selection

Restriction‐site‐associated DNA sequencing (RAD‐seq) and related methods are revolutionizing the field of population genomics in nonmodel organisms as they allow generating an unprecedented number of single nucleotide polymorphisms (SNPs) even when no genomic information is available. Yet, RAD‐seq data analyses rely on assumptions on nature and number of nucleotide variants present in a single locus, the choice of which may lead to an under‐ or overestimated number of SNPs and/or to incorrectly called genotypes. Using the Atlantic mackerel (Scomber scombrus L.) and a close relative, the Atlantic chub mackerel (Scomber colias), as case study, here we explore the sensitivity of population structure inferences to two crucial aspects in RAD‐seq data analysis: the maximum number of mismatches allowed to merge reads into a locus and the relatedness of the individuals used for genotype calling and SNP selection. Our study resolves the population structure of the Atlantic mackerel, but, most importantly, provides insights into the effects of alternative RAD‐seq data analysis strategies on population structure inferences that are directly applicable to other species.     

Genome scans for the contemporary response to selection in quantitative traits

Genome scans have been an important approach for discovering historical signatures of selection in both model and nonmodel species. An exciting new experimental design for genome scans is to measure the change in allele frequency before and after contemporary selection within a generation, from a single population. The most widely‐used methods, however, have two major limitations: they are based on testing one locus at a time, and they only have power to uncover loci that have evolved under relatively strong selection. On the other hand, complex quantitative traits are common in nature and are caused by several loci of small effect. Selection on a quantitative trait at the phenotypic level is predicted to be accompanied by subtle allele frequency changes in many loci that covary (a polygenic soft sweep), rather than a large, single‐effect allele (a selective sweep). In this issue of Molecular Ecology, Bourret et al. (2014) measure the contemporary response to natural selection across the genome in multiple cohorts of Atlantic salmon during their first year at sea. They introduce a multilocus framework based on groups of markers that covary in their genotypic distribution. While the traditional, single‐locus approach did not find evidence for repeated patterns of selection, the multivariate approach found that a group of covarying SNPs was selected for in different cohorts at one site. Their multilocus framework has potential to be a more fruitful approach for uncovering the genomic basis of adaptation in quantitative traits, although caution should be applied as the framework has yet to be validated with simulated data.     

Demographic inferences using short‐read genomic data in an approximate Bayesian computation framework: in silico evaluation of power,biases and proof of concept in Atlantic walrus

Approximate Bayesian computation (ABC) is a powerful tool for model‐based inference of demographic histories from large genetic data sets. For most organisms, its implementation has been hampered by the lack of sufficient genetic data. Genotyping‐by‐sequencing (GBS) provides cheap genome‐scale data to fill this gap, but its potential has not fully been exploited. Here, we explored power, precision and biases of a coalescent‐based ABC approach where GBS data were modelled with either a population mutation parameter (θ) or a fixed site (FS) approach, allowing single or several segregating sites per locus. With simulated data ranging from 500 to 50 000 loci, a variety of demographic models could be reliably inferred across a range of timescales and migration scenarios. Posterior estimates were informative with 1000 loci for migration and split time in simple population divergence models. In more complex models, posterior distributions were wide and almost reverted to the uninformative prior even with 50 000 loci. ABC parameter estimates, however, were generally more accurate than an alternative composite‐likelihood method. Bottleneck scenarios proved particularly difficult, and only recent bottlenecks without recovery could be reliably detected and dated. Notably, minor‐allele‐frequency filters – usual practice for GBS data – negatively affected nearly all estimates. With this in mind, we used a combination of FS and θ approaches on empirical GBS data generated from the Atlantic walrus (Odobenus rosmarus rosmarus), collectively providing support for a population split before the last glacial maximum followed by asymmetrical migration and a high Arctic bottleneck. Overall, this study evaluates the potential and limitations of GBS data in an ABC‐coalescence framework and proposes a best‐practice approach.