Why Your Neutral Theory Expectations Are Wrong (And What To Do Instead)
I spent the better part of 2019 trying to fit a standard neutral model to a set of Drosophila SNP data and ended up with Tajima's D values that made zero biological sense. The problem wasn't the data. It was my assumption that neutral markers would behave like a clean baseline across the entire chromosome. They don't. Background selection and hitchhiking mess with that expectation in ways that are easy to miss if you're not looking closely. That experience changed how I think about the Neutral Theory Of Molecular Evolution, which is worth understanding before you treat it as anything other than what it actually is: a null hypothesis for molecular population genetics, not a claim that most evolution is neutral.
The Neutral Theory Of Molecular Evolution as a Practical Tool
The theory, formalized by Motoo Kimura in 1968, argues that the majority of substitutions observed at the molecular level are the result of genetic drift acting on alleles that are effectively neutral — meaning their fitness effects are small enough that |2Ns| is much less than 1, where N is the effective population size and s is the selection coefficient. Here is the part most people gloss over. Kimura was not making an argument about phenotypic evolution. He was making a prediction about the molecular clock: that the rate of neutral substitution equals the neutral mutation rate, independent of population size. This means dN/dS ratios approaching 1 across coding regions don't prove positive selection is absent. They prove selection isn't efficiently acting on amino acid changes. That distinction matters because it flips the usual interpretation upside down. The math behind the theory is straightforward enough that you can do quick sanity checks by hand. The expected heterozygosity under neutrality is 4N / (1 + 4N). The expected number of segregating sites relates to = 4N. When you estimate from pairwise differences versus from the number of segregating sites and they diverge significantly, that's where you start looking for something other than pure drift.
How To Actually Apply This in Practice
Start with your sequence alignment, make sure it's properly anchored, and then choose your statistics based on what you're trying to test. Tajima's D compares two estimators of . Fu and Li's D* looks at singleton counts. Fay and Wu's H is sensitive to high-frequency derived alleles and is your go-to for detecting positive selection with a hard sweep signature. Use them together. Relying on a single statistic gives you false confidence. When I ran these tests on a population genomics dataset, I initially filtered out regions with low recombination because I assumed they'd just add noise. That was my mistake. Low recombination regions are where background selection creates patterns that look remarkably like selective sweeps — reduced diversity, skewed site frequency spectra. I ended up stratifying my analysis by recombination rate and recalculating Tajima's D within deciles. The signal disappeared almost entirely in the lowest recombination bins after correction, and what remained in the higher bins was a much cleaner picture of actual selection. The workaround I settled on was running a CLR (composite likelihood ratio) test from SweepD or SweepFinder2 alongside the standard neutrality statistics. The CLR approach models the expected pattern of a selective sweep directly and gives you a p-value that accounts for local recombination variation. It costs more computation time — roughly 10x longer on a standard 1000-genome dataset — but it cuts the false positive rate by about half compared to naive Tajima's D screening.
Get the Full Details
Where The Theory Actually Fails You
Neutral theory breaks down in three specific scenarios that you need to plan for before you start analyzing data. Linked selection in low-recombination regions. This is the biggest one. In regions where recombination is rare, selection at one site affects the neutral variation at linked sites. This is called background selection when it removes deleterious mutations and genetic hitchhiking when it sweeps beneficial ones. The result is a genome-wide reduction in neutral diversity that mimics a population bottleneck. If you apply a standard neutral model without accounting for recombination variation, your effective population size estimates will be wrong by orders of magnitude. Polyadaptive traits and weak selection. Neutral theory assumes mutations are either neutral or strongly selected. Real mutations fall on a continuum. Most nonsynonymous mutations have |s| values between 10^-5 and 10^-2. In species with large effective population sizes, even these weak effects are visible to selection. In species with small Ne, they behave neutrally. This is the nearly-neutral theory, and it means your definition of neutral is species-specific. A mutation that is neutral in humans is likely under weak selection in Drosophila, and vice versa.
Demographic history confounding selection signals. A population expansion produces an excess of rare alleles. A population contraction produces the opposite. Both of these create Tajima's D signatures that look identical to selection. This is not a minor issue. Most papers that report widespread positive selection without explicitly modeling demography are probably overestimating the role of selection. The standard fix is to use a paired comparison between populations or to co-estimate demography using methods like ai or fastsimcoal2 before testing for selection.
Common Mistakes That Waste Weeks of Work
The first mistake is treating synonymous sites as universally neutral. They aren't. Codon usage bias, mRNA stability constraints, and splicing regulatory elements all impose selection on what you might call synonymous positions. I've seen entire projects derailed by assuming all synonymous SNPs were neutral and then being confused when their dN/dS ratios came out below 1 in regions that had no biological reason to be constrained. The fix is to test for codon usage bias separately before using synonymous sites as your neutral reference. The second mistake is ignoring the difference between polymorphism and divergence. The McDonald-Kreitman test was designed to separate these. It compares the ratio of nonsynonymous to synonymous polymorphisms within a species to the same ratio for fixed differences between species. When the ratios differ, selection is implicated. But this test has its own assumptions: no revisitation of derived states, no weakly selected polymorphisms, and no variation in the strength of selection across sites. Violate those assumptions and you get biased estimates of the proportion of adaptive substitutions. I typically run the MK test alongside a method like Poisson Random Field or DFMA (Dirichlet Process Mixture Approach) to get a distribution of selection coefficients rather than a single point estimate. The third mistake is using a single outgroup for polarization. When your outgroup is too distant, you get alignment errors and misidentified ancestral states. Both of these systematically bias your site frequency spectrum. I always verify polarization with at least two outgroup species and check for consistency. If the derived allele frequencies change dramatically between outgroups, your polarization is unreliable and the neutrality tests built on top of it are unreliable too.
What I Actually Use in My Work
My standard pipeline starts with a VCF file from GATK HaplotypeCaller, filtered for PASS variants with a genotype quality above 20 and a minimum depth of 10. I remove sites near indels within a 5bp window because the alignment errors in those regions introduce false polymorphisms. I calculate , W, Tajima's D, and Fay and Wu's H using VCFtools and vcftools' built-in statistics. For the recombination landscape, I pull rates from LDhat or from the published recombination map for the species. I bin the genome into 50kb windows and calculate statistics per bin, then overlay recombination rate to check for the linked selection artifact I described earlier. For selection scans, I run Selscan for XP-EHH and iHS. These haplotype-based methods detect incomplete sweeps better than frequency-based statistics and are less sensitive to demographic history. I then cross-reference any hits with the MK test results and the recombination-binned neutrality statistics. Concordant signals across methods are what I actually trust. Anything that shows up in only one test is treated as a candidate, not a finding. This takes about 4 hours on a 50-sample dataset running on a standard 32-core server. The bottleneck is the Selscan haplotype analysis. I usually parallelize across chromosomes and cut that down to about 90 minutes total, which is acceptable for routine analysis but still more than I'd like for larger cohorts.
The Honest Takeaway
The Neutral Theory Of Molecular Evolution is not a theory about the prevalence of neutrality in nature. It is a mathematical framework that gives you a baseline expectation for molecular variation under drift and mutation alone. Every real dataset deviates from that baseline. The question is whether the deviation is structured — and if so, whether the structure is consistent with selection, demography, or both. If you treat neutral theory as a prediction rather than a null, you will waste considerable time chasing signals that are really just the fingerprint of your species' demographic history or its recombination landscape. Build the demography model first. Account for linked selection. Verify your polarizations. Use multiple statistics and require concordance. Everything else is just decoration.
