Information

Measuring genetic distance: $F_{ST}$ vs. Nei's distance

Measuring genetic distance: $F_{ST}$ vs. Nei's distance



We are searching data for your request:

Forums and discussions:
Manuals and reference books:
Data from registers:
Wait the end of the search in all databases.
Upon completion, a link will appear to access the found materials.

As far as I'm aware, Nei's genetic distance is quite old compared to $F_{ST}$. However, I have recently read more papers that utilized Nei's genetic distance alongside with $F_{ST}$. As I'm not very familiar with Nei, what are some advantages it has over $F_{ST}$?

Does Nei's genetic distance suffer from ascertainment bias?


Relationship between three measures of genetic differentiation GST, DEST and G’ST: how wrong have we been?

Table S1 Studies included in this meta-analysis.

Please note: Wiley-Blackwell are not responsible for the content or functionality of any supporting materials supplied by the authors. Any queries (other than missing material) should be directed to the corresponding author for the article.

Filename Description
MEC_4185_sm_TableS1.doc113 KB Supporting info item

Please note: The publisher is not responsible for the content or functionality of any supporting information supplied by the authors. Any queries (other than missing content) should be directed to the corresponding author for the article.


A SIMPLE TEST ON ALLELE SIZE INFORMATION CONTENT

The test indicates whether allele sizes provide information on population differentiation given a data set, that is, whether shifts in allele sizes resulting from step-wise-like mutations contribute to population differentiation. Contribution of stepwise-like mutations to genetic differentiation requires (1) that the mutation process is at least partially SMM-like and (2) that the mutation rate, μ, is large enough relative to the effect of drift and migration (e.g., μ≥ m otherwise new mutations are quickly spread beyond their native population by migration). Table 2 outlines the null hypotheses that can be tested, presenting a general null hypothesis as well as specific null hypotheses holding under particular prior assumptions.

The principle of the test is based on obtaining a distribution of a statistic under the null hypothesis (H0) that differences in allele sizes do not contribute to population differentiation. Therefore, we use a randomization procedure whereby the different allele sizes observed at a locus for a given data set are randomly permuted among allelic states. To better figure out the procedure, one may dissociate allelic state, identified, for example, by a letter (e.g., a, b, c, d, and e if there are five different alleles), and allele size, identified by a number (e.g., 4, 5, 7, 8, and 11, each representing the number of sequence repeats), given that there is a one-to-one correspondence between allelic state and allele size. Before randomization, the allele size attributed to each allelic state is the actual allele size (e.g., a, 4 b, 5 c, 7 d, 8 and e, 11). Throughout the randomization procedure, genotypes are defined in terms of allelic states and are not modified, but allele sizes are randomly reassigned among allelic states (e.g., a, 7 b, 4 c, 11 d, 5 and e, 8). After such a randomization, any two genes originally having the same allele size remain identical, although it can be for another allele size, whereas any two genes originally bearing different alleles of small size difference may bear alleles of large size difference, or reciprocally. Hence, the allele identity information is kept intact but not the allele size information. Under the null hypothesis (Table 2, case 1), the randomization procedure should not affect the expectation of a measure of differentiation such as RST. On the contrary, if allele sizes contribute to genetic differentiation, the RST computed after allele size permutation (hereafter called pRST) would depend solely on allele identity/nonidentity and hence have a smaller expectation than the value computed before randomization. The test can thus be designed by comparing the observed RST value (before randomization) to the distribution of pRST values obtained for all possible configurations of allele size permutations (or a representative subset of them, as the total number of different configurations quickly becomes enormous when the number of alleles exceeds 7 or 8). From this comparison, a probability that the null hypothesis holds can be estimated as the proportion of pRST values larger than the observed RST (one-tailed test). Note that the mean pRST should equal in expectation the FST computed on the same data (not accounting for potential statistical bias), as is confirmed later.

On a single locus, such a test can be applied only if a sufficient number of different alleles (n) are in the data set, as the number of different permutation configurations is equal to n!. Hence, five alleles (120 different configurations) appear to be a minimum to carry out such test at a type I error rate criterion of 5 or 1%. On a multilocus RST estimate, the test can be carried out by permuting allele sizes within each locus. It is noteworthy that the test makes no assumptions on the mutation model: A significant result (RST significantly >pRST) suggests that mutations contributed to genetic differentiation (e.g., because μ≥ m in an island model) and that the mutation process follows at least partially a SMM (the test remains valid under deviations from the SMM). Neutrality with respect to natural selection is, however, assumed. When the test is significant, FST is likely to provide a biased estimate of gene flow parameters, but it cannot be concluded a priori that RST would necessarily perform better given its larger variance (which is even more pronounced when mutations of more than one step can occur Z hivotovsky and F eldman 1995) and given the bias it may suffer when the mutation process deviates from the assumptions of the GSM (E stoup and A ngers 1998). A nonsignificant result (RST not significantly different from pRST) would suggest that allele size is not informative for population differentiation, because the mutation process is not step-wise-like and/or because mutations had not contributed to differentiation (e.g., because μ⪡ m in an island model). In this case, FST should surely be preferred over RST (although it would not ensure that FST provides a correct estimate of gene flow given the many other sources of bias related to population models W hitlock and M c C auley 1999).

Hypotheses tested by allele size permutations applied on RST

Which hypotheses can be tested and with which statistics? Simulations permit validation of the allele size permutation test and assess its power. But it is first necessary to insist on what can be tested (Table 2).

Randomizing allele sizes creates replicates of a data set for a mutation process following a KAM (or IAM) because, under this model, allele size is irrelevant and interchanging them is like replicating the past mutation processes leading to the present data set but with other randomly chosen alleles after each mutational event. Hence, one possible application of the allele size randomization procedure is to test whether the mutation process follows a KAM (Table 2, case 3). For this purpose, randomizing allele sizes can be applied on any statistic based on allele size, not only R-statistics but also various genetic distances for stepwise mutation models such as (δμ) 2 (e.g., G oldstein et al. 1995b S hriver et al. 1995), or simply on the total variance in allele size. It is, however, already well established that the large majority of microsatellite loci do not conform to a KAM, and the interesting question about the mutation process of microsatellites is rather how it deviates from an ideal SMM (E stoup and A ngers 1998). Therefore, using the allele size permutation procedure to test for the KAM is not discussed further.

A second application of the allele size permutation procedure, here assuming a priori that mutations follow at least partially a SMM-like process, is to test whether mutation has contributed to population divergence (Table 2, case 2). In other words, we can test whether the migration rate (m) among populations, or the reciprocal of the number of generations (t) since population divergence, is large compared to the mutation rates (μ⪡ m or μ⪡ 1/t, respectively Table 2, cases 2a and 2b). The allele size permutation test is the most interesting to address this question, because there is enough evidence that most microsatellites follow a SMM-like process (e.g., E llegren 2000 X u et al. 2000 Z hu et al. 2000 R enwick et al. 2001). However, for this purpose, allele size permutation cannot be applied to any statistic based on allele size: It performs well on R-statistics, which are ratios of allele size variance components, but not on genetic distances such as the G oldstein et al. (1995a) (δμ) 2 statistic, which is a between-populations component of allele size variance. The reason is that random permutations of allele sizes not only remove the within-population covariance between allele sizes for different alleles, but also modify the allele size variance under SMM or GSM, because the expected frequency distribution of allele sizes is not uniform (D onnelly 1999). Statistics expressing a component of allele size variance, such as the (δμ) 2 statistic, will always be affected by a change of the allele size variance, no matter whether or not mutations contributed to differentiation. On the contrary, statistics based on a ratio of variance components, such as RST, will not be affected if the within- and among-populations components of variance are multiplied by factors having the same expectations. The simulations presented hereafter show that this is what occurs when there is no within-population covariance between allele sizes for different alleles (i.e., differentiation due to drift and not stepwise mutations).

To show that the allele size permutation test is adequate for the RST statistic but not the (δμ) 2 statistic when testing m ⪢ μ or 1/t ⪢ μ (under the a priori assumption that the mutation process is stepwise-like Table 2, cases 2), we simulated a random-mating population of diploid individuals (population size N = 1000 individuals) at mutation-drift equilibrium (μ= 0.001) under the SMM. The allele size permutation test (1000 randomizations) was then applied on RST and (δμ) 2 computed between two independent samples (sample size n = 100 individuals) from that population for each of 200 simulated loci (the two samples thus represent undifferentiated subpopulations). The computer programs used for simulations and computations are described below. We report the percentage of loci for which the tests were significant (%RHo) according to the type I error rate criterion (α, the probability of rejecting the null hypothesis when it is true). Because the null hypothesis to be tested (1/t ⪢μ) is met by simulations, a valid testing procedure must ensure that %RHo =α otherwise it means that the procedure is not adequate to test this null hypothesis. Figure 1 shows that the allele size randomization testing procedure is indeed valid when applied on RST but not on (δμ) 2 .

Power of the test under SMM: To investigate the power of the test when testing if mutations contributed to population differentiation under the SMM (Table 2, cases 2), we checked the procedure on artificial data sets with realistic sample sizes derived from Monte Carlo simulations of populations made of diploid hermaphrodites. Three sets of demographic situations were simulated: (1) an island model at drift-migration-mutation equilibrium, (2) a model of two isolated populations having diverged from a common ancestral population at mutation-drift equilibrium, and (3) a linear stepping-stone model (gene flow restricted to adjacent populations) at drift-migration-mutation equilibrium. The island model was composed of 10 populations, consisting of 100 individuals each, and new generations were obtained by drawing genes at random from the population with probability 1 - m or from the other populations with probability m. The isolated population model was composed of two random-mating populations, consisting of 500 individuals each, and having diverged for t generations. The stepping-stone model was composed of 30 aligned populations, consisting of 50 individuals each, and new generations were obtained by drawing genes at random from the population with probability 1 - m or from the two adjacent populations with probability m.

—Control of the validity of the allele size permutation test when applied on RST (□) or (δμ) 2 (▵) statistics computed between two samples from a population at mutation-drift equilibrium under the SMM. The percentage of loci with the null hypothesis rejected (%RHo) is shown as a function of the type I error rate criterion (α), and the dashed line shows the %RHo =α relationship expected under the null hypothesis for a valid testing procedure. The null hypothesis of interest is whether the mutation rate is negligible, given that the mutation process is stepwise-like (Table 2, case 2). Results show that the allele size permutation procedure applied on (δμ) 2 is not suited to test this hypothesis.

The genetic parameters simulated were the following: At the initial stage all populations were fixed for one allele 10 loci were simulated with mutations following a SMM and μ= 10 -3 at all loci without size constraints. Simulations were run for a sufficient time to reach a steady state for total- and within-population gene diversity parameters, and then a sample of individuals representative of common experimental studies was extracted and analyzed. To obtain accurate estimates, 200 replicates were run for each set of conditions. Simulations were carried out using the software EASYPOP ver. 1.7.4 (B alloux 2001). Allele size permutation tests (with 1000 randomizations) and computations of FST and RST on the samples extracted were done with the program SPAGeDi (H ardy and V ekemans 2002). Single-locus and multilocus FST and RST were estimated following W eir and C ockerham (1984) and M ichalakis and E xcoffier (1996), respectively. It should be noted that this RST (an estimator of the parameter called ρST by R ousset 1996) differs somewhat from S latkin ’s (1995) original definition (M ichalakis and E xcoffier 1996) but is better suited for comparison with the FST estimator of W eir and C ockerham (1984) (called θ by these authors) and for demographic parameter estimations (R ousset 1996). Both these FST and RST estimators proceed by a standard hierarchical ANOVA where the observed variance (σ 2 ) of allele identity per locus and per allele (FST), or the variance of allele size per locus (RST), is partitioned into three components (random effects): among populations ( σ a 2 ) , among individuals within population ( σ b 2 ) , and between genes within individual within population ( σ c 2 ) . FST and RST are then estimated as σ a 2 ∕ ( σ a 2 + σ b 2 + σ c 2 ) (single-locus RST) or Σ σ a 2 ∕ Σ ( σ a 2 + σ b 2 + σ c 2 ) , where the summations apply over all loci (multilocus RST), all alleles of a locus (single-locus FST), or all alleles and loci (multilocus FST E xcoffier 2001).

For the island model, simulations were run for 5000 generations with migration rates among populations varying from 10 -4 to 10 -1 (i.e., m = 0.1-100μ) according to the runs. Global RST, FST, and pRST (for 1000 randomizations) were computed on a total sample of 300 individuals (30 individuals from each population). For the isolated populations model, a single population of 1000 individuals was simulated for 5000 generations, and then it was divided into two isolated subpopulations of 500 individuals that were run for 30-10,000 additional generations (i.e., 1/t = 0.1-33μ). RST, FST, and pRST (for 1000 randomizations) were computed on a total sample of 100 individuals (50 individuals from each subpopulation). For the stepping-stone model, 10,000 generations were simulated with a migration rate of 0.1 (0.05 between any two adjacent populations). Analyses were carried out on a sample of 20 individuals from each of the 30 populations (total sample size of 600 individuals). Pairwise FST/(1 - FST) and RST/(1 - RST) ratios were computed for each pair of populations, and these values were averaged over all pairs separated by 1, 2, 3. 20 steps (20 distance classes). Allele size permutation tests were applied on averaged pairwise RST/(1 - RST) ratios per distance class to provide pRST/(1 - pRST) values per distance class (1000 permutations). Here, pairwise FST/(1 - FST) and RST/(1 - RST) ratios were computed because theory predicts an approximate linear relationship with the linear distance between populations in one-dimensional isolation-by-distance models (R ousset 1997).

The validity of some of the simulation results could be verified by comparing them to theoretical expectations. For example, after 5000 generations of simulation of a single population of N = 1000 individuals (for the isolated population model), the average heterozygosity and average variance of allele size were equal to He = 0.68 and V = 1.96, respectively, with a mean number of alleles per locus of 5.8 (range, 3-11 alleles). These values are close to their expectations at mutation-drift equilibrium (E stoup and C ornuet 1999): Under strict SMM, He = 1 - (1 + 8Nμ) -0.5 = 0.67 and V = 2Nμ= 2. In the island model with 10 populations of 100 individuals each (d = 10, N = 100), average RST values were equal to 0.019, 0.197, 0.677, and 0.924 for m = 10 -1 , 10 -2 , 10 -3 , and 10 -4 , respectively (Figure 2A), in agreement with the expected values approximately equal to 1/(1 + 4Nm d/(d - 1)) = 0.022, 0.184, 0.692, and 0.957, respectively (R ousset 1996). In the isolated populations model (N = 500), divergence time t can be estimated from the relationship RST/(1 - RST) = t/2N (S latkin 1995 R ousset 1996), giving estimates of t = 97, 1132, and 11,301 for actual values of 100, 1000, and 10,000 generations, respectively. Finally, in the linear stepping-stone model (N = 50, m = 0.1), pairwise RST/(1 - RST) values increased linearly with the distance between populations (Figure 2C), giving a regression slope equal to 0.054, in agreement with the approximate expected value 1/(4Nm) = 0.050 for the linear stepping-stone model (R ousset 1997).

Results from all simulations confirm that mean pRST values (i.e., mean value computed after random permutations of allele size) are very close, though not exactly equal, to the FST values (Figure 2). For example, in the island model, the mean and standard deviation of the difference between FST and mean pRST values per locus were equal to 0.003 ± 0.007, 0.008 ± 0.012, and 0.010 ± 0.110 for m = 10 -2 , 10 -3 , and 10 -4 , respectively. Hence, mean pRST values were on average slightly lower than FST values although, for a given locus, the difference between the two could be quite substantial, especially under very low migration rates. For the other simulations, mean pRST values were generally slightly higher than FST (Figure 2, B and C). We also observed that the discrepancy between FST and mean pRST was much lower for multilocus than for single-locus estimates.

As expected, RST values are similar to FST values whenever m ⪢μ= 0.001 (island model), 1/t ⪢μ (diverging populations model), or populations are close (stepping-stone model with m ⪢μ). On the contrary, RST becomes considerably larger than FST when m ≤μ (island model), 1/t ≤μ (diverging populations model), or when populations are separated by more than five steps (stepping-stone model Figure 2).

To assess the power of the allele size permutation test, we present in Figure 2 (graphs on the right) the percentage of statistically significant tests (%RHo) among 200 simulation replicates (using α= 5%) according to (1) the migration rate m (island model), (2) the divergence time t in number of generations since isolation (isolated two-population model), and (3) the distance d in number of steps between populations (stepping-stone model). This is done for tests applied to each locus as well as to a multilocus estimate based on 10 loci.

—Simulation results for (1) an island model with migration rate m (A), (2) a two-population model isolated for t generations (B), and (3) a linear stepping-stone model of 30 populations (C). Graphs on the left show RST (□), FST (○), and mean pRST (⋄) values (mean multilocus estimates based on 10 loci and 200 replicates) according to m (A), t (B), or the number of steps separating populations (C). In C, averaged pairwise RST/(1 - RST), FST/(1 - FST), and mean pRST/(1 - pRST) ratios over all pairs separated by given numbers of steps are represented. Graphs on the right illustrate the power of the allele size permutation tests by giving the percentages of significant tests (%RHo) on RST estimates [or average pairwise RST/(1 - RST) ratios] based on a single locus (×) or 10 loci (▵) (i.e., multilocus estimate) and considering a type I error rate criterion α of 5% (dotted line). The symbols (× and ▵) on the horizontal axes of graphs A and B show the values at which the mean square errors of FST and RST are approximately equal.

In the island model, %RHo approaches α for relatively high migration rates (i.e., m = 10 -1 -10 -2 = 10-100μ), in accordance with our a priori expectation that we should not detect a significant effect when m ⪢μ (Figure 2A). On the contrary, for lower migration rates, mutation is no longer negligible compared to migration and the proportion of significant tests increases above α, reaching 88 and 100% when m = 10 -4 (m = 0.1μ) for tests on a single locus or 10 loci, respectively (Figure 2A). Tests based on 10 loci seem actually quite powerful for typical sample sizes encountered in experimental studies (300 individuals here), as 100% of the tests were significant when m =μ and already 24% when m = 10μ. Results of the two isolated population models are very similar to those of the island model if m is replaced by 1/t (Figure 2B). Here, however, tests seem less powerful than in the simulated island model (e.g., for 10 loci, %RHo > 50% when 1/t ≤μ in the isolated population model, and m ≤ 0.3μ in the island model), which is likely due to the smaller sample size (100 vs. 300 individuals) and the lower number of populations sampled (2 vs. 10). B alloux and G oudet (2002) showed indeed that the variance of RST increases substantially with fewer populations sampled. In the stepping-stone model, %RHo increases with the distance separating populations, but reaches a plateau beyond eight steps at ∼60% for estimates based on 10 loci and only 20% for single-locus estimates (Figure 2C). Surprisingly, %RHo is already significantly larger than α for populations separated by just one step and exchanging migrants at a high rate (m/2 = 0.05) relative to the mutation rate (μ= 0.001).

Usefulness of the test to determine the most appropriate statistics: To verify whether the test provides an adequate guideline to choose between RST and FST when assessing population differentiation, mean square errors (MSEs) of FST and RST were computed. The MSE is a synthetic measure of the efficiency of an estimator combining bias and variance (MSE = bias 2 + variance). It has already been used to compare the efficiency of FST and RST estimators (B alloux and G oudet 2002) or gene flow estimates based on FST or RST (G aggiotti et al. 1999). MSEs were computed as Σ(i - e) 2 /n, where i is the FST or RST estimate of the ith replicate, n is the number of replicates (n = 200), and e is the expected value given the demographic parameters. The expected value is e = 1/(1 + 4Nmd/(d - 1)) in the case of the island model (with N = 100 and d = 10), and e = t/(2N + t) in the case of the isolated population model (with N = 500). These are the values expected for RST under SMM and for FST under IAM (or KAM) and a low mutation rate (S latkin 1995 R ousset 1996). Note that e is not the expected FST under the conditions of the simulations (relatively high SMM and μ), but only a good approximation when mutation can be neglected.

For the island model and μ= 0.001 (SMM), with migration rate varying from 0.0001 to 0.1, the ratio MSE(RST)/MSE(FST) varied, respectively, from 0.06 to 2.1 for single-locus estimates and from 0.02 to 2.3 for multilocus estimates based on 10 loci. The migration rate at which MSE(RST) = MSE(FST) was between m = 0.001 and 0.002 for single-locus estimates and between m = 0.003 and 0.005 for multilocus estimates. As can be observed in Figure 2A, these migration rate limits under which RST performs better than FST, and above which the reverse occurs, closely match the migration rate under which the allele size permutation test becomes often significant (i.e., %RHo ≥ 30%). The same pattern is observed for the isolated populations model: For t varying from 30 to 10,000 generations, MSE(RST)/MSE(FST) varied from 2.37 to 0.41 and from 4.00 to 0.01 for single-locus and multilocus estimates, respectively, and MSE(RST) = MSE(FST) for t = 2000 (i.e., 2/μ) and t = 500 (i.e., 0.5/μ) for single-locus and multilocus estimates, respectively. Hence, the test becomes frequently significant when MSE(RST) is close to MSE(FST) (Figure 2B).

These results strongly suggest that the allele size permutation test is well suited to determine which of FST or RST is the most adequate for demographic parameters inferences, at least on the basis of the lowest MSE criterion. However, it must be pointed out that the statistic with lowest MSE is not necessarily the statistic that will provide the lowest MSE in the demographic estimate, because demographic estimates are usually not linear functions of FST or RST. For example, in the isolated population model, the τ= t/N estimates that can be derived using τF = 2FST/(1 - FST) and τR = 2RST/(1 - RST) give MSE(τR) > MSE(τF) for all simulated divergence time with single-locus estimates [τF can also be estimated as -ln(1 - FST) (R eynolds et al. 1983), but this leads essentially to the same results]. This occurs because whenever FST or RST approaches 1, the inferred τ quickly takes enormous values, so that the impact of the larger variance of RST relative to FST is greatly amplified in the inferred τ, although τR is much less biased than τF for τ≥ 1. The good news is that for multilocus estimates we obtained MSE(τR) = MSE(τF) for t = 500 and MSE(τR) < MSE(τF) for t > 500, as previously found for MSE(RST) = MSE(FST). Similarly, for the island model, where Nm can be estimated as NmF = (1/FST - 1)/4 and NmR = (1/RST - 1)/4, the m values corresponding to MSE(NmF) = MSE(NmR) were exactly equal to these obtained for MSE(RST) = MSE(FST) for both single- and multilocus estimates. Thus, the usefulness of the allele size permutation test to determine which of FST or RST is the most adequate for inferential purposes seems to be quite general, except probably with low sample size and/or low number of loci, when inferences are in any case doubtful because associated variances are too large.

Application examples: To illustrate the utility and power of the allele size permutation test with real data we present three examples of published data sets that we reanalyzed. These data were collected to assess population differentiation and check for isolation by distance in three different organisms. We computed global or pairwise FST and RST statistics as described above and applied the allele size permutation tests to obtain pRST values. These analyses were performed with SPAGeDi.

Biomphalaria pfeifferi, a selfing snail recently introduced in Madagascar: Biomphalaria pfeifferi, an intermediate host of a parasitic trematode causing intestinal bilharziasis, is a hermaphroditic freshwater snail distributed over most of Africa, the Middle East, and Madagascar. Madagascar was relatively recently invaded by this snail, probably as a result of human occupation a few hundred years ago (C harbonnel et al. 2002a). Moreover, according to a broad-scale survey of microsatellite variation throughout Madagascar, bottleneck (C ornuet and L uikart 1996) and admixture (B ertolle and E xcoffier 1998) tests suggest that at least three independent introductions from genetically differentiated sources occurred (C harbonnel et al. 2002a). A small-scale study of microsatellite variation also reveals that populations experienced recurrent bottlenecks and that migration has been frequent within watersheds but rare among them (C harbonnel et al. 2002b). This population dynamic and the high selfing rate experienced by this snail explain the high genetic differentiation among populations observed in Madagascar: FST = 0.80 and 0.58 for broad and small scales, respectively (C harbonnel et al. 2002a,b).

In this particular context, we can formulate a hypothesis regarding the information content that microsatellite allele sizes could bear. Given the postulated recent introductions of this snail in Madagascar, we expect that mutation has not contributed to differentiation among populations originating from the same introduction but has contributed to differentiation among populations originating from different introductions (at least if the source populations had diverged over enough time). The places and timing of the introductions are not known, but populations from a single watershed are likely to originate from a single introduction or, if genotypes from different introductions mixed in a watershed, migration within the watershed is likely to have prevented the buildup of a phylogeographical pattern at this scale. Therefore, we can expect RST to be close to FST for populations belonging to the same watershed and significantly larger than FST for populations from different watersheds when the latter were originally colonized by individuals from independent introductions.

Differentiation among populations of Biomphalaria pfeifferi at different scales

To test this hypothesis, we reanalyzed data from small-scale and large-scale studies by C harbonnel et al. (2002a,b). Global RST and FST values as well as pairwise RST and FST values between populations were computed. Distinguishing pairs of populations within or among watersheds, pairwise values were regressed on spatial distances (Mantel tests were used to assess the significance of the regression slopes), and average pairwise values were computed for a set of distance classes (defined in such a way that each contained ∼33 pairs of populations). One thousand random permutations of the allele sizes provided a distribution of pRST values, 95% confidence intervals covering the 25th to the 975th ordered values, and P values testing if RST > pRST.

Multilocus RST values are significantly higher than mean pRST at a broad scale but not at a local scale (Table 3). Applied to each locus, these tests were also significant for four out of eight loci at the broad scale but for none at the local scale.

The analysis of average pairwise multilocus FST and RST values per distance class at the broad scale shows the following (Figure 3):

Differentiation between populations occupying the same watershed is much lower than that between populations from different watersheds, even for populations separated by the same spatial distance. This is in line with the higher migration rate detected within watersheds than among them (C harbonnel 2002b).

A pattern of isolation by distance is detected within watersheds for both FST and RST (Mantel tests: P = 0.007 and 0.021, respectively). Among watersheds, such a pattern is not detected for FST but is for RST (Mantel tests: P = 0.18 and 0.002, respectively).

Within watersheds, RST’s are not significantly higher than pRST’s, whereas among watersheds, RST’s are significantly higher than pRST’s for all distance classes but the first one.

Average pairwise pRST values are always somewhat lower than pairwise FST values but they follow closely their pattern of variation with spatial distance.

In conclusion, at a local scale, RST values are close to FST values, and allele size permutation tests do not reveal any significant contribution of stepwise mutations to population differentiation. On the contrary, at a large scale, RST values are substantially higher than FST values and allele size permutation tests demonstrate that shifts in average allele sizes contribute significantly to population differentiation. Significant tests on RST values are expected if populations had diverged for a sufficiently long time and/or if populations exchanged migrants at a rate similar or inferior to the mutation rate. The results are thus very consistent with a priori expectations given that (1) at a large scale, both these conditions are probably met because populations far apart in Madagascar probably originated from relatively recent and independent introductions from source continental populations isolated for a long time, and migration rate is low among watersheds, and (2) at a local scale, particularly within watersheds, none of these conditions are likely to be met.

—Average pairwise FST (○ and •), RST (□ and ▪), and mean pRST (⋄ and ♦) values among populations of Biomphalaria pfeifferi throughout Madagascar for a set of distance classes, distinguishing comparisons between populations within watersheds (•, ▪, ♦) and among watersheds (○, □, ⋄). The dotted lines represent the range of the 95% central ordered pRST values (i.e., after allele size randomization). Each distance class contains 32-35 pairs of populations.

Fraxinus excelsior, a widespread European tree: Fraxinus excelsior (Oleaceae, common ash) is a widespread European wind-pollinated tree species found mostly in floodplain locations and with a scattered distribution within natural forests. The distribution of chloroplastic DNA (cpDNA) haplotypes throughout Europe suggests that F. excelsior was located in at least three different refuges during the last ice age, one putative refuge being the Balkan area (G. G. V endramin , unpublished data). H euertz et al. (2001) analyzed microsatellite polymorphism in 10 Bulgarian populations (Balkan area) from three regions (321 individuals). Populations were separated by 0.5-22 km within regions and 120-300 km among regions.

In the absence of evidence of long-term divergence between Bulgarian populations (no evidence of different refuges), and given that gene flow should be relatively extended in a wind-pollinated species, we may expect that stepwise-like mutations have not contributed significantly to population differentiation in Bulgaria. The data set of H euertz et al. (2001) was thus reanalyzed to compare average pairwise FST and RST values between populations, distinguishing pairs within and among Bulgarian regions, and testing RST values by allele size permutations (1000 randomizations).

Mean pairwise multilocus estimates were equal to FST = 0.074, RST = 0.091 within regions and FST = 0.097, RST = 0.180 among regions (Figure 4). Hence, whereas differentiation increases slightly from small to large geographical scales according to FST, it nearly doubles according to RST. Moreover, average pairwise RST is much larger than FST among regions, but only slightly larger than FST within regions. Within regions, observed RST’s are always within the 95% range of central pRST, but among regions, the multilocus RST estimate as well as the estimate for locus FEM19 is larger than the 95% range of pRST (Figure 4), demonstrating that stepwise-like mutations contributed to population differentiation at the large geographical scale for at least one locus.

Several causes may account for the significant allele size effect on population differentiation among regions in Bulgaria, for example:

The pattern may reflect isolation by distance. However, it seems unlikely that migration rate among regions is weak compared to the mutation rate given that pollen is wind dispersed.

The pattern may be due to postglacial recolonization from different refuges. There is, however, no evi-dence of different refuges from the maternally inherited cytoplasmic DNA as the same unique haplotype occurs in all three regions (M. H euertz , unpublished data).

The pattern may reflect human-mediated introduction of Fraxinus from remote regions.

The pattern may reflect locally occurring hybridization between F. excelsior and a related species such as F. angustifolia or F. pallisiae. Given that a total of four ash species (the former three and F. ornus) are found in Bulgaria and that different species occur in the same forests (M. H euertz , personal observation), this latter hypothesis merits further investigation. In any case, the observation that a significant effect of stepwise-like mutations is observed on a large scale but not on a small one remains very consistent with a priori expectations, as nearby populations should exchange genes at a relatively high rate.

—Mean pairwise RST, mean pRST, and FST values between Bulgarian populations of Fraxinus excelsior for populations belonging to the same region (A) or different regions (B). Values are given for each locus and the multilocus estimates. Bars of pRST indicate the mean pRST values over 1000 allele size permutations, and the corresponding intervals give the range of the 95% central pRST values.

Centaurea corymbosa, a rare and narrow-ranged cliff-dwelling herb: Centaurea corymbosa (Asteraceae) is a short-lived perennial herb species distributed over a very narrow range (within a 3-km 2 area of a calcareous massif along the French Mediterranean coast), where it occurs in only six small populations (C olas et al. 1997). It has specialized into an extreme habitat: the top of limestone cliffs where few other plant species survive. On more fertile ground, C. corymbosa is outcompeted, so that suitable habitat is highly fragmented, appearing as small islands dispersed in the landscape. Given that the species occupies only a small fraction of these “islands” (the whole massif extends over 50 km 2 ), colonization ability must be very limited, probably as a consequence of limited seed dispersal ability and the self-incompatibility system that prevents a potential newcomer from founding a new population on its own (C olas et al. 1997 F réville et al. 2001). Patterns of isozyme (C olas et al. 1997) and microsatellite (F réville et al. 2001) variation show high levels of differentiation among populations, with FST = 0.35 and 0.23, respectively, despite the narrow range of the species (2.3 km between the two most distant populations). High differentiation at such a small scale cannot be attributed to the mating system as the species is self-incompatible. It most likely results from small population sizes and low gene flow among populations. It might also be a consequence of more or less recurrent bottlenecks when new populations are founded (although the turnover should be relatively slow, given that no population extinction or foundation has been observed since 1994, when C. corymbosa populations began to be closely surveyed, and herbarium data show that five of the six populations were known >100 years ago).

In this context it is interesting to question whether gene flow among populations is sufficiently low to permit divergence by mutations. The higher observed FST value at allozyme loci than at microsatellite loci could indeed be caused by high mutation rates of microsatellites, provided that μ≥ m. F réville et al. (2001) pointed out that this hypothesis was also supported by the fact that FST values at the two most polymorphic microsatellite loci (12B1 and 21D9, Table 4), the ones likely to have the highest mutation rates, were lower than those for the two loci with intermediate levels of polymorphism (13D10 and 28A7, Table 4).

The allele size randomization procedure is adequate to address this question. Therefore, global RST, pRST, and FST were computed for microsatellite loci as described above, and RST was compared against the distribution of 1000 pRST values. Permutation tests did not detect any RST value significantly >pRST (Table 4). This suggests thus that differentiation is caused mainly by drift and that gene flow, m, and/or the reciprocal of divergence time, 1/t, are large compared to the mutation rate, μ. This result also implies that FST should be a better estimator than RST of population differentiation for this species. Actually, given the small population sizes (C olas et al. 1997, 2001), drift is expected to be high. For example, if populations had effective sizes of ∼100 individuals (there is actually much variance among populations) and conformed to an island model (there are actually some isolation-by-distance effects), a value of m = 0.006 would account for the observed FST, a value larger than typical microsatellite mutation rates (10 -3 -10 -4 ). Assuming that these populations have been in place for a sufficiently long time to potentially permit differentiation by mutations (shifting allele sizes), the absence of such mutation-driven differentiation also suggests that the migration rate is larger than the mutation rate, so that new mutation variants spread over all populations.

Differentiation among populations of Centaurea corymbosa, estimated by global RST, mean pRST, and FST values per locus and for a multilocus average

Nonsignificant tests could also be due to a lack of power, so the test should be applied to additional microsatellite loci to confirm these results (presently, only four out of six loci had a sufficient number of alleles to carry out permutation tests). Deviation from a SMM at some loci could also reduce the power of the test. For example, the dinucleotide locus 28A7 has six alleles with sizes following a sequence of one repeat step plus one allele at least six repeats smaller than the other ones. Although this pattern is not necessarily incompatible with a pure SMM (e.g., D onnelly 1999), it might suggest that a mutation of large effect created the outsider allele.


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 ith 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 ith 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.

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):

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.

FRP13:

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 PP2d. 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.


3. Materials and methods

(i) Seed collections and germination protocol

In summer 2004, we collected fruits from plants of nine populations along Tunisian coasts (Table 1). The size of these populations varies from almost 60 individuals to much more than 500 individual plants. Ten seeds per family were sown in a Petri dish on moistened filter paper. Germination was carried out under laboratory conditions (in the dark at 20–25°C). Four days later, four randomly chosen seedlings per family were planted in separate pots and randomized in an unheated greenhouse. Ten families per population (i.e. 40 seedlings per population) were cultivated under uniform environmental conditions and used for quantitative genetic study. Out of each group of 40 seedlings, 30 were used for allozyme analysis.

Table 1. Cakile maritima Tunisian populations. Their climate region and their alphanumeric code

(ii) Seed buoyancy and viability

The trait buoyancy (floating capacity) indicates the potential of species to be dispersed by water. It is given as the proportion of seeds floating after a defined time period. For this purpose, 100 seeds of C. maritima were placed in a closed bottle containing 500 ml of sea water in October 2005. The number of seeds that sank was recorded over time. In February 2006, ten of the floating seeds were sown on filtered sand and irrigated with tap water to test their viability.

(iii) Morphological measurements

Two sets of morphological traits were considered discriminator traits between populations: vegetative traits (leaf morphotype (leaf form (LFF)) and leaf length (LFL)) and reproductive traits (flower date emergence (FED), petal length (PTLT), petal width (PTW), petal colour (PTC), pistil length (PSL), fruit type (FRT), fruit sizes (length of the upper segment (LUP) and length of the lower segment (LLW) and seed number (SN)). Leaf characteristics were graded from 1 to 10, with 1 being leaves with entire margins, 10 being leaves with the most deeply pinnatified margins, and intermediates classified by graduations between these extremes (Fig. 1). Leaves selected for morphotype determination were also used for length measurements. Fruit size was measured separately for the upper and lower segments and was classified into three discriminate modals: unhorned, intermediate and horned quantified as a series of 0, 1 and 2, respectively (Fig. 2). Leaf length, petal length and width, pistil length and fruit size were measured by an electronic calliper. Colour and size measurements were performed on one petal chosen haphazardly from each individual.

Fig. 1. Leaf morphotypes in Cakile arranged as a standard sequence, 1–10, for field assessment of leaf variation.

Fig. 2. Siliculas modal of Cakile (0: unhorned fruit 1: intermediate fruit 2: horned fruit).

(iv) Data analysis

A nested analysis of variance (ANOVA) including population and family (nested within population) as random effects was conducted for each quantitative trait. The level of population differentiation in the quantitative traits was measured with Qst (Spitze, Reference Spitze 1993), which is analogous to Fst measured using allozyme marker loci. In order to estimate Qst, δb 2 is obtained directly from the population variance δp 2 , that is, (δb 2 =δp 2 ), whereas the family variance δf 2 has to be converted into δw 2 by multiplication with a coefficient (c) that depends on the relationship of individuals within families (δw 2 =cδf 2 ). For half-sibs, full-sibs and cloned individuals, c is 4, 2 and 1 (under the assumption of no dominance and epistasis), respectively. When populations are in Hardy–Weinberg disequilibrium (Fis≠0), the level of differentiation in quantitative traits could be expressed as: (Yang et al., Reference Yang, Yeh and Yanchukt 1996) where Qst has the expression (Spitze, Reference Spitze 1993) and others (Lande, Reference Lande 1992 Latta, Reference Latta 2004).

V W is estimated to be four times the between-family component of variance V fam under the assumption that maternal effects were weak and that within open-pollinated families were related as half-sibs (Yang et al., Reference Yang, Yeh and Yanchukt 1996), while V B is simply the among-population variance component. Approximate standard error (SE) values of Qst were obtained by the ‘delta’ technique (Stuart & Ord, Reference Stuart and Ord 1987). This method, used by Podolsky & Holtsford ( Reference Podolsky and Holtsford 1995), provides reliable SE estimates of genetic variance components for setting confidence limits (CLs) to genetic parameters (Hohls, Reference Hohls 1996).

(v) Allozyme experiment

For the allozyme analysis, an electrophoresis survey was used to estimate genetic variability within and among C. maritima populations. Approximately 200 mg leaf tissue were collected from each plant (1 month old), ground under liquid nitrogen and mixed with 100 μl of extraction buffer (PVP–potassium phosphate grinding buffer, pH 7), as described by Thrall et al. ( Reference Thrall, Andrew and Burdon 2000), and were centrifuged at 19 000 g for 20 min. Ground material was absorbed on Whatman 3MM filter paper wicks and stored in an ultra-cold freezer (−70°C) until analysis. Horizontal starch-gel electrophoresis was performed for seven enzyme systems revealing a minimum of 13 loci: peroxidase (Px EC 1.11.1.7), isocitrate dehydrogenase (IDH EC 1.1.1.42), glutamate oxaloacetate transaminase (GOT EC 2.6.1.1), shikimate dehydrogenase (SDH EC 1.1.1.25), leucine aminopeptidase (LAP EC 3.4.11.1), 6-phosphogluconate dehydrogenase (6-PGD EC 1.1.1.44) and malate dehydrogenase (MDH EC 1.1.1.37). The compositions of gel and electrode buffers were described in Soltis et al. ( Reference Soltis, Haufler, Darrow and Gastony 1983) and the methods used to stain allozyme bands were described in Michaud et al. ( Reference Michaud, Lumaret and Romane 1992) for Px and in Cardy et al. ( Reference Cardy, Stuber and Goodman 1980) for IDH, GOT, LAP, SDH, 6-PGD and MDH. For phosphatase acid (ACPH EC 3.1.3.2), vertical zoned polyacrylamide gels were prepared following Laemmli ( Reference Laemmli 1970) and were stained according to Selander et al. ( Reference Selander, Smith, Yang, Johnson and Gentry 1971). Loci were numbered sequentially with the most anodally migrating locus designated as locus 1. Genetic interpretation for all loci was straightforward.

(vi) Data analysis

At each of the 13 loci studied in C. maritima, genotypic and allelic frequencies were assessed from a survey of gel phenograms. Three coefficients, measuring genetic variation, were computed using BIOSYS-1: (i) the percentage of polymorphic loci (P) using 0·95 criterion (a locus is considered polymorphic if the most common allele has a frequency of less than 95% in all the populations analysed) (ii) the mean number of alleles per locus (A) and (iii) the expected heterozygosity (He) under Hardy–Weinberg equilibrium (Nei, Reference Nei 1978). These three coefficients were estimated for each of the sampled populations. The mean and standard deviation of the previously mentioned coefficients were then calculated over all populations sampled. The genetic structure within and among populations was also evaluated using Wright's ( Reference Wright 1965) F-statistics Fit, Fis and Fst. Heterozygote frequency for each polymorphic locus in each population was tested for significant deviation from the Hardy–Weinberg expectations with a χ 2 analysis (Li & Horvitz, Reference Li and Horvitz 1953) (BIOSYS-1). Estimate of Nm (the number of migrants per generation) was based on Fst (Wright, Reference Wright 1951) (Genetix 4.02). To test if populations' genetic differences followed the model of isolation by distance, we established the relationship between the genetic difference Fst values between pairs of populations and their geographical distance using the Mantel test (Genetix 4.02). In order to examine the genetic distance and the geographic distance, the UPGMA clustering analyses were conducted (BIOSYS-1).

(vii) Matrix correlation

The association between morphological, genetic, sea current speed and geographical distance matrices was examined with the Mantel test (1967). The Mantel test uses random permutations of matrix rows and columns to test if correlations between distance matrices are greater than expected by chance (Sokal & Rohlf, Reference Sokal and Rohlf 1995). We compared genetic and morphological distance matrices to each other with geographical matrices individually and then we compared Nm estimated from Fst and Qst to sea current speed. Pairwise population of sea current speed was measured by using ship-drift. Ship-drift measurement of surface current velocity consists of the vector difference between the velocity of a ship determined from two position fixes and the average estimated velocity of the ship through the water during the same time interval, usually 12–24 h. The vector difference is considered to be due to a surface current. The absence of connectedness was marked as 0.

(viii) Assignment tests

To test whether the marine currents, as described in Fig. 5, drive the direction of gene flow, it will be interesting to use the Assignment tests. These methods are strong tools to detect recent immigration events even when the overall population differentiation is low (Rannala & Mountain, Reference Rannala and Mountain 1997 Waser & Strobeck, Reference Waser and Strobeck 1998 Castric & Bernatchez, Reference Castric and Bernatchez 2004). In fact, assignment tests are particularly useful to trace the contemporary dynamics of natural populations without requiring equilibrium assumptions based on long-term genetic processes (Manel et al., Reference Manel, Gaggiotti and Waples 2005). As the pattern of misclassification of individuals within a population can be used to determine the direction of migration, assignment tests might be a useful method to examine the influence of sea currents on gene flow. The probability that an individual sampled within a given population is a migrant originating from another population was determined by using a Bayesian method (Geneclass2 software) (Piry et al., Reference Piry, Alapetite, Cornuet, Paetkau, Baudouin and Estoup 2004).


Discussion

Clustering approaches allow the partition of a sample of individuals into genetically distinct groups without an a priori definition of these groups. Most of the recent advances in clustering methodology have been made using Bayesian statistical models [3, 20, 5, 21, 22]. Bayesian methods assign individuals to groups based on their genotypes and the assumption that the markers are in Hardy-Weinberg and linkage equilibrium within each subpopulation.

In this study, a new method was used to infer the hidden structure in a population, based on the maximisation of the genetic distance and not making any assumption on HWE and LE, and we show that it yields a good performance under different simulated scenarios and with a real data set. Therefore, it could be a useful tool to determine genetically homogeneous groups, especially in those situations where the number of clusters is high, with complex population structure and where HWD and/or LD are present.

The simulation results indicate that the BAPS method is the least precise since it needed a large number of genotyped markers to reach the correct partition, especially when the population had reached the mutation-migration-drift equilibrium. For the original/basic scenarios, the performances of MGD and STRUCTURE were similar (good) whatever the parameter of comparison, although the new method presented a slight advantage (see Table 3 and Figure 2).

We have shown that departures from the implicit assumptions in the Bayesian methods about the Hardy-Weinberg and linkage equilibrium within populations affect their accuracy, especially for BAPS, leading to an overestimated number of clusters and a reduced proportion of correct groupings. These observations are in agreement with Kaeuffer et al. [35] who have shown that a high LD correlation coefficient value increases the probability of detecting spurious clustering with STRUCTURE. The randomisation of alleles (and also the randomisation of genotypes and haplotypes to some extent) re-establishes both HWE and LE. In these situations, the two methods evaluate correctly the number of clusters and give an increased proportion of correct groupings. On the contrary, MGD is more precise in disequilibrium situations and its performance does not change significantly after the randomisation, demonstrating the independence of the novel method from the existence or not of HWE and LE. From the results presented here, an alternative to test the accuracy of the results from any clustering method would be to compare the results obtained after the randomisation of the molecular information within each pre-defined subpopulation when this information is available.

The precision of all three methods is excellent for F STas low as 0.03. This is in agreement with the results of Latch et al. [10], who have proven that STRUCTURE and BAPS discern the population substructure extremely well at F ST= 0.02 - 0.03. However, in our simulations only STRUCTURE determines the correct number of clusters at F ST= 0.01. Notwithstanding, there is a controversy about the minimum differentiation level necessary for a population to be considered as genetically structured. Waples and Gaggiotti [36] have suggested that if F STis too reduced (e.g. F ST= 0.01) then it probably cannot be associated with statistically significant evidence for departures from panmixia. In these situations, it is not clear if the most appropriate solution for MGD (and also the other clustering methodologies) is to separate different subpopulations or to maintain the subpopulations as an undifferentiated population.

The simulated scenarios taking into account different selfing rates indicated both an increase in differentiation between subpopulations (i.e. higher F STvalues) and an increase in Hardy-Weinberg disequilibrium (F ISmoves from 0.01 to 0.81). However, the increase in F STvalues (from 0.27 to 0.42) are are not as great as that of the F ISvalues indicating that the Hardy-Weinberg disequilibrium can not be masked by the effect of the differentiation level. In addition, the increase in F STvalues should help to distinguish the different clusters and, therefore, the HWD should reach at least the lowest limit of its effect.

Our results obtained with the MGD method from the human data set are, in general, similar to those obtained with STRUCTURE [34] and also in concordance with a more recent study of 525910 SNP [37], although some discrepancies exist with the results of Li et al. [38] using 650000 SNP. Rosemberg et al. [34] have indicated multiple clustering solutions for K = 7 with STRUCTURE. However, the results obtained with MGD for K = 7 are in complete agreement with the seven geographical regions. A careful inspection of the results detects clusters where grouped individuals have multiple sources of ancestry, especially those in the Middle East and Central-South Asia. This situation (i.e. the estimated mixed ancestry) could be due either to recent admixture or to shared ancestry before the divergence of two populations but without subsequent gene flow between them. It has been indicated that global human genetic variation is greatly influenced by geography [39–41]. In addition, Serre and Pääbo [42] have indicated that the clusters obtained by Rosenberg et al. [34] have been generated by heterogeneous sampling and that these would disappear if more populations were analysed.

In this study, a simple island model with constant population sizes and invariant symmetrical migration has been considered, which are unlikely in natural systems. The performance of STRUCTURE has been recently evaluated [23] by simulating various dispersal scenarios and it seems to perform well with more complex population structures than the finite island model (hierarchical island model, contact zone model). In this study, the performance of the MGD method was better than that of the Bayesian approaches in the simulated scenarios with a higher number of clusters and a more complex population structure. However, further investigations are required to determine the capacity of the MGD method to deal with other kinds of population structure.

Computation time may be a limitation of the new method, especially when dealing with large amounts of markers. However, it should be noted that clustering analysis is not performed very often and the results are not usually needed urgently. Therefore, it may be worthwhile to wait for the results obtained with the most accurate method.

If the genetic distance calculated from the molecular coancestry has been evaluated as an alternative, then the use of other genetic distances previously published in the literature [24] could be investigated as the parameter to maximise both for codominant and dominant molecular markers. Moreover, the Nei minimum distance [25] could be inappropriate when working with various markers, for example when mixing data obtained with markers with different heterozygosis levels (e.g. mixing microsatellite and SNP data). In addition, a weighting procedure [43, 44] could also be implemented taking into account the subpopulation size, the number of loci or the number of alleles. Notwithstanding, the nature of the new method (i.e. the maximisation of the genetic distance) allows for the use of any measure which could better fit the available molecular data, beyond the Nei distance.

The informativity of the markers has a clear effect on the efficiency of the clustering methods, especially for BAPS. Increasing the number of markers (scenario 1 vs. 2, 3 vs. 4, 5 vs. 6 and 7 vs. 8) almost always yields better results: the correct number of clusters is estimated in more cases and the percentage of correct groupings is higher. In parallel, when comparing a similar number of markers but with different degrees of polymorphism (scenario 2 vs. 5, microsatellites vs. SNP) the biallelic markers yield worse performances. Notwithstanding, when using a reasonable number of markers (50 microsatelites and 300 SNP) MGD and STRUCTURE, at least, provide a high accuracy. However, when comparing results obtained with STRUCTURE, it is surprising that this method showed less accuracy with 10 microsatellites than with 50 microsatellites.

Although in the present work the method has been developed for co-dominant markers, whatever the approach (molecular coancestry or allelic frequencies), the methodology can also be easily extended to dominant molecular markers by replacing the molecular coancestry matrix with a matrix of any available measure of similarity for dominant markers [45] or estimating the allelic frequencies from recessives (see [46] and references therein) and then using the typical genetic distances.

The present formulation of the method does not explicitly account for the presence of admixed individuals. To do so, a different set of probabilities should be given to each locus in each individual (in the allelic frequencies approach) allowing for each locus to be assigned to different clusters. The increase in computation time and the ability of the optimisation algorithm to deal with a larger space of solutions deserve further investigations.

A compiled file of the code used to infer the number of clusters and the assignment of the individuals to each cluster in a given sample from the molecular coancestry matrix or the allele frequencies will be available on the web site http://www.uvigo.es/webs/c03/webc03/XENETICA/XB2/Jesus/Fernandez.htm.


RESULTS

Detecting hybridization with population genetic vs. community ecology diversity metrics

Of the four population genetic metrics and three species diversity indices, all calculated using a multilocus approach, only one (HN) detected a change in genetic diversity in hybrid populations ( Fig. 1 ). Relative to populations of wild sea beet, hybrid populations exhibited statistically significantly larger HN (Mann–Whitney U = 16, z = 𢄢.09, P = 0.037). However, we did not detect differences in other multilocus metrics of genetic diversity, including percent polymorphic loci (U = 40, z = 0.04, P = 0.10), number of alleles per polymorphic locus (U = 40, z = 0.04, P = 0.10), total number of alleles (U = 30, z = 𢄠.84, P = 0.40), Shannon–Weiner’s H (U = 23, z = 𢄡.47, P = 0.14), Simpson’s D (U = 26, z = 𢄡.20, P = 0.23), or McIntosh’s E (U = 47, z = 𢄡.41, P = 0.16).

In contrast to the multilocus results, single locus comparisons using species diversity metrics were more informative. Relative to wild sea beet populations, hybrid populations exhibited larger Shannon–Weiner diversity (H, at 11 of 12 loci, sign test: P = 0.003, e.g., Fig. 2 ), larger Simpson’s D (at 10 of 12 loci, P = 0.02), and larger McIntosh’s E values (at 11 of 12 loci, P = 0.003). Moreover, single locus comparisons of diversity using traditional population genetic metrics were less sensitive than community ecology metrics to changes in allele diversity or composition. Relative to wild sea beet populations, hybrid populations exhibited more alleles (at 10 of 12 loci, P = 0.02), but hybrid populations did not differ significantly from wild populations in percent polymorphic loci (decreased at nine of 12 loci, P = 0.073) or number of alleles per polymorphic locus (increased at eight of 12 loci, P = 0.19).

Mean single-locus Shannon–Weiner diversity estimates for 12 loci (separate lines) averaged across 10 wild or eight putative hybrid Beta vulgaris subsp. maritima populations. The 12 loci are represented with abbreviations here: AAT (aspartate amino transferase, E.C. 2.6.1.1), ACO (aconitase E.C. 4.2.1.3), GDH (glutamate dehydrogenase E.C. 1.4.1.2), LAP (leucine aminopeptidase E.C. 3.4.11.1), MDH1, MDH2 (NAD + malate dehydrogenase E.C. 1.1.1.37), PGM1, PGM2 (phosphoglucomutase E.C. 5.4.2.2), SKD (shikimate dehydrogenase E.C. 1.1.1.25), TPI1, TPI2 (triose phosphate isomerase E.C. 5.3.1.1), and UDP (uridine diphosphoglucose pyrophosphorylase E.C. 2.4.1.1). Error estimates not shown for clarity.

The consequences of hybridization for genetic diversity

Across the 12 loci, putative hybrid beet populations possessed only one-quarter of the rare alleles found in wild populations (χ 2 = 22.5, df = 1, P < 0.001 Fig. 3 ). Based on the Jaccard similarity coefficient, hybrid populations were 10.8% (±SE = 0.6%) more similar to chard cultivars than wild populations (one-sample t test: t = 18.66, df = 47, P < 0.001). Furthermore, hybrid populations were 3.7% (±SE = 1.1%) more similar to sugar beet cultivars than wild populations (t = 3.45, df = 47, P = 0.001). Finally, based on paired comparisons, hybrid populations were significantly more similar to chard than sugar beet cultivars (paired t test: t = 𢄦.62, df = 47, P < 0.001).

Mean number of alleles that do not contribute significantly to gene diversity (i.e., rare alleles, ±SE) per locus estimated for 12 allozyme loci averaged across 10 wild or eight putative hybrid Beta vulgaris subsp. maritima populations.


Results and discussion

To illustrate the applicability of our approach we apply it to two previously published data sets that were analyzed in [37] and [17], respectively.

Beetle Data

The first data set was used as part of a phylogeographic study of the beetle species Brachyderes rugatus rugatus on La Palma (Canary Islands) [37]. In this study 138 individual beetles were sampled. The 18 sampling locations are shown in Figure ​ Figure3. 3 . Using sequence data from the mitochondrial COII gene (for details see [37]), the 138 samples were subsequently grouped into 69 haplotypes, and a haplotype phylogeny based on the parsimony criterion was constructed using the TCS program [38]. This phylogeny is presented in Figure ​ Figure4 4 .

Sampling locations and regions for beetle data. A map of La Palma with sampling locations indicated by black dots [37]. Sampling locations where haplotypes from a particular phylogroup (cf. Figure 4) were found are depicted by the dashed curves. Note that the sampling location Altos de Jedey is the only one where haplotypes from two distinct phylogroups (namely 1 and 2) were found. The six groups of sampling locations corresponding to the six regions R1, R2, . R6 discussed in the text are also indicated.

Haplotype phylogeny for beetle data. The haplotype network presented in [37] for the haplotypes collected in La Palma. Note that all edges have length 1. The colored dots (black, red and green) represent the sampled haplotypes and the white dots hypothetical intermediates. Dashed boxes correspond to the three phylogroups, 1-3, identified in [37]. The haplotypes found in region R2 are highlighted in red, those found in R6 in green and those found in R3 are indicated by blue circles.

According to this phylogeny, the haplotypes were divided into 3 phylogroups, as indicated on the phylogeny and in Figure ​ Figure3. 3 . Based on these groupings it was concluded for Brachyderes rugatus rugatus that (i) there is a region of secondary contact, or melting-pot, in the South of the island at the overlap of regions 1 and 2, and (ii) that there is an ancestral region or hot-spot in the region containing the three sampling locations in the top right of region 2. Note that in [37] support for conclusion (i) was provided by performing the test given in [8] for detecting zones of secondary contact, which essentially involves calculation of the average distance between the geographic centers of clades at increasing nesting levels in a phylogeny on the haplotypes of interest.

To investigate whether our new method was supportive of conclusions (i) and (ii) or not, we first grouped the sampling locations into 6 regions R1, . R6 as shown in Figure ​ Figure3. 3 . We used these regions rather than the individual sampling locations, since the number of samples taken at each location was very small (between 2 and 8). When forming the groups, geographically close locations were grouped together. We also considered other groupings based on geographic proximity (data not shown) and the outcome was similar, though less pronounced when the number of groupings was reduced (smallest number of groupings used was 3). We then measured the diversity (using the measure PD) and haplotype connectivity for the haplotypes found in each region Rirelative to the phyletic distances given by the phylogeny in Figure ​ Figure4, 4 , as described in the Methods section.

The results for the 6 regions are summarized in Table ​ Table1. 1 . In this table, we present the size of the subset Y of haplotypes found in the region (column 2), the values PD(Y), PDmin(|Y|), PDmax(|Y|) (columns 3-5), and the normalized diversity score PD*(Y) (column 6) as defined in the Methods section. Similarly, we present the values HC(Y), HCmin(|Y|), HCmax(|Y|) and HC*(Y) (columns 7-10).

Table 1

RegionNumber of Haplotypes in regionDiversityHaplotype connectivity
PDPDminPDmaxPD*HCHCminHCmaxHC*
R6214725870.35143250.50
R3112810670.32161270.58
R2183320810.2173250.18
R47146550.1651270.15
R5182920810.1553250.09
R15104480.1471280.22

Diversity and haplotype connectivity scores for the geographic regions on La Palma indicated in Figure 3, ranked according to normalized phylogenetic diversity scores, PD*, as defined in the main text. The columns labeled with PDmin, PDmax, HCmin and HCmax contain the minimum/maximum score over all subsets containing the same number of haplotypes as found in the region.

As can be seen in Table ​ Table1, 1 , the two regions with the highest PD*score are R6 and R3, which also have a much higher HC* score than any of the other four regions. This is supportive of conclusion (i), i.e. that R6 is probably a melting-pot. Indeed, in Figure ​ Figure4 4 the haplotypes found in region R6 are highlighted in green, and it can be seen that they clump together into two groups. This also indicates why we obtained a high HC* score for this region. Similarly, the high PD* and HC* scores for region R3 suggests that this region is a melting-pot as well, a conclusion that is consistent with the findings in [37] where it is suggested that in R3 range expansions toward the South and the Northwest partially overlapped.

Concerning conclusion (ii), we see that amongst the remaining regions R2 clearly has the highest PD* score and a much lower HC* score than R6 and R3. This pattern of scores, i.e. relatively high diversity and low haplotype connectivity, is more supportive of a hot-spot scenario rather than a melting-pot scenario, in agreement with conclusion (ii). Examining Figure ​ Figure4, 4 , we see that the haplotypes in R2 (highlighted in red) are relatively spread out over the haplotype phylogeny, hence the low haplotype connectivity score.

Pine Data

The second data set that we consider formed part of a study of the phylogeographic history of the species Pinus pinaster around the Mediterranean [17]. Samples were taken from 10 locations as indicated in Figure ​ Figure5. 5 . Sequence data consisting of nine chloroplast simple sequence repeat markers gave rise to 34 different haplotypes (for details see [17]). For these 34 haplotypes a distance matrix was computed using the pairwise haplotypic difference (that is, for any two haplotypes, the sum of the difference between the allele size over the nine loci).

Sampling locations for pine data. Sampling locations for the data set in [17].

To understand the phylogeographic structure of this data, in [17] the frequency distribution of the pairwise distances between haplotypes, sometimes also called the genetic diversity spectrum (GDS) [12], was computed. We have recomputed this and depict the result in Figure ​ Figure6. 6 . In particular, based on considerations - such as the shape of the GDS for the Landes and Pantelleria locations - it was hypothesized that Landes and Pantelleria are hot-spots, although it was also stated that the hypothesis that they are melting-pots could not be excluded [[17], p.462]. Indeed, in a more recent extended phylogeographic study of Pinus pinaster [39] it was concluded that Landes was more likely to be a melting-pot.

Genetic diversity spectrum. The genetic diversity spectrum (GDS) for (a) the Landes location and (b) the Pantelleria location in Figure 5. For every possible distance, the number of pairs of haplotypes that are that distance apart is depicted.

Using the same distance matrix, we computed diversity and haplotype connectivity scores for each of the 10 sampling locations as explained in the Methods section (using the measure AD for diversity). These are presented in Table ​ Table2. 2 . Note that, in contrast to [17], our scores do not take into account how often a haplotype was found in a particular location but rather which haplotypes were found.

Table 2

Sampling locationNumber of Haplotypes in regionDiversityHaplotype connectivity
ADADminADmaxAD*HCHCminHCmaxHC*
Landes62.450.337.140.3161100.56
Pantelleria91.670.375.660.2531100.22
Leiria80.730.366.060.0611100.00
Sardinia90.700.375.660.0621100.11
Morocco80.690.366.060.0611100.00
Corsica80.680.366.060.0611100.00
Liguria50.640.318.060.0421110.10
Moncao60.330.337.140.0011100.00
Tuscany50.310.318.060.0011110.00
Alcacier50.310.318.060.0011110.00

Diversity and haplotype connectivity scores for the sampling locations pictured in Figure 5, ranked according to normalized average square-distance diversity score (AD*). The columns labeled with ADmin, ADmax, HCmin and HCmax contain the minimum/maximum score over all subsets containing the same number of haplotypes as found in the region.

As can be seen in Table ​ Table2, 2 , the two locations with highest AD* diversity scores are Landes and Pantelleria. In view of the HC* scores for these locations, this supports the melting-pot scenario, especially for the Landes location. Note that the bimodality of the GDS for the Landes location is also indicative of two clusters of haplotypes having low internal distances and high between cluster distances, which could also be regarded as a signature supporting a melting-pot scenario. However, the shape of the GDS for the Pantelleria location is somewhat less distinctive and so, in this case at least, the haplotype connectivity approach provides some useful additional information.


Concepts and Definitions

Let us start with a brief review on the species diversity (aka community diversity, biodiversity or ecological diversity) to explain the two essential elements of diversity concept in general, which should facilitate the introduction of our SNP diversity and similarity measures below. Species diversity refers to the ecological diversity of species in an ecological community, but diversity concept is equally applicable to genetic diversity (e.g. Nei 1973, Wehenkel et al., Bergmann et al.) 13,23,24 or other entities such as metagenome diversity (Ma and Li) 20 . Conceptually, diversity possesses two essential elements: the variety and the variability of varieties (Gaston Chao et al.) 10,25 . For example, the two elements of species diversity are species (variety) and the variability of species abundances. To quantify the concept of species diversity, one surveys a community (usually by sampling), counts the abundances of each species in the community, and obtains pi = (the relative abundance of species i) = (the number of individuals of species i)/(the total individuals of all species in the community), and also counts the number of species in the community (S). The dataset from such a survey (sampling) is a vector of species abundance in the form of (p1, p2, …, pi, …ps). For such a vector of relative abundances (frequencies), one approach to characterizing it is to fit a statistical distribution, which is known as species abundance distribution (SAD) in community ecology. The most widely used SADs include log-series, log-normal, and power law distributions a common property of SADs is that they are highly skewed, long tail distributions, but rarely follow the normal distribution or uniform distribution. Instead, the SAD is highly aggregated (skewed or non-random), just as the non-random SNP distribution previously mentioned in the introduction section. Although SAD fully describes the species abundance frequency and therefore adequately captures the full characteristics of species diversity, using a SAD to measure diversity fails to present intuitive measures to synthesize the two elements of diversity (i.e., variety and variability). An alternative approach to fitting SAD is to use various diversity metrics (also known as measures or indexes). Numerous diversity metrics for measuring species diversity have been proposed, with Shannon’s entropy being the most well known.

Diversity metrics belong to the so-termed aggregate functions, which combine several values into a single value (Beliakov et al., James) 6,7 . The arithmetic mean (average) is the most commonly utilized aggregation function, but it is a rather poor metric for measuring diversity due to the highly non-random distribution of species abundances. Instead, entropy-based aggregation function is suitable for measuring diversity. The first and also still one of the most widely utilized entropy-based diversity metric is Shannon entropy, which was attributed to Claude Shannon, the co-founder of information theory (Shannon, Shannon & Weaver) 8,26 , but Shannon had never studied biodiversity himself. What happened was that ecologists borrowed the idea from Shannon’s information theory, in which Shannon’s entropy measures the content of information or uncertainty in communication systems. Of course, Shannon’s entropy is indeed sufficiently general for measuring biodiversity because diversity is essentially heterogeneity, and heterogeneity and uncertainty both can be measured by the change of information, i.e., information lowers uncertainty.

Using Shannon entropy as example, species diversity (H), more accurately species evenness, can be computed with the following formula,

where S is the number of species in the community, and pi is the relative abundance of each species in the community. In terms of the variety-variability notion for defining diversity, the variety is the species and variability is the species abundance obviously. In fact, the variety-variability notion can be utilized to define diversity for any systems (not even limited to biological systems) that can be abstracted as the two elements of variety and variability, including SNP diversity, as exposed below.

Definitions for SNP diversities

Using an analogy, a chromosome that has many loci is similar to an ecological community of many species, and each locus may have different number of SNPs. With variety-variability notion for defining diversity, the locus is the variety (similar to species in a community), and the number of SNPs at each locus is the variability (similar to species abundance in a community). Assuming S is the number of loci with any SNP, and pi is the relative abundance of SNPs at locus i (i.e., the number or abundance of SNPs at locus i divided by the total number of SNPs from all loci), then SNP diversity can be measured with Shannon entropy (Eq. 1). Strictly speaking, SNP may also be termed locus diversity, since locus is essentially the ‘habitat’ where SNPs reside. Figure 1 conceptually illustrated the distribution of SNPs on a chromosome specifically how pi is defined and computed.

A conceptual diagram showing the distribution of SNPs on a chromosome with reference to the reference chromosome: the chromosome is similar to an ecological community, and the number of SNPs on a gene locus is similar to the species abundance in an ecological community. For example, there are three SNPs on the locus of gene-1, assuming the total SNPs on the chromosome is N (or 10 displayed with the first 3 genes displayed), then the relative SNP abundance for gene-1 is equal to 3/N (or 3/10 = 0.3 with the 3 genes displayed). Similarly, p2, p3, … can be computed. When the relative abundances of SNPs are available, the diversity (Hill numbers) can be computed based on the diversity definitions [Eqs. (2–15)]. The R-codes computing alpha-diversity, beta-diversity (including similarity) profiles are provided in the OSI.

Although Shannon’s entropy has been widely used for measuring species diversity, a recent consensus among ecologists is that Hill numbers, which are based on Renyi’s general entropy, offer the most appropriate metrics for measuring alpha-diversity and for multiplicatively partitioning beta-diversity (Chao et al. 2012, 2014, Ellison 2010, Kaplinsky & Arnaout) 9,10,12,19 . Given the advantages of Hill numbers over other existing diversity indexes, we believe that the Hill numbers should also be a preferred choice for defining the SNP diversity.

SNP alpha-diversity

Hill numbers were derived by Hill (1973) based on Renyi’s (1961) general entropy 15,16 . Here we propose to apply it for defining the SNP alpha-diversity, i.e.,

where G is the number of gene loci with any SNP, pi is the relative abundance (i.e., the frequency of occurrence) of SNPs at locus i, q = 0, 1, 2, … is the order number of SNP diversity, q D is the SNP alpha-diversity at diversity order q, i.e., the Hill numbers of the q-th order.

The Hill number is undefined for q = 1, but its limit as q approaches to 1 exists in the following form:

The diversity order (q) determines the sensitivity of the Hill number to the relative abundance (i.e., the frequency of occurrence) of SNP. When q = 0, the SNP frequency does not count at all and 0 D = G, i.e., the SNP richness, similar to the species richness in species diversity concept. When q = 1, 1 D equals the exponential of Shannon entropy, and is interpreted as the number of SNPs with typical or common frequencies. Hence, Shannon index is essentially a special case of Hill numbers at diversity order q = 1. When q = 2, 2 D equals the reciprocal of Simpson index, i.e.,

which is interpreted as the number of dominant or very frequently occurred SNPs. Therefore, two most widely used diversity indexes, Shannon index and Simpson index are the special cases, and more accurately, the functions of the Hill numbers.

In general, we need to specify an entity (unit or scope) for defining and measuring SNP diversity. For demonstrative purpose in this article, we choose individual chromosome as the entity for defining SNP diversity, similar to using community for defining species diversity. The general interpretation of diversity of order q is that the chromosome contains q D = x loci with equal SNP frequency. Note that the entity for defining SNP diversity can be other appropriate units such as the whole genome of an organism or segment of chromosome.

The above-defined SNP diversity measures the diversity of SNP on an individual genetic entity (such as chromosome or genome), similar to the concept of alpha diversity in community species diversity, and we term it SNP alpha-diversity. In the following, we define the counterparts of species beta-diversity and gamma-diversity in community ecology for SNPs, i.e., SNP beta-diversity and SNP gamma-diversity.

SNP gamma diversity

While the previously defined SNP alpha-diversity is aimed to measure the SNP diversity within a genetic entity (such as a chromosome or genome), the following SNP gamma-diversity is defined to measure the total SNP diversity of pooled, multiple (N) chromosomes from a population (cohort) of N different individuals, one from each individual but with the same chromosome numbering.

Assuming there are N individuals in a population (cohort), we define the SNP gamma-diversity with the following formula, similar to the species gamma-diversity in ecology (e.g., Chao et al. Chiu et al.) 9,10,27 ,

where (overline<

_>) is the SNP frequency on the i-th locus (i = 1, 2, …,G) in the pooled population of N individuals (termed N-population).

Comparing Eq. (5) for gamma diversity with Eq. (2) for alpha diversity reveals that the gamma-diversity is the Hill numbers based on the SNP frequency at i-th locus in the N-population. Similar to Chao et al. 9,10 Chiu et al. 27 , derivation for species gamma-diversity in ecological community, assuming yij is the SNP frequency at i-th locus of j-th individual, yi+ is the total value of SNP at i-th locus contained in the N individuals, y+j is the total SNP from j-th individual, y++ is the total SNP contained in N individuals, pij is the SNP frequency at i-th locus of j-th individual, wj is the weight of the j-th individual,

it can be easily derived that,

Plug Eq. (6) for (overline<

_>) into the definition of SNP gamma diversity [Eq. (5)], we obtain the following formulae for computing SNP gamma-diversity of N-population as follows:

SNP beta diversity

In community ecology, there are two schemes for defining beta-diversity: one is the additive partition and another is the multiplicative partition of gamma diversity into assumingly independent alpha-diversity and beta-diversity. Recent consensus (e.g., Jost Ellison Chao et al., Gotelli & Chao, Gotelli & Ellison) 9,10,11,12,28,29 recommended the use of multiplicative partition. Let ( (<>^D_) ) and ( (<>^D_) ) are alpha- and gamma-diversity measured with the Hill numbers, respectively, beta-diversity is defined as:

We adopt the exactly same multiplicative partition of the Hill numbers in species diversity for measuring SNP beta-diversity except that both alpha- and gamma- diversities are computed with SNP frequency (relative abundance), rather than with species abundances.

This SNP beta-diversity ( (<>^D_<eta >) ) derived from the above multiplicative partition takes the value of 1 if all communities are identical, and the value of N (the number of individuals in the population) when all individuals are completely different from each other (i.e., no shared SNPs).

Although Eq. (2) correctly defines the SNP alpha-diversity, it requires some adaptations to apply for the partition of gamma diversity in order to obtain beta-diversity with Eq. (9). Similar to the derivation for species alpha diversity as demonstrated in Chiu et al. 27 , we can derive the following formulae for SNP alpha diversity in N-population setting, i.e.,

The computation of SNP beta-diversity can then be accomplished with Eqs. (7–11), i.e., Eqs. (7 and 8) for gamma diversity, (9) for beta-diversity and (10–11) for alpha-diversity.

We define a series of the Hill numbers for SNP diversity at different diversity order q = 0, 1, 2, … as SNP diversity profile, that is, a series of Hill numbers corresponding to different non-linearity levels weighted differently with the SNP frequency distribution.


MANAGEMENT IMPLICATIONS

Our study demonstrates that repeated genetic sampling across biologically relevant time scales (e.g., 10 generation intervals or

5 years for the bilby) will allow stakeholders to assess whether established management practices are sufficient to maintain genetic diversity at levels comparable to ancestral populations. Where excessive loss of heterozygosity is a concern (i.e., genetic diversity has dropped to a level significantly lower than that of the founding group[s]), the translocation of individuals from genetically divergent populations at a rate of 1–10 migrants per 1–2 generations should be sufficient to mitigate the worst effects of inbreeding, while ensuring that gene flow does not fully disrupt genetic distinctiveness among the individual management units. Where practical, translocation programs should introduce individuals from genetically divergent populations that occupy comparable environments and climatic zones.

When the number of founders is high, it may be advantageous to split captive or conservation-fenced metapopulations into several discrete management units or subpopulations. Although this can lead to a short-term loss of heterozygosity at the level of the individual management units, the crossing of representatives from different sub-populations during or within 1 generation of translocation should be sufficient to restore genetic diversity to pre-fragmentation levels, preserving the adaptive potential of the species as a whole. Maintaining the metapopulation as series of discrete management units will also ensure that existing genetic reserves can be used for future translocations, minimizing the likelihood of bottlenecking events due to the repeated movement of individuals from the same source populations.


Watch the video: Naturfag DNA (August 2022).