Controlling type-I error of the McDonald-Kreitman test in genomewide scans for selection on noncoding DNA

Abstract

Departures from the assumption of homogenously interdigitated neutral and putatively selected sites in the McDonald-Kreitman test can lead to false rejections of the neutral model in the presence of intermediate levels of recombination. This problem is exacerbated by small sample sizes, nonequilibrium demography, recombination rate variation, and in comparisons involving more recently diverged species. I propose that establishing significance levels by coalescent simulation with recombination can improve the fidelity of the test in genomewide scans for selection on noncoding DNA.

Figure 1.—

The MK test applied to genomic regions with different spatial organization of neutral (n) and selected (X) polymorphic and divergent mutations. D shows the spatial organization of neutral and selected sites in the study of Jeong et al. (2008). The tan gene cis-regulatory element, tMSE, regulates the expression of abdominal pigmentation in D. yakuba and is inferred to have an excess of divergence relative to polymorphism when compared to synonymous sites from two neighboring genes.

Figure 2.—

Performance of the MK test when comparing neutral sites to closely linked putatively selected sites (Figure 1B). Tests were only performed on MK tables where each marginal sum exceeded five counts. ρ/θ is the ratio of the population recombination rate ρ = 4N_erL to the population mutation rate, θ = 4N_eμL, where N_e is the species effective population size, L is the locus length in base pairs, and r and μ are the recombination and mutation rates per site per generation. The standard neutral model is simulated with θ = 1.5 (light shading), θ = 3 (shaded), θ = 6 (dark shading), and θ = 12 (solid). Sample size, n = 12 and the divergence time, T, were set to 2 (in units of 4N_e generations) in each case. All points are based on 10,000 replicates of the neutral coalescent with recombination implemented using the program ms, except θ = 1.5 (100,000 replicates).

Figure 3.—

The effect of various factors on the performance of the MK test in the context of Figure 1B. Tests were only performed on MK tables where each marginal sum exceeded five counts. (A) The effect of sample size. For open squares, θ = 3 (shaded line) and θ = 12 (solid line) are plotted for ρ/θ = 1. For solid squares, θ is adjusted with sample size such that the mean observed number of polymorphisms, E(S), is the same as for n = 12. (B) The effect of demography. The sample size, n = 12. θ and T are varied such that observed average number of polymorphic and divergent sites [E(S) ∼ 4.5 and E(D) ∼ 12, respectively] is similar for all models and corresponds to expectations for θ = 1.5 and T = 2 under the SNM. Models: SNM, the standard neutral model (light shading, ms command line –t 1.5); PS, a two-deme population structure model with all individuals sampled from one deme (shaded, ms command line –t 0.825–I 3 12 0 1–ma × 1 0 1 × 0 0 0 × –ej 4 2 1 –ej 4 3 1); BN, a bottleneck model (dark shading, ms command line –t 4.3 -en 1 0.0041 0.03 -en 1 0.019 1 –ej 0.45 2 1 –ej 0.45 3 1) based on parameters estimated for D. melanogaster by Thornton and Andolfatto (2006); BN + PS, a combined population structure and bottleneck model (solid, with parameters as above, except –t 2.4 –ej 0.7 2 1 –ej 0.7 3 1). (C) The effect of divergence depth. All parameters are as in B, except species divergence time is varied. For SNM, T = 2 and for BN + PS, T = 0.7. (D) The effect of an intron or a recombination hotspot. The two demographic models shown are the SNM with θ = 6 (squares) and the BN + PS model (circles) from B. For the intron model (solid lines), a 1-kb intron is added in the middle of a 1200-bp sequence. Note that recombination rate is uniform and the intron is not included in the MK test but changes the cumulative genetic distance between the tested regions. For the hotspot model (dark shaded lines), a 100-bp hotspot with 10-fold intensity relative to the background in the middle of the surveyed region was modeled (implemented with msHOT (Hellenthal and Stephens 2007), with locus length = 1200 bp and command line: -v 1 550 650 10). In the hotspot model, background levels of recombination were set such that the cumulative genetic distance of the surveyed region is the same in the presence and absence of the hotspot. For each demographic model, results for “no intron” and “no hotspot” are indicated with light shaded lines for comparison.