## Abstract

The advent of high throughput sequencing and genotyping technologies enables the comparison of patterns of polymorphisms at a very large number of markers. While the characterization of genetic structure from individual sequencing data remains expensive for many nonmodel species, it has been shown that sequencing pools of individual DNAs (Pool-seq) represents an attractive and cost-effective alternative. However, analyzing sequence read counts from a DNA pool instead of individual genotypes raises statistical challenges in deriving correct estimates of genetic differentiation. In this article, we provide a method-of-moments estimator of for Pool-seq data, based on an analysis-of-variance framework. We show, by means of simulations, that this new estimator is unbiased and outperforms previously proposed estimators. We evaluate the robustness of our estimator to model misspecification, such as sequencing errors and uneven contributions of individual DNAs to the pools. Finally, by reanalyzing published Pool-seq data of different ecotypes of the prickly sculpin *Cottus asper*, we show how the use of an unbiased estimator may question the interpretation of population structure inferred from previous analyses.

IT has long been recognized that the subdivision of species into subpopulations, social groups, and families fosters genetic differentiation (Wahlund 1928; Wright 1931). Characterizing genetic differentiation as a means to infer unknown population structure is therefore fundamental to population genetics and finds applications in multiple domains, including conservation biology, invasion biology, association mapping, and forensics, among many others. In the late 1940s and early 1950s, Malécot (1948) and Wright (1951) introduced *F*-statistics to partition genetic variation within and between groups of individuals (Holsinger and Weir 2009; Bhatia *et al.* 2013). Since then, the estimation of *F*-statistics has become standard practice (see, *e.g.*, Weir 1996, 2012; Weir and Hill 2002) and the most commonly used estimators of have been developed in an analysis-of-variance framework (Cockerham 1969, 1973; Weir and Cockerham 1984), which can be recast in terms of probabilities of identity of pairs of homologous genes (Cockerham and Weir 1987; Rousset 2007; Weir and Goudet 2017).

Assuming that molecular markers are neutral, estimates of are typically used to quantify genetic structure in natural populations, which is then interpreted as the result of demographic history (Holsinger and Weir 2009): large values are expected for small populations among which dispersal is limited (Wright 1951), or between populations that have long diverged in isolation from each other (Reynolds *et al.* 1983). When dispersal is spatially restricted, a positive relationship between and the geographical distance for pairs of populations generally holds (Slatkin 1993; Rousset 1997). It has also been proposed to characterize the heterogeneity of estimates across markers for identifying loci that are targeted by selection (Cavalli-Sforza 1966; Lewontin and Krakauer 1973; Beaumont and Nichols 1996; Vitalis *et al.* 2001; Akey *et al.* 2002; Beaumont 2005; Weir *et al.* 2005; Lotterhos and Whitlock 2014, 2015; Whitlock and Lotterhos 2015).

Next-generation sequencing (NGS) technologies provide unprecedented amounts of polymorphism data in both model and nonmodel species (Ellegren 2014). Although the sequencing strategy initially involved individually tagged samples in humans (The International HapMap Consortium 2005), whole-genome sequencing of pools of individuals (Pool-seq) is being increasingly used for population genomic studies (Schlötterer *et al.* 2014). Because it consists of sequencing libraries of pooled DNA samples and does not require individual tagging of sequences, Pool-seq provides genome-wide polymorphism data at considerably lower cost than sequencing of individuals (Schlötterer *et al.* 2014). However, non-equimolar amounts of DNA from all individuals in a pool and stochastic variation in the amplification efficiency of individual DNAs have raised concerns with respect to the accuracy of the so-obtained allele frequency estimates, particularly at low sequencing depth and with small pool sizes (Cutler and Jensen 2010; Anderson *et al.* 2014; Ellegren 2014). Nonetheless, it has been shown that, at equal sequencing efforts, Pool-seq provides similar, if not more accurate, allele frequency estimates than individual-based analyses (Futschik and Schlötterer 2010; Gautier *et al.* 2013). The problem is different for diversity and differentiation parameters, which depend on second moments of allele frequencies or, equivalently, on pairwise measures of genetic identity: with Pool-seq data, it is indeed impossible to distinguish pairs of reads that are identical because they were sequenced from a single gene from pairs of reads that are identical because they were sequenced from two distinct genes that are identical in state (IIS) (Ferretti *et al.* 2013).

Appropriate estimators of diversity and differentiation parameters must therefore be sought to account for both the sampling of individual genes from the pool and the sampling of reads from these genes. There has been several attempts to define estimators for the parameter for Pool-seq data (Kofler *et al.* 2011; Ferretti *et al.* 2013), from ratios of heterozygosities (or from probabilities of genetic identity between pairs of reads) within and between pools. In the following, we will argue that these estimators are biased (*i.e.*, they do not converge toward the expected value of the parameter) and that some of them have undesired statistical properties (*i.e.*, the bias depends on sample size and coverage). Here, following Cockerham (1969, 1973), Weir and Cockerham (1984), Weir (1996), Weir and Hill (2002), and Rousset (2007), we define a method-of-moments estimator of the parameter using an analysis-of-variance framework. We then evaluate the accuracy and precision of this estimator, based on the analysis of simulated data sets, and compare it to estimates defined in the software package PoPoolation2 (Kofler *et al.* 2011) and in Ferretti *et al.* (2013). Furthermore, we test the robustness of our estimators to model misspecifications (including unequal contributions of individuals in pools and sequencing errors). Finally, we reanalyze the prickly sculpin (*Cottus asper*) Pool-seq data (published by Dennenmoser *et al.* 2017), and show how the use of biased estimators in previous analyses may challenge the interpretation of population structure.

Note that throughout this article, we use the term “gene” to designate a segregating genetic unit (in the sense of the “Mendelian gene” from Orgogozo *et al.* 2016). We further use the term “read” in a narrow sense, as a sequenced copy of a gene. For the sake of simplicity, we will use the term “Ind-seq” to refer to analyses based on individual data, for which we further assume that individual genotypes are called without error.

## Model

*F*-statistics may be described as intraclass correlations for the IIS probability of pairs of genes (Cockerham and Weir 1987; Rousset 1996, 2007). is best defined as:(1)where is the IIS probability for genes sampled within subpopulations, and is the IIS probability for genes sampled between subpopulations. In the following, we develop an estimator of for Pool-seq data by decomposing the total variance of read frequencies in an analysis-of-variance framework. A complete derivation of the model is provided in the Supplemental Material, File S1.

For the sake of clarity, the notation used throughout this article is given in Table 1. We first derive our model for a single locus and eventually provide a multilocus estimator of *F*_{ST}. Consider a sample of subpopulations, each of which is made of genes sequenced in pools (hence is the haploid sample size of the *i*th pool). We define as the number of reads sequenced from gene *j* in subpopulation *i* at the locus considered. Note that is a latent variable that cannot be directly observed from the data. Let be an indicator variable for read *r* from gene *j* in subpopulation *i*, such that if the *r*th read from the *j*th gene in the *i*th deme is of type *k*, and otherwise. In the following, we use standard dot notation for sample averages, *i.e.*: and The analysis-of-variance is based on the computation of sums of squares, as follows:(2)As is shown in File S1, the expected sums of squares depend on the expectation of the allele frequency over all replicate populations sharing the same evolutionary history, as well as on the IIS probability that two genes in the same pool are both of type *k*, and the IIS probability that two genes from different pools are both of type *k*. Taking expectations (see the detailed computations in File S1), one has:(3)for reads within individual genes, since we assume that there is no sequencing error, *i.e.*, all the reads sequenced from a single gene are identical and for all *r*. For reads between genes within pools, we get:(4)where is the total number of reads in the full sample (total coverage), is the coverage of the *i*th pool, and arises from the assumption that the distribution of the read counts is multinomial (*i.e.*, that all genes contribute equally to the pool of reads; see Equation A15 in File S1). For reads between genes from different pools, we have:(5)where and (see Equation A16 in File S1). Rearranging Equation 4 and Equation 5 and summing over alleles, we get:(6)and(7)where Let and Then, using the definition of from Equation 1, we have:(8)which yields the method-of-moments estimator(9)where(10)and(11)(see Equations A25 and A26 in File S1). In Equation 10 and Equation 11, is the average frequency of reads of type *k* within the *i*th pool, and is the average frequency of reads of type *k* in the full sample. Note that from the definition of is the weighted average of the sample frequencies with weights equal to the pool coverage. This is equivalent to the weighted analysis-of-variance in Cockerham (1973) (see also Weir and Cockerham 1984; Weir 1996; Weir and Hill 2002; Rousset 2007; Weir and Goudet 2017). Finally, the full expression of in terms of sample frequencies develops as:If we take the limit case where each gene is sequenced exactly once, we recover the Ind-seq model: assuming for all then and Therefore, and Equation 9 reduces exactly to the estimator of for haploids: see Weir (1996), p. 182, and Rousset (2007), p. 977.

As in Reynolds *et al.* (1983), Weir and Cockerham (1984), Weir (1996), and Rousset (2007), a multilocus estimate is derived as the sum of locus-specific numerators over the sum of locus-specific denominators:(13)where and are subscripted with *l* to denote the *l*th locus. For Ind-seq data, Bhatia *et al.* (2013) refer to this multilocus estimate as a “ratio of averages” as opposed to an “average of ratios,” which would consist of averaging single-locus over loci. This approach is justified in the appendix of Weir and Cockerham (1984) and in Bhatia *et al.* (2013), who analyzed both estimates by means of coalescent simulations. Note that Equation 13 assumes that the pool size is equal across loci. Also note that the construction of the estimator in Equation 13 is different from Weir and Cockerham’s (1984). These authors defined their multilocus estimator as a ratio of sums of components of variance (*a*, *b*, and *c* in their notation) over loci, which give the same weight to all loci whatever the number of sampled genes at each locus. Equation 13 follows GENEPOP’s rationale (Rousset 2008) instead, which gives more weight to loci that are more intensively covered.

## Materials and Methods

### Simulation study

#### Generating individual genotypes:

We first generated individual genotypes using ms (Hudson 2002), assuming an island model of population structure (Wright 1931). For each simulated scenario, we considered eight demes, each made of haploid individuals. The migration rate (*m*) was fixed to achieve the desired value of (0.05 or 0.2), using equation 6 in Rousset (1996) leading to, *e.g.*, for and for The mutation rate was set at giving We considered either fixed or variable sample sizes across demes. In the latter case, the haploid sample size *n* was drawn independently for each deme from a Gaussian distribution with mean 100 and SD 30; this number was rounded up to the nearest integer, with a minimum of 20 and maximum of 300 haploids per deme. We generated a very large number of sequences for each scenario and sampled independent single nucleotide polymorphisms (SNPs) from sequences with a single segregating site. Each scenario was replicated 50 times (500 times for Figure 3 and Figure S2).

#### Pool sequencing:

For each ms simulated data set, we generated Pool-seq data by drawing reads from a binomial distribution (Gautier *et al.* 2013). More precisely, we assume that for each SNP, the number of reads of allelic type *k* in pool *i* follows:(14)where is the number of genes of type *k* in the *i*th pool, is the total number of genes in pool *i* (haploid pool size), and is the simulated total coverage for pool *i*. In the following, we either consider a fixed coverage, with for all pools and loci, or a varying coverage across pools and loci, with

#### Sequencing error:

We simulated sequencing errors occurring at rate which is typical of Illumina sequencers (Glenn 2011; Ross *et al.* 2013). We assumed that each sequencing error modifies the allelic type of a read to one of three other possible states with equal probability (there are therefore four allelic types in total, corresponding to four nucleotides). Note that only biallelic markers are retained in the final data sets. Also note that, since we initiated this procedure with polymorphic markers only, we neglect sequencing errors that would create spurious SNPs from monomorphic sites. However, such SNPs should be rare in real data sets, since markers with a low minimum read count (MRC) are generally filtered out.

#### Experimental error:

Nonequimolar amounts of DNA from all individuals in a pool and stochastic variation in the amplification efficiency of individual DNAs are sources of experimental errors in Pool-seq. To simulate experimental errors, we used the model derived by Gautier *et al.* (2013). In this model, it is assumed that the contribution of each gene *j* to the total coverage of the *i*th pool follows a Dirichlet distribution:(15)where the parameter *ρ* controls the dispersion of gene contributions around the value which is expected if all genes contributed equally to the pool of reads. For convenience, we define the experimental error *ϵ* as the coefficient of variation of *i.e.*, (see Gautier *et al.* 2013). When *ϵ* tends toward 0 (or equivalently, when *ρ* tends to infinity), all individuals contribute equally to the pool and there is no experimental error. We tested the robustness of our estimates to values of *ϵ* between 0.05 and 0.5. The case could correspond, for example, to a situation where (for ) five individuals contribute more reads than the other five individuals.

### Other estimators

For the sake of clarity, a summary of the notation of the estimators used throughout this article is given in Table 2.

#### PP2_{d}:

This estimator of is implemented by default in the software package PoPoolation2 (Kofler *et al.* 2011). It is based on a definition of the parameter as the overall reduction in average heterozygosity relative to the total combined population (see, *e.g.*, Nei and Chesser 1983):(16)where is the average heterozygosity within subpopulations, and is the average heterozygosity in the total population (obtained by pooling together all subpopulations to form a single virtual unit). In PoPoolation2, is the unweighted average of within-subpopulation heterozygosities:(17)(using the notation from Table 1). Note that in PoPoolation2, is restricted to the case of two subpopulations only (). The two ratios in the right-hand side of Equation 17 are presumably borrowed from Nei (1978) to provide an unbiased estimate, although we found no formal justification for the expression in Equation 17 for Pool-seq data. The total heterozygosity is computed as (using the notation from Table 1):

#### PP2_{a}:

This is the alternative estimator of provided in the software package PoPoolation2. It is based on an interpretation by Kofler *et al.* (2011) of Karlsson *et al.*’s (2007) estimator of , as:(19)where and are the frequencies of identical pairs of reads within and between pools, respectively, computed by simple counting of IIS pairs. These are estimates of the IIS probability for two reads in the same pool (whether they are sequenced from the same gene or not), and the IIS probability for two reads in different pools. Note that the IIS probability is different from in Equation 1, which, from our definition, represents the IIS probability between distinct genes in the same pool. This approach therefore confounds pairs of reads within pools that are identical because they were sequenced from a single gene from pairs of reads that are identical because they were sequenced from distinct, yet IIS genes.

#### FRP_{13}:

This estimator of was developed by Ferretti *et al.* (2013) (see their equations 3, 10, 11, 12, and 13). Ferretti *et al.* (2013) use the same definition of as in Equation 16 above, although they estimate heterozygosities within and between pools as “average pairwise nucleotide diversities,” which, from their definitions, are formally equivalent to IIS probabilities. In particular, they estimate the average heterozygosity within pools as (using the notation from Table 1):(20)and the total heterozygosity among the populations as:

### Analyses of Ind-seq data

For the comparison of Ind-seq and Pool-seq data sets, we computed on subsamples of 5000 loci. These subsamples were defined so that only those loci that were polymorphic in all coverage conditions were retained, and the same loci were used for the analysis of the corresponding Ind-seq data. For the latter, we used either the Nei and Chesser’s (1983) estimator based on a ratio of heterozygosity (see Equation 16 above), hereafter denoted by or the analysis-of-variance estimator developed by Weir and Cockerham (1984), hereafter denoted by

All the estimators were computed using custom functions in the R software environment for statistical computing, version 3.3.1 (R Core Team 2017). All of these functions were carefully checked against available software packages to ensure that they provided strictly identical estimates.

### Application example: *C. asper*

Dennenmoser *et al.* (2017) investigated the genomic basis of adaption to osmotic conditions in the prickly sculpin (*C. asper*), an abundant euryhaline fish in northwestern North America. To do so, they sequenced the whole genome of pools of individuals from two estuarine populations (Capilano River Estuary, CR; Fraser River Estuary, FE) and two freshwater populations (Pitt Lake, PI; Hatzic Lake, HZ) in southern British Columbia (Canada). We downloaded the four corresponding BAM files from the Dryad Digital Repository (http://dx.doi.org/10.5061/dryad.2qg01) and combined them into a single mpileup file using SAMtools version 0.1.19 (Li *et al.* 2009) with default options, except the maximum depth per BAM that was set to 5000 reads. The resulting file was further processed using a custom awk script to call SNPs and compute read counts, after discarding bases with a base alignment quality (BAQ) score <25. A position was then considered a SNP if: (1) only two different nucleotides with a read count >1 were observed (nucleotides with read being considered as a sequencing error); (2) the coverage was between 10 and 300 in each of the four alignment files; (3) the minor allele frequency, as computed from read counts, was in the four populations. The final data set consisted of 608,879 SNPs.

Our aim here was to compare the population structure inferred from pairwise estimates of using the estimator (Equation 12) with that of PP2_{d}. To determine which of the two estimators performs better, we then compared the population structure inferred from and to that inferred from the Bayesian hierarchical model implemented in the software package BayPass (Gautier 2015). BayPass allows the robust estimation of the scaled covariance matrix of allele frequencies across populations for Pool-seq data, which is known to be informative about population history (Pickrell and Pritchard 2012). The elements of the estimated matrix can be interpreted as pairwise and population-specific estimates of differentiation (Coop *et al.* 2010) and therefore provide a comprehensive description of population structure that makes full use of the available data.

### Data availability

An R package called poolfstat, which implements estimates for Pool-seq data, is available at the Comprehensive R Archive Network (CRAN): https://cran.r-project.org/web/packages/poolfstat/index.html.

The authors state that all data necessary for confirming the conclusions presented in this article are fully represented within the article, figures, and tables. Supplemental material (including Figures S1–S4, Tables S1–S3, and a complete derivation of the model in File S1) available at Figshare: https://doi.org/10.25386/genetics.6856781.

## Results

### Comparing Ind-seq and Pool-seq estimates of

Single-locus estimates of are highly correlated with the classical estimates of (Weir and Cockerham 1984) computed on the individual data that were used to generate the pools in our simulations (see Figure 1). The variance of across independent replicates decreases as the coverage increases. The correlation between and is stronger for multilocus estimates (see Figure S1A).

### Comparing Pool-seq estimators of

We found that our estimator has extremely low bias (<0.5% over all scenarios tested: see Table 3 and Tables S1–S3). In other words, the average estimates across multiple loci and replicates closely equal the expected value of the parameter, as given by equation 6 in Rousset (1996), which is based on the computation of IIS probabilities in an island model of population structure. In all the situations examined, the bias does not depend on the sample size (*i.e.*, the size of each pool) or on the coverage (see Figure 2). Only the variance of the estimator across independent replicates decreases as the sample size increases and/or as the coverage increases. At high coverage, the mean and root mean squared error (RMSE) of over independent replicates are virtually indistinguishable from that of the estimator (see Table S1).

Figure 3 shows the RMSE of estimates for a wide range of pool sizes and coverages. The RMSE decreases as the pool size and/or the coverage increases. The estimates are more precise and accurate when differentiation is low. Figure 3 provides some clues to evaluate the pool size and the coverage that is necessary to achieve the same RMSE as for Ind-seq data. Consider, for example, the case of samples of haploids. For (in the conditions of our simulations), the RMSE of estimates based on Pool-seq data tends to the RMSE of estimates based on Ind-seq data either by sequencing pools of ∼200 haploids at 20×, or by sequencing pools of 20 haploids at ∼200×. However, the same precision and accuracy are achieved by sequencing ∼50 haploids at ∼50×.

Conversely, we found that (the default estimator of implemented in the software package PoPoolation2) is biased when compared to the expected value of the parameter. We observed that the bias depends on both the sample size and the coverage (see Figure 2). We note that, as the coverage and the sample size increase, converges to the estimator (Nei and Chesser 1983) computed from individual data (see Figure S1B). This argument was used by Kofler *et al.* (2011) to validate their approach, even though the estimates of depart from the true value of the parameter (Figure S1, B and C).

The second of the two estimators of implemented in PoPoolation2, which we refer to as is also biased (see Figure 2). We note that the bias decreases as the sample size increases. However, the bias does not depend on the coverage (only the variance over independent replicates depends on coverage). The estimator developed by Ferretti *et al.* (2013), which we refer to as is also biased (see Figure 2). However, the bias does not depend on the pool size or on the coverage (only the variance over independent replicates depends on coverage). converges to the estimator computed from individual data (see Figure 2). At high coverage, the mean and RMSE over independent replicates are virtually indistinguishable from that of the estimator.

Lastly, we stress that our estimator provides estimates for multiple populations and is therefore not restricted to pairwise analyses, contrary to PoPoolation2’s estimators. We show that, even at low sample size and low coverage, Pool-seq estimates of differentiation are virtually indistinguishable from classical estimates for Ind-seq data (see Table 3).

### Robustness to unbalanced pool sizes and variable sequencing coverage

We evaluated the accuracy and the precision of the estimator when sample sizes differ across pools and when the coverage varies across pools and loci (see Figure 4). We found that, at low coverage, unequal sampling or variable coverage causes a negligible departure from the median of estimates computed on individual data, which vanishes as the coverage increases. At 100× coverage, the distribution of estimates is almost indistinguishable from that of (see Figure 4 and Tables S2 and S3).

### Robustness to sequencing and experimental errors

Figure 5 shows that sequencing errors cause a negligible negative bias for estimates. Filtering (using an MRC of 4) improves estimation slightly, but only at high coverage (Figure 6B). It must be noted, however, that filtering increases the bias in the absence of sequencing error, especially at low coverage (Figure 6A). With experimental error, *i.e.*, when individuals do not contribute evenly to the final set of reads, we observed a positive bias for estimates (Figure 5). We note that the bias decreases as the size of the pools increases. Figure S2 shows the RMSE of estimates for a wider range of pool sizes, coverage, and experimental error rate (*ϵ*). For increasing the coverage cannot improve the quality of the inference if the pool size is too small. When Pool-seq experiments are prone to large experimental error rates, increasing the size of pools is the only way to improve the estimation of Filtering (using an MRC of 4) does not improve estimation (Figure 6C).

### Application example

The reanalysis of the prickly sculpin data revealed larger pairwise estimates of multilocus using the estimator, as compared to (see Figure 7A). Furthermore, we found that estimates are smaller for within-ecotype pairwise comparisons as compared to between-ecotype comparisons. Therefore, the inferred relationships between samples based on pairwise estimates show a clear-cut structure, separating the two estuarine samples from the freshwater ones (see Figure 7C). We did not recover the same structure using estimates (see Figure 7B). Additionally, the scaled covariance matrix of allele frequencies across samples is consistent with the structure inferred from estimates (see Figure 7D).

## Discussion

Whole-genome sequencing of pools of individuals is increasingly popular for population genomic research on both model and nonmodel species (Schlötterer *et al.* 2014). The development of dedicated software packages (reviewed in Schlötterer *et al.* 2014) undoubtedly has something to do with the breadth of research questions that have been tackled using Pool-seq. However, the analysis of population structure from Pool-seq data are complicated by the double sampling process of genes from the pool and sequence reads from those genes (Ferretti *et al.* 2013).

The naive approach that consists of computing from read counts as if they were allele counts (*e.g.*, as in Chen *et al.* 2016) ignores the extra variance brought by the random sampling of reads from the gene pool during Pool-seq experiments. Furthermore, such computation fails to consider the actual number of lineages in the pool (haploid pool size). Altogether, these limits may result in severely biased estimates of differentiation when the pool size is low (see Figure S3). A possible alternative is to compute from allele counts imputed from read counts using a maximum-likelihood approach conditional on the haploid size of the pools (*e.g.*, as in Smadja *et al.* 2012; Leblois *et al.* 2018), or from allele frequencies estimated using a model-based method which accounts for the sampling effects and the sequencing error probabilities inherent to pooled NGS experiments (see Fariello *et al.* 2017). However, these latter approaches may only be accurate in situations where the coverage is much larger than pool size, allowing for a reduction of the sampling variance of reads (see Figure S3). We therefore developed a new estimator of the parameter for Pool-seq data in an analysis-of-variance framework (Cockerham 1969, 1973). The accuracy of this estimator is barely distinguishable from that of the Weir and Cockerham’s (1984) estimator for individual data. Furthermore, it does not depend on the pool size or on the coverage, and it is robust to unequal pool sizes and varying coverage across demes and loci.

In our analysis, the frequency of reads within pools is a weighted average of the sample frequencies, with weights equal to the pool coverage. Therefore, our approach follows Cockerham’s (1973) one, which he referred to as a weighted analysis-of-variance (see also Weir and Cockerham 1984; Weir 1996; Weir and Hill 2002; Weir and Goudet 2017). With unequal pool sizes, weighted and unweighted analyses differ. As discussed recently in Weir and Goudet (2017), the unweighted approach seems appropriate when the between component exceeds the within component, *i.e.*, when is large (Tukey 1957). It turns out that optimal weighting depends upon the parameter to be estimated (Cockerham 1973) and is only efficient at lower levels of differentiation (Robertson 1962). In a likelihood analysis of the island model, Rousset (2007) derived asymptotically efficient weights that are proportional to for the sum of squares of different samples (see also Robertson 1962). To the best of our knowledge, such optimal weighting has never been considered in the literature.

### Analysis-of-variance and probabilities of identity

In the analysis-of-variance framework, is defined in Equation 1 as an intraclass correlation for the probability of IIS (Cockerham and Weir 1987; Rousset 1996). Extensive statistical literature is available on estimators of intraclass correlations. Beside analysis-of-variance estimators, introduced in population genetics by Cockerham (1969, 1973), estimators based on the computation of probabilities of identical response within and between groups have been proposed (see, *e.g.*, Fleiss 1971; Fleiss and Cuzick 1979; Mak 1988; Ridout *et al.* 1999; Wu *et al.* 2012), which were originally referred to as kappa-type statistics (Fleiss 1971; Landis and Koch 1977). These estimators have later been endorsed in population genetics, where the “probability of identical response” was then interpreted as the frequency with which the genes are alike (Cockerham 1973; Cockerham and Weir 1987; Weir 1996; Rousset 2007; Weir and Goudet 2017).

This suggests that, with Pool-seq data, another strategy could consist of computing from IIS probabilities between (unobserved) pairs of genes, which requires that unbiased estimates of such quantities are derived from read count data. We have done this in the second section of File S1 and we provide alternative estimators of for Pool-seq data (see Equations A44 and A48 in File S1). These estimators (denoted by and ) have exactly the same form as the analysis-of-variance estimator if the pools all have the same size and if the number of reads per pool is constant (Equation A33 in File S1). This echoes the derivations by Rousset (2007) for Ind-seq data, who showed that the analysis-of-variance approach (Weir and Cockerham 1984) and the simple strategy of estimating IIS probabilities by counting identical pairs of genes provide identical estimates when sample sizes are equal (see Equation A28 in File S1 and also Cockerham and Weir 1987; Weir 1996; Karlsson *et al.* 2007). With unbalanced samples, we found that analysis-of-variance estimates have better precision and accuracy than IIS-based estimates, particularly for low levels of differentiation (see Figure S4). Interestingly, we found that IIS-based estimates of for Pool-seq data have generally lower bias and variance if the overall estimates of IIS probabilities within and between pools are computed as unweighted averages of population-specific or pairwise estimates (see Equations A39 and A43 in File S1), as compared to weighted averages (Equations A46 and A47 in File S1). Equation A28 in File S1 further shows that our estimator may be rewritten as a function close to except that it also depends on the sum in both the numerator and the denominator. This suggests that if the ’s differ among subpopulations, then our estimator provides an estimate of an average of population-specific (Weir and Hill 2002; Weir and Goudet 2017).

It follows from the derivations in File S1 that the estimator (Equation 19) is biased because the IIS probability between pairs of reads within a pool is a biased estimator of the IIS probability between pairs of distinct genes in that pool (see Equations A34–A36 in File S1). This is the case because the former confounds pairs of reads that are identical because they were sequenced from a single gene from pairs of reads that are identical because they were sequenced from distinct, yet IIS genes.

A more justified estimator of has been proposed by Ferretti *et al.* (2013), based on previous developments by Futschik and Schlötterer (2010). Note that, although they defined as a ratio of functions of heterozygosities, they actually worked with IIS probabilities (see Equation 20 and Equation 21). However, although Equation 20 is strictly identical to Equation A39 in File S1, we note that they computed the total heterozygosity by integrating over pairs of genes sampled both within and between subpopulations (compare Equation 21 with Equation A43 in File S1), which may explain the observed bias (see Figure 2).

### Comparison with alternative estimators

An alternative framework to Weir and Cockerham’s (1984) analysis-of-variance has been developed by Masatoshi Nei and coworkers to estimate from gene diversities (Nei 1973, 1977, 1986; Nei and Chesser 1983). The estimator (see Equation 16, Equation 17, and Equation 18) implemented in the software package PoPoolation2 (Kofler *et al.* 2011) follows this logic. However, it has long been recognized that both frameworks are fundamentally different in that the analysis-of-variance approach considers both statistical and genetic (or evolutionary) sampling, whereas Nei and coworkers’ approach do not (Weir and Cockerham 1984; Excoffier 2007; Holsinger and Weir 2009). Furthermore, the expectation of Nei and coworkers’ estimators depend on the number of sampled populations, with a larger bias for lower numbers of sampled populations (Goudet 1993; Excoffier 2007; Weir and Goudet 2017). This is the case because the computation of the total diversity in Equation 18 and Equation 21 includes the comparison of pairs of genes from the same subpopulation, whereas the computation of IIS probabilities between subpopulations do not (see, *e.g.*, Excoffier 2007). Therefore, we do not recommend using the estimator implemented in the software package PoPoolation2 (Kofler *et al.* 2011).

### Applications in evolutionary ecology studies

Pool-seq is being increasingly used in many application domains (Schlötterer *et al.* 2014), such as conservation genetics (see, *e.g.*, Fuentes-Pardo and Ruzzente 2017), invasion biology (see, *e.g.*, Dexter *et al.* 2018), and evolutionary biology in a broader sense (see, *e.g.*, Collet *et al.* 2016). These studies use a large range of methods, which aim at characterizing fine-scaled population structure (see, *e.g.*, Fischer *et al.* 2017), reconstructing past demography (see, *e.g.*, Chen *et al.* 2016; Leblois *et al.* 2018), or identifying footprints of natural or artificial selection (see, *e.g.*, Chen *et al.* 2016; Fariello *et al.* 2017; Leblois *et al.* 2018).

Here, we reanalyzed the Pool-seq data produced by Dennenmoser *et al.* (2017), who investigated the adaptive genomic divergence between freshwater and brackish-water ecotypes of the prickly sculpin *C. asper*, an abundant euryhaline fish in northwestern North America. Measuring pairwise genetic differentiation between samples using , we found a clear-cut structure separating the freshwater from the brackish-water ecotypes. Such genetic structure supports the hypothesis that populations are locally adapted to osmotic conditions in these two contrasted habitats, as discussed in Dennenmoser *et al.* (2017). This structure, which is at odds with that inferred from estimates, is not only supported by the scaled covariance matrix of allele frequencies, but also by previous microsatellite-based studies, which showed that populations were genetically more differentiated between ecotypes than within ecotypes (Dennenmoser *et al.* 2014, 2015).

### Limits of the model and perspectives

We have shown that the stronger source of bias for the estimate is unequal contributions of individuals in pools. This is because we assume in our model that the read counts are multinomially distributed, which supposes that all genes contribute equally to the pool of reads (Gautier *et al.* 2013), *i.e.*, that there is no variation in DNA yield across individuals and that all genes have equal sequencing coverage (Rode *et al.* 2018). Because the effect of unequal contribution is expected to be stronger with small pool sizes, it has been recommended to use Pool-seq with at least 50 diploid individuals per pool (Lynch *et al.* 2014; Schlötterer *et al.* 2014). However, this limit may be overly conservative for allele frequency estimates (Rode *et al.* 2018) and we have shown here that we can achieve very good precision and accuracy of estimates with smaller pool sizes. Furthermore, because genotypic information is lost during Pool-seq experiments, we assume in our derivations that pools are haploid (and therefore that is nil). Analyzing nonrandom mating populations (*e.g.*, in selfing species) is therefore problematic.

Finally, our model, as in Weir and Cockerham (1984), formally assumes that all populations provide independent replicates of some evolutionary process (Excoffier 2007; Holsinger and Weir 2009). This may be unrealistic in many natural populations, which motivated Weir and Hill (2002) to derive a population-specific estimator of for Ind-seq data (see also Vitalis *et al.* 2001). Even though the use of Weir and Hill’s (2002) estimator is still scarce in the literature (but see Weir *et al.* 2005; Vitalis 2012), Weir and Goudet (2017) recently proposed a reinterpretation of population-specific estimates of in terms of allelic matching proportions, which are strictly equivalent to IIS probabilities between pairs of genes. It is therefore straightforward to extend Weir and Goudet’s (2017) estimator of population-specific for the analysis of Pool-seq data, using the unbiased estimates of IIS probabilities provided in File S1.

## Acknowledgments

We thank Alexandre Dehne-Garcia for his assistance in using computer farms. We thank two anonymous reviewers for their positive comments and suggestions. Analyses were performed on the GenoToul bioinformatics platform Toulouse Midi-Pyrénées (http://bioinfo.genotoul.fr) and the High Performance Computational platform of the Centre de Biologie pour la Gestion des Populations. This work is part of V.H.’s Ph.D.; V.H. was supported by a grant from the Institut National de la Recherche Agronomique’s Plant Health and Environment (SPE) Division and by the BiodivERsA project EXOTIC (ANR-13-EBID-0001). Part of this work was supported by the project SWING (ANR-16-CE02-0015) of the French National Research Agency, and by the CORBAM project of the French region Hauts-de-France.

## Footnotes

Supplemental material available at Figshare: https://doi.org/10.25386/genetics.6856781.

*Communicating editor: M. Beaumont*

- Received March 9, 2018.
- Accepted July 21, 2018.

- Copyright © 2018 by the Genetics Society of America