Bottlenecks and genetic diversity

population genetics

ecology tutorial

ggplot2

Push a simulated population through a bottleneck in R and watch which measures of diversity notice: allelic richness, heterozygosity and the M ratio disagree.

Author

Tidy Ecology

Published

2026-08-12

A population that spends five generations at ten individuals and then recovers carries a genetic record of those five generations for a long time. Which record, exactly, depends on what you measure. Allelic richness, heterozygosity and the ratio of allele number to allele size range respond at different speeds and to different features of the crash, and a monitoring programme that tracks only one of them can miss a bottleneck or misdate it.

This tutorial builds a stepwise-mutation population, sends replicate copies through bottlenecks of different depths, and measures all three.

A population with alleles worth losing

Microsatellite-style markers mutate by adding or removing a repeat unit, so alleles are integers and the mutation model is a random walk. The simulation below tracks allele sizes for every gene copy at 60 loci in 300 diploids, with the mutation rate set so that the scaled diversity theta = 4Nu is 10, which is the range typical of real microsatellite panels.

library(ggplot2)

te_pal <- list(forest = "#275139", green = "#2f8f63", sage = "#93a87f",
               clay = "#b5534e", gold = "#cda23f", line = "#dad9ca",
               ink = "#16241d")

theme_te <- function() {
  theme_minimal(base_size = 12) +
    theme(panel.grid.minor = element_blank(),
          panel.grid.major = element_line(colour = "#e7e6dc"),
          plot.background = element_rect(fill = "white", colour = NA),
          panel.background = element_rect(fill = "white", colour = NA),
          plot.title = element_text(face = "bold", colour = te_pal$ink),
          axis.title = element_text(colour = "#2c3a31"))
}

step_gen <- function(pop, n_next, u) {
  n_now <- nrow(pop) / 2
  n_loci <- ncol(pop)
  pick <- function(par) {
    rows <- 2 * par - 1 + matrix(sample(0:1, n_next * n_loci, replace = TRUE),
                                 n_next, n_loci)
    matrix(pop[cbind(as.vector(rows), rep(seq_len(n_loci), each = n_next))],
           n_next, n_loci)
  }
  h1 <- pick(sample(n_now, n_next, replace = TRUE))
  h2 <- pick(sample(n_now, n_next, replace = TRUE))
  out <- matrix(0L, 2 * n_next, n_loci)
  out[seq(1, 2 * n_next, by = 2), ] <- h1
  out[seq(2, 2 * n_next, by = 2), ] <- h2
  mut <- runif(length(out)) < u
  out[mut] <- out[mut] + sample(c(-1L, 1L), sum(mut), replace = TRUE)
  out
}

stats_pop <- function(pop) {
  k <- apply(pop, 2, function(col) length(unique(col)))
  spread <- apply(pop, 2, function(col) diff(range(col)) + 1)
  het <- apply(pop, 2, function(col) 1 - sum((table(col) / length(col))^2))
  c(richness = mean(k), spread = mean(spread), m_ratio = mean(k / spread),
    het = mean(het))
}

set.seed(467)
n_anc <- 300
n_loci <- 60
theta <- 10
u <- theta / (4 * n_anc)
pop <- matrix(as.integer(round(rnorm(2 * n_anc * n_loci, 0, sqrt(theta / 2)))),
              nrow = 2 * n_anc)
for (g in seq_len(150)) pop <- step_gen(pop, n_anc, u)
print(round(stats_pop(pop), 4))

richness   spread  m_ratio      het 
  9.8333  10.1000   0.9740   0.8013

The ancestral population carries 9.83 alleles per locus spread over 10.1 repeat units, so the ratio of allele number to size range is 0.974, and gene diversity is 0.8013. That is the reference state. Every scenario below starts from this same population, which is why the comparisons are between paths rather than between independent simulations.

The bottleneck itself

Five generations at ten individuals, then ten generations back at 300. The control path stays at 300 for all fifteen generations, so the difference between them is the bottleneck and nothing else. Forty replicate paths of each.

run_path <- function(pop, sizes, u) {
  for (n_next in sizes) pop <- step_gen(pop, n_next, u)
  pop
}

reps <- 40
bott_sizes <- c(rep(10, 5), rep(300, 10))
ctrl_sizes <- rep(300, 15)

set.seed(468)
bott <- t(replicate(reps, stats_pop(run_path(pop, bott_sizes, u))))
ctrl <- t(replicate(reps, stats_pop(run_path(pop, ctrl_sizes, u))))

summ <- function(x, nm) {
  data.frame(path = nm, richness = round(mean(x[, "richness"]), 3),
             m_ratio = round(mean(x[, "m_ratio"]), 4),
             het = round(mean(x[, "het"]), 4))
}
print(rbind(summ(ctrl, "no bottleneck"), summ(bott, "bottleneck to 10")),
      row.names = FALSE)

             path richness m_ratio    het
    no bottleneck    9.698  0.9694 0.7996
 bottleneck to 10    7.270  0.9346 0.6633

cat("relative loss against the control: richness",
    round(1 - mean(bott[, "richness"]) / mean(ctrl[, "richness"]), 4),
    " M ratio", round(1 - mean(bott[, "m_ratio"]) / mean(ctrl[, "m_ratio"]), 4),
    " heterozygosity",
    round(1 - mean(bott[, "het"]) / mean(ctrl[, "het"]), 4), "\n")

relative loss against the control: richness 0.2504  M ratio 0.0359  heterozygosity 0.1705

cat("expected heterozygosity loss from 5 generations at N = 10:",
    round(1 - (1 - 1 / (2 * 10))^5, 4), "\n")

expected heterozygosity loss from 5 generations at N = 10: 0.2262

Allelic richness falls by 25.0 percent, heterozygosity by 17.1 percent, and the M ratio by only 3.6 percent. The ordering of the first two is the classic result and the reason allelic measures are recommended for detecting recent bottlenecks: richness counts alleles, most alleles are rare, and rare alleles are the ones a small sample of parents fails to pass on. Heterozygosity weights alleles by frequency, so losing a rare allele barely moves it.

The heterozygosity number also has a prediction attached. Five generations at ten individuals should cost 1 - (1 - 1/(2N))^5 of the gene diversity, which is 22.6 percent, and the measurement is 17.1 percent. The gap is the ten recovery generations: at a mutation rate this high, a large population starts rebuilding diversity immediately, and part of the loss is already repaid by the time anybody samples.

set.seed(471)
track <- function(pop, sizes, u) {
  out <- rbind(stats_pop(pop))
  for (n_next in sizes) {
    pop <- step_gen(pop, n_next, u)
    out <- rbind(out, stats_pop(pop))
  }
  as.data.frame(out)
}
traj_b <- track(pop, bott_sizes, u)
traj_c <- track(pop, ctrl_sizes, u)
print(round(traj_b[c(1, 3, 6, 11, 16), c("richness", "m_ratio", "het")], 4))

   richness m_ratio    het
1    9.8333  0.9740 0.8013
3    5.3833  0.7634 0.7110
6    4.3167  0.7004 0.6309
11   7.0833  0.9084 0.6487
16   7.5333  0.9286 0.6593

Generation by generation the picture is sharper than the endpoint comparison suggests. Richness falls from 9.8333 to 5.3833 within two generations and bottoms at 4.3167 at the end of the bottleneck, because 10 diploids can carry at most 20 copies at a locus. Ten generations later it has climbed back to 7.0833 and after fifteen to 7.5333, all of that from new mutation rather than recovery of anything lost. Gene diversity goes 0.8013, then 0.6309 at the bottom, then 0.6487 and 0.6593: it loses less and regains less. The M ratio is the most striking: 0.9740 before, 0.7004 at the bottom, and 0.9286 by generation 15. Measured during the crash it is unmistakable; measured fifteen generations later it is nearly gone.

mk <- function(df, nm) {
  data.frame(generation = 0:(nrow(df) - 1),
             richness = df$richness / df$richness[1],
             het = df$het / df$het[1], path = nm)
}
traj_df <- rbind(mk(traj_b, "bottleneck"), mk(traj_c, "control"))
long_df <- rbind(
  data.frame(traj_df[c("generation", "path")], value = traj_df$richness,
             measure = "allelic richness"),
  data.frame(traj_df[c("generation", "path")], value = traj_df$het,
             measure = "gene diversity"))

ggplot(long_df, aes(generation, value, colour = measure, linetype = path)) +
  annotate("rect", xmin = 0, xmax = 5, ymin = -Inf, ymax = Inf,
           fill = te_pal$sage, alpha = 0.18) +
  geom_line(linewidth = 0.9) +
  scale_colour_manual(values = c(`allelic richness` = te_pal$clay,
                                 `gene diversity` = te_pal$forest),
                      name = NULL) +
  scale_linetype_manual(values = c(bottleneck = "solid", control = "dashed"),
                        name = NULL) +
  labs(title = "A deep crash and a partial, mutation-fed recovery",
       subtitle = "shaded: five generations at N = 10, then ten generations back at N = 300",
       x = "generation", y = "value relative to the start") +
  theme_te()

Two solid lines dropping sharply inside a shaded band at the left, richness falling further than gene diversity, then richness climbing part of the way back while gene diversity stays flat; the dashed control lines stay near one. — Figure 1: Allelic richness and gene diversity relative to their starting values along the bottleneck path (five generations at ten individuals, shaded) and the control path.

Which alleles actually went

Averages hide the mechanism, so count alleles by frequency class before and after.

spectrum_pop <- function(pop) {
  cls <- c(0, 0.01, 0.05, 0.10, 1)
  out <- sapply(seq_len(ncol(pop)), function(j) {
    fr <- table(pop[, j]) / nrow(pop)
    table(cut(fr, cls, include.lowest = TRUE))
  })
  rowMeans(out)
}

set.seed(470)
pop_ctrl <- run_path(pop, ctrl_sizes, u)
pop_bott <- run_path(pop, bott_sizes, u)
spec <- rbind(control = spectrum_pop(pop_ctrl),
              bottleneck = spectrum_pop(pop_bott))
print(round(spec, 3))

           [0,0.01] (0.01,0.05] (0.05,0.1] (0.1,1]
control       1.367       2.467        2.1   3.800
bottleneck    1.800       1.667        0.8   2.967

The damage is concentrated in the middle of the spectrum. Alleles between five and ten percent frequency drop from 2.1 per locus to 0.8, and the one to five percent class from 2.467 to 1.667, while alleles above ten percent are largely intact (3.8 against 2.967). The rarest class actually rises, from 1.367 to 1.8, because ten generations at 300 individuals with this mutation rate produce fresh singletons. So a bottlenecked population is not simply a smaller version of the original: it has lost its intermediate-frequency alleles and is refilling the bottom of the spectrum with new mutations.

spec_df <- data.frame(
  class = rep(colnames(spec), 2),
  alleles = c(spec["control", ], spec["bottleneck", ]),
  path = rep(c("control", "bottleneck"), each = ncol(spec)))
spec_df$class <- factor(spec_df$class, levels = colnames(spec))

ggplot(spec_df, aes(class, alleles, fill = path)) +
  geom_col(position = position_dodge(width = 0.75), width = 0.65,
           colour = "white") +
  scale_fill_manual(values = c(control = te_pal$sage,
                               bottleneck = te_pal$clay), name = NULL) +
  labs(title = "The middle of the frequency spectrum is what disappears",
       subtitle = "alleles per locus by frequency class, 60 loci",
       x = "allele frequency class", y = "alleles per locus") +
  theme_te()

Paired bars in four frequency classes; the bottlenecked bars are shorter in the middle two classes and in the commonest class, and taller in the rarest class. — Figure 2: Alleles per locus by frequency class in the control and bottlenecked populations.

How deep does a bottleneck have to be?

The three measures also differ in how sharply they respond to depth. Repeat the five-generation bottleneck at four sizes.

set.seed(469)
sweep_b <- c(5, 10, 25, 50)
sw <- t(sapply(sweep_b, function(nb) {
  b <- t(replicate(25, stats_pop(run_path(pop, c(rep(nb, 5), rep(300, 10)),
                                          u))))
  c(bottleneck_N = nb, richness = mean(b[, "richness"]),
    m_ratio = mean(b[, "m_ratio"]), het = mean(b[, "het"]))
}))
print(round(sw, 4))

     bottleneck_N richness m_ratio    het
[1,]            5   6.0740  0.9354 0.5510
[2,]           10   7.3200  0.9327 0.6664
[3,]           25   8.5247  0.9505 0.7469
[4,]           50   9.0093  0.9563 0.7720

z_shift <- function(stat) {
  (mean(ctrl[, stat]) - mean(bott[, stat])) / sd(ctrl[, stat])
}
cat("shift in control standard deviations: richness",
    round(z_shift("richness"), 1), " M ratio", round(z_shift("m_ratio"), 1),
    " heterozygosity", round(z_shift("het"), 1), "\n")

shift in control standard deviations: richness 22.3  M ratio 6.6  heterozygosity 39.6

cat("control standard deviation: richness",
    round(sd(ctrl[, "richness"]), 4),
    " M ratio", round(sd(ctrl[, "m_ratio"]), 4),
    " heterozygosity", round(sd(ctrl[, "het"]), 4), "\n")

control standard deviation: richness 0.1088  M ratio 0.0053  heterozygosity 0.0034

A bottleneck to 50 costs 0.7 alleles per locus and 0.03 of gene diversity: real, and easy to miss with a small marker panel. A bottleneck to five costs 3.6 alleles and 0.25 of diversity. The M ratio moves least across the whole range (0.9354 to 0.9563 against a control 0.9694), which is the honest limit of that statistic in this simulation: five generations is not long enough to erode the allele size range, and the ratio was designed for bottlenecks lasting tens of generations.

The last two lines correct a folk claim worth stating carefully. Measured in units of the control variation, the heterozygosity shift is the largest of the three (39.6 standard deviations against 22.3 for richness), because with 60 loci the mean gene diversity is estimated very precisely: its control standard deviation is 0.0034 while richness carries 0.1088. So allelic richness loses a bigger fraction, and heterozygosity is measured more sharply. Which one detects a given bottleneck first depends on the panel size, and the sensible answer is to report both.

What this does and does not license

Three limits belong with any bottleneck analysis.

A single sample cannot separate a recent bottleneck from a long history of small size. Both leave low diversity; only the relationship between allelic and heterozygosity measures, or a second sample from the past, distinguishes them. That relationship is exactly what the heterozygosity-excess and M-ratio tests use, and it fades within a few dozen generations.

Recovery of census size is not recovery of diversity. In the simulation the population is back at 300 within ten generations and its diversity is not, because the only source of new alleles is mutation. On the timescales conservation cares about, alleles lost in a bottleneck are lost, and the only fast route back is immigration.

Diversity is not fitness. The simulated loci are neutral markers and say nothing directly about inbreeding depression, adaptive potential or population growth. They measure how much drift the population has experienced, which is a useful proxy precisely because it is a proxy.

Where to go next

The cluster now has an effective size, two ways to estimate it, and a demonstration of what a crash does to the markers. The last tutorial turns all of this into checks: which Ne an estimate refers to, over what window, what migration does to it, and how a fluctuating history quietly breaks the assumptions the estimators are built on.

References

Nei M, Maruyama T, Chakraborty R 1975 Evolution 29(1):1-10 (10.1111/j.1558-5646.1975.tb00807.x)

Cornuet JM, Luikart G 1996 Genetics 144(4):2001-2014 (10.1093/genetics/144.4.2001)

Garza JC, Williamson EG 2001 Molecular Ecology 10(2):305-318 (10.1046/j.1365-294X.2001.01190.x)

Allendorf FW 1986 Zoo Biology 5(2):181-190 (10.1002/zoo.1430050212)

--- title: "Bottlenecks and genetic diversity" description: "Push a simulated population through a bottleneck in R and watch which measures of diversity notice: allelic richness, heterozygosity and the M ratio disagree." date: "2026-08-12 11:00" categories: [R, population genetics, ecology tutorial, ggplot2] image: thumbnail.png image-alt: "Allelic richness collapsing during five bottleneck generations while heterozygosity declines gently, both relative to their starting values." --- A population that spends five generations at ten individuals and then recovers carries a genetic record of those five generations for a long time. Which record, exactly, depends on what you measure. Allelic richness, heterozygosity and the ratio of allele number to allele size range respond at different speeds and to different features of the crash, and a monitoring programme that tracks only one of them can miss a bottleneck or misdate it. This tutorial builds a stepwise-mutation population, sends replicate copies through bottlenecks of different depths, and measures all three. ## A population with alleles worth losing Microsatellite-style markers mutate by adding or removing a repeat unit, so alleles are integers and the mutation model is a random walk. The simulation below tracks allele sizes for every gene copy at 60 loci in 300 diploids, with the mutation rate set so that the scaled diversity `theta = 4Nu` is 10, which is the range typical of real microsatellite panels. ```{r setup} #| message: false #| warning: false library(ggplot2) te_pal <- list(forest = "#275139", green = "#2f8f63", sage = "#93a87f", clay = "#b5534e", gold = "#cda23f", line = "#dad9ca", ink = "#16241d") theme_te <- function() { theme_minimal(base_size = 12) + theme(panel.grid.minor = element_blank(), panel.grid.major = element_line(colour = "#e7e6dc"), plot.background = element_rect(fill = "white", colour = NA), panel.background = element_rect(fill = "white", colour = NA), plot.title = element_text(face = "bold", colour = te_pal$ink), axis.title = element_text(colour = "#2c3a31")) } ``` ```{r engine} step_gen <- function(pop, n_next, u) { n_now <- nrow(pop) / 2 n_loci <- ncol(pop) pick <- function(par) { rows <- 2 * par - 1 + matrix(sample(0:1, n_next * n_loci, replace = TRUE), n_next, n_loci) matrix(pop[cbind(as.vector(rows), rep(seq_len(n_loci), each = n_next))], n_next, n_loci) } h1 <- pick(sample(n_now, n_next, replace = TRUE)) h2 <- pick(sample(n_now, n_next, replace = TRUE)) out <- matrix(0L, 2 * n_next, n_loci) out[seq(1, 2 * n_next, by = 2), ] <- h1 out[seq(2, 2 * n_next, by = 2), ] <- h2 mut <- runif(length(out)) < u out[mut] <- out[mut] + sample(c(-1L, 1L), sum(mut), replace = TRUE) out } stats_pop <- function(pop) { k <- apply(pop, 2, function(col) length(unique(col))) spread <- apply(pop, 2, function(col) diff(range(col)) + 1) het <- apply(pop, 2, function(col) 1 - sum((table(col) / length(col))^2)) c(richness = mean(k), spread = mean(spread), m_ratio = mean(k / spread), het = mean(het)) } set.seed(467) n_anc <- 300 n_loci <- 60 theta <- 10 u <- theta / (4 * n_anc) pop <- matrix(as.integer(round(rnorm(2 * n_anc * n_loci, 0, sqrt(theta / 2)))), nrow = 2 * n_anc) for (g in seq_len(150)) pop <- step_gen(pop, n_anc, u) print(round(stats_pop(pop), 4)) ``` The ancestral population carries 9.83 alleles per locus spread over 10.1 repeat units, so the ratio of allele number to size range is 0.974, and gene diversity is 0.8013. That is the reference state. Every scenario below starts from this same population, which is why the comparisons are between paths rather than between independent simulations. ## The bottleneck itself Five generations at ten individuals, then ten generations back at 300. The control path stays at 300 for all fifteen generations, so the difference between them is the bottleneck and nothing else. Forty replicate paths of each. ```{r paths} run_path <- function(pop, sizes, u) { for (n_next in sizes) pop <- step_gen(pop, n_next, u) pop } reps <- 40 bott_sizes <- c(rep(10, 5), rep(300, 10)) ctrl_sizes <- rep(300, 15) set.seed(468) bott <- t(replicate(reps, stats_pop(run_path(pop, bott_sizes, u)))) ctrl <- t(replicate(reps, stats_pop(run_path(pop, ctrl_sizes, u)))) summ <- function(x, nm) { data.frame(path = nm, richness = round(mean(x[, "richness"]), 3), m_ratio = round(mean(x[, "m_ratio"]), 4), het = round(mean(x[, "het"]), 4)) } print(rbind(summ(ctrl, "no bottleneck"), summ(bott, "bottleneck to 10")), row.names = FALSE) cat("relative loss against the control: richness", round(1 - mean(bott[, "richness"]) / mean(ctrl[, "richness"]), 4), " M ratio", round(1 - mean(bott[, "m_ratio"]) / mean(ctrl[, "m_ratio"]), 4), " heterozygosity", round(1 - mean(bott[, "het"]) / mean(ctrl[, "het"]), 4), "\n") cat("expected heterozygosity loss from 5 generations at N = 10:", round(1 - (1 - 1 / (2 * 10))^5, 4), "\n") ``` Allelic richness falls by 25.0 percent, heterozygosity by 17.1 percent, and the M ratio by only 3.6 percent. The ordering of the first two is the classic result and the reason allelic measures are recommended for detecting recent bottlenecks: richness counts alleles, most alleles are rare, and rare alleles are the ones a small sample of parents fails to pass on. Heterozygosity weights alleles by frequency, so losing a rare allele barely moves it. The heterozygosity number also has a prediction attached. Five generations at ten individuals should cost `1 - (1 - 1/(2N))^5` of the gene diversity, which is 22.6 percent, and the measurement is 17.1 percent. The gap is the ten recovery generations: at a mutation rate this high, a large population starts rebuilding diversity immediately, and part of the loss is already repaid by the time anybody samples. ```{r trajectory} set.seed(471) track <- function(pop, sizes, u) { out <- rbind(stats_pop(pop)) for (n_next in sizes) { pop <- step_gen(pop, n_next, u) out <- rbind(out, stats_pop(pop)) } as.data.frame(out) } traj_b <- track(pop, bott_sizes, u) traj_c <- track(pop, ctrl_sizes, u) print(round(traj_b[c(1, 3, 6, 11, 16), c("richness", "m_ratio", "het")], 4)) ``` Generation by generation the picture is sharper than the endpoint comparison suggests. Richness falls from 9.8333 to 5.3833 within two generations and bottoms at 4.3167 at the end of the bottleneck, because 10 diploids can carry at most 20 copies at a locus. Ten generations later it has climbed back to 7.0833 and after fifteen to 7.5333, all of that from new mutation rather than recovery of anything lost. Gene diversity goes 0.8013, then 0.6309 at the bottom, then 0.6487 and 0.6593: it loses less and regains less. The M ratio is the most striking: 0.9740 before, 0.7004 at the bottom, and 0.9286 by generation 15. Measured during the crash it is unmistakable; measured fifteen generations later it is nearly gone. ```{r fig-trajectory} #| fig-width: 7 #| fig-height: 4.2 #| fig-cap: "Allelic richness and gene diversity relative to their starting values along the bottleneck path (five generations at ten individuals, shaded) and the control path." #| fig-alt: "Two solid lines dropping sharply inside a shaded band at the left, richness falling further than gene diversity, then richness climbing part of the way back while gene diversity stays flat; the dashed control lines stay near one." mk <- function(df, nm) { data.frame(generation = 0:(nrow(df) - 1), richness = df$richness / df$richness[1], het = df$het / df$het[1], path = nm) } traj_df <- rbind(mk(traj_b, "bottleneck"), mk(traj_c, "control")) long_df <- rbind( data.frame(traj_df[c("generation", "path")], value = traj_df$richness, measure = "allelic richness"), data.frame(traj_df[c("generation", "path")], value = traj_df$het, measure = "gene diversity")) ggplot(long_df, aes(generation, value, colour = measure, linetype = path)) + annotate("rect", xmin = 0, xmax = 5, ymin = -Inf, ymax = Inf, fill = te_pal$sage, alpha = 0.18) + geom_line(linewidth = 0.9) + scale_colour_manual(values = c(`allelic richness` = te_pal$clay, `gene diversity` = te_pal$forest), name = NULL) + scale_linetype_manual(values = c(bottleneck = "solid", control = "dashed"), name = NULL) + labs(title = "A deep crash and a partial, mutation-fed recovery", subtitle = "shaded: five generations at N = 10, then ten generations back at N = 300", x = "generation", y = "value relative to the start") + theme_te() ``` ## Which alleles actually went Averages hide the mechanism, so count alleles by frequency class before and after. ```{r spectrum} spectrum_pop <- function(pop) { cls <- c(0, 0.01, 0.05, 0.10, 1) out <- sapply(seq_len(ncol(pop)), function(j) { fr <- table(pop[, j]) / nrow(pop) table(cut(fr, cls, include.lowest = TRUE)) }) rowMeans(out) } set.seed(470) pop_ctrl <- run_path(pop, ctrl_sizes, u) pop_bott <- run_path(pop, bott_sizes, u) spec <- rbind(control = spectrum_pop(pop_ctrl), bottleneck = spectrum_pop(pop_bott)) print(round(spec, 3)) ``` The damage is concentrated in the middle of the spectrum. Alleles between five and ten percent frequency drop from 2.1 per locus to 0.8, and the one to five percent class from 2.467 to 1.667, while alleles above ten percent are largely intact (3.8 against 2.967). The rarest class actually rises, from 1.367 to 1.8, because ten generations at 300 individuals with this mutation rate produce fresh singletons. So a bottlenecked population is not simply a smaller version of the original: it has lost its intermediate-frequency alleles and is refilling the bottom of the spectrum with new mutations. ```{r fig-spectrum} #| fig-width: 7 #| fig-height: 4 #| fig-cap: "Alleles per locus by frequency class in the control and bottlenecked populations." #| fig-alt: "Paired bars in four frequency classes; the bottlenecked bars are shorter in the middle two classes and in the commonest class, and taller in the rarest class." spec_df <- data.frame( class = rep(colnames(spec), 2), alleles = c(spec["control", ], spec["bottleneck", ]), path = rep(c("control", "bottleneck"), each = ncol(spec))) spec_df$class <- factor(spec_df$class, levels = colnames(spec)) ggplot(spec_df, aes(class, alleles, fill = path)) + geom_col(position = position_dodge(width = 0.75), width = 0.65, colour = "white") + scale_fill_manual(values = c(control = te_pal$sage, bottleneck = te_pal$clay), name = NULL) + labs(title = "The middle of the frequency spectrum is what disappears", subtitle = "alleles per locus by frequency class, 60 loci", x = "allele frequency class", y = "alleles per locus") + theme_te() ``` ## How deep does a bottleneck have to be? The three measures also differ in how sharply they respond to depth. Repeat the five-generation bottleneck at four sizes. ```{r severity} set.seed(469) sweep_b <- c(5, 10, 25, 50) sw <- t(sapply(sweep_b, function(nb) { b <- t(replicate(25, stats_pop(run_path(pop, c(rep(nb, 5), rep(300, 10)), u)))) c(bottleneck_N = nb, richness = mean(b[, "richness"]), m_ratio = mean(b[, "m_ratio"]), het = mean(b[, "het"])) })) print(round(sw, 4)) z_shift <- function(stat) { (mean(ctrl[, stat]) - mean(bott[, stat])) / sd(ctrl[, stat]) } cat("shift in control standard deviations: richness", round(z_shift("richness"), 1), " M ratio", round(z_shift("m_ratio"), 1), " heterozygosity", round(z_shift("het"), 1), "\n") cat("control standard deviation: richness", round(sd(ctrl[, "richness"]), 4), " M ratio", round(sd(ctrl[, "m_ratio"]), 4), " heterozygosity", round(sd(ctrl[, "het"]), 4), "\n") ``` A bottleneck to 50 costs 0.7 alleles per locus and 0.03 of gene diversity: real, and easy to miss with a small marker panel. A bottleneck to five costs 3.6 alleles and 0.25 of diversity. The M ratio moves least across the whole range (0.9354 to 0.9563 against a control 0.9694), which is the honest limit of that statistic in this simulation: five generations is not long enough to erode the allele size range, and the ratio was designed for bottlenecks lasting tens of generations. The last two lines correct a folk claim worth stating carefully. Measured in units of the control variation, the heterozygosity shift is the largest of the three (39.6 standard deviations against 22.3 for richness), because with 60 loci the mean gene diversity is estimated very precisely: its control standard deviation is 0.0034 while richness carries 0.1088. So allelic richness loses a bigger fraction, and heterozygosity is measured more sharply. Which one detects a given bottleneck first depends on the panel size, and the sensible answer is to report both. ## What this does and does not license Three limits belong with any bottleneck analysis. A single sample cannot separate a recent bottleneck from a long history of small size. Both leave low diversity; only the relationship between allelic and heterozygosity measures, or a second sample from the past, distinguishes them. That relationship is exactly what the heterozygosity-excess and M-ratio tests use, and it fades within a few dozen generations. Recovery of census size is not recovery of diversity. In the simulation the population is back at 300 within ten generations and its diversity is not, because the only source of new alleles is mutation. On the timescales conservation cares about, alleles lost in a bottleneck are lost, and the only fast route back is immigration. Diversity is not fitness. The simulated loci are neutral markers and say nothing directly about inbreeding depression, adaptive potential or population growth. They measure how much drift the population has experienced, which is a useful proxy precisely because it is a proxy. ## Where to go next The cluster now has an effective size, two ways to estimate it, and a demonstration of what a crash does to the markers. The last tutorial turns all of this into checks: which `Ne` an estimate refers to, over what window, what migration does to it, and how a fluctuating history quietly breaks the assumptions the estimators are built on. ## References Nei M, Maruyama T, Chakraborty R 1975 Evolution 29(1):1-10 (10.1111/j.1558-5646.1975.tb00807.x) Cornuet JM, Luikart G 1996 Genetics 144(4):2001-2014 (10.1093/genetics/144.4.2001) Garza JC, Williamson EG 2001 Molecular Ecology 10(2):305-318 (10.1046/j.1365-294X.2001.01190.x) Allendorf FW 1986 Zoo Biology 5(2):181-190 (10.1002/zoo.1430050212) ## Related tutorials - [Effective population size](../effective-population-size/) - [Estimating Ne from genetic data](../estimating-effective-size-from-data/) - [Checking an effective size estimate](../checking-an-effective-size-estimate/) - [Hardy-Weinberg expectations](../hardy-weinberg-and-genotype-frequencies/)