Hacdivsel: a Program Implementing a New Haplotype-f St Composite Method for the Detection of Divergent Selection in Pairs of Populations of Non-model Species

Running title: HacDivSel: detection of divergent selection. CC-BY-NC-ND 4.0 International license not peer-reviewed) is the author/funder. It is made available under a The copyright holder for this preprint (which was. Abstract In this work two complementary methods for detection of divergent selection between pairs of populations connected by migration are introduced. The new statistics do not require knowledge on the ancestral or derived allelic state and are robust to false positives. Additionally, it is not necessary to perform neutral simulations to obtain critical cutoff values for the identification of candidates. The utility of the methods is demonstrated in simulations. First method, called nvdF ST , combines information from the haplotype patterns with inter-population differences in allelic frequency. Remarkably, this is not a F ST outlier test because highest F ST s are not a priori selected to identify loci. On the contrary, candidate loci are chosen based upon haplotype allelic class metric and then the F ST for these loci are estimated and compared to the overall F ST. Evidence of divergent selection is concluded only when both the haplotype pattern and the F ST value indicate that outcome. It is shown that powers ranging from 70-90% are achieved in many of the scenarios assayed while the false positive rate is controlled below the desired nominal level ( = 0.05). Additionally, the method is also robust to demographic scenarios involving population bottleneck and expansion. The second method is developed for cases with independently segregating markers, where the information contained in the haplotypes vanishes. In this case, the power to detect selection is attained by developing a new F ST extreme-outlier test based on a k-means clustering algorithm. Both kinds of strategies have been implemented in the program HacDivSel.. CC-BY-NC-ND 4.0 International license not peer-reviewed) is the author/funder. It is made available under a The copyright holder for this preprint (which was .


Introduction
Current population genetics has a main focus in the detection of the signatures of selection at the molecular level.In previous years, the main effort was focused in human and other model organisms but now, the increasing amount of information on genomes of non-model species also permits to focus the search in many other situations of interest.One of these includes the sought for local adaptation and selection in structured populations.By nonmodel species we mean a species for which we lack a priori information on candidate loci with known function that could be potentially adaptive.As well there is not outgroup so that the allelic state, ancestral or derived, is unknown.There are several methods aiming to detect selection in genome regions undergoing local adaptation.Some of them are based on finding outlier loci when measuring genetic differentiation between populations.From its original formulation (Lewontin& Krakauer 1973) this technique has been both, questioned and improved in many different ways (Akey 2009;Akey et al. 2002;Bonhomme et al. 2010;Chen et al. 2010;Duforet-Frebourg et al. 2014;Excoffier et al. 2009;Fariello et al. 2013;Foll& Gaggiotti 2008).
Different problems have been encountered regarding to the use of F ST measures as a guide in the search for selection.For example, under the infinite island model, the effect of gene flow and the corresponding correlations in the gene frequencies among local subpopulations could inflate the neutral variance in F ST , leading to high rate of false positives ( Bonhomme et al. 2010).Moreover, several non-local processes, as background selection, species-wide selective sweeps or bottleneck and population expansions scenarios, can also produce F ST outliers (Bierne et al. 2013;Maruki et al. 2012).Finally, F ST methods are not designed to work-well when there are many loci under selection (Bierne et al. 2013;Li et al. 2012).All of the above provokes that, still under optimal conditions, some of the methods tend to produce many false positives (Perez-Figueroa et al. 2010).
Recently, some promising alternatives have appeared to deal with these issues (Bonhomme et al. 2010;Duforet-Frebourg et al. 2014;Fariello et al. 2013;Frichot et al. 2013) allowing for more accurate identification of loci under local divergent selection (De Villemereuil et al. 2014;Lotterhos& Whitlock 2014).However, besides the fact that these methods come with stronger computational cost, they still present some caveats that difficult their use in exploratory studies with non-model organisms.For example, it is necessary to known the allelic state, ancestral or derived, or having information on the diploid genotype, or the methods are dependent on some a priori assumptions on the outlier model and/or need an outgroup to account for the population structure.
Furthermore, if two populations connected by migration are under divergent selection, the selected alleles will be at more or less intermediate frequencies.Depending on the interplay between selection, migration and drift, a critical migration threshold will exist to maintain polymorphism (Yeaman& Otto 2011;Yeaman& Whitlock 2011).In this case, the detection of selection would be difficult using frequency spectrum methods but see (Achaz 2009), it would be also difficult under most of the F ST -based existing approaches, but see hapFLK (Fariello et al. 2013).However, there is no program specialized to detect divergent selection at genomic or sub-genomic level, while working without information on the structure of the population tree, the state, ancestral or derived, of the alleles and on the genotypes.
The goal of the present work is to develop such a program, for detecting ongoing divergent selection in pairs of populations with gene flow.This new method should be adequate for working with non-model species.It would also desirable that selection could be detected without simulating any neutral demographic scenario.Moreover, it should work at the genomic level with automatic decision-making for the window size and must be protected from false positives.The latter is one of the main problems associated to the genome-wide selection detection methods.The algorithm we will describe achieves these goals, at least on the several scenarios assayed; being powerful for detecting single or polygenic selection when two populations connected by migration undergo divergent selection.
The design of the work is as follows: In the first part, the model is developed, which includes the computation of a normalized variance difference, nvd, for detecting haplotype patterns under selection.The nvd value is then used to perform an F ST index measure.Additionally, the algorithm for a conservative, extreme positive outlier test, to deal with the case of fully unlinked SNPs, is also developed.
These methods are implemented in the computer program HacDivSel.
Second, the simulation setting for testing the methods is explained.
Finally, the results concerning power and false positive rate (declared significant versus true nulls) are given and discussed.

Generalized HAC variance difference
For a given sample, let define the major-allele-reference haplotype (MARH) as the haplotype carrying only major frequency alleles of its constituting SNPs (Hussin et al. 2010).Define the mutational distance between any haplotype and MARH as the number of sites (SNPs) in the haplotype carrying a non-major allele.An haplotype allelic class (HAC) groups the haplotypes having the same mutational distance.Therefore, the HAC of a given haplotype corresponds to the number of non-major (i.e.minor) alleles it carries, so that every haplotype having the same number of minor alleles belongs to the same HAC class.
Given the definitions above, consider a sample of n haplotypes of length L SNPs.For each evaluated SNP i (i∈[1,L]) we can perform a partition of the HAC classes into P 1 , the subset of HACs for the haplotypes carrying the most frequent or major allele at the SNP i under evaluation and P 2 the subset with the remaining haplotypes carrying the minor allele at i.
That is, let '0' to be the major allele for the SNP i and accordingly '1' is the minor allele.Then, P 1 includes every haplotype carrying the allele '0' for the SNP i and P 2 the remaining haplotypes carrying '1' for that SNP.In P 1 we have different HAC values depending on the distance of each haplotype from MARH and similarly in P 2 .In each subset we can compute the variance of the HACs.That is, in P 1 we have the variance v 1i and correspondingly variance v 2i in P 2 .The rationale of the HAC-based methods is that if the SNP i is under ongoing selection then the variance in the partition 1 will tend to be zero because the allele at higher frequency (i.e. in the partition 1) should be the favored one and the sweeping effect will make the HAC values in this partition to be lower (close to the MARH) and consequently provoking lower variance values.The variance in the second partition should not be affected by the sweeping effect because it does not carry the favored allele.So, the difference v 2i -v 1i would be highly positive in the presence of selection and not so otherwise (Hussin et al. 2010).For a window size of L SNPs, the variance difference between P 2 and P 1 can be computed to obtain a summary statistic called Svd (Hussin et al. 2010) that can be immediately generalized to Where, f is the frequency of the derived allele, and the parameters a and b permit to give different weights depending on if it is desired to detect higher frequencies (a = 0) or more intermediate ones (a > 0) of the derived allele.If a = 0 and b = 1 the statistic corresponds to the original Svd and if a =1 and b = 4 it corresponds to the variant called SvdM (Rivas et al. 2015).Note that when taking a = 1 it is not necessary to distinguish between ancestral and derived alleles because f and 1-f are interchangeable.
A drawback in the gSvd statistic is its dependence on the window size L as has already been reported for the original Svd (Hussin et al. 2010;Rivas et al. 2015).Although gSvd is normalized by L, the effect of the window size on the computation of variances is quadratic (see Appendix for details) which explains why the normalization is not effective in avoiding a systematic increase of the statistic under larger window sizes.This impact of the window size is important because the two different partitions may experience different scaling effects, which would increase the noise in the estimation.Additionally, the change in the scale due to the window size will be dependent on the recombination rate and on the effect of selection.Thus, it is desirable to develop a HAC-based statistic not dependent on the window size.In what follows, the between partition variance difference will be reworked in order to develop a new normalized HAC-based statistic, specially focused on detecting divergent selection in local adaptation scenarios with migration.
Note that, for any given sample of size n, the corresponding means and variances at each partition are related via the general mean and variance in that sample.Consider m, m 1 , m 2 the mean HAC distances in the sample and in the partitions P 1 and P 2 respectively, for any given candidate SNP.We have the following relationships for the mean m and sample variance S 2 values (the subscripts identify the partition, see Appendix for details) with  ̅ = ( 1 −1) 1 2 +( 2 −1) 2 2  −1 ; n 1 and n 2 are the sample sizes (n 1 > n 2 by definition) and Using the relationships in (1), the between partitions variance difference can be recomputed and some non-informative term discarded (see details in the Appendix) to finally obtain a new statistic for the variance difference Note that (2) will augment with decreasing S 1 and increasing S 2 (because the latter increases S).Therefore, if selection is favoring the major allele of the SNP i, then the variance S 1 2 will tend to zero (Hussin et al. 2010;Rivas et al. 2015), the sample variance S 2 will be a function of the variance S 2 2 and the value in (2) will be positive.Because we are interested in detecting intermediate allele frequencies (see below), the parameters a and b from gSvd has been substituted by a = 1 and b = 4 as these are the values that permit to ignore the allelic state while maximizing (2) for intermediate frequencies.

Variance upper bound and normalized variance difference
Note that HAC values vary in the range [0, L] which provokes that the sample variance S 2 has an upper bound at nL 2 / [4(n-1)].Then the maximum variance difference occurs when f i = 0.5, S 1 2 = 0, n 2 = n/2 and by substituting in (2) we get an upper bound If we divide (2) by the right side in (3) we have a normalized variance difference The quantity from (4) can be computed for each SNP in a sample of sequences of a given length L and the SNP giving the maximum nvd can be considered as a candidate for selection.Furthermore, it is possible to compute (4) for each population or to combine the two populations in a unique sample.The latter is better for our purpose of looking for divergent selection in populations undergoing gene flow.When pooling both populations the frequencies tend to be intermediate in the selective sites.Therefore, we compute the normalized variance difference in (4) for the data obtained by merging the shared SNPs from the two population samples.Note however that the reference haplotype (MARH) is computed just from one of the populations (by default population 1).
Recall that ( 4) is already normalized by the square of the window size L.However, the problem remains, as with other genome-wide estimators based on sliding-window approaches, on how to choose an optimal window size when computing the quantities of interest.A solution to this problem is to automate the choice of the window size by selecting the size which gives the maximum value for the statistic (Rivas et al. 2015).Therefore, we focus on the candidate having maximum nvd for every SNP and window size.
At this point we already have a HAC-based statistic, nvd, that is independent of the window size and that should produce higher positive values for pairs of populations undergoing divergent selection.However, if selection is not at work, the maximum nvd value would be a false positive.Unfortunately, we ignore the distribution of the statistic and cannot decide if a given nvd value is supporting the hypothesis of selection or not.As well we might not have enough information on the species to simulate its evolution under a given neutral demography.Therefore, we still need to identify if the value obtained for a given sample is due to the effect of selection.In doing so, we will compute two more measures before giving a diagnostic about the presence of divergent selection.The first is a sign test based on the lower bound of (4), the second is the F ST of the SNP having the maximum nvd compared with the global F ST .

Sign test
From a lower bound of (4) we derive the quantity (see Appendix for details) where hac 1i are the HAC values measured at each haplotype i in the partition 1 and the sum is over the n 1 sequences in that partition.A negative sign in (5) suggests that the value of nvd is not the result of divergent selection.Indeed, we require (5) to be positive to count a given candidate as significant.
Combined method: nvdF ST The sign test defined above is a good strategy for discarding some false candidates.However, we still lack a method for obtaining p-values associated to the sites chosen by the nvd algorithm.Performing neutral simulations to assess significance in a given data set could be challenging especially with non-model organisms because the demographic history and related evolutionary parameters are not well known.As an alternative, it is possible to combine the information on candidate SNPs as given by nvd with the interpopulation differentiation measure at that site.The significance of the obtained quantity is far easier to asses.The joint use of these methods produces the combined measure nvdF ST .The rationale of the approach is that if divergent selection acts on an specific site then the F ST at that site will be higher compared to the overall F ST .To obtain a p-value an obvious strategy would be to perform an LK test (Lewontin& Krakauer 1973).However, LK and its derivatives have several problems and a tendency to produce high rates of false positives (De Mita et al. 2013;De Villemereuil et al. 2014;Lotterhos& Whitlock 2014).Instead, we proceed as follows, let i to be a candidate site chosen because it has the maximum nvd value, then we measure the index I i = F STi -F ST comparing the candidate with the overall F ST .To get the pvalue for a given I i , the data is resampled several times to generate an empirical distribution.
In doing so, the mean frequency for every SNP in the pooled populations is considered as this is the expectation under the homogenizing effect of migration provided that Nm > 1 (Crow& Kimura 1970).Then, for any iteration, the probability of a given allele at each population is obtained from a binomial B(p,n), where p is the mean allelic frequency at that site and n the local population sample size.The p-values correspond to the proportion of times that the resampled indexes were larger than I i .Thus, a low p-value indicates that the observed I i was larger than the resampled ones most of the times.Note that, for each site, the resampling procedure has variance pqn which will be larger at intermediate frequencies.
For a candidate with a high mean frequency (i.e. both populations has high frequency) the value of I i is low and the p-value will be high.If the mean frequency is low the situation is similar, I i is low and the p-value high.When the mean frequency is intermediate two situations are possible, first, each population has similar intermediate frequencies which again imply high p-values or alternatively, the frequencies can be extreme and opposite at each population.In the latter, I i is high and its p-value low.Recall that we are seeking for selection in populations connected by migration and working only with SNPs shared between them.Thus, the SNPs that are fixed in one of the populations are not considered.
Note that if Nm < 1 it is likely that the shared SNPs are scarce.Therefore, by considering only shared SNPs, we avoid to test data with Nm < 1 that could incur in false positives.
. The F ST values were computed following the algorithm in Ferretti et al (Ferretti et al. 2013).
The number of resamplings was set to 500 times.

Multiple test correction by the effective number of independent SNPs, sign test and significance
We have computed nvd and the F ST index and got a candidate site with its p-value.Since nvd was obtained after testing a number of positions on a given window size, it is desirable to apply a multiple test correction for the number of independent SNPs in the window.To roughly estimate the number of independent SNPs, we calculate the linkage disequilibrium measure D' (Devlin& Risch 1995;Lewontin 1988) at each pair of consecutive sites and then store the quantity r' = 1 -|D'| for each pair.The effective number of independent SNPs (M effs ) between site w ini and w end is then obtained as one plus the summation of the r' values in the interval [w ini , w end ).
The Šidák correction (Cheverud 2001;Sidak 1967) can now be applied to get the corrected significance level c = 1 -(1 -) 1/Meffs with nominal level  = 0.05.Thus, the algorithm nvdF ST ends with a candidate being considered significant only when its p-value (computed as explained in the previous section) is lower than c and the sign in (5) is positive.

The k-means extreme positive outlier (EPO) test
The nvdF ST method assumes the existence of a dense map of linked genetic markers.If the data consists mostly in independent markers this would provoke the failure to detect selection by the nvdF ST method because the HAC-based information will not be accessible.
To deal with this situation, a second method was implemented consisting in a heuristic twostep procedure that performs a conservative test for identifying extreme outliers.
We intend our method to be conservative because, as mentioned above, the variance of the F ST distribution is quite unpredictable under a variety of scenarios.This provokes high rates of false positives associated with the F ST outlier tests.Therefore, our test takes advantage on the fact that the involved regions under divergent selection may produce a kind of extreme outliers that would be clustered apart from the neutral ones.Only when this kind of outliers is detected the subsequent LK test is performed.This approach is heuristic, i.e. a problem solving method that finds a satisfactory (but not necessarily optimal) solution to a difficult problem.
The rationale of the algorithm is as follows: The first step consists in computing the extreme positive outliers (EPO) in the sense of Tukey i.e. those sites having a F ST value higher than 3IQR where IQR is the interquartile range (Tukey 1977).The second step identifies different classes inside the extreme outlier set.This is done by a k-means algorithm which is a well-known method widely used in clustering and digital image processing (Schubert et al. 2012;Vattani 2011).Here, a k-modal distribution is assumed and all the elements of the EPO set are classified in one of the k classes.The class with lower values is discarded and only the elements, if any, in the upper classes having values higher than the adjacent inter-mode middle point are maintained in the set.By default k = 2 and two modes {0, 1} were used.Finally, for each of the candidates remaining in the EPO set a LK test (Lewontin& Krakauer 1973) is performed to compute its p-value.The Šidák correction (Cheverud 2001;Sidak 1967) for the number of remaining outliers in the set is applied to get the significance level.

Software description
Both nvdF ST and the EPO test have been implemented in the program HacDivSel.Complete details of the software can be found in the accompanying manual.We here just mention that the input program files must be in MS-format (Hudson 2002)

Simulations and analysis
There are several examples of adaptation to divergent environments connected by migration such as the intertidal marine snail L. saxattilis (Rolan-Alvarez 2007), wild populations of S. salar (Bourret et al. 2013), lake whitefish species (Renaut et al. 2011) and so on.To perform simulations as real as possible, a model resembling the most favorable conditions for the formation of ecotypes under local adaptation with gene flow was implemented.Some relevant demographic information from L. saxatilis, such as migration rates and population size as estimated from field data (Rolan-Alvarez 2007), was used.Concerning selection intensities, we considered moderate selection pressures and few loci with large effects (Sadedin et al. 2009;Thibert-Plante & Gavrilets 2013).Therefore, the simulation design includes both a single selective locus model plus some simulations under a polygenic architecture with 5 selective loci.Two populations of facultative hermaphrodites were simulated under divergent selection and migration.Each individual consisted of a diploid chromosome of length 1Mb.The contribution of each selective locus to the fitness was 1-hs with h = 0.5 in the heterozygote or h = 1 otherwise (Table 1).In the polygenic case the fitness was obtained by multiplying the contribution at each locus.In both populations the most frequent initial allele was the ancestral.The selection coefficient for the ancestral allele was always s = 0 while s = ± 0.15 for the derived.That is, in population 1 the favored allele was the derived (negative s, i.e. 1 + h|s| in the derived) which was at initial frequency of 10 -3 while in the other population the favored was the ancestral (positive s, i.e. 1 -h|s| in the derived) and was initially fixed.
Table 1.Fitness Model.The ancestral allele is noted with uppercase A and the derived as a.

AA Aa aa
In the single locus model the selective site was located at different relative positions 0, 0.01, 0.1, 0.25 and 0.5.In the polygenic model the positions of the five sites were 4×10 -6 , 0.2, 0.5, 0.7 and 0.9.Under both architectures, the overall selection pressure corresponded to  = 4Ns = 600 with N = 1000.Simulations were run in long term scenarios during 5,000 and 10,000 generations and in short-term scenarios during 500 generations.Some extra cases with weaker selection  = 140 (s = ± 0.07, N = 500) in the long-term (5,000 generations) and stronger selection,  = 6000 (s = ± 0.15, N = 10,000) in the short-term were also run.
The mating was random within each population.The between population migration was Nm In most scenarios, the number of SNPs in the data ranged between 100 and 500 per Mb.
However, only the SNPs shared between populations were considered thus giving numbers between 60-300 SNPs per Mb i.e. medium to high density SNP maps.
The interplay between divergent selection, drift and migration (Yeaman& Otto 2011) under the given simulation setting should permit that the adaptive divergence among demes persists despite the homogeneity effects of migration (see Critical migration threshold section in the Appendix).

Combined method (nvdF ST )
Under the single locus architecture, the power of nvdF ST vary between 80-90% for both medium (60 SNPs/Mb) and high density (250 SNPs/Mb) maps.(Table 2).These results can be compared with previously published analysis (Rivas et al. 2015) of the same files corresponding to rows 1-6 in Table 2.In the previous study the methods Svd, SvdM and OmegaPlus (Alachiotis et al. 2012) were evaluated with best powers of 63-79% obtained by Svd and SvdM in cases with high mutation and recombination (Table 2 in Rivas et al. 2015).
When these methods were applied considering the two populations merged then the best powers were attained by SvdM ranging from 42 to 94% (Table 3 in Rivas et al. 2015).Recall that the methods Svd, SvdM and OmegaPlus oblige the user to perform simulations of a neutral demography to obtain the p-values for the tests.As it can be appreciated from rows 1 to 6 in Table 2, nvdF ST performs quite well with powers from 81 to 93% without the need of performing additional neutral simulations.The given results are for 10,000 generations; the results for 5,000 generations were quite similar and are therefore omitted.
Under the polygenic architecture (n = 5 in Table 2) the power is also good (80-99%).At least one candidate was found 99% of the times and more than one were found 80% of the time.
However, the number of correctly identified sites was quite variable ranging between 1 and 3.
An exception in the power attained by nvdF ST occurs when all SNPs segregate independently (last row in Table 2).In this case, the method fails to detect selection which is not surprising because the information from the haplotype allelic classes is absent under linkage equilibrium; the adequate patterns are not found which provokes both a negative sign in the test and a candidate with low F ST index measure.

Extreme Positive Outlier (EPO) test
The EPO test is very conservative as can be observed in Table 3 where the false positive rate is 0 in every case.Its power increases with the density and the independence of the markers reaching 60% of detection in the case of independent SNPs and maps with 250-300 SNPs/Mb.As expected for an outlier test, the power undergoes a breakdown under a polygenic setting (row with n = 5 in Table 3).Therefore, the EPO test is complementary to nvdF ST having its maximum power when the latter has its minimum and viceversa.

Bottleneck-expansion scenarios
The robustness of the methods was additionally tested under a bottleneck-expansion scenario that is known to leave signatures that mimic the effect of positive selection.Thus, no selection was applied to this scenario in order to test its effect on the false positive rate of the methods.As a result, for nvdF ST the false positive rate is maintained below the nominal level (4.8%) and for the EPO test it continues to be 0%.Therefore, both methods seem to be robust to these kinds of confounding scenarios.

Short-term strong and long-term weak selection scenarios
The performance under the strong selection scenario is presented in Table 6.Not surprisingly, the number of segregating sites is considerably reduced.In fact the minimum window size allowed by the program had to be shortened from 51 to 25 to perform the analyses.The power of detection range between 46-66% with 0 false positive rate.These results can be compared with gSvd methods, Svd and SvdM, in Table 6 (t = 500 generations) from Rivas et al. (Rivas et al. 2015).The results in Rivas et al. were more dependent on the recombination rate having low powers (14-28%) under full linkage and great power (70-96%) with high recombination.Recall however that to asses significance with these methods the exact neutral demography was simulated by Rivas and coworkers.
Concerning very weak selection in long-term scenarios (Table 6, α = 140) the power varied between 47-52% and localization between 2-32kb from the real selective position.

Discussion
The goal in this work was to develop a method for detection of divergent selection in pairs of populations connected by migration with the requisite of being protected from false positives which is a known concern for differentiation-based methods (De Mita et al. 2013;De Villemereuil et al. 2014;Lotterhos& Whitlock 2014).Additionally, the model should be useful for non-model species and it should not be necessary to perform neutral simulations to obtain critical cut-off values for the candidates.
It has been shown that combining an haplotype-based method with the F ST differentiation measure, the so-called nvdF ST , is quite powerful strategy for detecting divergent selection.
However, when the whole set of markers is segregating independently, the information in the haplotype vanishes.Therefore, nvdF ST is complemented with a new F ST outlier test called EPO that was developed as a conservative test because the mentioned tendency of such methods to produce false positives.Its conservativeness is attained by assuming that among the detected outliers, the ones which are due to divergent selection would be clustered apart from neutral ones.The implemented extreme positive outlier test behaves acceptably well when markers are independent, reaching powers between 60-80% while maintaining a false positive rate of virtually zero.
All the developed methods work without any a priori knowledge about candidate positions neither any information about the state, ancestral or derived, of the alleles.Consequently, they seem an interesting option for exploratory studies in non-model species.
In general, the F ST -based methods are affected by the presence of polygenic scenarios (Bierne et al. 2013;De Villemereuil et al. 2014) because those tests are specifically designed for finding larger than average F ST values which are difficult to discover if the overall mean is high.The nvdF ST performs even better in this scenario because the distributed selective signal facilitates the discovery of the corresponding patterns.Since only the F ST of the specific site is compared with the overall F ST and the null distribution is obtained using inter-population mean frequencies, the nvdF ST maintains high power under the polygenic setting.
Nevertheless, the localization of various selective sites is not easy at least under the sample size assayed.
Besides the detection of the signal of selection in the sequences, the accuracy in inferring the location of the selected allele was also studied.When the target site is located at the centre of the studied region and the overall recombination rate is at least 0.3 cM/Mb (ρ ≥ 12), the nvdF ST method performs quite well under weak selection, with the inferred location within 32 kb of distance from the true location in the worst case.Under strong selection, the localization is worst, 77 kb, but this could be due to the lower number of segregating sites (only 62 in Table 6).The performance is poorer when the target is not at the center.In this case, with recombination of 1.5 cM/Mb, the inferred location changes from an almost perfect localization (<1 kb) to distances of 10-120 kb as the target is shifted away.Fortunately, the complementary EPO test works well, at least when the overall recombination rate is above 1.5 cM/Mb, giving detection powers near 80% and seems not to be so affected for such positional effect, resulting in distances from the true selective position between <1 -36 kb.
Clearly, inferring the selective location is still a pending issue for many of the selection detection methods.There is also plenty of room for improvement under the nvdF ST and EPO methods in this regard, for example, trying to further explore the relationship between recombination and the window sizes producing the highest scores.Indeed, the interplay among divergent selection, recombination, drift and migration should be considered for further improving the efficiency of the methods.
. As a conclusion, nvdF ST combines haplotype and population differentiation information for the detection of divergent selection and seems to work well when knowledge of the haplotype phase is at hand.The complementary EPO method is a conservative alternative useful when the full set of SNPs is unlinked.Both strategies can be applied without the need of performing neutral simulations and have false positive rates below the desired nominal level.
Crecimiento, GPC2013-011) and fondos FEDER.The author declares to have no conflict of interest.
higher or lower than (A-8) depending on the HAC values of the first partition.In any case, a negative value in (A-9) may be caused by m 1 being equal or higher than m 2 and suggests, whether it be n 1 = n/2 or not, that the value of nvd is not the result of divergent selection.
Indeed, we require (A-9) to be positive to count a given candidate as significant.

Critical migration threshold
Our simulation model is a particular case (with symmetric migration and intermediate dominance) of the model in Yeaman and Otto (2011) that these authors develop to study the interplay between drift, divergent selection and migration on the maintenance of polymorphism between the interconnected populations.Thus, we can compute the critical migration threshold below which adaptive divergence among demes is likely to persist.By rearranging terms in equation ( 11) from Yeaman and Otto (2011) after substituting the fitness relationships from our system, we obtain the critical migration threshold: (A-10) where  = 4Ns .For each selective pressure , we can therefore compute the critical number of migrants (Nm crit ) below which the selective polymorphism should be present in the data.
The weaker the selection the lower the threshold so, for  = 140 the minimum critical number of migrants is 355 individuals.Thus our highest migration Nm = 50 is far below the threshold.This means that both scenarios Nm = 10 and 50, would tend to maintain the locally adaptive allele for every selective scenario assayed (weak, intermediate and strong) despite the homogeneity effects of migration.
and should contain sequence samples from two populations.A typical command line for calling the program to analyze a file named sel.txtcontaining 50 sequences from each population would be HacDivSel -name sel.txt -sample 50 -candidates 10 -SL 0.05 -output anyname Where -sample is the sample size for each population and the label -candidates 10, indicates that the ten highest nvd values must be included in the output.The program would analyze the file and produce as output the highest 10 values and its significance at the 0.05 level for different window sizes after the nvdF ST test.It also performs the EPO test and gives the candidate outliers, if any, and their significance.Only the SNPs shared by the two populations are considered.Which imply that there are at least 4 copies of each SNP in the metapopulation.

=
10 plus some cases with Nm = 0 or Nm = 50 in a short-term scenario.Recombination ranged from complete linkage between pairs of adjacent SNPs (no recombination, ρ = 0), intermediate values ρ = 4Nr = {4, 12, 60} and fully independent SNPs.A bottleneck-expansion scenario was also studied consisting in a neutral case with equal mutation and recombination rates, θ = ρ = 60, and a reduction to N = 10 in one of the populations in the generation 5,000 with the subsequent expansion following a logistic growth with rate 2 and K max = 1000.For every selective case, 1000 runs of the corresponding neutral model were simulated.To study the false positive rate (FPR) produced by the selection detection tests, the significant results obtained in the neutral cases were counted.The simulations were implemented using the last version of the program GenomePop2 (Carvajal-Rodriguez 2008).

Table 2 .
Performance of the combined method (nvdF ST ) with n true selective sites.Selection was  = 600 and Nm = 10.Mean localization is given in distance kb from the real selective position.

Table 3 .
Performance of the extreme outlier test (EPO) with n true selective sites.Selection was  = 600 and Nm = 10.Mean localization is given in distance kb from the real selective position.

Table 5 .
Performance of the combined method (nvdF ST ) with a single selective site in the short term (500 generations).Selection was  = 600 and Nm = 50.Mean localization is given in distance kb from the real selective position.Mean number of shared SNPs per Mb.θ: Mutation rate.ρ: Recombination rate.FPR: false positive rate. ∑:

Table 6 .
Performance of the combined method (nvdF ST ) with a single selective site in the short-term strong ( = 6000) and the long-term weak ( = 140) selection scenarios.Nm was 10.Mean localization is given in distance kb from the real selective position.
∑: Mean number of shared SNPs per Mb.θ: Mutation rate.ρ: Recombination rate.t: number of generations.FPR: false positive rate.*: only 40 runs having a minimum of 25 SNPs.