Detecting Selection in Sequences

Detecting selection from sequence data is fundamentally a problem of distinguishing signal from a well-characterized statistical null. The null model is the neutral theory of molecular evolution: most substitutions observed between species, and most polymorphisms observed within a species, are not driven by selection but by the joint action of mutation and genetic drift acting on selectively neutral or nearly neutral variants (building on Nearly Neutral Theory). Under this null, the rate of nucleotide substitution is expected to equal the neutral mutation rate, independent of population size, and the frequency distribution of neutral alleles within a population follows a specific, calculable shape. Selection breaks these expectations in directionally diagnostic ways. Purifying (negative) selection removes deleterious mutations before they can drift to appreciable frequency, so it depresses substitution rates and skews allele frequencies toward rare classes. Positive (directional) selection drives an advantageous variant to fixation faster than drift would, elevating substitution rates at the selected site and dragging linked neutral variation along with it — a phenomenon called a selective sweep, which leaves a localized trough of genetic diversity around the selected locus. Balancing selection, including heterozygote advantage and negative frequency-dependence (building on Balancing Selection and Heterozygote Advantage), maintains multiple alleles at intermediate frequency for long periods, elevating diversity and skewing the frequency spectrum toward intermediate-frequency classes rather than rare or fixed ones. The three test families surveyed here — dN/dS, McDonald-Kreitman, and site-frequency-based tests — each exploit a different observable consequence of these processes, operate over different evolutionary timescales (dN/dS looks across millions of years of divergence; SFS-based tests look within a single population over thousands of generations), and are vulnerable to different confounders (demographic history, recombination rate variation, mutation rate heterogeneity, and linked selection). Understanding all three, and how they complement each other, is essential for making a credible claim that a particular gene, or a particular nucleotide, has been shaped by natural selection.

dN/dS: comparing rates of amino-acid-changing and silent substitution

The logic of dN/dS begins with the distinction between synonymous (silent) substitutions, which change a codon's nucleotide sequence without changing the encoded amino acid, and nonsynonymous (replacement) substitutions, which do change the amino acid. Synonymous substitutions are, to a first approximation, invisible to natural selection at the protein level: they do not alter the phenotype of the protein product, so their fixation rate should track the underlying neutral mutation rate reasonably closely (though codon usage bias, splicing signals, and mRNA stability can impose weak selection even here — a caveat worth flagging). Nonsynonymous substitutions, in contrast, change the protein and are therefore exposed to selection. Comparing the two rates isolates the effect of selection on protein function from the background mutation rate, because both classes of site experience the same mutation process but only one is masked from selection's view.

The procedure: 1. Align homologous protein-coding sequences from two or more species. 2. For each codon, count the number of synonymous sites (S) and nonsynonymous sites (N) — fractional counts, since a single codon can offer both types of mutational path depending on which nucleotide changes. 3. Count the observed synonymous differences (Sd) and nonsynonymous differences (Nd) between the sequences. 4. Calculate dS = Sd / S (synonymous substitutions per synonymous site) and dN = Nd / N (nonsynonymous substitutions per nonsynonymous site). 5. Compute ω = dN / dS. Interpretation of ω:

Because ω is a ratio between two counts made from the same sequence, it automatically controls for local mutation rate, sequence length, and time since divergence — this is what makes it powerful for comparing genes across a genome or across lineages. Most protein-coding genes, when averaged over their entire length, show ω well below 1 (typically 0.1–0.3), reflecting the fact that most positions in most proteins are structurally or functionally constrained. Genes or codon positions with ω > 1 are rare and biologically striking: they mark sites actively being reshaped by selection, often in genes involved in immune recognition, sperm-egg interaction, or host-pathogen arms races. A crucial refinement is that ω averaged over an entire gene can mask localized positive selection at a few codons buried within an otherwise strongly constrained protein. Maximum-likelihood methods (implemented in software such as PAML's codeml) allow ω to vary across codon sites and across lineages on a phylogeny, dramatically increasing power to detect selection acting on just a handful of residues — exactly the pattern expected when, say, only the antigen-binding cleft of an immune receptor is under diversifying selection while the rest of the protein remains functionally constrained.

McDonald-Kreitman: polymorphism versus divergence

The dN/dS test compares two substitution rates between species but says nothing about what is happening within a species — and it can be fooled by weakly deleterious mutations that occasionally fix by drift, inflating dN without any positive selection at all. The McDonald-Kreitman (MK) test, published in 1991, solves this by adding a second axis of comparison: polymorphism within species versus fixed differences (divergence) between species. The logic: if a class of mutations (synonymous or nonsynonymous) is truly neutral, the ratio of polymorphic-to-divergent counts should be the same whether you look at synonymous or nonsynonymous sites — because both within-population segregation and between-species fixation are governed by the same neutral parameters (mutation rate and effective population size cancel out of the ratio). Build a 2×2 contingency table:

Under strict neutrality, Pn/Ps = Dn/Ds. A statistically significant departure from this equality — tested with a chi-square or Fisher's exact test — indicates selection, and the direction of the departure tells you which kind:

• Dn/Ds > Pn/Ps (an excess of nonsynonymous fixed differences relative to nonsynonymous polymorphism): adaptive substitutions are being fixed by positive selection faster than neutral theory predicts, because advantageous nonsynonymous mutations sweep to fixation without lingering as polymorphism.
• Pn/Ps > Dn/Ds (an excess of nonsynonymous polymorphism relative to fixed differences): purifying selection is keeping mildly deleterious nonsynonymous variants segregating at low frequency without letting them ever reach fixation, or balancing selection is preserving nonsynonymous polymorphism.

A useful derived quantity is the neutrality index (NI) = (Pn/Ps) / (Dn/Ds). NI < 1 indicates an excess of adaptive fixations (positive selection); NI > 1 indicates an excess of segregating deleterious variation (purifying selection dominating the polymorphism pool). From NI one can estimate α, the proportion of nonsynonymous substitutions fixed by positive selection:

α = 1 − NI = 1 − (Pn·Ds)/(Ps·Dn) This single number, applied genome-wide, has become one of the most cited statistics in molecular evolution for estimating how much of a genome's adaptive protein evolution is driven by natural selection versus drift.

Site-frequency spectrum (SFS) tests

The third family looks only within a single species, using the site-frequency spectrum: for every polymorphic site in a sample of sequences, record how many copies of the derived (or minor) allele are present, then tabulate this into a histogram across all polymorphic sites. Under the standard neutral model (constant population size, no selection, free recombination), this spectrum has a precisely predictable shape: an excess of rare-frequency variants, following roughly a 1/i distribution (where i is allele count), because most mutations are young and have not had time to drift to higher frequency. Two classic summary statistics compare features of this spectrum:

Watterson's θw, based purely on the number of segregating sites, and π (nucleotide diversity, or "pairwise"), based on the average number of pairwise differences between sequences, weighting each site by its allele frequency. Under strict neutrality with constant population size, both statistics are unbiased estimators of the same underlying parameter, θ = 4Neμ (for diploids), and should agree.

Tajima's D, introduced by Fumio Tajima in 1989, formalizes the comparison:

D = (π − θw) / √(variance of the difference)

• D ≈ 0: consistent with neutral equilibrium.
• D < 0: an excess of rare (low-frequency) variants relative to neutral expectation. This is the signature of a recent selective sweep (a beneficial mutation dragging linked neutral variants to fixation, then new mutations quickly restoring rare variants) or of purifying selection removing intermediate/high-frequency deleterious variants, or of recent population expansion — a crucial demographic confound.
• D > 0: an excess of intermediate-frequency variants. This is the signature of balancing selection maintaining multiple alleles at appreciable frequency for long periods, or of a recent population bottleneck — again, a demographic confound that must be ruled out independently.

Because Tajima's D cannot distinguish selection from demography using a single locus, robust inference requires comparing the candidate locus against the genome-wide background (which reflects the shared demographic history of the whole population) — an outlier locus relative to that background is much better evidence of selection than an outlier relative to the neutral theoretical expectation alone. A complementary, spatially explicit signature of a selective sweep is a valley of reduced diversity flanking the selected site, because a beneficial mutation rising rapidly in frequency drags along the haplotype background on which it arose, eliminating variation at linked neutral sites through genetic hitchhiking. The width of this valley depends on the strength of selection and the local recombination rate — strong selection and low recombination produce wide swept regions, which is why genomic scans for selection often combine diversity troughs with haplotype-based statistics (such as extended haplotype homozygosity, EHH) that detect unusually long, unbroken haplotypes around a recently and rapidly fixed allele.

McDonald and Kreitman (1991) established the polymorphism/divergence framework using the Adh (alcohol dehydrogenase) locus in three species of Drosophila: D. melanogaster, D. simulans, and D. yakuba. They sequenced Adh alleles within and between species and tabulated the now-classic 2×2 table of synonymous versus nonsynonymous, polymorphic versus fixed changes. Their actual data showed:

A neutral model predicts these ratios should match: 7/17 (fixed) should equal roughly 2/42 (polymorphic). Instead, nonsynonymous fixed differences were dramatically overrepresented relative to nonsynonymous polymorphism — a highly significant departure (Fisher's exact test) that could only be explained by adaptive fixation of amino-acid-changing mutations at Adh. This was one of the first pieces of direct molecular evidence for recurrent positive selection acting on a specific, well-characterized protein-coding gene, and it did so without requiring any external calibration of the mutation rate, because the test is entirely internal to the polymorphism/divergence comparison.

Hudson, Kreitman, and Aguadé (1987), in the paper that immediately preceded and helped motivate the MK framework, first formalized the logic of comparing observed polymorphism to a neutral expectation calibrated from divergence (the HKA test), applying it also to Adh and showing reduced polymorphism relative to divergence at that locus in a pattern consistent with a recent selective sweep — an early demonstration that within-species variation, not just between-species substitution, carries a testable neutral expectation.

Tajima (1989) derived the statistical distribution of D under the standard neutral coalescent model, providing the null distribution against which real data could be judged, and demonstrated using both simulated and real human mitochondrial DNA and Drosophila data that significant deviations could be reliably detected with modest sample sizes. This turned the site-frequency spectrum from a descriptive curiosity into a rigorous hypothesis-testing framework, and it remains the entry point for essentially every modern genome-wide scan for selection (building directly toward the coalescent-based reasoning used in more sophisticated modern tests). A landmark applied confirmation came from studies of the human sickle-cell allele (HbS) at the β-globin locus. Anthony Allison's fieldwork in East Africa in the 1950s first established the ecological correlation between malaria prevalence and sickle-cell heterozygote frequency, and decades later, molecular surveys of the haplotype background surrounding HbS found strikingly reduced diversity and a single dominant haplotype extending over a wide genomic region — the textbook signature of a rapid selective sweep, later reinforced by the observation that at least four to five independent HbS haplotypes arose on different genetic backgrounds in different regions of Africa and Asia, each showing its own localized sweep signature. This combination of a classical fitness experiment (heterozygote survival advantage against malaria) with a molecular sweep signature (haplotype homozygosity and reduced flanking diversity) is one of the cleanest examples in human genetics of converging phenotypic and molecular evidence for balancing selection maintaining a polymorphism that is, at the same time, the product of repeated local positive selection events.

Example 1: A McDonald-Kreitman table with asymmetric counts. Suppose a population genomics survey of a stickleback immune gene yields the following counts comparing a freshwater and a marine population as the "within-species" polymorphism sample, and a distantly related species as the outgroup for fixed differences:

Under strict neutrality, Dn/Ds should equal Pn/Ps. Dn/Ds = 24/12 = 2.0 Pn/Ps = 9/41 ≈ 0.22 These are wildly different — nonsynonymous fixed differences are about 9-fold overrepresented relative to what neutral expectation predicts from the polymorphism data. Compute the neutrality index: NI = (Pn/Ps) / (Dn/Ds) = 0.22 / 2.0 ≈ 0.11 Because NI is far below 1, this indicates a strong excess of adaptive nonsynonymous fixation. Estimate α, the fraction of nonsynonymous substitutions driven to fixation by positive selection: α = 1 − NI = 1 − 0.11 = 0.89 Interpretation: roughly 89% of the amino-acid-changing substitutions that became fixed between these lineages are estimated to have been fixed by positive selection rather than drift — consistent with an immune gene locked in an evolutionary arms race with pathogens, where new resistance alleles are repeatedly favored and swept to fixation, leaving little standing nonsynonymous polymorphism behind (small Pn) while accumulating many nonsynonymous fixed differences (large Dn). A Fisher's exact test on the 2×2 table (values 24, 12, 9, 41) would return a highly significant p-value, confirming this is not a chance fluctuation.

Example 2: Tajima's D at a candidate sweep locus, checked against genomic background. A study of a plant population resistant to a fungal pathogen sequences a 5-kb window around a candidate resistance gene in 40 chromosomes and finds:

• Number of segregating sites, S = 18
• Watterson's estimator, θw ≈ 4.9 (using S divided by the harmonic-number correction factor for n = 40)
• Nucleotide diversity, π ≈ 1.6

Tajima's D is strongly negative here because π « θw: there are many segregating sites (so θw, which only counts how many sites vary, is moderate), but most of the variants are at very low frequency, so pairwise diversity π (which weights by allele frequency and is dominated by common variants) is much smaller than θw predicts under neutrality. A formal calculation gives D ≈ −2.1, which falls below the 5th-percentile critical value for n = 40 under the standard neutral model (roughly −1.7 to −2.0 depending on θ), suggesting statistical significance. But before concluding "selective sweep," compare against genome-wide background: if the same population shows a genome-wide average Tajima's D of −1.3 (because the population recently expanded, which also produces excess rare variants genome-wide), the resistance-locus value of −2.1 is only modestly more extreme than the background, and a formal outlier test (e.g., comparing the empirical percentile of this window against thousands of other windows simulated or observed under the same demographic history) might find it is not significant after all — or might still be significant if it sits far in the tail. This example illustrates why single-locus Tajima's D can never distinguish selection from demography on its own; the diagnostic power comes from treating the genome-wide distribution as the empirical null and asking whether a specific locus is an outlier relative to that, ideally combined with an independent signal such as a haplotype-homozygosity peak or an MK-test excess of fixed nonsynonymous differences at the same gene.

John H. McDonald and Martin Kreitman brought together two threads that had previously been separate: the codon-based logic of counting synonymous and nonsynonymous changes (which had existed since the 1970s work of Motoo Kimura and later Miyata and Yasunaga), and the population-genetic idea of comparing within-species polymorphism to between-species divergence (which Kreitman had helped pioneer with Hudson and Aguadé in 1987). Kreitman's earlier 1983 study, sequencing the complete Adh gene in eleven Drosophila melanogaster lines, had itself been a landmark: it was one of the first direct DNA sequencing surveys of natural population variation, revealing far more silent polymorphism than the classical allozyme methods of the day had detected, and setting the stage for finer-grained tests of neutrality once multi-species sequence comparisons became feasible. The 1991 McDonald-Kreitman paper condensed this lineage of work into an elegant, internally-controlled statistical test that required no assumptions about absolute mutation rates — a design choice that is a major reason the test remains in routine use, essentially unchanged, more than three decades later.

Genome-wide scans built on these three test families now underpin practical work across medicine, agriculture, and conservation. In human genetics, dN/dS and haplotype-based sweep scans identified the malaria-resistance signatures around the HbS, G6PD, and DARC (Duffy antigen) loci, and site-frequency and haplotype methods mapped the EPAS1 high-altitude adaptation in Tibetan populations, a variant later shown to have introgressed from Denisovans. In cancer genomics, dN/dS-style analyses distinguish driver mutations (an excess of nonsynonymous changes relative to neutral expectation within tumor sequencing cohorts, indicating positive selection for cancer-promoting changes) from passenger mutations that accumulate neutrally alongside them, directly informing which mutations are prioritized for targeted therapy development. In infectious disease surveillance, codon-based ω models applied to viral phylogenies (influenza hemagglutinin, HIV envelope, and SARS-CoV-2 spike protein) flag antigenic sites under diversifying selection from host immune pressure, guiding vaccine strain selection. In agriculture, MK-style and SFS-based scans applied to crop and livestock genomes identify loci under selection during domestication, distinguishing genuine targets of artificial selection from neutral hitchhikers, which sharpens marker-assisted breeding programs. In conservation genetics, genome-wide Tajima's D and diversity scans help assess whether small, bottlenecked populations retain adaptive genetic variation or have lost it to drift, informing decisions about genetic rescue and captive breeding management. All of these applications rest on the same foundational logic developed for Adh and human hemoglobin: compare what selection predicts to what neutral drift alone would produce, and let the sequence data adjudicate — a logic that, when extended to genealogies of whole populations, leads directly into the machinery of The Coalescent: An Introduction.

Ask the tutor: How do biologists distinguish a selective sweep from a population bottleneck when both produce a negative Tajima's D across the whole genome?