An overview of current population genomics methods for the analysis of whole‐genome resequencing data in eukaryotes

Characterizing the population history of a species and identifying loci underlying local adaptation is crucial in functional ecology, evolutionary biology, conservation and agronomy. The constant improvement of high‐throughput sequencing techniques has facilitated the production of whole genome data in a wide range of species. Population genomics now provides tools to better integrate selection into a historical framework, and take into account selection when reconstructing demographic history. However, this improvement has come with a profusion of analytical tools that can confuse and discourage users. Such confusion limits the amount of information effectively retrieved from complex genomic data sets, and impairs the diffusion of the most recent analytical tools into fields such as conservation biology. It may also lead to redundancy among methods. To address these isssues, we propose an overview of more than 100 state‐of‐the‐art methods that can deal with whole genome data. We summarize the strategies they use to infer demographic history and selection, and discuss some of their limitations. A website listing these methods is available at www.methodspopgen.com.


| INTRODUC TI ON
Comprehensive analyses of species history and selection contribute to our understanding of causation in biology, an effort that has included genetics, developmental science and ecology (Laland et al., 2011). The number of population genomic studies aimed at elucidating the history of natural populations has increased enormously in the last 10 years. A few examples include an improved understanding of the history of human migrations, admixture and adaptation (e.g., Abi-Rached et al., 2011;Li & Durbin, 2011;Sabeti et al., 2002), the origin of domesticated species (e.g., Axelsson et al., 2013;Cubry et al., 2018;Schubert et al., 2014) and the genetic basis of local adaptation (e.g., Kolaczkowski et al., 2011;Kubota et al., 2015;Legrand et al., 2009;Roux et al., 2013). Developments in whole-genome resequencing have continually improved the throughput of genetic data, while reducing the time and cost of their production. Increased data production has been accompanied by a drive to develop efficient computational methods to interpret patterns of genetic variation at the genomic scale. These interconnected developments have allowed species histories to be inferred even when little preliminary knowledge is available. Investigating variation across multiple genomes sampled across populations or closely related species is now a common task for teams studying evolutionary processes, who can rely on a diverse array of methods to infer demography and selection. Such progress has confirmed the value of population genomics in understanding biological diversity, beyond the initial handful of model species upon which most of the field was built (Abzhanov et al., 2008;Ellegren et al., 2012;Jenner & Wills, 2007;Mandoli & Olmstead, 2000;Poelstra et al., 2014;Weber et al., 2013;White et al., 2010). Such advances are needed to broaden our view of evolutionary processes and improve sampling of distant clades.
Ultimately, this should provide a more balanced picture than the one brought by the study of a few model species (Abzhanov et al., 2008). From an applied perspective, genomic approaches also have the potential to improve conservation genetic inference by scaling up the amount of data available (Shafer et al., 2015), understanding the past history of species (Leitwein et al., 2020), and identifying loci and alleles important for local adaptation, which can then be used to define relevant conservation units (Fraser & Bernatchez, 2001).
Much effort has recently been made in facilitating the dissemination of sometimes complex, state-of-the-art methods. Nevertheless, the last comprehensive review of methods for population genetics was performed more than 10 years ago (Excoffier & Heckel, 2006). Recent methodological advances have brought increased analytical complexity to the field, and an inflation in the number of methods covering any one topic. The widespread use of sophisticated analytical tools is made difficult by the lack of communication between fields (Shafer et al., 2015), little user-friendliness of software, inflation of data formats (Lischer & Excoffier, 2012) and the ever-increasing number of methods made available. As a consequence, it has become increasingly difficult for all potential users (and also developers) to follow developments and be sure of selecting the most appropriate method for the question and data at hand.
Combining approaches is one of the current grand challenges in evolutionary biology (Cushman, 2014). While large-scale collaborations and sharing of skills between researchers allow for detailed analyses, a global summary of methods that can handle whole-genome resequencing data would be valuable for smaller research teams, so they can quickly start new projects and evaluate their experimental design. It would also facilitate communication between different subfields of evolutionary biology, by providing a common resource that can be used to identify methodological convergence and possible synergies. It may also avoid situations where similar methods are developed in parallel. Furthermore, the issue of anthropogenic environmental change and decline in biodiversity is pressing, and merits enhanced efforts to disseminate methods that can leverage genomic data, ultimately improving our understanding of the response of biodiversity to environmental change. Many conservation practitioners are receptive to using genetic tools, but do not always have access to the relevant expertise (Taylor et al., 2017). A freely accessible methodological summary may be useful in this context.
In this review, we assume that the reader is already familiar with the main concepts and current issues in population genomics, but needs an overview of the different methods associated with these concepts.
We promote the idea that multiple approaches must be used and compared in any population genomics project. This has several benefits: it gives the investigator a better idea of the robustness of results and may reveal issues in raw data processing. In addition, different methods aim to detect slightly different signals, and their combination may provide a more comprehensive overview of the processes acting. We aim at providing a resource that, if not fully comprehensive, can act as an efficient starting point for researchers investigating whole-genome variation in the next few years. This article can be used in combination with other recent methodological reviews on selection (Haasl & Payseur, 2016;Koropoulis et al., 2020), demographic inference or simulation-based approaches Smith & Flaxman, 2020).
For the sake of simplicity, we divide our review into two sections ( Figure 1): (i) methods devoted to the study of population structure and history (Tables 1 and 2), and (ii) detecting signatures of evolutionary processes along the genome (Tables 3 and 4). We end this review by outlining how different analyses can be combined, and present future directions that may be taken by the field of population genomics. We particularly insist on the interest-but also the challenges-of model-based approaches to test specific hypotheses, benchmark different methods and incorporate intrinsic properties of genomes (Table 4). The tables and a summary of the methods discussed in this paper will be kept updated to follow improvements, and are available at www.metho dspop gen.com. Contributions are of course welcome, and can be sent to the following email address: methodspopgen@gmail.com.

| P OPUL ATION S TRUC TURE AND HIS TORY
Genetic diversity and its genome-wide variance are directly impacted by variation in many factors including effective population sizes, population structure, inbreeding and migration. Moreover, the effects of selection on diversity at linked sites depends directly on local variation in the recombination rate. All these factors are important to characterize in any study of genome-wide variation. In this section, we describe methods aiming to quantify these aspects (see also Tables 1 and 2).

| Estimating familial relationships and reconstructing pedigrees
Understanding relatedness and structure both within and between populations is an important starting point for any study making inferences of selection or demographic history. Methods for estimating the relatedness of individuals are suited to studies relying on pedigree information (for example in quantitative genetics studies), or if there are reasons to suspect that familial relationships and inbreeding can play a major role in shaping the genetic structure of the population(s) considered. The most powerful methods in this category are likelihood-based and make use of heterozygous sites in each individual (e.g., colony in Wang, 2019, see also the detailed review in Huisman, 2017). Each pedigree configuration can be assigned a likelihood at each locus which depends on the probability of observing a given genotype conditional on the genotypes of assigned parents. Assuming independent loci, a composite likelihood can then be derived for a set of unlinked single nucleotide polymorphisms (SNPs). Information about pedigrees is important in order to filter out related individuals before carrying other population genetics analyses. Furthermore, mendelian constraints provide important information about haplotypes that can be used by phasing programs. Including related individuals can be useful when attempting to phase genotypes and generate a reference panel for further phasing in unrelated samples. The popular phasing algorithm Shapeit (Delaneau et al., 2019;Williams et al., 2012) can include familial information when reconstructing phased haplotypes.

| Using nonsupervised models to estimate relatedness and population structure
An elegant and efficient class of methods relies on using multivariate approaches such as principal component analysis (PCA) to infer relatedness between individuals and populations without a priori knowledge. These methods apply a dimension reduction procedure to matrices of individual genotypes, projecting genotypic variability along several axes of variation ). These approaches have been especially useful to study the consistency between geographical and genetic structure in human populations of Europe (Novembre et al., 2008). Procrustes rotation (Novembre et al., 2008) can be used to match geographical coordinates with PCA axes, showing how isolation by distance has shaped genetic structure. Since these methods do not have underlying assumptions based on (diploid) population genetics, they are suitable for analysing species displaying polyploidy or mixed-ploidy (Dufresne et al., 2014).
They go beyond a mere description of data, since projections of individuals on PCA axes can be used to infer admixture proportions, and contain information about demographic processes shaping genetic diversity (McVean, 2009). PCA can be used as a summary of genetic F I G U R E 1 Graphical summary of this review. This work is divided into two main sections: the first section covers methods that generally assume neutrality and are generally used for demographic inference. The second section covers methods that aim at identifying loci under the direct and indirect effects of selection. The table listing the relevant methods is indicated at the bottom right of each box. The linear structure of this review does not necessarily reflect the network of possible comparisons between the results of different methods. These possible comparisons are indicated by double arrows. Results obtained from methods listed in different sections can be used to inform the next steps of an analysis (single arrows). We acknowledge that there is no one-size-fits-all pipeline, and elements of this general framework may be entirely omitted from an analysis depending on the research question

| Model-based inference of population structure
Unlike the previous set of "algorithmic" approaches (see taxonomy

| Heterogeneous structure in space: Landscape genomics
Some methods can explicitly use spatial information to inform clustering, allowing improved consideration of the effect of landscape heterogeneity on selection against migrants and drift (e.g., spacemix, tess3, Table 1).
This spatial perspective can be useful to visualize the location and shape of hybrid zones (Guedj & Guillot, 2011 There exist methods that complement the approaches described above, by identifying which combination of geographical and ecological distance limits dispersal. A good example is bedassle, which uses the deviation of allele frequencies at unlinked sites in local populations from the global average, and estimates genetic covariance between all pairs of populations. It then uses a spatial model to estimate the strength of association of covariance with environmental features, assuming a negative relationship between genetic and environmental distances. However, disentangling these effects has proved to be complex. A deeper analysis of genes more strongly impacted by either geography or ecology may be more informative when it comes to the proximate causes of reduced dispersion and differentiation, such as biased dispersal (Bolnick & Otto, 2013;Edelaar & Bolnick, 2012) or selection against migrants (Hendry, 2004). Landscape genomics now extends its focus to adaptive genetic variation, and benefits from new methods targeting signatures of selection (see below).

| Inferring phylogenetic relationships
Recent advances in molecular phylogenetic methods, and the employment of different types of next-generation sequencing (NGS) data is well beyond the scope of this review (see, e.g., Moriarty Lemmon & Lemmon 2013;Cruaud et al. 2014;Wen et al. 2015). In this respect, both maximum likelihood and Bayesian approaches have become popular to investigate evolutionary relationships between individuals from different populations, even when divergence is very recent (e.g., Wagner et al., 2013).
These methods are implemented in software such as raxml (Stamatakis, 2014) and beast2 (Drummond & Rambaut, 2007). Ultimately, all molecular phylogenies reconstruct the genealogy of the genes with which they have been constructed. A problem when applying phylogenetic methods, especially in the context of recent divergence, is the assumption that gene trees are representative of lineage history. This assumption is likely to be violated at the population level, since the influences of gene flow and incomplete lineage sorting are strong at this scale, and may cause gene trees to deviate from population history. Fortunately, many recent coalescent-based methods in phylogenomics explicitly model gene trees to fit inside a speciation framework.
When using genome-wide data at the population level, methods specifically dedicated to reconstructing multiple species coalescent (MSCs) models such as *beast (starbeast) may be preferred over concatenation (Edwards et al., 2016), since they accommodate fluctuations in genealogical history across the genome, allowing discordance between species trees and individual gene trees to be identified.
However, in the presence of strong gene flow, MSC models can underestimate divergence times and overestimate effective population sizes (Leaché et al., 2014), because they attempt to explained the observed diversity with a strict isolation model. This issue is partially tackled by methods such as phrapl (Jackson et al., 2017), which estimates the likelihood of complex histories by examining genealogies at multiple genes and comparing them with coalescent simulations.
Such integration is particularly needed for species and populations that are in the "grey zone of speciation" (Roux et al., 2016).
While useful to infer topologies, caution is advised when using branches lengths obtained from SNP-only data sets, foe example to calculate divergence times between different groups or species (Leaché et al., 2015). For this purpose, it might be more straightforward and reliable to extract from the data both variant and invariant sites at several genes (e.g., coding or conserved sequences), and analyse the whole sequences in software like beast2. Such analyses can also be performed in two steps: first, estimate the phylogenetic relationships between samples, then apply a molecular clock model to obtain times since divergence. This style of approach is implemented in the Bayesian method MCMCTree in the paml package (Yang, 2007).

| Inferring demographic history with likelihood methods based on the allele frequency spectrum
The allele frequency spectrum (AFS) is the distribution of allele frequencies at polymorphic loci in one or several populations (called in that case joint or multipopulation spectrum). Different patterns of gene flow and demographic events all shape the AFS in specific ways (e.g., alleles are likely to occur at more similar frequencies if divergence is recent or if populations are highly connected). Several methods use the AFS to infer the demographic events explaining current genomic diversity. Two of the most popular methods (∂a∂i and fastsimcoal2; Gutenkunst et al., 2009;Excoffier et al., 2013) fit population genetics model specified by the user to the observed spectrum using a maximum-likelihood approach. The AFS expected under a given scenario is obtained through simulation, either using a diffusion approach (∂a∂i), or coalescent simulations (fastsimcoal2).
These approaches quickly estimate parameters using composite likelihoods, but do not explicitly take into account correlations induced by linkage disequilibrium (LD) between physically linked markers (but see able; Beeravolu et al., 2018). This might limit power to detect recent demographic events (e.g., migration, Jenkins et al., 2012).
Including SNPs that are physically close together should not strongly bias parameter estimation. However, such an approach prevents direct comparisons of likelihoods from different models. Therefore, physically independent SNPs should be used to consider composite likelihoods as quasi likelihoods for model comparison (Excoffier et al., 2013). Using allele frequencies estimated from pooled data sets is also feasible, as illustrated by a recent study on hybridization in Populus species where AFS was estimated from pooled whole genome resequencing data (Christe et al., 2016). The same applies for low-depth-sequencing data, with software such as angsd or atlas that are able to extract the most likely AFS and other relevant summary statistics. Such approaches are particularly promising to analyse whole-genome data from species with large genomes, ancient DNA samples, or when sequencing costs would otherwise be too prohibitive (Box 1).

Methods have been developed to infer variation in population sizes
with time using the whole genome of one or several diploid individuals.
Briefly, these methods model successive genealogies along the genome sequence as a Markov process: the genealogy at one locus only depends on the genealogy at the previous locus. Changes in the topology

BOX 2 Efficiently simulate whole-genome data
Simulations of whole-genome data are poised to become a standard tool for researchers, and recent initiatives such as stdpopsim, an open library of population genetics simulation models for multiple species, might help design reproducible simulations (Adrion, Cole, et al., 2020). More than 145 genetic simulators are currently available, but not all can handle genome-sized data (see https://surve illan ce.cancer. gov/genet ic-simul ation -resou rces/). Simulated data can be used to define significance thresholds for summary statistics when trying to scan the genome for regions under selection. Simulations are also at the core of simulation-based algorithms such as ABC or supervised machine-learning.
By comparing simulations with observed data, these methods can identify the processes that underlie diversity in any given genomic region.
There are two main categories of simulators, those based on coalescent ("backward in time" simulators), and forward-in-time simulators. The ms software, with its extensions (such as msms, Moreover, despite the abundance of species that practice self-fertilization and asexual reproduction, only facsexcoalescent is able to model coalescence in facultatively sexual species. Forward-in-time simulators such as slim3 ( are sequenced as a single library, can be an option to reduce costs. Summary statistics along the genome and allele frequency spectra can then be extracted for each population (e.g., using methods such as popoolation; Kofler, Orozco-terWengel, et al., 2011;. Since individual information is not available, variation in LD across individuals cannot be fully exploited, but methods such as ∂a∂i can still be used to test complex demographic scenarios (Gutenkunst et al., 2009). Shallow shotgun sequencing (1-5× per individual) is another approach that gives access to individual information for a similar cost (Buerkle & Gompert, 2013), but might prevent using methods requiring accurate phasing and unbiased individual genotypes.
Nevertheless, recent methods such as those implemented in the packages angsd  or atlas (Link et al., 2017) are promising. For example, atlas includes an approach to reconstruct past demographic histories by applying the pairwise sequentially Markovian coalescent (PSMC) to low-depth ancient DNA samples. angsd comes with several methods that estimate relatedness in lowdepth samples (ngsrelatev2; Hanghøj et al., 2019), and can estimate allele frequency spectra that can be used for demographic and selection inference. Among the most powerful methods available, recent versions of arg-weaver are promising since they can take into account genotype quality when reconstructing genealogies along the genome, and can therefore be applied to "low-quality" samples. One of the main drawbacks is that such analyses take time, making arg-weaver more suited to investigating genealogies in a limited set of genomic regions of interest.
are due to recombination events reconnecting branches in the tree.
The whole genealogy is usually not estimated, however, which results in drastic gains in speed. Such methods have the advantage of requiring only a small number of individuals (1-10), no a priori knowledge of population history, and permitting time-varying gene flow to be incorporated (see msmc-im). One general drawback, however, is that they are limited to rather simple scenarios, and do not handle more than two populations as yet (but see dical2, Table 2). While powerful, they are sensitive to confounding factors such as population structure (Orozco-terWengel, 2016) that lead to false signatures of expansion or bottlenecks. These methods also do not allow extremely recent demographic events to be investigated, since the coalescence of two alleles from a single individual in the recent past (a few tens to hundreds of generations) is infrequent. Moreover, most of these methods require the data to be phased (but see smc++; Terhorst et al., 2016), for example with fastphase (Scheet & Stephens, 2006) or beagle (Browning & Browning, 2011). In addition, phasing errors can lead to strong biases in parameter estimates for recent times (Terhorst et al., 2016). An extension of these methods takes into account population structure and aims to identify the number of islands contributing to a single genome, assuming it is sampled from a Wright n-island metapopulation (Mazet et al., 2015;Rodríguez et al., 2018). Such developments should improve the amount of information retrieved from only a few genomes.

Methods based on tracts of identity-by-descent (IBD, Palamara
& Pe'er, 2013) constitute an interesting alternative for more complex model testing when whole genomes are available in large number. Such methods allow recent demographic events to be inferred with relative precision. They are used to predict the length of haplotypes shared by two individuals that are inherited from a common ancestor without recombination. However, IBD detection requires large cohorts and accurate phasing, and therefore application of these methods has been largely restricted to human populations so far (Browning & Browning, 2011;Palamara & Pe'er, 2013). Another approach has used tracts of identity-by-state (IBS) to perform demographic inference over a range of timescales (Harris & Nielsen, 2013). IBS tracts are directly observable since they are simply the intervals between pairwise differences in an alignment of sequences and do not require any assumption about coancestry to be defined. The method predicts the length distribution of IBS tracts for pairs of haplotypes under a range of demographic parameters. These predicted spectra are then compared to empirical data under a likelihood framework, as with methods based on the AFS.

| Selection, introgression and their impact on sequence variation
While demographic forces such as drift and migration will affect the whole genome, selection in the presence of recombination is expected to be specific to particular portions of the genome, and therefore yield discrepancies with genome-wide polymorphism (Lewontin & Krakauer, 1973; but see section 3.9). Both positive and negative selection have long-distance effects on sites that are adjacent to those under selection, an effect often put under the umbrella term of "linked selection" (Cruickshank & Hahn, 2014;Ravinet et al., 2017). These effects are stronger in regions of low recombination, and may explain the correlations observed between nucleotide diversity, divergence metrics and recombination rates that are observed across many clades (Charlesworth et al., 1997;Cruickshank & Hahn, 2014). Using wholegenome resequencing data, it is possible to estimate the effective recombination rate along the genome (see the Recombination class of methods in Table 1). Such estimates are particularly useful in the absence of pre-existing genetic maps to assess how recombination and linked selection may bias estimates of diversity statistics or scans for selection. It also provides a way to determine a suitable window size to compute "independent" statistics along genomic windows.
In the sections that follow, we describe different methods aiming at identifying regions under selection by contrasting local patterns of diversity and divergence with genome-wide patterns. We begin with approaches focusing on single populations, and then summarize those focused on multiple populations.

| Quantifying positive and purifying selection on coding regions
The ratio between the number of nonsynonymous and synonymous mutations (also called dn/ds, K A /K S or ω) is often used to detect whether a specific gene is undergoing negative (ω < 1) or positive can be compared and tested in polydfe (Tataru & Bataillon, 2019).
Note that a very detailed tutorial with scripts is available in Tataru and Bataillon (2020) for the latter method.

| Detecting selective sweeps (recent positive selection)
Selective sweeps reduce diversity in genomic regions flanking the selected site(s). This leads to local deviations in the shape of the AFS that can be captured by several summary statistics computed over genomic windows, such as π, the nucleotide diversity (Nei & Li, 1979), Tajima's D (Tajima, 1989), and Fay and Wu's H (Fay & Wu, 2000). Using a combination of these statistics allows targets of selection to be identified with greater precision, and minimizes the confounding effects of demography. However, defining a threshold beyond which the values of a set of statistics supports selection is nontrivial. Recent developments in machine learning and Approximate Bayesian Computation (ABC) may assist in this regard (see section 3.9 below). Directly contrasting genome-wide with local AFS is another option that does not require combining results from multiple summary statistics. This approach has been used to develop composite tests, such as the composite likelihood ratio (CLR) test (Degiorgio et al., 2016;Stamatakis et al., 2013) that aims to detect recent selective sweeps by maximizing the likelihood of a model with selection in a genomic window, and comparing it to a model built on SNPs sampled from the genomic background.
In regions near to a selected allele, it is expected that LD is increased and diversity is decreased, especially after recent positive selection. A class of methods are aimed at targeting those regions that display an excess of long homozygous haplotypes, such as the extended haplotype homozygosity (EHH) test (Sabeti et al., 2002). It is also possible to compare haplotype extension across populations, with the Cross Population Extended Haplotype Homozygosity test (XP-EHH; Sabeti et al., 2007) or Rsb (the standardized ratio of EHH at a given SNP site; Tang et al., 2007). These methods require data to be phased in order to reconstruct haplotypes, which can make them susceptible to switch-errors. Nevertheless, methods based on LD may be more sensitive to selection on standing variation or on multiple alleles that leave a more subtle signature (so called soft sweeps).
Statistics dedicated to the detection of soft sweeps include the nSL statistics (Ferrer-Admetlla et al., 2014) in selscan or the H2/H1 statistics (Garud et al., 2015). These statistics usually examine the distri-  (Schrider et al., 2015), as well as their relative importance (Jensen, 2014;Messer & Petrov, 2013). Even hard selective sweeps can be challenging to detect with LD-based statistics especially under unstable demography and weak selection (Jensen, 2014). It is advisable to combine several approaches to improve confidence when pinpointing candidate genes for selection. Methods based on LD alone can sometimes miss the actual variants under selection due to the impact of recombination on local polymorphism that can mimic soft or ongoing hard sweeps (Schrider et al., 2015).

| Detecting long-term balancing selection
Unlike directional selection, balancing selection can lead to the maintenance of polymorphism at selected loci over long periods of time. This type of selection is extremely relevant for evolutionary biologists (Sellis et al., 2011), since it is at the core of strong coevolutionary dynamics such as host-parasite interactions (Ebert & Fields, 2020). Despite its importance, balancing selection has often been overlooked. This is mostly due to its narrow effects, particularly in the case of long-term balancing selection where recombination erodes association between loci under selection and neutral neighbouring loci. Nevertheless, the emergence of whole-genome resequencing data has facilitated the investigation of these narrow signals. Several recent methods and summary statistics (see Table 3, "Detecting balancing selection") have been specifically developed to detect this type of selection (Bitarello et al., 2018;DeGiorgio et al., 2014;Rasmussen et al., 2014;Siewert & Voight, 2017 (Charlesworth, 2006), and be detected by methods aimed at detecting long haplotypes and low diversity.  (Song et al., 2011).

| Detecting introgressed genomic regions
Summary statistics can be useful to obtain a first set of candidates for introgression and selection. One may, for example, plot the distribution of a differentiation measure such as F ST (Weir & Cockerham, 1984) between populations, estimates of effective recombination rates and nucleotide diversity along the genome. Such an approach has been used in Darwin's finches, which uncovered genomic islands of divergence with low recombination rates resisting gene flow (Han et al., 2017). Other approaches, such as chromosome painting (Table 1)

| Identifying highly differentiated loci and associations between allele frequencies and environmental features
When an allele is under positive selection in a population, its frequency tends to rise to fixation, unless gene flow from other populations or strong drift prevents this from happening (Charlesworth et al., 1997). It is therefore possible to contrast patterns of differentiation between populations adapted to their local environment to detect loci under divergent selection (e.g., displaying a high F ST ).
However, it is essential to control for population structure, as it may strongly affect the distribution of differentiation measures and produce high rates of false positives. Modern methods based on this principle (Table 3) correct for relatedness across populations, and can test association between allele frequencies and environmental features (see the extensive review by François et al., 2015). Methods such as baypass (Gautier, 2015) are convenient in both describing population structure and providing preliminary insights into the proportion of loci that do not follow neutral expectations. When this proportion is not too high, outliers can be removed to avoid bias  and the remaining loci can be used for demographic inference and model-testing. These estimated parameters can then be used to simulate sequences or independent SNPs and generate a neutral expectation. Loci that are more likely to be neutral can be used to further calibrate tests for selection (Lotterhos & Whitlock, 2014).
Detecting an association between environment and allele frequencies does not necessarily imply a role for local adaptation. For example, in the case of secondary contact, intrinsic genetic incompatibilities can lead to the emergence of tension zones that may shift until they reach an environmental barrier where they can be trapped (Bierne et al., 2011). In addition, the effects of selection at linked sites might generate false positives. The sampling strategy must take into account the particular historical and demographic features of the species investigated to gain power (Nielsen et al., 2007). The sequencing strategy must also be carefully considered to control for spatial autocorrelation of genotypes due to IBD and shared demographic history. For example, localized range expansion may produce a spurious association between environmental features and allele frequencies due to repeated founder effects and allele surfing (Excoffier & Ray, 2008). Including samples from populations not affected by such an expansion may avoid reaching biased conclusions by examining signatures of association at a broader scale.

| Identifying significant genotype-phenotype associations and and epistatic interactions between variants
The methods described above focus on allele frequencies at the population scale, but do not test association with traits that vary be- Uncovering the genetic basis of complex, polygenic traits remains challenging, even in model species (Pritchard & Di Rienzo, 2010;Rockman, 2012). It may be unavoidable as a first step to focus only on traits that are under relatively simple genetic determinism. This can, however, lead to the overrepresentation of loci of major phenotypic effect, a fact that should be acknowledged when discussing the impact of selection on genome variation. The fact that loci of major effect are the easiest to target does not imply that they are necessarily the main substrate of selection (Rockman, 2012). Association methods may help to target variants undergoing soft sweeps, weak selection or those involved in polygenic control of traits .
In such cases, signatures of selection may be subtle and sometimes difficult to retrieve from allele frequency data. Nevertheless, recent tools may have a higher sensitivity to polygenic selection (see section below on Ancestral recombination graphs), and a recent method uses genome-wide patterns of LD between a candidate gene (the "bait") and loci along the genome to detect candidate genes that may be involved in epistatic interactions (Boyrie et al., 2020). Such developments hold great promise in addressing the issue of nonadditive genetic effects.

| Inferring differences in history along the genome with ancestral recombination graphs
Ancestral recombination graphs (ARGs) are a generalization of the coalescent and describe the sequence of genealogies along a sample of recombining sequence. Genealogies are estimated for each nonrecombining block, and recombination between adjacent blocks is described by breaking the branch leading to the recombining haplotype and allowing it to recoalesce to the rest of the tree. This succession of local trees joined by recombination events provides a full description of the genealogical history of the data and is therefore a promising approach to characterize all modes of selection, introgression and demography while taking into account variation in recombination and mutation rate. Methods that are able to estimate or approximate these genealogies have long been computationally intensive and unapplicable to whole genomes.
Fortunately, recent improvements make their application to whole genomes feasible (Table 4). A good example is argweaver (Rasmussen et al., 2014), which has allowed candidate genes for long-term balancing selection to be recovered from human data, and has recently been used in combination with machine learning methods to study speciation in capuchino seedeater birds (Hejase et al., 2020) and introgression in humans (Kuhlwilm et al., 2016). However, argweaver remains slow when analysing more that 50 diploid genomes, and is even slower for low-depth or unphased data. Another promising method is implemented in the relate software (Speidel et al., 2019). relate is able to estimate genome-wide genealogies, and uses this information to reconstruct past demographic trajectories, changes in mutation rate over time, identify loci under positive selection, and estimate when selection acted on candidate mutations. relate can handle thousands of genomes in a manageable amount of time, but requires an outgroup sequence to polarize derived alleles, and a recombination map. relate comes with two add-on methods, clues and palm, which can use the relate output to estimate the strength of selection for single loci and a set of candidate loci identified by GWAS respectively. The latter method therefore provides a way to quantify polygenic selection and adaptive introgression at multiple loci, and constitutes a major advance in the field of population genomics.

| Jointly inferring demographic history and selection using ABC and supervised machine learning
It has recently become clear that the interactions between mutation and recombination rate, introgression, demography, selective sweeps and background selection have to be integrated into analyses of genetic variation (Andrew et al., 2013;Li et al., 2012;Ravinet et al., 2017). Simulation-based approaches hold great promise for incorporating this complexity. Two nonexclusive methodologies are promising: ABC and supervised machine learning approaches. Both rely on simulated data simulated under a range of parameters which are used to identify the combination of parameter values that are most likely to have generated the observed data (see Box 2 for a discussion on simulators). Both methodologies are flexible and powerful, and have become increasingly popular in population genomics (Csilléry et al., 2010;. These methods can accomodate any type of marker and arbitrarily complex models. By measuring the distance between carefully chosen summary statistics describing each simulation with those from the observed data set, it is possible to infer which combination of selective and demographic parameters best explains the data. These methods enable the rate of false positives to be estimated, for example by estimating how many times simulations of neutral sequences are classified as selected. However, using summary statistics leads to the loss of potentially useful information (Robert et al., 2011). Machine learning presents three advantages with respect to this problem. First, whereas ABC is prone to show lower performance as the number of statistics summarizing simulations increases, machine learning tends to display better performance. Second, ABC usually excludes many simulations by retaining only those closest to the observed data, whereas machine learning makes use of all simulations to form a model. Third, machine learning algorithms can identify the set of summary statistics that is the most useful for inference, and some neural networks algorithms may also be trained directly on images of aligned sequences (e.g., imagene, Table 4). An example of application of ABC and deep learning is provided in a study of African populations of Drosophila melanogaster (Sheehan & Song, 2016), in which the authors use these approaches to identify genomic regions under balancing and positive selection, and infer past demography.
The flexibility of ABC and supervised machine learning means that researchers can adapt the pipeline to their need, for example by using the simulator that generate simulations that are as close to the specifics of their study system as possible. However, this flexibility also means a steeper training curve for researchers learning these methods, due to the lack of a clearly unified analytical pipeline. Moreover, recent discussions on the power of machine learning pipelines to detect selective sweeps have highlighted that a careful consideration of demographic null models and methodological limitations is imperative (Harris et al., 2018).

| CON CLUS ION
As illustrated by sections 3.8 and 3.9, the field of population genomics is now moving towards both better integrating the demographic framework in inferences of selection, and, conversely, taking into account selection when reconstructing demographic history. The joint inference of loci under selection and quantification of demographic dynamics is of crucial importance in fields such as landscape genomics or the study of ongoing speciation. It might provide insights into the role of selection, recombination and gene flow in promoting or impairing local adaptation to new habitats. The growing availability of genome-wide data for nonmodel species is therefore promising, but requires caution and high stringency in our interpretation of observed patterns. With the decreasing cost of sequencing, it has been suggested that NGS will rapidly broaden our perspective on complex evolutionary processes, from biogeography (Lexer et al., 2013) to the genetic basis of traits (Hohenlohe, 2014) and the maintenance of polymorphisms (Hedrick, 2006). However, the study of DNA sequence variation is already challenging in its own right, and prone to storytelling (Pavlidis et al., 2012). In order to be informative about processes such as selection and demography, population genomics should ultimately be combined with other disciplines such as ecology and functional analyses (Habel et al., 2015). This can be achieved, for example, by assessing the function of selected genes, the consistency of demographic history with information retrieved from the fossil record or geological history, and the broader integration of population genomics with other fields and methods whenever possible, such as niche modelling, common garden experiments or the study of macro-evolutionary patterns of selection and diversification.

ACK N OWLED G EM ENTS
The University of Basel, New York University Abu Dhabi and the

AUTH O R CO NTR I B UTI O N
YB wrote the first draft of the manuscript and maintained the website. YB and BW agreed on the structure and revised the initial draft.

DATA AVA I L A B I L I T Y S TAT E M E N T
No data to provide. The tables shown in this article are available at www.methodspopgen.com