Accounting for population structure and data quality in demographic inference with linkage disequilibrium methods

Abstract

Linkage disequilibrium methods for demographic inference usually rely on panmictic population models. However, the structure of natural populations is generally complex and the quality of the genotyping data is often suboptimal. We present two software tools that implement theoretical developments to estimate the effective population size (N_e): GONE2, for inferring recent changes in N_e when a genetic map is available, and currentNe2, which estimates contemporary N_e even in the absence of genetic maps. These tools operate on SNP data from a single sample of individuals, and provide insights into population structure, including the F_ST index, migration rate, and subpopulation number. GONE2 can also handle haploid data, genotyping errors, and low sequencing depth data. Results from simulations and laboratory populations of Drosophila melanogaster validated the tools in different demographic scenarios, and analysis were extended to populations of several species. These results highlight that ignoring population subdivision often leads to N_e underestimation.

Fig. 1. currentNe2 estimates in different simulated scenarios.

The top plot in each panel shows estimates of contemporary effective population size (N_e) for the entire metapopulation (blue dots), the sum of the effective sizes of the subpopulations (N_T, red dots), and N_e assuming panmixia (green squares; green arrows indicate values below the scale limit). Note that the labels ‘N’ on the x-axes here refer generically to the three estimates. The three lower plots show estimates for migration rates (m), genetic differentiation indices (F_ST), and the number of subpopulations (s). All estimates were calculated using the currentNe2 under the assumption of population subdivision (option ‘-x’), except for N_e under panmixia. Each estimate represents the geometric mean of ten simulated metapopulations, with 95% confidence intervals. True metapopulation values are indicated by grey dashed lines (N_e and N_T lines often overlap, with N_e typically higher). A Effect of varying migration rates in a metapopulation with two subpopulations of 1000 individuals each. B Effect of increasing the number of equally sized subpopulations while maintaining a total metapopulation size of 2000 individuals (m = 0.001). C Effect of unequal sampling proportions between two subpopulations of 1000 individuals each (m = 0.001). D Effect of unequal subpopulation sizes in a metapopulation with two subpopulations totalling 2000 individuals (m = 0.001). All estimates are based on random samples of 100 individuals from the entire metapopulation, except in (C), where sampling proportions vary. Source data are provided as a Source Data file.

Fig. 2. GONE2 estimates of the historical N_e in different simulated scenarios.

Four different profiles of demographic change are represented in the columns by the grey-shaded areas: constant population size, recent decline, recent expansion, and bottleneck. The scenarios are also arranged in rows in four different situations. First row: Estimates from panmictic populations using the default option of the program (green lines) and the ‘-x’ option, which incorrectly assumes that the population is divided in two subpopulations (red lines). Second row: Estimates from metapopulations composed of two subpopulations of equal size N and migration rate m = 0.001 per generation. The analyses consider either a metapopulation divided into two subpopulations (option ‘-x’, red lines) or the incorrect assumption of a single panmictic population (default program options, green lines). Third row: Estimates from low depth of sequencing data (average depth 2x) of panmictic populations. Green lines show estimates not corrected for low sequencing depth (default program options), while red lines show estimates corrected for low sequencing depth (option ‘-g 3’ to be used for any sequencing depth). Fourth row: Estimates from sequencing data with 5% base calling errors. Green lines represent unadjusted estimates not corrected for errors (program default option), while red lines represent estimates corrected for the true 5% error rate (option ‘-b 0.05’). Simulation results are the geometric mean of 20 replicates. The shaded areas around the estimate lines represent the 95% interval of the distribution of the estimates. Source data are provided as a Source Data file.

Fig. 3. Historical N_e estimates for different species using GONE2.

Estimates of N_e trajectories over the last 50 generations for the same samples analyzed in Table 1. Each panel corresponds to a different species, with estimates obtained using GONE2 under two assumptions: the default option (green lines), which assumes panmixia, and the ‘-x’ option (red lines), which takes population structure into account. The salmon panel combines the analyses of two tributaries of the river Dee. Seabream panel combines the analysis from two parallel breeding programmes. GONE2 did not produce results for cattle for more than five generations, probably due to the complexity of the mixing of 134 highly differentiated breeds. Source data are provided as a Source Data file.

Fig. 4. Estimates of the historical N_e using GONE2 with samples from 17 males of Drosophila melanogaster.

The experimental populations were maintained at constant size (A), underwent a recent drastic reduction in size (B), recovered from a previous bottleneck (C), and underwent further expansion after a bottleneck (D), as shown by the grey shaded areas. The green lines show the estimates of N_e from autosomes using the modifier ‘-g 0’ for unphased diploid genotypes. The red lines show the estimates of N_e from X chromosomes using the haploid option ‘-g 1’. The N_e estimates from X chromosomes have been transformed to the autosomal scale (see Methods) to compare sex ratios. Source data are provided as a Source Data file.