Quantifying the relationship between genetic diversity and population size suggests natural selection cannot explain Lewontin’s Paradox

  1. Vince Buffalo  Is a corresponding author
  1. Institute for Ecology and Evolution, University of Oregon, United States


Neutral theory predicts that genetic diversity increases with population size, yet observed levels of diversity across metazoans vary only two orders of magnitude while population sizes vary over several. This unexpectedly narrow range of diversity is known as Lewontin’s Paradox of Variation (1974). While some have suggested selection constrains diversity, tests of this hypothesis seem to fall short. Here, I revisit Lewontin’s Paradox to assess whether current models of linked selection are capable of reducing diversity to this extent. To quantify the discrepancy between pairwise diversity and census population sizes across species, I combine previously-published estimates of pairwise diversity from 172 metazoan taxa with newly derived estimates of census sizes. Using phylogenetic comparative methods, I show this relationship is significant accounting for phylogeny, but with high phylogenetic signal and evidence that some lineages experience shifts in the evolutionary rate of diversity deep in the past. Additionally, I find a negative relationship between recombination map length and census size, suggesting abundant species have less recombination and experience greater reductions in diversity due to linked selection. However, I show that even assuming strong and abundant selection, models of linked selection are unlikely to explain the observed relationship between diversity and census sizes across species.


A longstanding mystery in evolutionary genetics is that the observed levels of genetic variation across sexual species span an unexpectedly narrow range. Under neutral theory, the average number of nucleotide differences between sequences (pairwise diversity, π) is determined by the balance of new mutations and their loss by genetic drift (Kimura and Crow, 1964; Malécot, 1948; Wright, 1931). In particular, expected pairwise diversity at neutral sites in a panmictic population of Nc diploids is π4Ncμ, where μ is the per basepair per generation mutation rate. Given that metazoan germline mutation rates only differ 10-fold (10−8–10−9, Kondrashov and Kondrashov, 2010; Lynch, 2010), and census sizes vary over several orders of magnitude, under neutral theory one would expect that pairwise diversity also vary over several orders of magnitude. However, early allozyme surveys revealed that diversity levels across a wide range of species varied just an order of magnitude (Lewontin, 1974, p. 208); this is known as Lewontin’s ‘‘Paradox of Variation’. With modern sequencing-based estimates of π across taxa ranging over only three orders of magnitude (0.01–10%, Leffler et al., 2012), Lewontin’s paradox remains unresolved through the genomics era.

Early on, explanations for Lewontin’s Paradox have been framed in terms of the neutralist–selectionist controversy (Lewontin, 1974; Kimura, 1984; Gillespie, 1991; Gillespie, 2001). The neutralist view is that beneficial alleles are sufficiently rare and deleterious alleles are removed sufficiently quickly, that levels of genetic diversity are shaped predominantly by genetic drift and mutation (Kimura, 1984). Specifically, non-selective processes decouple the effective population size implied by observed levels of diversity π^, N~e=π^/4μ, from the census size, Nc. By contrast, the selectionist view is that direct selection and the indirect effects of selection on linked neutral diversity suppress diversity levels across taxa, specifically because the impact of linked selection is greater in large populations. Undoubtedly, these opposing views represent a false dichotomy, as population genomic studies have uncovered evidence for the substantial impact of both demographic history (e.g. Zhao et al., 2013; Palkopoulou et al., 2015) and linked selection on genome-wide diversity (e.g. Elyashiv et al., 2016; Begun and Aquadro, 1992; Aguade et al., 1989; McVicker et al., 2009).

Possible resolutions of Lewontin’s Paradox

A resolution of Lewontin’s Paradox would involve a mechanistic description and quantification of the evolutionary processes that prevent diversity from scaling with census sizes across species. This would necessarily connect to the broader literature on the empirical relationship between diversity and population size (Frankham, 1996; Nei and Graur, 1984; Soulé, 1976; Leroy et al., 2021), and the ecological and life history correlates of genetic diversity (Nevo, 1978; Powell, 1975; Nevo et al., 1984). Three categories of processes stand out as potentially capable of decoupling census sizes from diversity: non-equilibrium demography, variance and skew in reproductive success, and selective processes.

It has long been appreciated that effective population sizes are typically less than census population sizes, tracing back to early debates between R.A. Fisher and Sewall Wright (Fisher and Ford, 1947; Wright, 1948). Possible causes of this divergence between effective and census population sizes include demographic history (e.g. population bottlenecks), extinction and recolonization dynamics, or the breeding structure of populations (e.g. the variance in reproductive success and population substructure). Early explanations for Lewontin’s Paradox suggested bottlenecks during the last glacial maximum severely reduced population sizes (Kimura, 1984; Ohta and Kimura, 1973; Nei and Graur, 1984), and emphasized that large populations recover to equilibrium diversity levels more slowly (Nei and Graur, 1984, Kimura, 1984 p. 203–204). Another explanation is that cosmopolitan species repeatedly endure extinction and recolonization events, which reduces effective population size (Maruyama and Kimura, 1980; Slatkin, 1977).

While intermittent demographic events like bottlenecks and recent expansions have long-term impacts on diversity (since mutation-drift equilibrium is reached on the order of size of the population), characteristics of the breeding structure such as high variance (Vw) or skew in reproductive success continuously suppress diversity below the levels predicted by the census size (Wright, 1938). For example, in many marine animals, females are highly fecund, and dispersing larvae face extremely low survivorship, leading to high variance in reproductive success (Waples et al., 2018; Waples et al., 2013; Hedgecock and Pudovkin, 2011; Hauser and Carvalho, 2008). Such ‘‘sweepstakes’ reproductive systems can lead to remarkably small ratios of effective to census population size (e.g. Ne/Nc can range from 10–6–10–2), since Ne/N1/Vw(Hedgecock, 1994; Wright, 1938; Nunney, 1993), and require multiple-merger coalescent processes to describe their genealogies (Eldon and Wakeley, 2006). Overall, these reproductive systems diminish the diversity in some species, but seem unlikely to explain Lewontin’s Paradox broadly across metazoans.

Alternatively, selective processes, and in particular the indirect effects of selection on linked neutral variation, could potentially explain the observed narrow range of diversity. The earliest mathematical model of hitchhiking was proffered as a explanation of Lewontin’s Paradox (Smith and Haigh, 1974). Since, linked selection has been shown to impact diversity levels in a variety of species, as evidenced by the correlation between recombination and diversity (Aguade et al., 1989; Begun and Aquadro, 1992; Cutter and Payseur, 2003; Stephan and Langley, 1998; Cai et al., 2009). Theoretic work to explain this pattern has considered the impact of a steady influx of beneficial mutations (recurrent hitchhiking; Stephan et al., 1992; Stephan, 1995), and purifying selection against deleterious mutations (background selection, BGS; Charlesworth et al., 1993; Nordborg et al., 1996; Hudson and Kaplan, 1994). Indeed, empirical work indicates background selection diminishes diversity around genic regions in a variety of species (McVicker et al., 2009; Hernandez et al., 2011; Charlesworth, 1996), and now efforts have shifted towards teasing apart the effects of positive and negative selection on genomic diversity (Elyashiv et al., 2016).

A class of models that are of particular interest in the context of Lewontin’s Paradox are recurrent hitchhiking models that decouple diversity from the census population size. These models predict diversity levels when strongly selected beneficial mutations regularly enter and sweep through the population, trapping lineages and forcing them to coalesce (Kaplan et al., 1989; Gillespie, 2000). In general, decoupling occurs under these hitchhiking models when the rate of coalescence due to selection is much greater than the rate of neutral coalescence (e.g. Coop and Ralph, 2012, Equation 22). In contrast, under other linked selection models, the resulting effective population size is proportional to population size; these models cannot decouple diversity, all else equal. For example, models of background selection and polygenic fitness variation predict diversity is proportional to population size, mediated by the total recombination map length and the deleterious mutation rate or fitness variation (Charlesworth et al., 1993; Nicolaisen and Desai, 2012; Nordborg et al., 1996; Robertson, 1961; Santiago and Caballero, 1995).

Recent approaches towards resolving Lewontin’s Paradox

Recently, Corbett-Detig et al., 2015 used population genomic data to estimate the reduction in diversity due to background selection and hitchhiking across 40 species, and showed that the impact of selection increases with two proxies of census population size, species range and the inverse of body size. Based on this evidence, they argued that selection could explain Lewontin’s Paradox; however, in a re-analysis, Coop, 2016 demonstrated that the observed magnitude of these reductions is insufficient to explain the orders-of-magnitude shortfall between observed and expected levels of diversity across species. Other recent work has found that life history characteristics related to parental investment, such as propagule size, are good predictors diversity in animals (Romiguier et al., 2014; Chen et al., 2017). Nevertheless, while these diversity correlates are important clues, they do not propose a mechanism by which these traits act to constrain diversity within a few orders of magnitude.

Here, I revisit Lewontin’s Paradox by integrating several data sets in order to compare the observed relationship between diversity and census size with the predicted relationship under different selection models. Prior surveys of genetic diversity either lacked census population size estimates, used allozyme-based measures of heterozygosity, or included fewer species. To address these shortcomings, I first estimate census sizes by combining predictions of population density based on body size with ranges estimated from geographic occurrence data. Using these estimates, I quantify the relationship between census size and previously-published genomic diversity estimates across 172 metazoan taxa within nine phyla, thus characterizing the relationship between π and Nc that underlies Lewontin’s Paradox.

Past work looking at the relationship between π and Nc has been unable to fully account for phylogenetic non-independence across taxa (Felsenstein, 1985). To address this, I use phylogenetic comparative methods (PCMs) with a synthetic time-calibrated phylogeny to account for shared phylogenetic history. Moreover, it is disputed whether considering phylogenetic non-independence is necessary in population genetics, given that coalescent times within species are much less than divergence times (Whitney and Garland, 2010; Lynch, 2011). Using PCMs, I address this by estimating the degree of phylogenetic signal in the diversity census size relationship, and investigating how these traits evolve along the phylogeny.

Finally, I explore whether the predicted reductions of diversity under background selection and recurrent hitchhiking are sufficiently strong to resolve Lewontin’s Paradox. I do so using selection parameters from Drosophila melanogaster, a species known to be strongly affected by linked selection. Given the effects of linked selection are mediated by recombination map length, I also investigate how map lengths vary with census population size using data from a previously-published survey (Stapley et al., 2017). I find map lengths are typically shorter in large–census-size species, increasing the effects of linked selection in these species, which could further decouple diversity from census size. Still, I find the combined impact of these modes of linked selection fall short in explaining Lewontin’s Paradox, and discuss future avenues through which the Paradox of Variation could be fully resolved.


Estimates of census population size

An impediment in resolving Lewontin’s Paradox is characterizing the relationship between diversity and census population sizes. This is difficult because census population sizes are unavailable for many taxa, especially for extremely abundant, cosmopolitan species that define the upper limit of ranges. Previous work has surveyed the literature for census size estimates (Nei and Graur, 1984; Soulé, 1976; Frankham, 1996), or used range, body size, or qualitative categories as proxies for census size (Corbett-Detig et al., 2015; Leffler et al., 2012). To quantify the relationship between genomic estimates of diversity and census population sizes, I first approximate census population sizes for 172 metazoan taxa (Figure 1). I estimate population densities based on an empirical linear relationship between body sizes and density that holds across metazoans (see Figure 1—figure supplement 1; Damuth, 1981; Damuth, 1987). Then, from geographic occurrence data, I estimate range sizes. Finally, I estimate population size as the product of these predicted densities and range estimates (see Materials and methods: Macroecological Estimates of Population Size). Note that the relationship between population density and body size is driven by energy budgets, and thus reflects macroecological equilibria (Damuth, 1987). Consequently, population sizes are underestimated for taxa like humans and their domesticated species, and overestimated for species with anthropogenically reduced densities or fragmented ranges. For example, the population size of Lynx lynx is likely around 50,000 (IUCN, 2020) which is around two orders of magnitude smaller than my estimate. Additionally, the range size estimates do not consider whether an area has unsuitable habitat, and thus may be overestimated for species with particular niches or patchy habitats. While my approach produces approximate and sometimes crude estimates, it has the advantage that it can be efficiently calculated for numerous taxa, which is sufficient to estimate the magnitude of Lewontin’s Paradox (see Population Size Validation for more on validation based on biomass and other approaches).

Figure 1 with 5 supplements see all
The distribution of approximate census population sizes estimated by this study.

Some phyla containing few species were excluded for clarity.

Figure 1—source data 1

The population size estimates for 172 metazoan taxa.


Characterizing the Diversity–Census-size Relationship

To determine which ecological or evolutionary processes could decouple diversity from census population size, we first need to quantify this relationship across a wide variety of taxa. Previous work has found there is a significant relationship between heterozygosity and the logarithm of population size or range size, but these studies relied on heterozygosity measured from allozyme data (Soulé, 1976; Frankham, 1996; Nei and Graur, 1984). I confirm these findings using pairwise diversity estimates from genomic sequence data and the estimated census sizes (Figure 2). The pairwise diversity estimates are from three sources: Leffler et al., 2012, Corbett-Detig et al., 2015, and Romiguier et al., 2014, and are predominantly from either synonymous or non-coding DNA (see Methods and Materials: 4.1 Diversity and Map Length Data). Overall, an ordinary least squares (OLS) relationship on a log-log scale fits the data well (Figure 2, gray dashed line). The OLS slope estimate is significant and implies a 13% percent increase in differences per basepair for every order of magnitude census size grows (95% confidence interval [12%, 14%], adjusted R2=0.26; see also the OLS fit per-phyla, Figure 2—figure supplement 2).

Figure 2 with 4 supplements see all
A visualization of Lewontin’s Paradox of Variation.

Pairwise diversity (data from Leffler et al., 2012, Corbett-Detig et al., 2015, and Romiguier et al., 2014), which varies over three orders of magnitude, shows a weak relationship with approximate population size, which varies over 12 orders of magnitude. The shaded curve shows the range of expected neutral diversity if Ne were to equal Nc under the four-alleles model, log10(π)=log10(θ)log10(1+4θ/3) where θ=4Ncμ, for two mutation rates, μ=10-8 and μ=10-9, and the light gray dashed line represents the maximum pairwise diversity under the four alleles model. The dark gray dashed line is the OLS regression fit, and the blue dashed line is the regression fit using a phylogenetic mixed-effects model. Points are colored by phylum. The species Equus ferus przewalskii (Nc103 and π=3.6×10-3) was an outlier and excluded from this figure for visual clarity.

Figure 2—source data 1

The diversity and population size dataset for 172 metazoan taxa.


Notably, this relationship has few outliers and is relatively homoscedastic. This is in part because of the log-log scale, in contrast to previous work (Nei and Graur, 1984; Soulé, 1976); see Figure 2—figure supplement 1 for a version on a log-linear scale. However, it is noteworthy that few taxa have diversity estimates below 10−3.5 differences per basepair. Those that do, lynx (Lynx Lynx), wolverine (Gulo gulo), and Massasauga rattlesnake (Sistrurus catenatus) face habitat loss and declining population sizes. These three species are all in the IUCN Red List, but are listed as least concern (though their presence in the Red List indicates they are of conservation interest). In Appendix D, Appendix D Diversity and IUCN Red List Status, I explore the relationships between IUCN Red List status, diversity, and population size.

Phylogenetic non-independence and the population size diversity relationship

One limitation of using ordinary least squares is that shared phylogenetic history can create correlation structure in the residuals, which violates an assumption of the regression model and can lead to bias (Felsenstein, 1985; Revell, 2010). To address this shortcoming, I fit the diversity–census-size relationship using a phylogenetic mixed-effects model, investigated whether there is a signal of phylogenetic non-independence, estimated the continuous trait values on the phylogeny, and explored how diversity and population size evolve. Prior population genetic comparative studies have lacked time-calibrated phylogenies and assumed unit branch lengths (Whitney and Garland, 2010), a shortcoming that has drawn criticism (Lynch, 2011). I use a synthetic time-calibrated phylogeny created from the DateLife project (O’Meara et al., 2020) to account for shared phylogenetic history (see Materials and methods: Phylogenetic Comparative Methods).

Using a phylogenetic mixed-effects model (Lynch, 1991; Hadfield and Nakagawa, 2010; de Villemereuil and Nakagawa, 2014) implemented in Stan (Carpenter et al., 2017; Stan Development Team, 2020), I estimated the linear relationship between diversity and population size (on a log-log scale) accounting for phylogeny, for the 166 taxa without missing data and present in the synthetic chronogram. This type of model is needed because closely-related species may differ from the average trend between Nc and π in similar ways due to shared phylogenetic history, similar life history traits, etc., and thus do not represent independent observations as is assumed by the standard regression model. This is a form of phylogenetic pseudoreplication, and can be accounted for with a phylogenetic mixed-effects model. The phylogenetic mixed-effects model does not assume that there is phylogenetic structure in either Nc or π (which itself is not a violation of the standard regression model, Revell, 2010 and Uyeda et al., 2018), but rather accounts for phylogenetic correlation structure in the residuals if any is present. Importantly, phylogenetic mixed-effects models simultaneously estimate the degree of phylogenetic structure in the residuals while fitting the relationship between Nc and π. If the residuals are distributed independently, the estimated relationship would be similar to that found by ordinary least squares, and the estimated phylogenetic signal would be zero. Overall, this approach is conservative, making no assumptions about the source of the phylogenetic signal while accounting for violations of the regression model due to dependence among the residuals if present (see Revell, 2010 for a discussion of this).

As with the linear regression, I find this relationship is positive and significant (95% credible interval 0.03, 0.11), though somewhat attenuated compared to the OLS estimates (Figure 3B). Since the population size estimates are based on range and body mass, they are essentially a composite trait; fitting phylogenetic mixed-effects models separately on body mass and range indicates these have significant positive and negative effects, respectively (Figure 2—figure supplement 3; see also Figure 2—figure supplement 4 for the relationship between diversity and the range categories of Leffler et al., 2012).

Figure 3 with 3 supplements see all
Phylogenetic comparative models of diversity and population size.

(A) The ancestral continuous trait estimates for the population size and diversity (differences per bp, log scaled) across the phylogeny of 166 taxa. The phyla of the tips are indicated by the color bar in the center. (B) The posterior distributions of the intercept, slope, and phylogenetic signal (λ, de Villemereuil and Nakagawa, 2014) of the phylogenetic mixed-effects model of diversity and population size (log scaled). Also shown are the 90% credible interval (light blue shading), posterior mean (blue line), OLS estimate (gray solid line), and bootstrap OLS confidence intervals (light gray shading). (C) The node-height tests of diversity, population size, and the two components of the population size estimates, body mass, and range (all traits on log scale before contrast was calculated). Each point shows the standardized phylogenetic independent contrast and branching time for a pair of lineages. Red lines are robust regression estimates (and are only shown for statistically significant relationships at the α=0.05 level). Note that some outlier pairs with very high phylogenetic independent contrasts were excluded (in all cases, these outliers were in the genus Drosophila).

Since the phylogenetic mixed-effects model simultaneously estimates the variance of the phylogenetic effect (σp2) and the residual variance (σr2), these can be used to estimate the phylogenetic signal, λ=σp2/(σp2+σr2) (Lynch, 1991; de Villemereuil and Nakagawa, 2014; see Freckleton et al., 2002 for a comparison to Pagel’s λ). When residuals are free of correlations due to shared phylogenetic history, then λ=0 and all the variance could be explained by evolution or noise on the tips. In the relationship between population size and diversity, the posterior mean of λ=0.67 (90% credible interval [0.58,0.75]) indicates a majority of the variance perhaps might be due to shared phylogenetic history (Figure 3B).

This high degree of phylogenetic signal substantiates Gillespie's concern (Gillespie, 1991) that the π–Nc relationship may be driven by chordate-arthropod differences. A visual inspection of the estimated ancestral continuous values for diversity and population size on the phylogeny indicates the high phylogenetic signal seems to be driven in part by chordates having low diversity and small population sizes compared to non-chordates (Figure 3A). This problem resembles Felsenstein’s worst-case scenario (Felsenstein, 1985; Uyeda et al., 2018), where a singular event on a lineage separating two clades generates a spurious association between two traits.

To investigate whether clade-level differences dominated the relationship between diversity and population size, I fit phylogenetic mixed-effects models to phyla-level subsets of the data for clades with sufficient sample sizes (see Methods: 4.4 Phylogenetic Comparative Methods). This analysis shows a significant positive relationship between diversity and population size in arthropods, and positive weak relationships in molluscs and chordates (Figure 3—figure supplement 1). Each of the 90% credible intervals for slope overlap, suggesting the relationship between π and Nc is similar across these clades.

Additionally, I have explored the rate of trait change through time using node-height tests (Freckleton and Harvey, 2006). Node-height tests regress the absolute values of the standardized contrasts between lineages against the branching time (since present) of these lineages. Under Brownian Motion (BM), standardized contrasts are estimates of the rate of character evolution (Felsenstein, 1985); if a trait evolves under constant rate BM, this relationship should be flat. For both diversity and population size, node-height tests indicate a significant increase in the rate of evolution towards the present (robust regression p-values 0.023 and 0.00018 respectively; Figure 3C). Considering the constituents of the population size estimate, range and body mass, separately, the rate of evolution of range but not body mass shows a significant increase (p-value 1.03 × 10−7) towards the present.

Interestingly, the diversity node-height test reveals two rate shifts at deeper splits (Figure 3C, top left) around 570 Mya. These nodes represent the branches between tunicates and vertebrates in chordates, and cephalopods and pleistomollusca (bivalves and gastropods) in molluscs. While the cephalopod-pleistomollusca split outlier may be an artifact of having a single cephalopod (Sepia officinalis) in the phylogeny, the tunicate-vertebrate split outlier is driven by the low diversity of vertebrates and the previously-documented exceptionally high diversity of tunicates (sea squirts; Nydam and Harrison, 2010; Small et al., 2007). This deep node representing a rate shift in diversity could reflect a change in either effective population size or mutation rate, and there is some evidence of both in this genus Ciona (Small et al., 2007; Tsagkogeorga et al., 2012). Neither of these deep rate shifts in diversity is mirrored in the population size node-height test (Figure 3C, top right). Rather, it appears a trait impacting diversity but not census size (e.g. mutation rate or offspring distributions) has experienced a shift on the lineage separating tunicates and vertebrates. At nearly 600 Mya, these deep nodes illustrate that expected effective population sizes (and thus coalescence times) can share phylogenetic history, due to phylogenetic inertia in some combination of population size, reproductive system, and mutation rates.

Finally, an important caveat is the increase in rate towards the tips could be caused by measurement noise, or possibly uncertainty or bias in the divergence time estimates deep in the tree. Inspecting the lineage pairs that lead to this increase in rate towards the tips indicates these represent plausible rate shifts, e.g. between cosmopolitan and endemic sister species like Drosophila simulans and Drosophila sechellia; however, ruling out measurement noise entirely as an explanation would involve modeling the uncertainty of diversity and population size estimates.

Assessing the impact of linked selection on diversity across taxa

The above analyses reemphasize the drastic shortfall of diversity levels as compared to census sizes. Linked selection has been proposed as the mechanism that acts to reduce diversity levels from what we would expect given census sizes (Smith and Haigh, 1974; Gillespie, 2000; Corbett-Detig et al., 2015). Here, I test this hypothesis by estimating the scale of diversity reductions expected under background selection and recurrent hitchhiking, and comparing these to the observed relationship between π and Nc.

I quantify the effect of linked selection on diversity as the ratio of observed diversity (π) to the estimated diversity in the absence of linked selection (π0), R=π/π0. Here, π0 would reflect only demographic history and non-heritable variation in reproductive success. There are two difficulties in evaluating whether linked selection could resolve Lewontin’s Paradox. The first difficulty is that π0 is unobserved. Previous work has estimated π0 using methods that exploit the spatial heterogeneity in recombination and functional density across the genome to fit linked selection models that incorporate both hitchhiking and background selection (Elyashiv et al., 2016; Corbett-Detig et al., 2015). The second difficulty is understanding how R varies across taxa, since we lack estimates of critical model parameters for most species. Still, I can address a key question: if diversity levels were determined by census sizes (π0=4Ncμ), would the combined effects of background selection and recurrent hitchhiking be sufficient to reduce diversity to observed levels? Furthermore, does the relationship between census size and predicted diversity under linked selection across species, πBGS+HH=Rπ0, match the observed relationship in Figure 2?

Since we lack estimates of selection parameters across species, I parameterize the hitchhiking and BGS models using estimates from Drosophila melanogaster, a species known to be strongly affected by linked selection (Sella et al., 2009). Under a generalized model of hitchhiking and background selection (Elyashiv et al., 2016; Coop and Ralph, 2012) and assuming Ne=Nc, the expected diversity is

(1) πBGS+HHθ1/B(U,L)+2NcS(γ,J,L)

where θ=4Ncμ, B(U,L) is the effect of background selection, and S(γ,J,L) is the rate of coalescence caused by sweeps (Elyashiv et al., 2016, Equation 1, Coop and Ralph, 2012, Equation 20). Under background selection models with recombination, the reduction is B(U,L)=exp(U/L) where U is the per diploid genome per generation deleterious mutation rate, and L is the recombination map length in Morgans (Hudson and Kaplan, 1994; Nordborg et al., 1996). This BGS model is similar to models of effective population size under polygenic fitness variation, and can account for other modes of linked selection (Robertson, 1961; Santiago and Caballero, 1995; Santiago and Caballero, 1998, see Appendix 2, Background Selection and Polygenic Fitness Models). The coalescence rate due to sweeps is S(γ,J,L)=γLJ, where γ is the number of adaptive substitutions per generation, and J is the probability a lineage is trapped by sweeps as they occur across the genome (J2,2 in Equation 15 of Coop and Ralph, 2012).

Parameterizing the model this way, I then set the key parameters that determine the impact of recurrent hitchhiking and background selection (γ, J, and U) to strong selection values estimated for Drosophila melanogaster by Elyashiv et al., 2016. My estimate of the adaptive substitutions per generation (γDmel2.3×10-3) based Elyashiv et al. implies a rate of sweeps per basepair of νBP,Dmel2.34×10-11, which is close to other estimates from D. melanogaster (see Figure 4—figure supplement 5A). The rate of deleterious mutations per diploid genome, per generation is parameterized using the estimate from Elyashiv et al., UDmel=1.6, which is slightly greater than previous estimates based on Bateman-Mukai approaches (Mukai, 1985; Mukai, 1988; Charlesworth, 1987). Finally, the probability that a lineage is trapped in a sweep, JDmel4.5×10-4, is calculated from the estimated genome-wide average coalescence rate due to sweeps from Elyashiv et al. (see Figure 4—figure supplement 5B and Materials and methods: Predicted Reductions in Diversity for more details on parameter estimates). Using these parameters, I then explore how the predicted range of diversity levels varies across species with recombination map length (L) and census population size (Nc).

Previous work has found that the impact of linked selection increases with Nc (Corbett-Detig et al., 2015; see Figure 4—figure supplement 4A), and it is often thought that this is driven by higher rates of adaptive substitutions in larger populations (Ohta, 1992), despite equivocal evidence (Galtier, 2016). However, there is another mechanism by which species with larger population sizes might experience a greater impact of linked selection: recombination map length, L, is known to correlate with body mass (Burt and Bell, 1987) and thus varies inversely with population size. As this is a critical parameter that determines the genome-wide impact of both hitchhiking and background selection, I examine the relationship between recombination map length (L) and census population size (Nc) across taxa, using available estimates of map lengths across species (Stapley et al., 2017; Corbett-Detig et al., 2015). I find a significant non-linear relationship using phylogenetic mixed-effects models (Figure 4A; see Methods and materials: 4.4 Phylogenetic Comparative Methods). There is also a correlation between map length and genome size (Figure 4—figure supplement 2) and genome size and population size (Figure 4—figure supplement 1). These findings are consistent with the hypothesis that non-adaptive processes increase genome size in small-Ne species (Lynch and Conery, 2003) which in turn could increase map lengths, as well as the hypothesis that map lengths are adaptively longer to more efficiently select against deleterious alleles in smaller populations (Roze, 2021). Overall, the negative relationship between map length and census size indicates linked selection is expected to be stronger in species with short map lengths, which are high-Nc species.

Figure 4 with 5 supplements see all
Predicting the impact of linked selection on diversity.

(A) The observed relationship between recombination map length (L) and census size (Nc) across 136 species with complete data and known phylogeny. Triangle points indicate six social taxa excluded from the model fitting since these have adaptively higher recombination map lengths (Wilfert et al., 2007). The dark gray line is the estimated relationship under a phylogenetic mixed-effects model, and the gray interval is the 95% posterior average. (B) Points indicate the observed π–Nc relationship across taxa shown in Figure 2, and the blue ribbon is the range of predicted diversity were Ne=Nc for μ=10-810-9, and after accounting for the expected reduction in diversity due to background selection and recurrent hitchhiking under Drosophila melanogaster parameters. In both plots, point color indicates phylum.

Figure 4—source data 1

The map length, population size, and linked selection estimates for 136 metazoan taxa.


Then, I predict the expected diversity (πBGS+HH) under background selection and hitchhiking, assuming Ne=Nc and that all species had the rate of sweeps and strength of BGS as D. melanogaster. Since neutral mutation rates μ are unknown and vary across species, I calculate the range of predicted πBGS+HH estimates for µ = 10−9–10−8 (using the four-alleles model, Tajima, 1996), and compare this to the observed relationship between π and Nc in Figure 4B. Under these parameters and assumptions, linked selection begins to appreciably constrain diversity for Nc107, since S(γDmel,JDmel,L)108107 and linked selection dominates drift when S(γ,J,L)>1/2N. Overall, this reveals two problems for the hypothesis that linked selection could solve Lewontin’s Paradox. First, low to mid-Nc species (census sizes between 104–107) have sufficiently long map lengths that their diversity levels are only moderately reduced by linked selection, leading to a wide gap between predicted and observed diversity levels. For this not to be the case, the rate of adaptive mutations or the deleterious mutation rate would need to be orders of magnitude higher for species within this range than in Drosophila melanogaster, which is incompatible with the rate of adaptive protein substitutions across species (Galtier, 2016) and overall mutation rates (Lynch, 2010). Furthermore, linked selection has been quantified in humans, which fall in this census size range, and has been found to be relatively weak (McVicker et al., 2009; Hernandez et al., 2011; Hellmann et al., 2008; Cai et al., 2009; Boyko et al., 2008). Second, while hitchhiking and BGS can reduce predicted diversity levels for high-Nc species (Nc>1012) to observed levels, this would imply available estimates of π0 are underestimated by several orders of magnitude in Drosophila (Figure 4—figure supplement 4B). The high reductions in π predicted here (compared to those of Elyashiv et al., 2016) are a result of using Nc, rather than Ne=π0/4μ in the denominator of Equation (1), which leads to a very high rate of sweeps in the population. I do not consider selective interference, though the saturation of adaptive substitutions per Morgan would only act to limit the reduction in diversity (Weissman and Barton, 2012), and thus these results are conservative.

Finally, the poor fit between observed and predicted levels of diversity across species is not remedied by stronger selection parameters. In Figure 4—figure supplement 3B, I increase both selection parameters U and γ ten-fold each, and find the same qualitative pattern: on a log-log scale the relationship between Nc and π is linear, while the predicted diversity under linked selection is non-linear with Nc. Under this ten-fold higher selection regime, there is more overlap between observed and predicted levels of diversity, but diversity is severely under-predicted for high-Nc species. Additionally, this would imply that selection in low-to-mid-Nc species is ten-folder higher than estimated in Drosophila melanogaster, which is implausible. Overall, this suggests that present models of linked selection, even with very strong selection across species, are qualitatively incapable of matching the observed relationship between Nc and π and thus cannot explain Lewontin’s Paradox.


Nearly fifty years after Lewontin’s description of the Paradox of Variation, how evolutionary, life history, and ecological processes interact to constrain diversity across taxa to a narrow range remains a mystery. I revisit Lewontin’s Paradox by first characterizing the relationship between genomic estimates of pairwise diversity and approximate census population size across 172 metazoan species. Previous surveys have used allozyme-based estimates, fewer taxa, or proxies of population size. My estimates of census population sizes are rough approximates, since they use body size to predict density. An improved estimate might account for vagility (as Soulé, 1976 did), though this is harder to do systematically across many taxa. Future work might also use other ecological information, such as total biomass, or species distribution modeling to improve census size estimates (Bar-On et al., 2018; Mora et al., 2011). Still, it seems more accurate estimates would be unlikely to change the qualitative findings here, which resemble those of early surveys (Nei and Graur, 1984; Soulé, 1976).

One limitation of this study is that diversity estimates are collated from a variety of sources rather than estimated with a single bioinformatic pipeline. This leads to technical noise across diversity estimates; perhaps the relationship between π and Nc found here could be tighter with a standardized bioinformatic pipeline. In addition, there might be systematic bioinformatic sources of bias: for example high-diversity sequences may fail to align to the reference genome and end up unaccounted for, leading to a downward bias. Alternatively, a high-diversity sequences might map to the reference genome, but adjacent mis-matching SNPs might be mistaken for a short insertion or deletion. While these issues might affect estimates in high-diversity species, it is unlikely to change the qualitative relationship between π–Nc.

Macroevolution and Across-Taxa population genomics

Lewontin’s Paradox arises from a comparison of diversity across species, yet it has been disputed whether such comparisons require phylogenetic comparative methods. Extending previous work that has accounted for phylogeny in particular clades (Leffler et al., 2012), or using taxonomical-level averages (Romiguier et al., 2014), I show that the positive relationship between diversity and census size is significant using a mixed-effects model with a time-calibrated phylogeny. Additionally, I find a high degree of phylogenetic signal, evidence of deep shifts in the rate of evolution of genetic diversity, and that arthropods and chordates form clusters. Overall, this suggests that previous concerns about phylogenetic non-independence in comparative population genetic studies were warranted (Gillespie, 1991; Whitney and Garland, 2010). Notably, Lynch, 2011 has argued that PCMs for pairwise diversity are unnecessary, since mutation rate evolution is fast and thus free of phylogenetic inertia, sampling variance should exceed the variance due to phylogenetic shared history, and coalescence times are much shorter than divergence times. Since my findings suggest PCMs are necessary in some cases, it is worthwhile to address these points.

First, Lynch has correctly pointed out that while coalescence times are much less than divergence times and should be free of phylogenetic shared history, the factors that determine coalescence times (e.g. mutation rates and effective population size) may not be (Lynch, 2011). In other words, coalescence times are free from phylogenetic shared history were we to condition on these causal factors that could be affected by shared phylogenetic history. My estimates of phylogenetic signal in the residuals, by contrast, are not conditioned on these factors. Importantly, even "correcting for" phylogeny implicitly favors certain causal interpretations over others (Westoby et al., 1995; Uyeda et al., 2018). Future work could try to untangle what causal factors determine coalescence times across species, as well as how these factors evolve across macroevolutionary timescales. Second, it is a misconception that a fast rate of trait evolution necessarily reduces phylogenetic signal (Revell et al., 2008), and that if either or both variables in a regression are free of phylogenetic signal, PCMs are unnecessary (Revell, 2010; Uyeda et al., 2018). The evidence of high phylogenetic signal found in this study suggests PCMs are necessary when fitting the relationship between Nc and π in order to account for correlated residuals among closely-related species, and to avoid spurious results from phylogenetic pseudoreplication.

Finally, beyond just accounting for phylogenetic non-independence, macroevolution and phylogenetic comparative methods are a promising way to approach across-species population genomic questions. For example, one could imagine that diversification processes could contribute to Lewontin’s Paradox. If large-Nc species were to have a rate of speciation that is greater than the rate at which mutation and drift reach equilibrium (which is indeed slower for large Nc species), this could act to decouple diversity from census population size. That is to say, even if the rate of random demographic bottlenecks were constant across taxa, lineage-specific diversification processes could lead certain clades to be systematically further from demographic equilibrium, and thus have lower diversity than expected for their census population size.

How could selection still explain Lewontin’s Paradox?

Even assuming selection parameters estimated from Drosophila melanogaster, where the effects of linked selection are thought to be especially strong, the predicted patterns of diversity under linked selection poorly fit observed patterns of diversity across species. My results support the analysis by Coop, 2016 showing that levels of π0 estimated by Corbett-Detig et al., 2015 are not decoupled from genome-wide average π, as would occur if linked selection were to explain Lewontin’s Paradox. Additionally, my analysis goes a step further, showing that current linked selection models under a wide range of selection parameters are incapable of explaining the observed relationship between census size and diversity. This is in part because mid-Nc species have sufficiently long recombination map lengths to diminish the effects of even strong selection. Overall, while this suggests hitchhiking and background selection seem unlikely to explain patterns of diversity across taxa, there are three major potential limitations of my approach that need further evaluation.

First, I approximate the reduction in diversity using homogeneous background selection and recurrent hitchhiking models (Kaplan et al., 1989; Hudson and Kaplan, 1995; Coop and Ralph, 2012), when in reality, there is genome-wide heterogeneity in functional density, recombination rates, and the adaptive substitutions across species. Each of these factors mediate how strongly linked selection impacts diversity across the genome. Despite these model simplifications, the predicted reduction in diversity in Drosophila melanogaster is 85% (when using Ne, not Nc), which is reasonably close to the estimated 77% from the more realistic model of Elyashiv et al. that accounts for the actual position of substitutions, annotation features, and recombination rate heterogeneity (though it should be noted that these both use the same parameter estimates). Furthermore, even though my model fails to capture the heterogeneity of functionality density and recombination rate in real genomes, it is still conservative, likely overestimating the effects of linked selection to see if it could be capable of decoupling diversity from census size and explain Lewontin’s Paradox. This is in part because the strong selection parameter estimates from Drosophila melanogaster used, but also because I assume that the effective population size is equal to the census size. Even then, this decoupling only occurs in very high–census-size species, and implies that the diversity in the absence of linked selection, π0, is currently underestimated by several orders of magnitude. Moreover, the study of Corbett-Detig et al., 2015 did consider recombination rate and functional density heterogeneity in estimating the reduction due to linked selection across species, yet their predicted reductions are orders of magnitude weaker than those considered here by assuming that Ne=Nc (Figure 4—figure supplement 4B). Overall, given the effects estimated under more realistic inference models are still orders of magnitude weaker than those used in this study, current models of linked selection seem fundamentally unable to fit the diversity–census-size relationship.

Second, my model here only considers hard sweeps, and ignores the contribution of soft sweeps (e.g. from standing variation or recurrent mutations; Hermisson and Pennings, 2005; Pennings and Hermisson, 2006), partial sweeps (e.g those that do not reach fixation), and the interaction of sweeps and spatial processes. While future work exploring these alternative types of sweeps is needed, the predicted reductions in diversity found here under the simplified sweep model are likely relatively robust to these other modes of sweeps for a few reasons. First, the shape of the diversity–recombination curve is equivalent under models of partial sweeps and hard sweeps, though these imply different rates of sweeps (Coop and Ralph, 2012). Second, in the limit where most fitness variation is due to weak soft sweeps from standing variation scattered across the genome (i.e. due to polygenic fitness variation), levels of diversity are well approximated by quantitative genetic linked selection models (Robertson, 1961; Santiago and Caballero, 1995). The reduction in diversity under these models is nearly identical to that under background selection models, in part because deleterious alleles at mutation-selection balance constitute a considerable component of fitness variation (see Appendix Section B; Charlesworth and Hughes, 2000; Charlesworth, 2015). Third, the parameters from Elyashiv et al., 2016 could reflect a mixture of types of sweeps (Elyashiv et al., 2016 p. 14 and p. 19 of their Supplementary Online Materials). Finally, I also disregarded the interaction of sweeps and spatial processes. For populations spread over wide ranges, limited dispersal slows the spread of sweeps, allowing for new beneficial alleles to arise, spread, and compete against other segregating beneficial variants (Ralph and Coop, 2015; Ralph and Coop, 2010). Through limited dispersal should act to ‘‘soften sweeps’ and not impact my findings for the reasons described above, future work could investigate how these processes impact diversity in ways not captured by hard sweep models.

Third, other selective processes, such as fluctuating selection or hard selective events (i.e. selection resulting in a reduction in the population size), could reduce diversity in ways not captured by the background selection and hitchhiking models. Since frequency-independent fluctuating selection reduces diversity under most conditions (Novak and Barton, 2017), this could lead seasonality and other sources of temporal heterogeneity to reduce diversity in large-Nc species with short generation times more than longer-lived species with smaller population sizes. Future work could consider the impact of fluctuating selection on diversity under simple models (Barton, 2000) if estimates of key parameters governing the rate of such fluctuations were known across taxa. Additionally, another mode of selection that could severely reduce diversity across taxa, yet remains unaccounted for in this study, is periodic hard selective events. These selective events could occur regularly in a species’ history yet be indistinguishable from demographic bottlenecks with just population genomic data.

Spatial and demographic processes

One limitation of this study is the inability to quantify the impact of spatial and demographic population genetic processes on the relationship between diversity and census population sizes across taxa. The genomic diversity estimates collated in this study unfortunately lack details about the sampling process and spatial data, which can have a profound impact on population genomic summary statistics (Battey et al., 2020). These issues could systematically bias species-wide diversity estimates; for example, if diversity estimates from a cosmopolitan species were primarily from a single region or subpopulation, diversity would be an underestimate relative to the entire population. However, biased spatial sampling alone seems incapable of explaining the π-Nc divergence in high-Nc taxa. In the extreme scenario in which only one subpopulation was sampled, FST would need to be close to one for population subdivision alone to sufficiently reduce the total population heterozygosity to explain the orders-of-magnitude shortfall between predicted and observed diversity levels. This can be seen by rearranging the expression for FST as HS=(1-FST)HT, where HS and HT are the subpopulation and total population heterozygosities; if HT=4Ncμ, then only FST1 can reduce HS several orders of magnitude. Yet, across-taxa surveys indicate that FST is almost never this high within species (Roux et al., 2016). Future work could quantify the extent to which more realistic spatial processes contribute to Lewontin’s Paradox. For example, high-Nc taxa usually experience range expansions, with repeated founder effects and local extinction/recolonization dynamics that depress diversity (Slatkin, 1977). In particular, with the appropriate data, one could estimate the empirical relationship between dispersal distance, range size, and coalescent effective population size across taxa.

In this study, I have focused entirely on assessing the role of linked selection, rather than demography, in reducing diversity across taxa. In contrast to demographic models, models of linked selection have comparatively fewer parameters and more readily permit rough estimates of diversity reductions across taxa. Given that I find that models of linked selection are incapable of explaining the observed relationship between Nc and π, this supports the hypothesis the diversity across species are shaped primarily by past demographic fluctuations. Still, a full resolution of Lewontin’s Paradox would require understanding how the demographic processes across taxa with incredibly heterogeneous ecologies and life histories transform Nc into Ne. With population genomic data becoming available for more species, this could involve systematically inferring the demographic histories of tens of species and looking for correlations in the frequency and size of bottlenecks with Nc across species.

Measures of effective population size, Timescales, and Lewontin’s Paradox

Lewontin’s Paradox describes the extent to which the effective population sizes implied by diversity, N~e, diverge from census population sizes. However, there are a variety other effective population size estimators calculable from different data and summary statistics (Wang et al., 2016; Caballero, 1994; Galtier and Rousselle, 2020). These include estimators based on the site frequency spectrum, observed decay in linkage disequilibrium, or temporal estimators that use the variance in allele frequency change through time. These various estimators capture different summaries of effective population size on shorter timescales than coalescent-based estimators (see Wang, 2005 for a review), and thus could be used to tease apart processes that impact the Ne-Nc relationship in the more recent past.

Temporal Ne estimators already play an important role in understanding another summary of the Ne-Nc relationship: the ratio Ne/Nc, which is an important quantity in conservation genetics (Frankham, 1995; Mace and Lande, 1991) and in understanding evolution in highly fecund marine species. Surveys of the short-term Ne/Nc relationship across taxa indicate mean Ne/Nc is on order of 0.1 (Frankham, 1995; Palstra and Ruzzante, 2008; Palstra and Fraser, 2012), though the uncertainty in these estimates is high, and some species with sweepstakes reproduction systems like Pacific Oyster (Crassostrea gigas) can have Ne/Nc106 (Hedgecock, 1994). Estimates of the Ne/Nc ratio may be an important, yet under appreciated piece of solving Lewontin’s Paradox. For example, if Ne is estimated from the allele frequency change across a single generation (i.e. Waples, 1989), Ne/Nc constrains estimates of the variance in reproductive success (Wright, 1938; Nunney, 1993; Nunney, 1996). This implies that apart from species with sweepstakes reproductive systems, the variance in reproductive success each generation (whether heritable or non-heritable) is likely insufficient to significantly contribute to constraining N~e for most taxa. Still, further work is needed to characterize (1) how Ne/Nc varies with Nc across taxa (though see Palstra and Fraser, 2012, Figure 2), and (2) the variance of Ne/Nc over longer time spans (i.e. how periodic sweepstakes reproductive events act to constrain Ne). Overall, characterizing how Ne/Nc varies across taxa and correlates with ecology and life history traits could provide clues into the mechanisms that leads propagule size and survivorship curves to be predictive of diversity levels across taxa (Romiguier et al., 2014; Hallatschek, 2018; Barry et al., 2020).

Finally, short-term temporal Ne estimators may play an important role in resolving Lewontin’s Paradox. These estimators, along with short-term estimates of the impact of linked selection (Buffalo and Coop, 2019; Buffalo and Coop, 2020), can inform us how much diversity is depressed by selection on shorter timescales, free from the rare strong selective events or severe bottlenecks that impact pairwise diversity. It could be that in any one generation, selection contributes more to the variance of allele frequency changes than drift, yet across-taxa patterns in diversity are better explained processes acting sporadically on longer timescales, such as colonization, founder effects, and bottlenecks. Thus, the pairwise diversity may not give us the best picture of the generation to generation evolutionary processes acting in a population to change allele frequencies. Furthermore, certain observed adaptations occur at a pace that is inexplicable given small effective population sizes implied by diversity, and are only possible if short-term effective population sizes are orders of magnitude larger (Karasov et al., 2010; Barton, 2010).


In Building a Science of Population Biology (Lewontin et al., 2004), Lewontin laments the difficulty of uniting population genetics and population ecology into a cohesive discipline of population biology. Lewontin’s Paradox of Variation remains a major unsolved problem at the nexus of these two different disciplines: we fail to understand the processes that connect a central parameter of population ecology, census size, to a central parameter of population genetics, effective population size across species. Given that selection seems to fall short in resolving Lewontin’s Paradox, a full resolution will require a mechanistic understanding the ecological, life history, and macroevolutionary processes that connect Nc to Ne across taxa. While I have focused exclusively on metazoan taxa since their population densities are more readily approximated from body mass, a full resolution must also include plant species (with the added difficulties of variation in selfing rates, different dispersal strategies, pollination, etc.).

Looking at Lewontin’s Paradox through an macroecological and macroevolutionary lens begets interesting questions outside of the traditional realm of population genetics. Here, I have found that diversity and Nc have a consistent relationship without many outliers, despite the wildly disparate ecologies, life histories, and evolutionary histories of the taxa included. Furthermore, taxa with very large census sizes have surprisingly low diversity. Is this explained by macroevolutionary processes, such as different rates of speciation for large-Nc taxa? Or, are the levels of diversity we observe today an artifact of our timing relative to the last glacial maximum, or the last major extinction? Did large-Nc prehistoric animal populations living in other geological eras have higher levels of diversity than our present taxa? Or, does ecological competition occur on shorter timescales such that strong population size contractions transpire and depress diversity, even if a species is undisturbed by climatic shifts or mass extinctions? Overall, patterns of diversity across taxa are determined by many overlaid evolutionary and ecological processes occurring on vastly different timescales. Lewontin’s Paradox of Variation may persist unresolved for some time because the explanation requires synthesis and model building at the intersection of all these disciplines.

Materials and methods

Diversity and map length data

Request a detailed protocol

The data used in this study are collated from a variety of previously published surveys. Of the 172 taxa with diversity estimates, 14 are from Corbett-Detig et al., 2015, 96 are from Leffler et al., 2012, and 62 are from Romiguier et al., 2014. The Corbett-Detig et al. data is estimated from four-fold degenerate sites, the Romiguier et al. data is synonymous sites, and the Leffler et al. data is estimated predominantly from silent, intronic, and non-coding sites. All types of diversity estimates from Leffler et al., 2012 were included to maximize the taxa in the study, since the variability of diversity across functional categories is much less than the diversity across taxa. Multiple diversity estimates per taxa were averaged. The total recombination map length data were from both (Stapley et al., 2017; 127 taxa), and (Corbett-Detig et al., 2015; 9 taxa). Both studies used sex-averaged recombination maps estimated with cross-based approaches; in some cases errors in the original data were found, documented, and corrected. These studies also included genome size estimates used to create Figure 4—figure supplement 2 and Figure 4—figure supplement 1.

Macroecological estimates of population size

Request a detailed protocol

A rough approximation for total population size (census size) is Nc=DR, where D is the population density in individuals per km2 and R is the range size in km2. Since population density estimates are not available for many taxa included in this study, I used the macroecological abundance-body size relationship to predict population density from body size. Since body length measurements are more readily available than body mass, I collated body length data from various sources (see https://github.com/vsbuffalo/paradox_variation; copy archived at swh:1:rev:8fa6b5834f6536319b1e5cd9722ca02d317183df, Buffalo, 2021); body lengths were averaged across sexes for sexually dimorphic species, and if only a range of lengths was available, the midpoint was used.

Then, I re-estimated the relationship between body mass and population density using the data in the appendix table of Damuth, 1987, which includes 696 taxa with body mass and population density measurements across mammals, fish, reptiles, amphibians, aquatic invertebrates, and terrestrial arthropods. Though the abundance-body size relationship can be noisy at small spatial or phylogenetic scales (Chapter 5, Gaston and Blackburn, 2008), across deeply diverged taxa such as those included in this study and Damuth, 1987, the relationship is linear and homoscedastic (see Figure 1—figure supplement 1). Using Stan (Stan Development Team, 2020), I jointly estimated the relationship between body mass from body length using the Romiguier et al., 2014 taxa, and used this relationship to predict body mass for the taxa in this study. These body masses were then used to predict population density simultaneously, using the Damuth, 1981 relationship. The code of this routine (pred_popsize_missing_centered.stan) is available in the GitHub repository (https://github.com/vsbuffalo/paradox_variation/).

To estimate range, I first downloaded occurrence records from Global Biodiversity Information Facility (Global Biodiversity Information Facility, 2020) using the rgbif R package (Chamberlain et al., 2014; Chamberlain and Boettiger, 2017). Using the occurrence locations, I inferred whether a species was marine or terrestrial, based on whether the majority of their recorded occurrences overlapped a continent using rnaturalearth and the sf packages (South, 2017; Pebesma, 2018). For each taxon, I estimated its range by finding the minimum α-shape containing these occurrences. The α parameters were set more permissive for marine species since occurrence data for marine taxa were sparser. Then, I intersected the inferred ranges for terrestrial taxa with continental polygons, so their ranges did not overrun landmasses (and likewise with marine taxa and oceans). I inspected diagnostic plots for each taxa for quality control (all of these plots are available in paradox_variation GitHub repository), and in some cases, I manually adjusted the α parameter or manually corrected the range based on known range maps (these changes are documented in the code data/species_ranges.r and data/species_range_fixes.r). The range of C. elegans was conservatively approximated as the area of the Western US and Western Europe based on the map in Frézal and Félix, 2015. Drosophila species ranges are from the Drosophila Speciation Patterns website, (Yukilevich, 2012; Yukilevich, 2017). To further validate these range estimates, I have compared these to the qualitative range descriptions Leffler et al., 2012 (Figure 1—figure supplement 4) and compared my α-shape method to a subset of taxa with range estimates from IUCN Red List (Chamberlain, 2020; IUCN, 2020; Figure 1—figure supplement 3). Each census population size is then estimated as the product of range and density.

Population size validation

Request a detailed protocol

I validated the approximate census sizes by comparing the implied biomass of these estimates to estimates of the total carbon biomass on earth by phylum (Bar-On et al., 2018). For species i with wet body mass mi and census size Ni, the implied biomass is miNi. For all species in a phylum S, this total sample biomass is bS=iSmiNi. I then compare this wet biomass to the carbon biomasses by phylum by Bar-On et al., 2018. Across animal species, the ratio of dry to wet body mass, and carbon body mass dry body mass varies little. In their study, Bar-On et al. assume wet body mass has a 70% water content, and 50% of dry body mass is carbon mass, leading to a wet body mass to carbon mass factor of 10.7/0.5=0.15. I use this factor to convert the total wet biomass to carbon biomass per phylum.

First, I compared the relative carbon biomass in this study to the relative carbon biomass on earth per phylum. This shows that this study’s sample over represents chordate biomass (by a factor of ∼3), and under represents in arthropod biomass (by a factor of 0.02) relative to the proportion of carbon biomass of these phyla on earth (see column eight of Table 1). Second, to check whether the carbon biomass per phylum in the sample was broadly consistent with the total on earth by phylum (BS for phylum S), I calculated the expected sample biomass if species were sampled randomly from the total species in a phylum, (BS×nS/TS, where nS is the total number of species in the sample in phylum S, TS is the total number of species in phylum S on earth). The fraction of total species on earth included in the sample in this study is depicted in Figure 1—figure supplement 2.

Table 1
How the total carbon biomass estimates by phylum from Bar-On et al., 2018 compare to the implied biomass estimates from this study.

All biomass estimates are carbon biomass, and the proportions are of total biomass with respect to the study. The proportion of biomass in this study compared to the Bar-On et al. estimates Bar-On et al., 2018 indicates chordates are overrepresented and arthropods are underrepresented in the present study; the factor that each phylum is overrepresented is given in the eighth column. Total species by phylum estimates are from Reaka-Kudla et al., 1996; Nicol, 1969; Zhang, 2013; Chapman, 2009. The ratio column is the ratio of total biomass implied by the Nc estimates of each species in a phylum to the actual biomass of that phylum.

Bar-On et al.Present study
phylumtotal species (T)biomass (B)prop. biomassbiomass (b)prop. biomassnum. species (n)factor overrepresentedprop. total species (f=n/T)factor (b/f⁢B)
Arthropoda1.26 × 1061.200.46352.80 × 10−40.0102680.025.41 × 10−54.31
Chordata5.41 × 1040.870.33572.67 × 10−20.9715682.891.26 × 10−324.40
Annelida1.70 × 1040.200.07721.23 × 10−50.000430.011.76 × 10−40.35
Mollusca9.54 × 1040.200.07724.56 × 10–40.0166130.211.36 × 10−416.70
Cnidaria1.60 × 1040.100.03863.07 × 10−50.001120.031.25 × 10−42.45
Nematoda2.50 × 1040.020.00774.03 × 10−60.000110.024.00 × 10−55.03

Next, I look at the ratio of sample biomass per phylum, bS, to this expected biomass per phylum (Table 1). The consistency is quite close for this rough approach and the non-random sample of taxa included in this study. The carbon biomass estimates for chordates implied by the census size estimates are ∼24-fold higher than expected, but is well within reasonable expectations given that the chordate sample includes many larger-bodied domesticated species (and is a biased sample in other ways). Similarly, the implied arthropod carbon biomass is quite close to what one would expect. Overall, these values indicate that the census size estimates here do not lead to implied biomasses per phylum that are outside the range of plausibility. For other population size consistency checks, see Appendix 3.

Phylogenetic comparative methods

Request a detailed protocol

Of the full dataset of 172 taxa with diversity and population size estimates, a synthetic calibrated phylogeny was created for 166 species that appear in phylogenies in DateLife project (O’Meara et al., 2020; Sanchez-Reyes and O’Meara, 2019). This calibrated synthetic phylogeny was then subset for the analyses based on what species had complete trait data. The diversity-population size relationship assessed by a linear phylogenetic mixed-effects model implemented in Stan (Stan Development Team, 2020), according to the methods described in de Villemereuil and Nakagawa, 2014, (see stan/phylo_mm_regression.stan in the GitHub repository). This same Stan model was used to estimate the same relationship between arthropod, chordate, and mollusc subsets of the data, though a reduced model was used for the chordate subset due to identifiability issues leading to poor MCMC convergence (Figure 3—figure supplement 1).

The relationship between recombination map length and the logarithm of population size is non-linear and heteroscedastic, and was fit using a lognormal phylogenetic mixed-effects model on the 130 species with complete data. Since social insects have longer recombination map lengths (Wilfert et al., 2007), social taxa were excluded when fitting this model. All Rhat (Vehtari et al., 2019) values were below 1.01 and the effective number of samples was over 1,000, consistent with good mixing; details about the model are available in the GitHub repository (phylo_mm_lognormal.stan). Continuous trait maps (Figure 3A, Figure 3—figure supplement 3, and Figure 3—figure supplement 2) were created using phytools (Revell, 2012). Node-height tests were implemented based on the methods in Geiger (Pennell et al., 2014; Harmon et al., 2008), and use robust regression to fit a linear relationship between phylogenetic independent contrasts and branching times.

Predicted reductions in diversity

Request a detailed protocol

The predicted reductions in diversity due to linked selection are approximated using selection and deleterious mutation parameters from Drosophila melanogaster, and the recombination map length estimates from Stapley et al., 2017 and Corbett-Detig et al., 2015. The mathematical details of the simplified sweep model are explained in the Appendix Section A. I use estimates of the number of substitutions, m, in genic regions between D. melanogaster and D. simulans from Hu et al., 2013. Following Elyashiv et al., 2016, only substitutions in UTRs and exons are included, since they found no evidence of sweeps in introns. Then, I average over annotation classes to estimate the mean proportion of substitutions that are beneficial, αDmel=0.42, which are consistent with the estimates of Elyashiv et al. and estimates from MacDonald–Kreitman test approaches (see Eyre-Walker, 2006, Table 1). Then, I use divergence time estimates between D. melanogaster and D. simulans of 4.2×106 and estimate of ten generations per year (Obbard et al., 2012), calculating there are γDmel=αm/2T=2.26×103 substitutions per generation. Given the length of the Drosophila autosomes, G, this implies that the rate of beneficial substitutions per basepair, per generation is νBP,Dmel=γDmel/G=2.34×1011. Finally, I estimate JDmel4.5×10-4 from the estimate of genome-wide average rate of sweeps from Elyashiv et al. (Supplementary Table S6) and assuming DrosophilaNe=106. These Drosophila melanogaster hitchhiking parameter estimates are close to other previously-published estimates (Figure 4—figure supplement 5). Finally, I use UDmel=1.6, from Elyashiv et al., 2016. With these parameter estimates from D. melanogaster, the recombination map lengths across species, and Equation (1), I estimate πBGS+HH (assuming Nc=Nc) across all species. This leads to a range of predicted diversity ranges across species corresponding to μ=109108; to visualize these, I take a convex hull of all diversity ranges and smooth this with R’s smooth.spline function.

Appendix 1

Simplified sweep effects model

I use a simplified model of the effects of recurrent hitchhiking and background selection (BGS) occurring uniformly along a genome. Expected diversity is given by

(2) E(π)=θθ+1/B+2NS
(3) θ1/B+2NS

(Equation 1 Elyashiv et al., 2016, Equation 4 of Kim and Stephan, 2000, and Equation 20 of Coop and Ralph, 2012). The BGS component is given by Hudson and Kaplan, 1995,

(4) B(U,L)=Neexp(-UL)

and the hitchhiking component is

(5) S=νBPrBPJ

(Coop and Ralph, 2012, Equation 20) where νBP and rBP are the substitutions and recombination per basepair respectively, J is the probability that two lineages coalesce down to one, given sweeps occur uniformly along the genome. Under this homogeneous sweep model, J is

(6) J=0Lqf(r)2𝑑r

where qf(r) is the approximate probability that a lineage is trapped by a sweep to frequency f when it is r recombination fraction away from this sweep (Coop and Ralph, 2012; Equation 15).

Since I use Drosophila melanogaster parameter estimates from Elyashiv et al., 2016, I now reconcile their model’s S term with the simple model above. They estimate S in Drosophila melanogaster using a composite likelihood model that considers hitchhiking and background selection simultaneously, using substitutions and stratifying by annotation. For a neutral position at site x, the coalescence rate due to sweeps is given by Elyashiv et al.’s Equation 3,

(7) S(x)=1TiSα(iS)ya(iS)exp(-r(x,y)τ(s,N))g(s|iS)𝑑s

where T is the length of the lineage (in generations) on which substitutions accrue, iS=1,,IS is the annotation class (e.g. exons, introns, UTRs), α(iS) is the fraction of substitutions in annotation class iS that are beneficial, a(iS) is the set of all substitutions in annotation class iS, τ(s,N) is the fixation time of a site with additive effect s, and g(s|iS) is the distribution of selection coefficients for annotation class iS.

Note, that we can recover the model of Coop and Ralph, 2012 from this expression. Suppose there is only one annotation class, and α fraction of substitutions are beneficial, and one selection coefficient s¯, (i.e. g(s)=δ0(s-s¯)), then

(8) S(x)=αTyaexp(-r(x,y)τ(s¯,N)).

Let the number of substitutions be m:=|a|, and imagine their positions are uniformly distributed on a segment of length G basepairs with the focal site is the middle at position x=0. Then, each substitution y is a random distance lyU(G/2,G/2) away from the focal site. Assuming the recombination rate is a constant rBP per basepair, and approximating the sum with an integral, we have,

(9) S=αTi=1mEli(exp(-rBPliτ(s¯,N)))
(10) =αTGi=1m0Gexp(-rBPτ(s¯,N))𝑑
(11) =αmTG0Gexp(-rBPτ(s¯,N))𝑑

Using u-substitution with r=rBP this simplifies to

(12) S=αmTGrBP0Lexp(-rτ(s¯,N))𝑑r

where L=GrBP.

To simplify this notation, note that the rate of adaptive substitutions per basepair per generation is νBP=αm/GT, so

(13) S=νBPrBP0Lexp(-rτ(s¯,N))𝑑r

This is analogous to the second term of Coop and Ralph, 2012, Equation 17, with k=i=2 and x=1 (e.g. conditioning on a sweep to fixation). Note that there appears to be a factor of two error in Elyashiv et al., 2016 compared to Coop and Ralph, 2012; here I include the factor of two. Then,

(14) S=νBPrBP0Lexp(-2rτ(s¯,N))𝑑rJ

where the integral is equal to J (J2,2 of Equation 15 in Coop and Ralph, 2012) since a simple model of qf(r)=fexp(-2rτ(s,N)) and if we condition on fixation, f=1. This expression is useful to generalize across species, since we know N and L. Additionally, we have estimates of α and m/T in Drosophila and other species. In Elyashiv et al, they consider the number of substitutions per generation in genic regions only; it should be noted that the number of coding basepairs varies little across species. For convenience, I define γ=αm/T as the number of adaptive substitutions per generation per entire genome, such that S(γ,J,L)=γLJ used in the main text. Using the estimates of m4.5×105, α0.42, and T8.4×107 from the Supplementary Material of Elyashiv et al., I arrive at γ0.00226 adaptive substitutions per generation, per genome. For a 100 megabase genome, this translates to a νBP2.34×10-11, which is close to previous estimates (Figure 4—figure supplement 5A). For J, I use an empirical estimate calculated from the genome-wide average of the rate of coalescent events due to sweeps, from Supplementary Table S6 of Elyashiv et al. (rs=2NS0.92; see Figure 4—figure supplement 5B). This implies J4.46×10-4. Alternatively, I have tried using the estimated distribution of selection coefficients from Elyashiv et al., but this led to a weaker estimate of J, since the adaptive substitutions considered tend to cluster around genic regions.

Appendix 2

Background selection and polygenic fitness models

Throughout the main text, I use recurrent hitchhiking and background selection models to estimate the reduction in diversity due to linked selection. Another class of linked selection models, which I refer to as quantitative genetic linked selection models (QGLS; Robertson, 1961; Santiago and Caballero, 1995), can also depress genome-wide diversity. Furthermore, these models may depress diversity at neutral sites unlinked to the regions containing fitness variation. While I did not explicitly incorporate these models into my estimates of the diversity reductions, their effect is implicit in background selection models because they are analytically nearly identical. Here, I briefly sketch out the connection between BGS and QGLS models.

Under the Santiago and Caballero, 1998 model, the effective population size is NeSC98=Nexp(C2/(1Z)L), where C2 is the standardized heritable fitness variation, 1-Z is the decay of genetic variance through time, and L is the recombination map length. This model can accommodate a variety of modes of selection such as selection on an infinitesimal trait (Santiago and Caballero, 1995, p. 1016), and the flux of either weakly advantageous or deleterious alleles (Santiago and Caballero, 1998, p. 2109). If the source of fitness variation is entirely the input of new deleterious mutations with heterozygous effect sh at rate U per diploid genome per generation, then under mutation-selection balance, the equilibrium relative variance in reproductive success C2=Ush (Crow and Kimura, 1970; Caballero, 2020, p. 167), and Z=1sh1/2Nc (Santiago and Caballero, 1998). Thus, if 1/2Nc<<sh<<1, then C2/(1Z)U and NeSC98Nexp(U/L), which is the BGS model used in the main text and is a result of many background selection models with similar assumptions (Hudson and Kaplan, 1994, Equation 15; Hudson and Kaplan, 1995, Equation 9; Nordborg et al., 1996, Equation 4; Barton, 1995, Equation 22b). Intuitively, the similarity of these models reflects the fact that a substantial proportion of heritable fitness variation is caused by the continual flux of deleterious alleles across the genome under mutation-selection balance (Charlesworth, 2015; Charlesworth and Hughes, 2000).

Appendix 3

Additional population size validation

In addition to the biomass-based validation described in the main text, I also conducted a few other consistency checks. First, note that the body-mass-based estimates of density for Drosophila are similar to previously used estimates in surveys of census size and diversity. Nei and Graur, 1984 suggested a maximum of 5 Drosophila per m2, including regions of the range that are not inhabitable. Across Drosophila, the body mass based estimates suggest 106.7–107.6 individuals per km2, or 4.5 – 36.3 individuals per m2, which are consistent with this previous estimate. Nei and Graur’s estimates of Drosophila pseudoobscura’s census size are four orders of magnitude smaller than mine, but their approach uses a speculated ratio of population sizes of different Drosophila species rather than range sizes (Nei and Graur, 1984, p. 81).

As another consistency check, I looked at the rank order of mammals by biomass. Whale species have the first and third highest biomass with 11.4 and 3.9 megatons of carbon biomass (for Balaenoptera bonaerensis and Eschrichtius robustus, respectively). While this seems high, a recent study shows that across whale species, pre-whaling carbon biomass was at the tens of megatons level (Pershing et al., 2010, Table 1 and Figure 1). Given that my census size estimates represent populations at a macroecological equilibrium, they would not reflect reduced density due to whaling or other anthropogenic causes. Humans had the second largest biomass, followed by wolf species (Canis lupus and C. latrans); as with whales, the population sizes for wolf species represent pre-anthropogenic densities and are overestimates compared to current population sizes, as expected.

Finally, there are other estimates of approximate population sizes for some species that I compared my estimates to. The United Nation’s FAOSTAT database estimates the total number of horses (Equus caballus) on earth as ∼60 million; the estimate in this study is close to 40 million. For other domesticated species like chicken (Gallus gallus), estimates range from 25 million to 19.6 billion (FAOSTAT statistics database, 2021; Robinson et al., 2014); the present study’s estimate lies in the middle at ∼175 million. Again, this is a known limitation of this method, as the range is estimated from occurrence data and does not consider species’ niches. This present study’s estimate of the number of king penguins (Aptenodytes patagonicus) is about 3 million; the population size was recently estimated as 2.23 million pairs (Shirihai, 2008).

Appendix 4

Diversity and IUCN Red List Status

I also investigated the relationship between species’ IUCN Red List categories (an ordinal scale of how threatened a species is) and both diversity and population size, finding that species categorized as more threatened have both smaller population sizes and reduced diversity, compared to non-threatened species (Appendix 4—figure 1) consistent with past work (Spielman et al., 2004). A linear model of diversity regressed on population size has lower AIC when the IUCN Red List categories are included, and the estimates of the effect of IUCN status are all negative on diversity, though not all are significant in part because some categories have three or fewer species (Appendix 4—table 1).

Appendix 4—figure 1
A version of Figure 2 with points colored by their IUCN Red List conservation status.

Margin boxplots show the diversity and population size ranges (thin lines) and interquartile ranges (thick lines) for each category. NA/DD indicates no IUCN Red List entry, or Red List status Data Deficient; LC is Least Concern, NT is Near Threatened, VU is Vulnerable, EN is Endangered, and CR is Critically Endangered.

Appendix 4—table 1
The regression estimates of full IUCN Red List population size model for diversity, log10(π)=β0+βLCLC+βNTNT+βVUVU+βENEN+βCRCR+βNclog10(Nc); df=165.

Using AIC to compare this full model to a reduced model of log10(π)=β0+βNclog10(Nc), AICfull=204.9, AICreduced=216.4.

Mean2.5 %97.5 %

Data availability

All primary datasets collated by this study, including new census size and range estimates, are available on Github at http://github.com/vsbuffalo/paradox_variation (copy archived at https://archive.softwareheritage.org/swh:1:rev:8fa6b5834f6536319b1e5cd9722ca02d317183df). An archived version of this repository is also available at Zenodo.

The following data sets were generated
    1. Vince B
    (2021) Zenodo
    vsbuffalo/paradox_variation: biorxiv v.1 with minor corrections.
The following previously published data sets were used
    1. Stapley J
    2. Feulner PGD
    3. Johnston SE
    4. Santure AW
    5. Smadja CM
    (2017) figshare
    Supplementary material from "Variation in recombination frequency and distribution across eukaryotes: patterns and processes".


    1. Barton NH
    (2000) Genetic hitchhiking
    Philosophical Transactions of the Royal Society of London. Series B: Biological Sciences 355:1553–1562.
  1. Software
    1. Buffalo V
    Code and Data for Why do species get a thin slice of π?, version swh:1:rev:8fa6b5834f6536319b1e5cd9722ca02d317183df https://archive.softwareheritage.org/swh:1:rev:8fa6b5834f6536319b1e5cd9722ca02d317183df
    Software Heritage.
  2. Book
    1. Caballero A
    Quantitative Genetics
    Cambridge University Press.
  3. Software
    1. Chamberlain S
    2. Ram K
    3. Barve V
    4. Mcglinn D
    rgbif: interface to the global biodiversity information facility API, version R package version 0. 7, 7
    rgbif: interface to the global biodiversity information facility API.
  4. Software
    1. Chamberlain S
    rredlist: ‘IUCN’ red list client
    rredlist: ‘IUCN’ red list client.
  5. Book
    1. Chapman AD
    Numbers of Living Species in Australia and the World
    Department of the Environment, Water, Heritage and the Arts Canberra.
  6. Book
    1. Charlesworth B
    Sexual Selection: Testing the Alternatives
  7. Book
    1. Charlesworth B
    2. Hughes KA
    The maintenance of genetic variation in Life-History traits
    In: Singh R. S, Krimbas C, editors. Evolutionary Genetics: From Molecules to Morphology, 1. Cambridge: University Press. pp. 369–392.
  8. Book
    1. Crow JF
    2. Kimura M
    An Introduction to Population Genetics Theory
    New York, Evanston and London: Harper and Row Publishers.
  9. Book
    1. de Villemereuil P
    2. Nakagawa S
    (2014) General quantitative genetic methods for comparative biology
    In: Garamszegi L. Z, editors. Modern Phylogenetic Comparative Methods and Their Application in Evolutionary Biology: Concepts and Practice. Berlin: Berlin, Heidelberg. pp. 287–303.
  10. Book
    1. Gaston K
    2. Blackburn T
    Pattern and Process in Macroecology
    John Wiley & Sons.
  11. Book
    1. Gillespie JH
    The Causes of Molecular Evolution
    Oxford: Oxford University Press Google Scholar.
    1. Hedgecock D
    Does variance in reproductive success limit effective population sizes of marine organisms
    Genetics and Evolution of Aquatic Organisms 122:122–134.
  12. Book
    1. Hudson RR
    2. Kaplan NL
    (1994) Gene trees with background selection
    In: Golding B, editors. Non-Neutral Evolution: Theories and Molecular Data. Boston: Springer. pp. 140–153.
  13. Book
    1. Kimura M
    The Neutral Theory of Molecular Evolution
    Cambridge University Press.
  14. Book
    1. Lewontin RC
    2. Singh RS
    3. Uyenoyama MK
    (2004) Building a science of population biology
    In: Singh RS, Uyenoyama MK, editors. The Evolution of Population Biology. Cambridge University Press. pp. 7–20.
  15. Book
    1. Malécot G
    Les mathématiques de l'hérédité
    Paris: Masson.
  16. Book
    1. Mukai T
    Experimental verification of the neutral theory
    In: Ohta T, Aoki K, editors. Population Genetics and Molecular Evolution. Berlin: Springer-Verlag. pp. 125–145.
  17. Conference
    1. Mukai T
    Genotype-environment interaction in relation to the maintenance of genetic variability in populations of Drosophila melanogaster
    Proceedings of the Second International Conference on Quantitative Genetics.
    1. Nei M
    2. Graur D
    Extent of protein polymorphism and the neutral mutation theory
    Evolutionary Biology 17:73–118.
  18. Book
    1. Nevo E
    2. Beiles A
    3. Ben-Shlomo R
    (1984) The evolutionary significance of genetic diversity: Ecological, demographic and life history correlates
    In: Mani GS, editors. Evolutionary Dynamics of Genetic Diversity. Heidelberg: Springer Berlin. pp. 13–213.
  19. Book
    1. Powell JR
    Protein variation in natural populations of animals
    In: Theodosius D, Hecht M, William C. S, editors. Evolutionary Biology, 8. New York: Plenum Press. pp. 79–199.
  20. Book
    1. Reaka-Kudla ML
    2. Wilson DE
    3. Wilson EO
    Biodiversity II: Understanding and Protecting Our Biological Resources
    Joseph Henry Press.
  21. Book
    1. Shirihai H
    The Complete Guide to Antarctic Wildlife: Birds and Marine Mammals of the Antarctic Continent and the Southern Ocean
    Princeton University Press.
  22. Book
    1. Soulé ME
    Allozyme variation, its determinants in space and time
    In: Ayala F. J, editors. Molecular Evolution. Sunderland, Massachusetts: Sinauer Associates. pp. 60–77.
  23. Book
    1. South A
    Rnaturalearth: World Map Data From Natural Earth
    Natural Earth.
  24. Book
    1. Stan Development Team
    Stan Modeling Language Users Guide and Reference Manual
    Stan Developer.
    1. Wang J
    (2005) Estimation of effective population sizes from data on genetic markers
    Philosophical Transactions of the Royal Society B: Biological Sciences 360:1395–1409.
    1. Wright S
    Size of population and breeding structure in relation to evolution
    Science 87:430–431.
  25. Book
    1. Zhang Z-Q
    (2013) Animal biodiversity: An update of classification and diversity in 2013
    In: Zhang Z. Q, editors. Animal Biodiversity: An Outline of Higher-Level Classification and Survey of Taxonomic Richness (Addenda 2013), 3703. Zootaxa. pp. 5–11.

Article and author information

Author details

  1. Vince Buffalo

    Institute for Ecology and Evolution, University of Oregon, Eugene, United States
    Conceptualization, Resources, Data curation, Software, Formal analysis, Validation, Investigation, Visualization, Methodology, Writing - original draft, Writing - review and editing
    For correspondence
    Competing interests
    No competing interests declared
    ORCID icon "This ORCID iD identifies the author of this article:" 0000-0003-4510-1609


National Institutes of Health (1R01GM117241)

  • Vince Buffalo

The funders had no role in study design, data collection and interpretation, or the decision to submit the work for publication.


I would like to thank Andy Kern and Peter Ralph for helpful discussions and supporting me during this work, and Graham Coop for inspiration and helpful feedback during socially distanced nature walks at Yolo Basin. I thank Jessica Stapley for kindly providing the recombination map length data, and Yaniv Brandvain, Amy Collins, Doc Edge, Tyler Kent, Chuck Langley, Matt Osmond, Sally Otto, Molly Przeworski, Jeff Ross-Ibarra, Aaron Stern, Anastasia Teterina, Michael Turelli, Margot Wood, and my Kern-Ralph labmates for helpful discussions. Sarah Friedman, Katherine Corn, and Josef Uyeda provided very useful advice about phylogenetic comparative methods; yet I take full responsibility for any shortcomings of my analysis. Finally, I am indebted to Guy Sella, Matt Pennell, and two other anonymous reviewers for helpful feedback. I would like to also thank UO librarian Dean Walton for helping me track down some rather difficult to find older papers. This work was supported by an NIH Grant (1R01GM117241) awarded to Andrew Kern.

Version history

  1. Received: February 12, 2021
  2. Accepted: August 16, 2021
  3. Accepted Manuscript published: August 19, 2021 (version 1)
  4. Accepted Manuscript updated: August 31, 2021 (version 2)
  5. Version of Record published: October 1, 2021 (version 3)


© 2021, Buffalo

This article is distributed under the terms of the Creative Commons Attribution License, which permits unrestricted use and redistribution provided that the original author and source are credited.


  • 9,366
  • 851
  • 75

Views, downloads and citations are aggregated across all versions of this paper published by eLife.

Download links

A two-part list of links to download the article, or parts of the article, in various formats.

Downloads (link to download the article as PDF)

Open citations (links to open the citations from this article in various online reference manager services)

Cite this article (links to download the citations from this article in formats compatible with various reference manager tools)

  1. Vince Buffalo
Quantifying the relationship between genetic diversity and population size suggests natural selection cannot explain Lewontin’s Paradox
eLife 10:e67509.

Share this article


Further reading

    1. Evolutionary Biology
    Raphael Aguillon, Mieka Rinsky ... Oren Levy
    Research Article

    The circadian clock enables anticipation of the day/night cycle in animals ranging from cnidarians to mammals. Circadian rhythms are generated through a transcription-translation feedback loop (TTFL or pacemaker) with CLOCK as a conserved positive factor in animals. However, CLOCK’s functional evolutionary origin and mechanism of action in basal animals are unknown. In the cnidarian Nematostella vectensis, pacemaker gene transcript levels, including NvClk (the Clock ortholog), appear arrhythmic under constant darkness, questioning the role of NvCLK. Utilizing CRISPR/Cas9, we generated a NvClk allele mutant (NvClkΔ), revealing circadian behavior loss under constant dark (DD) or light (LL), while maintaining a 24 hr rhythm under light-dark condition (LD). Transcriptomics analysis revealed distinct rhythmic genes in wild-type (WT) polypsunder LD compared to DD conditions. In LD, NvClkΔ/Δ polyps exhibited comparable numbers of rhythmic genes, but were reduced in DD. Furthermore, under LD, the NvClkΔ/Δ polyps showed alterations in temporal pacemaker gene expression, impacting their potential interactions. Additionally, differential expression of non-rhythmic genes associated with cell division and neuronal differentiation was observed. These findings revealed that a light-responsive pathway can partially compensate for circadian clock disruption, and that the Clock gene has evolved in cnidarians to synchronize rhythmic physiology and behavior with the diel rhythm of the earth’s biosphere.

    1. Computational and Systems Biology
    2. Evolutionary Biology
    Ryan T Bell, Harutyun Sahakyan ... Eugene V Koonin
    Research Article

    A comprehensive census of McrBC systems, among the most common forms of prokaryotic Type IV restriction systems, followed by phylogenetic analysis, reveals their enormous abundance in diverse prokaryotes and a plethora of genomic associations. We focus on a previously uncharacterized branch, which we denote coiled-coil nuclease tandems (CoCoNuTs) for their salient features: the presence of extensive coiled-coil structures and tandem nucleases. The CoCoNuTs alone show extraordinary variety, with three distinct types and multiple subtypes. All CoCoNuTs contain domains predicted to interact with translation system components, such as OB-folds resembling the SmpB protein that binds bacterial transfer-messenger RNA (tmRNA), YTH-like domains that might recognize methylated tmRNA, tRNA, or rRNA, and RNA-binding Hsp70 chaperone homologs, along with RNases, such as HEPN domains, all suggesting that the CoCoNuTs target RNA. Many CoCoNuTs might additionally target DNA, via McrC nuclease homologs. Additional restriction systems, such as Type I RM, BREX, and Druantia Type III, are frequently encoded in the same predicted superoperons. In many of these superoperons, CoCoNuTs are likely regulated by cyclic nucleotides, possibly, RNA fragments with cyclic termini, that bind associated CARF (CRISPR-Associated Rossmann Fold) domains. We hypothesize that the CoCoNuTs, together with the ancillary restriction factors, employ an echeloned defense strategy analogous to that of Type III CRISPR-Cas systems, in which an immune response eliminating virus DNA and/or RNA is launched first, but then, if it fails, an abortive infection response leading to PCD/dormancy via host RNA cleavage takes over.