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).
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/N≈1/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).
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. 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).
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).
 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 (Nc≈103 and π=3.6×10-3) was an outlier and excluded from this figure for visual clarity. 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.
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).
(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.
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 (γDmel≈2.3×10-3) based Elyashiv et al. implies a rate of sweeps per basepair of νBP,Dmel≈2.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, JDmel≈4.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.
(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-8–10-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. 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 Nc≳107, since S(γDmel,JDmel,L)≈10−8–10−7 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.
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.
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.
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 FST≈1 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.
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/Nc≈10−6 (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. 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.
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.
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=∑i∈SmiNi. 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 1−0.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.
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.
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.
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.
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×10−3 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×10−11. Finally, I estimate JDmel≈4.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 μ=10−9–10−8; to visualize these, I take a convex hull of all diversity ranges and smooth this with R’s smooth.spline function. 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 (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, and the hitchhiking component is (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 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)=1T∑iSα(iS)∑y∈a(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)=αT∑y∈aexp(-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 ly∼U(−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=αT∑i=1mEli(exp(-rBPliτ(s¯,N))) (10) =αTG∑i=1m∫0Gexp(-rBPℓτ(s¯,N))𝑑ℓ (11) =αmTG∫0Gexp(-rBPℓτ(s¯,N))𝑑ℓ Using u-substitution with r=ℓrBP this simplifies to (12) S=αmTGrBP∫0Lexp(-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=νBPrBP∫0Lexp(-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=νBPrBP∫0Lexp(-2rτ(s¯,N))𝑑r⏟J 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 m≈4.5×105, α≈0.42, and T≈8.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 νBP≈2.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=2NS≈0.92; see Figure 4—figure supplement 5B). This implies J≈4.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. 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/(1−Z)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=1−sh−1/2Nc (Santiago and Caballero, 1998). Thus, if 1/2Nc<<sh<<1, then C2/(1−Z)≈U and NeSC98≈Nexp(−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). 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). 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). 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.
Using AIC to compare this full model to a reduced model of log10(π)=β0+βNclog10(Nc), AICfull=204.9, AICreduced=216.4.
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.
Thank you for submitting your article "Why do species get a thin slice of π? Revisiting Lewontin's Paradox of Variation" for consideration by eLife. Your article has been reviewed by 4 peer reviewers, including Guy Sella as the Reviewing Editor and Reviewer #1, and the evaluation has been overseen by Detlef Weigel as the Senior Editor. The following individual involved in review of your submission has agreed to reveal their identity: Matthew Pennell (Reviewer #4). The reviewers have discussed their reviews with one another, and the Reviewing Editor has drafted this to help you prepare a revised submission. Essential revisions: The reviewers appreciated the scholarship and extensive work that went into this manuscript. At the same time, they had several issues with the novel analyses. After some discussion, they suggested that if the paper is revised into more of a review, then it could be of interest to the broad readership of eLife. This would imply substantial streamlining, down-weighting and even removing some analyses, and addressing the reviewers comments on others. Specifically: – Find ways to validate the census size estimates (Reviewers 1 and 2). – The reviewers had a hard time understanding the motivation for the phylogenetic analysis (e.g., Reviewers 1-3). One option is to clarify this motivation, as well as the interpretation of the results; another is to remove this analysis or parts of it. In addition, Reviewer 4, who is an expert on these analyses, had a number of concrete comments that should be addressed. – While the reviewers found the analysis of linked selection effects intriguing, they were unclear about its interpretation. Notably, to what extent is it simply affirming the conclusions of Coop (2016), as opposed to illustrating a much more dramatic effect (e.g., Reviewers 1-3). Reviewer #1 (Recommendations for the authors): 1) Substantial streamlining and editing would make the manuscript much clearer. 2) The abstract is way too long: not on the intended resolution. 3) If you aim for a broad readership, you may consider having the background be less of a historical review (without sacrificing scholarship). Also, no need to recap the historical review at the beginning of the Discussion. 4) You recap what you do at length several times (e.g., L 130-156), which is repetitive. 5) It would be clearer to use consistent terminology, e.g., instead of "enigma", "anomaly", "paradox" and "explanation", "resolution"… 6) The writing switches between acknowledging that several factors are plausibly at work and seeking "a solution". 7) It is claimed that "an ordinary least squares relationship on a log-log scale fits the data well" but I did not find a quantification to this effect. Namely, what proportion of the variance does it explain? Moreover, it is claimed that this relationship is homoscedastic. Can that be quantified as well? From looking at the figure it seems that the regression may explain ~1 of the ~3 orders of magnitude in diversity levels and the residuals explain ~2. It would also be helpful to say what we learn from this fit that we did not learn from staring at the plot. Does one (or more) of the potential explanations for Lewontin's paradox posit a log-log relationship? Reviewer #2 (Recommendations for the authors): I would love to see a revised version of this piece published. Reviewer #3 (Recommendations for the authors): My main suggestion is to expand a little more the arguments for why the author feels particular factors are unlikely to explain the patterns qualitatively, and to touch on some of the alternative explanations raised in the public review. [Editors' note: further revisions were suggested prior to acceptance, as described below.] Thank you for submitting your article "Quantifying Lewontin's Paradox Suggests Natural Selection is Unlikely to Explain the Narrow Range of Diversity Among Species" for consideration by eLife. Your article has been reviewed by 4 peer reviewers, including Guy Sella as the Reviewing Editor and Reviewer #1, and the evaluation has been overseen by Detlef Weigel as the Senior Editor. The following individual involved in review of your submission has agreed to reveal their identity: Matthew Pennell (Reviewer #4). The reviewers have discussed their reviews with one another, and the Reviewing Editor has drafted this to help you prepare a revised submission. All of us find the topic important, the manuscript very interesting and the revised version improved. At the same time, we also felt that the potential impact of the paper is diminished by not (or only partially) addressing key issues raised in the previous round of review. All the more so, if this is to be more of a research paper than a review (although there can be some range in between). Notably: 1) We found that the purpose and insights of the phylogenetic regression model were still not clear. While we had Matthew to answer some of our questions, most of the readers will not have such a resource. 2) We think that more precisely articulating how the linked selection analysis may have moved us forward (i.e., beyond just saying that it did not resolve the paradox entirely) would be helpful. For example, if you believe that it accounts for several orders of magnitude even if it does not account for the shape of the dependency of π on Nc for intermediate values of Nc, then that is important progress. 3) In that regard, we also think that exploring a somewhat wider range of selection parameters (along the lines suggested by reviewer 1 in the previous round) would be helpful in elucidating how far linked selection can and cannot take us. In summary, we think that the impact of your nice work would be enhanced substantially by addressing these issues, but we leave that to your discretion and will not require another round of reviews. Reviewer #1 (Recommendations for the authors): The manuscript has improved substantially in both substance and presentation. I believe that it would be all the more impactful with another thorough round of editing. I attach a pdf with comments/suggestions, but am not a native speaker and am also mildly dyslexic, which is to say that it may be worth getting additional feedback from better writers/editors. A few comments about substance: – The interpretation and necessity of controlling for phylogenetic non-independence. Having read the revised manuscript, your replies and the comments of Reviewer 4, I remain confused about the interpretation of the analysis. For example, I don't understand whether it is about controlling for phylogenetic relationships in errors or in factors that affect the relationship between π and Nc. Moreover, to the extent that it does the latter, what are the assumptions that go into the correction. In the discussion you say that (l. 415-16) "The evidence of high phylogenetic signal found in this study suggests PCMs are needed, in part to avoid spurious results from phylogenetic pseudo-replication." If I understand you correctly, you are suggesting that factors that modify the relationship between π and Nc, such as life history traits, are likely to be more similar in more closely related species, and consequently, the (pi, Nc) points corresponding to closely related species may be considered replication in that they reflect the same processes. If this is what you meant, I think it should be spelled out. Moreover, it would be good to do so early on, in the results, such that the reader can better understand what the PCM analysis is plausibly about. In the Results section you say that (l. 220-223): "If the relationship between diversity and population size was free of shared phylogenetic history, \λ = 0 and all the variances could be explained by evolution on the tips; this is analogous to Lynch's conjecture that coalescent times should be free of phylogenetic signal (2011)." Doesn't this contradict the interpretation discussed in the previous paragraph? Also, I still don't understand the claim that PCM corrects for "spurious pseudo-replication". Shouldn't the determination whether it is "pseudo-replication" depend on the notion of the "true" relationship that you are trying to estimate? Stated differently, maybe some of the phylogenetic signal arises from similarities in factors that affect the pi-Nc relationship and others that just affect π (e.g., mutation rates), a given notion of the "true" relationship would suggest that you want to correct for the latter but not the former. How can PCMs do that without specifying what it is that you are correcting for? I realize that some of these problems may reflect my ignorance about PCMs, but doubt that the general readership you are aiming is much more knowledgeable. – Relationship between genetic map length and Nc. You note that (l 333-5): "These findings are consistent with both the hypothesis that non-adaptive processes increase genome size in small-Ne species (Lynch and Conery, 2003) which in turn could increase map lengths…". I think the map length is largely explained (e.g., R^2=0.96) by the requirement of having one cross-over per chromosome or chromosomal arm (can't remember which). Specifically, I think that this relationship is much stronger than with genome size, and am not sure whether there is any residual effect of size after controlling for the number of chromosomes. – Discussion about different measures of Ne. I found several points in this section deep and insightful. As you point out, there are different N_e's for different quantities, and comparisons among them informs us about different processes. I think it would be helpful to emphasize that these are different quantities rather than estimators of the same quantity and that Lewontin's paradox is specifically about the one defined by diversity levels. A few general comments about presentation: – You often use several different terms for the same thing, e.g., N_e and expected coalescence rate; explain, solve and resolve; genetic/recombinational map length, heterozygosity/pairwise diversity. This is confusing to readers, who wonder whether there was a reason for using the specific term. Choosing one for each term and sticking with it would be clearer. – You sometimes use terms that seem to be private abbreviations, e.g., "strong selection parameters" and "quantifying… paradox". – Perhaps avoid hyphenations as abbreviations, e.g., "low-Nc species" and "across-taxa relationship". - The discussion would benefit from extensive editing. Very nice work! Reviewer #2 (Recommendations for the authors): I have mixed feelings about this revision. The author did not follow the reviewers' suggestion of shifting from a research to a review article, and rather tried to reinforce his original results. One problem I'm having is that I still do not see the point of the sophisticated phylogenetic analysis that is presented. It is clear from the existing literature that π has a strong phylogenetic inertia; I guess the interesting question is what causes this inertia; whether λ is 0.4, 0.6 or 0.8 is good to know but not a major achievement, in my opinion. I have a similar comment regarding the report that genetic map tends to be shorter in large census sized species, which is cool but not totally novel, and not a very strong effect – I do not understand why social insects are excluded from the model, by the way. I like much figure 4b, which substantiates Coop's 2016 argument via empirical data analysis – but again the added value is mainly in the illustration; the argument existed already. In my opinion this version has essentially the same strengths and weaknesses as the previous one: text book-like figures and excellent writing, but no breakthrough as far as the newly reported results are concerned. A missed opportunity, I would say. Reviewer #3 (Recommendations for the authors): In reviewing the author's response letter and the revised draft, the revisions have certainly streamlined the paper and clarified the author's justification for the analyses. However, I still have some concerns about how the work is motivated at the outset, and how the novelty of the results is explained: In the Introduction, the bottom of page 3 the author says that a limitation of past work is that others do not propose a mechanism by which traits act to constrain diversity within a few orders of magnitude. This seems to be a major motivation for doing a new study, but how does this study tackle this question? Did the author have an a priori reason to expect that a more explicit estimation of Nc would resolve Lewontin's paradox, and if so, why? Could the phylogenetic correction have been expected to have resolved Lewontin's paradox, and if so, why (see below)? Correcting for the effect of Nc on map length in tandem with the linked selection model fitting does seem to fall into this category, and with this in mind, I think the author could go a long way towards clarifying the novelty of the work by being more explicit about these results. As noted by reviewer 1 in the first round of reviews the 'glass half full' interpretation of this study for the genetic draft model is that the map length correction combined with linked selection model fitting CAN actually go a significant way towards shrinking down the range of diversity expected. While it can't go ALL the way, this is an important and novel result that is more important to clarify than to simply conclude that we are still in the exact same conclusion zone as before this work. The author has gone some way towards using linked selection to explain the variation, and this is worth clarifying, and perhaps quantifying in the discussion. I still find the motivation for the phylogenetic correction hard to parse in the introduction, and I would suggest front-loading some of the points the author makes in the response about the importance of phylogenetic correction even if coalescent times themselves are not constrained by phylogeny into the intro. Clearly, multiple reviewers found it difficult to understand why phylogenetic correction was needed and wasn't just eroding power, so this is important to front-load. Relatedly, I am still not convinced by the argument that for Lynch's conjecture to be true λ must be zero (page 7). A number of possible interpretations made by the author in the results provide possible mechanisms (like phylogenetic changes in mutation rate) that would still enable coalescent times themselves to be free of a direct phylogenetic effect. The revisions and response now provide a better explanation for the importance of a phylogenetic correction in any case, but I don't think any of the analyses have definitively told us that coalescent times have a phylogenetic signal- it could simply be the phylogenetic inertia of traits associated with Ne and mutation rate. Reviewer #4 (Recommendations for the authors): I liked this paper before and I like it even better now. You have done a thorough and thoughtful job at responding to reviewer comments and have addressed my major critiques. This is a really excellent study. https://doi.org/10.7554/eLife.67509.sa1
Thank you for the feedback. At present (and discussed in more detail in the replies below), the three novel analyses and findings of this study (see reply to Reviewer 2) make this outside the scope of a review, so I am submitting it as a research article.
Reviewer 1’s feedback helped me discover an error in my previous analysis. The census sizes are now a few orders of magnitude smaller for cosmopolitan arthropod species like Drosophila melanogaster. I have also conducted a consistency check by comparing the implied biomasses from my study to previously-published estimates across phyla.
I have significantly reworked how this analysis was framed, in part thanks to the feedback of Reviewer 4. I also have reworked some sections of the discussion to discuss why such analyses are needed.
I have reworked part of the discussion (lines 518-524) to discuss how the findings in this paper fits in, and builds off of the findings of Coop (2016). I discuss this in much more detail below.
Thank you, following this and other feedback, I have edited some unclear sections.
I have reduced the abstract significantly, to 200 words (the limit of eLife).
I have removed the recap at the beginning of the discussion. I have not cut down the historical context as, (1) I believe such a mini review (in the spirit of Coop and Ralph, 2012) is needed, and (2) other reviewers and readers did seem to like this aspect.
I am not sure I follow this comment — this is the first time I describe what I am doing in the manuscript, which seems like an important part of the introduction?
"… I have removed some synonyms for clarity.
I agree that “explanation” or “resolution” are better terms, and have made these changes.
The log-log relationship is used because of heteroscedasticity on a linear-log scale, and because both axes vary over several orders of magnitude. I do not think it is too surprising that π varies a few orders of magnitude, since Nc varies over so many (and Ne a fraction of that). Regarding the proportion of variance, I have included an R2 in the paper (the adjusted R2 is about 0.25). However, it should be noted that R2 has numerous statistical issues (see p. 181-182 of Cosma Schalizi’s The Truth about Linear Regression, http://www.stat.cmu.edu/~cshalizi/TALR/TALR.pdf). Likewise, I don’t think a formal test of heteroscedasticity is particularly helpful; a plot of the residuals versus x values is convincing enough (though I think unnecessary to include in the supplementary materials).
This is useful feedback that other reviewers have voiced as well. I have expanded on how my results connection to that of Coop (2016) on lines 517-524, as well as addressed some more points about the need for PCMs in the discussion (lines 440-473).
I have significantly edited this (lines 197-284) to more clearly explain the purpose of the phylogenetic mixed-effect model and better frame the statistical issues it addresses.
The main problem I show here is that assuming strong selection (i.e. the levels predicted under D. melanogaster), theoretic linked selection models do not match the observed data if Ne = Nc. Unfortunately, this analysis cannot shed light on how much linked selection impacts diversity across species because I use Drosophilamelanogaster parameters that are only suitable for qualitatively checking if linked selection could explain the observed pattern, not quantify the impact across species. Still, my findings are a serious qualitative mismatch between the observed and predicted relationship between diversity and population size. Furthermore, as the new analysis described below shows, this qualitative finding is not sensitive to parameter choice.
I have added a new supplementary figure (p. 47) looking at how a ten-fold increase in BGS (the U parameter) and HH (the γ parameter) would impact the predicted relationship between diversity and population size. This figure also depicts the predicted effect of BGS and HH on selection individually (subfigure A). Overall, this demonstrates that the poor fit found between observed and predicted levels of diversity is not something that can be remedied with stronger parameters. I discuss this on lines 377-386, which also addresses the above issue as well.
Thank you to all the reviewers for feedback.
I agree this was unclear in the way I presented the results; I have fixed this to make it clear I am dealing with phylogenetic correlation in the residuals.
I added a lot in the Results section on pseudo-replication (lines 197-284) to address this and more clearly explain why this correction is needed, and what it does and does not assume.
I agree this was incorrect and have removed this.
This gets at some deeper issues about PCMs, but I think the important part is that phylogenetic signal in either the X or Y variable itself is not a violation of the standard regression model, only error in the residuals (Revell, 2010). Even if the true relationship is driven by, say, a mutation rate change along a lineage, sister taxa that inherit this mutation rate along their parent lineage are not independent observations of this effect in a way that an entirely different clade with its own mutation rate change would be (Maddison and FitzJohn, 2014). The phylogenetic mixed-effect model is only looking for correlation in errors and down-weighting the samples as if they were correlated observations. This is a conservative approach to deal with the structure found in the residuals, and importantly, the findings are not qualitatively different than the OLS.
This is an interesting point, and I do have a figure showing the relationship between chromosome number and map length in the project’s GitHub repository. I’ve run a quick regression of map_length ~ chrom_number + genome_size and find both are significant. However, the direction of causality is unclear and this could be plagued by collider bias. Overall, I think the discussion of this is best left as is (the original sentence was added only to address a reviewer’s concerns in the last round of revisions).
I have tried to make this clear with “These various estimators capture different summaries of effective population size on shorter timescales than coalescent-based estimators”.
I have made these changes (though in some cases allelic heterozygosity, not per basepair diversity was measured by these earlier allozyme studies).
I have made these changes.
I have kept the low-Nc, high-Nc, etc as this saves a lot of text and makes sentences shorter.
I have made the majority of these changes, thank you for your feedback!
I think it is important to point out that Coop (2016) only refuted Corbett-Detig et al’s claim that their estimated impact of linked selection is strong enough to explain Lewontin’s Paradox. Since Corbett-Detig’s work does not have census population sizes, only proxies of population size, we have no concept of the magnitude of the effect to be explained, and whether linked selection can explain it. This can be seen in Coop’s (2016) Figure 1, where the dashed red line is a hypothetical relationship — in my work, by quantifying the π-Nc relationship, we can get a sense of the scale. Finally, my work goes a step further than Coop’s work, showing that not only do Corbett-Detig’s estimates not go far enough in showing linked selection can explain Lewontin’s Paradox, but that no estimates, not even those orders of magnitude larger, can explain Lewontin’s Paradox (see Figure 4-Supplement 4).
See the response to reviewer #2 above. Also, if the relationship between π and Nc were no longer significant after accounting for phylogeny, and within-clade regressions were similarly non-significant, this would not solve Lewontin’s Paradox. Rather, this would suggest that the relationship is even less dependent on Nc than we expected, and perhaps some macroecological or macroevolutionary process could explain this.
I have addressed this with the main comments. The issue here is that the findings here are qualitative; they show that linked selection models are incapable of explaining Lewontin’s Paradox under strong selection expected of D. melanogaster. Given that it would be erroneous to assume that linked selection in all species is as strong as it is in D. melanogaster, I cannot use these results to quantify the impact of linked selection across species.
I have clarified this in the PCM section.
I discuss this in the discussion a bit (see the mention of Uyeda et al. around lines 432-438). Overall, a trait that has phylogenetic signal is not a violation of the regression model, and not a problem; only if the residuals have correlated errors is this a problem (see Revell 2010 and Uyeda et al. 2018). https://doi.org/10.7554/eLife.67509.sa2 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.
© 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.
Article citation count generated by polling the highest count across the following sources: Crossref, PubMed Central, Scopus.
|
The following image depicts four selection scenarios on the trait of body size
Pos Terkait
Copyright © 2024 idkuu.com Inc.
