## Abstract

The effective population size (*N*_{e}) is one of the most fundamental parameters in population genetics. It is thought to vary across the genome as a consequence of differences in the rate of recombination and the density of selected sites due to the processes of genetic hitchhiking and background selection. Although it is known that there is intragenomic variation in the effective population size in some species, it is not known whether this is widespread or how much variation in the effective population size there is. Here, we test whether the effective population size varies across the genome, between protein-coding genes, in 10 eukaryotic species by considering whether there is significant variation in neutral diversity, taking into account differences in the mutation rate between loci by using the divergence between species. In most species we find significant evidence of variation. We investigate whether the variation in *N*_{e} is correlated to recombination rate and the density of selected sites in four species, for which these data are available. We find that *N*_{e} is positively correlated to recombination rate in one species, *Drosophila melanogaster*, and negatively correlated to a measure of the density of selected sites in two others, humans and *Arabidopsis thaliana*. However, much of the variation remains unexplained. We use a hierarchical Bayesian analysis to quantify the amount of variation in the effective population size and show that it is quite modest in all species—most genes have an *N*_{e} that is within a few fold of all other genes. Nonetheless we show that this modest variation in *N*_{e} is sufficient to cause significant differences in the efficiency of natural selection across the genome, by demonstrating that the ratio of the number of nonsynonymous to synonymous polymorphisms is significantly correlated to synonymous diversity and estimates of *N*_{e}, even taking into account the obvious nonindependence between these measures.

THE effective population size (*N*_{e}) is one of the most fundamental quantities in population genetics, evolutionary biology, and molecular ecology, since it determines the effectiveness of natural selection and the level of neutral genetic diversity that a population contains (Charlesworth 2009). Populations and regions of the genome with small *N*_{e} tend to have low levels of genetic diversity, to be susceptible to the accumulation of deleterious mutations through genetic drift, and to have potentially low rates of adaptive evolution (Charlesworth 2009).

The effective population size is expected to vary across the genome as a consequence of genetic hitchhiking (Smith and Haigh 1974) and background selection (Charlesworth *et al.* 1993). The action of both positive and negative natural selection, particularly in regions of the genome with low rates of recombination, is expected to reduce the effective population size leading to lower levels of genetic diversity and reduced effectiveness of selection. Hence variation in the rate of recombination and the density of selected sites is expected to generate variation in *N*_{e.}

The evidence that there is variation in *N*_{e} within a genome comes from three sources. First, it has been shown that levels of neutral genetic diversity are correlated to rates of recombination in *Drosophila* (Begun and Aquadro 1992), humans (Hellmann *et al.* 2003), and some plant species (Tenaillon *et al.* 2004; Roselius *et al.* 2005). This could be due to variation in the mutation rate since neutral genetic diversity is proportional to the effective population size multiplied by the mutation rate. However, the level of neutral sequence divergence between species, which should be proportional to the mutation rate, is not correlated to the rate of recombination in *Drosophila* (Begun and Aquadro 1992) and the plant species (Roselius *et al.* 2005) that have been investigated. Furthermore, although there is a correlation between neutral sequence divergence and recombination rate in humans, this correlation is not sufficient to explain the correlation between diversity and the recombination rate (Hellmann *et al.* 2005). It is has also been shown that the Y and W chromosomes, which have no recombination over most of their length, have substantially lower diversity than other chromosomes and that this cannot be attributed to differences in the mutation rate or the fact that there are fewer Y and W chromosomes than autosomes (Filatov *et al.* 2001; Montell *et al.* 2001; Bachtrog and Charlesworth 2002; Hellborg and Ellegren 2004). It thus seems that the effective population size varies across genomes and is positively correlated to the rate of recombination.

Second, under the neutral theory of molecular evolution it is expected that levels of diversity and divergence should be proportional to each other, since both depend on the neutral mutation rate. Deviations from this hypothesis, caused by variation in *N*_{e}, can be tested using the HKA test and derivatives of it (Hudson *et al.* 1987; Ingvarsson 2004; Wright and Charlesworth 2004; Innan 2006). Evidence for departures from the neutral hypothesis, based on the HKA test, comes from multiple multilocus surveys in plants (Roselius *et al.* 2005; Schmid *et al.* 2005), the chicken Z chromosome (Sundström *et al.* 2004), humans (Zhang *et al.* 2002), and *Drosophila* (Moriyama and Powell 1996; Machado *et al.* 2002).

Third, variation in the effective population size should manifest itself as variation in the effectiveness of selection and this has also been observed. In *Drosophila* it has been shown that codon usage bias is lower in the regions of the genome with very low rates of recombination (Hey and Kliman 2002; Kliman and Hey 2003; Marais *et al.* 2003). It has also been shown that the number of nonsynonymous polymorphisms (*P*_{n}) relative to the number of synonymous polymorphisms (*P*_{s}) is higher in the low-recombining parts of the *Drosophila melanogaster* genome (Presgraves 2005), that the rate of nonsynonymous (*d*_{N}) relative to the rate of synonymous (*d*_{S}) substitution is positively correlated to the frequency of recombination (Betancourt and Presgraves 2002), and that the overall efficiency of selection appears to be lower in the regions of the genome with low rates of recombination (Presgraves 2005; Larracuente *et al.* 2008). Likewise it has been shown that *d*_{N}/*d*_{S} is higher on the Y or W chromosome than on the other chromosomes in humans (Wyckoff *et al.* 2002) and birds (Berlin and Ellegren 2006) and on the fourth chromosome of *Drosophila* species (Arguello *et al.* 2010). In contrast, Bullaughey *et al.* (2008) found no correlation between *d*_{N}/*d*_{S} and the rate of recombination in primates.

It is thought that the correlation between *d*_{N}/*d*_{S} or *P*_{n}/*P*_{s} and the rate of recombination is due to regions of the genome with little or no recombination having low effective population size and hence reduced effectiveness of natural selection (Betancourt *et al.* 2009). *P*_{n}/*P*_{s} is negatively correlated to the rate of recombination because regions with low effective population size allow more slightly deleterious mutations to segregate for a longer time. In contrast, *d*_{N}/*d*_{S} can be either positively or negatively correlated to the rate of recombination depending on the prevalence of advantageous mutations. If advantageous mutations are common, then regions of the genome with high rates of recombination are expected to evolve faster because they have a higher effective mutation rate and because selection is effective on a greater proportion of mutations. In contrast, if advantageous mutations are rare, then regions of the genome with high rates of recombination may have low values of *d*_{N}/*d*_{S} because selection against slightly deleterious mutations is more effective.

Although it is well established that *N*_{e} varies across the genome in a few species, it is unclear whether this is true of all species and, more importantly, how much variation in *N*_{e} there is and whether this variation results in differences in the effectiveness of selection. Here we test whether there is variation in the effective population size by considering whether there is significant variation in neutral diversity, taking into account that this might be due to variation in the mutation rate by using the divergence between species to control for differences in the mutation rate. We also quantify the variation in *N*_{e}. We estimate *N*_{e} from the nucleotide diversity at putatively neutral sites, since this is expected to be equal to 4*N*_{e}μ in a diploid organism, where *N*_{e} is the effective population size and μ is the mutation rate per generation. We use the divergence between two species at neutral sites as an estimate of the mutation rate per generation. Note that since we are comparing loci within a genome, they all share the same generation time (unless they are on the sex chromosomes or in the mitochondrial DNA) and so this does not have to be explicitly taken into account. We can therefore estimate the effective population size for each locus. However, although each individual estimate is unbiased, the distribution of these values has a variance that is greater than the true variance because of sampling error; a locus might have a particularly low diversity just by chance and not because its effective population size is particularly low. To get around this problem we use a hierachical Bayesian framework to estimate the distribution of *N*_{e} across genes, taking into account the sampling error associated with both the polymorphism and the divergence data.

We test for and investigate the variation in the effective population size in 10 eukaryotic species including humans, *D. melanogaster*, *A. thaliana*, and *Saccharomyces paradoxus* (Table 1). We find that there is statistically significant variation in *N*_{e} across genes, but that it is rather modest in most of the organisms. We also investigate whether variation in *N*_{e} within a genome leads to variation in the proportion of effectively neutral mutations, by testing whether the ratio of the number of nonsynonymous to synonymous polymorphisms is correlated to the effective population size, in a way that circumnavigates the obvious nonindependence between the two variables. We find overall evidence for a correlation between these two parameters and hence conclude that even modest variation in the effective population size is sufficient to generate variation in the effectiveness of natural selection.

## Materials and Methods

### Sequence data

We obtained data from different plant species, mouse, fruitfly, and yeast using publicly available data from GenBank (http://www.ncbi.nlm.nih.gov/Genbank). Polymorphism data for *Homo sapiens* were downloaded from Enivironmental Genome Project (egp.gs.washington.edu) and Seattle SNPs (pga.gs.washington.edu) Web sites and for *Arabidopsis thaliana* from http://walnut.usc.edu/2010. The annotated protein-coding genome of *A. thaliana* was obtained from TAIR 8 (ftp://ftp.arabidopsis.org), and the annotated *A. lyrata* genome was obtained from JGI http://genome.jgi-psf.org/. The annotated protein-coding genomes of *Pan troglodytes*, *Macaca mulatta*, and *Rattus norvegicus* were obtained from Ensembl (http://www.ensembl.org/info/data/ftp/index.html). The *S. cerevisae* genome chromosome III was obtained from http://www.yeastgenome.org. We restricted our analysis of *D. melanogaster* to data from the Zimbabwe population, of the *S. paradoxus* data set to the European population, and of the human data set to African populations, since all of these represent the ancestral populations of the three species (Garrigan and Hammer 2006; Stephan and Li 2007; Liti *et al.* 2009). Qualitatively similar results were obtained in the three cases when using global data.

### Preparation of the data

The analysis was performed using protein-coding sequences. Coding regions were assigned using protein-coding genomic data or, if given, were taken from the GenBank input files. Sequences were aligned using Clustalw, using default parameter values (Thompson *et al.* 1994). The outgroup ortholog was assigned using the best BLAST (Altschul *et al.* 1990) hit or, if given, was taken from the polymorphism data set. We used only polymorphism data for which we could assign an outgroup sequence. For all analyses the number of synonymous substitutions and polymorphisms served as the neutral standard. For computational reasons all sites had to have been sampled in the same number of chromosomes within each species; because some loci had been sampled in more individuals than others and other loci had missing data, we reduced the data set to a common number of chromosomes by randomly sampling the polymorphisms at each site without replacement. The numbers of synonymous and nonsynonymous sites and substitutions were estimated by randomly selecting one allele from the polymorphism data and comparing it against the outgroup using the F3x4 model implemented in PAML (Yang 1997) in which codon frequencies are estimated from the nucleotide frequencies at the three codon positions. The proportion of sites estimated by PAML was also used to compute the numbers of synonymous and nonsynonymous sites for the polymorphism data. Although how we choose to define a site can be important in some circumstances (Bierne and Eyre-Walker 2003), this is not likely to be a problem in the current context because we use the same definition for both the divergence and the polymorphism data; as such, the number of sites effectively cancels out in most of our analyses (however, see discussion of selection on synonymous codon bias below). Statistics concerning numbers of loci and numbers of sites as well as polymorphic sites are shown in Table 1.

### Testing for variation in diversity and the effective population size

We investigated whether there was significant variation in the level of diversity across the genome, using two tests. If we assume there is free recombination within and between loci (or no recombination within and between loci), then variation in diversity can be tested using a simple (2 × *k*) χ^{2}-test of independence across the *k* loci within each species, where for each locus we have the number of sites with a polymorphism and the number of sites without a polymorphism. Note that this test is valid only when the same numbers of chromosomes have been sampled across all loci. However, some of the variation in diversity between loci might be due to variation in the genealogy if there is limited or no recombination between loci. We therefore applied a variant of the classic HKA test, but we removed the divergence information from the test. The test statistic *X*^{2} is set up as(1)where and are the expected value and variance of the number of segregating polymorphisms, *P*, in gene *i*,(2)(3)with *n* being the number of alleles, *M* the number of loci, θ = 4*N*_{e}μ, and *L _{i}* the number of sites in gene

*i*. Estimates of θ were obtained by minimizing the value of

*X*

^{2}.

*X*

^{2}is expected to be χ

^{2}-distributed with (

*L*− 1) d.f.

Any variation that we detect in diversity might be due to variation in the mutation rate or variation in the local effective population size. We therefore performed two further analyses to investigate whether there was variation in diversity that could not be explained by variation in the mutation rate, as measured by synonymous divergence between species. The first test was a second approximate (2 × *k*) χ^{2}-test of independence, performed as follows. For each locus we have the number of sites used to estimate the level of silent site divergence (*L*_{d}), the estimated number of substitutions (*D*), the number of sites used to estimate silent site diversity (*L*_{p}), and the number of sites with a polymorphism (*P*). Since *L*_{d} and *L*_{p} can be different, we reduced the divergence or polymorphism data set, whichever was larger, to the size of the other, resampling without replacement the numbers of substitutions or polymorphisms as appropriate; for example, if *L*_{d} was half *L*_{p}, we sampled *L*_{p} sites from the divergence data to generate a subsample of the substitutions (*D*′) over sites. We can then perform a (2 × *k*) χ^{2}-test where the cells for each gene are the number of sites with a substitution (*D*′) and the number of sites with a polymorphism (*P*′). Note that the data set will be reduced using this method, resulting in a loss of power. Furthermore, this test is only approximate because we assume that the number of substitutions is binomially distributed, whereas in fact it has a more complex distribution because of the correction for multiple hits. Some of the expected values can be very small in both χ^{2}-tests: we therefore checked the *P*-values from the χ^{2}-tests by generating the null distribution for the test. This was performed by randomly assigning polymorphisms and substitutions across the contingency table preserving the marginal totals. We then recalculated the statistic and performed this 1000 times. The *P*-value was the proportion of such randomly generated values that exceeded the observed value. Generally we found that the *P*-value from randomization and the *P*-value assuming our test statistics were χ^{2}-distributed were similar (Supporting Information, Table S1). We therefore present the results from the standard χ^{2}-test.

This test assumes free recombination between sites within loci and loci (or no recombination between sites and loci). A more conservative test is the classic HKA test that tests for heterogeneity in the ratio of diversity divided by divergence between loci assuming no recombination within loci, but free recombination between loci. We performed the multiple-locus HKA test using software provided by J. Hey (http://genfaculty.rutgers.edu/hey/software#HKA). To perform this test we had to exclude loci with zero divergence; for most species this constituted a small fraction of the total number of loci. However, we had to exclude *Sorghum bicolor* from the analysis because too many loci showed zero divergence.

### Recombination and density of selected sites

We obtained estimates of recombination rate variation along chromosomes for *A. thaliana* (Singer *et al.* 2006), *D. melanogaster* (Hey and Kliman 2002), *H. sapiens* (Kong *et al.* 2002), and *Mus musculus* (Dumont *et al.* 2011). Gene density was estimated as the proportion of coding sites in window sizes of 50 kb, 500 kb, and 5 Mb. Since results are qualitatively similar, we discuss only results for the window size of 500 kb. Conservation scores (Siepel *et al.* 2005) were obtained from the UCSC genome browser (http://genome.ucsc.edu/) for *D. melanogaster* across 15 species, *H. sapiens* across 17 species, and *M. musculus* across 30 species.

### Bayesian analysis

To estimate the distribution of *N*_{e} we used a hierachical Bayesian analysis in which we estimate the parameters of the distribution of *N*_{e} (Figure S1). If we assume that the population size is stationary, the expected number of polymorphisms segregating in a sample of *n* sequences, *P*, and the number of differences between the outgroup and a single sequence from the ingroup, *D*_{s}, are(4)(5)where *L*_{p} and *L*_{d} are the numbers of sites that can have a polymorphism or a substitution, respectively, μ is the nucleotide mutation rate per generation, and *t* is the time of divergence. We are interested in the distribution of *N*_{e}. To estimate this distribution we assume that *N*_{e} and μ follow a log-normal or a gamma distribution. Assuming free recombination and using Equations 4 and 5 above, we can write the likelihood of observing polymorphisms and substitutions,(6)where *X*(*S*, *S*(*x*)) is the Poisson distribution and is the probability density of the distribution of *N*_{e}, and *M*(μ|σ_{μ}) is the probability density of the distribution of the mutation rate; these distributions are parameterized such that the mean is fixed at unity, leaving us to estimate the shape parameter. If there is no recombination within a locus, then we can rewrite Equation 4 as(7)where τ is the length of the genealogy scaled such that *E*[τ] = 1. We can rewrite Equation 6, and the likelihood then becomes

To calculate the probability density distribution *M*(τ | *n*) of genealogy lengths we randomly simulated 10,000 genealogies, scaling them such that the average total length was unity. In theory it is possible to accommodate ancestral polymorphism into the method; however, we found that the method rarely gave stable estimates of , particularly in the no recombination model. We therefore concentrated on data sets in which the influence of ancestral polymorphism was likely to be minimal—*i.e.*, in which the average divergence was more than five times the average of θ_{W} (Table 1). If we assume that the ancestral *N*_{e} of a locus is correlated to the current *N*_{e}, we expect ancestral polymorphism to decrease the apparent variation in *N*_{e}.

To estimate the posterior distribution of the parameters and σ_{μ} we used a Monte Carlo Markov chain running the Metropolis–Hastings algorithm (Hastings 1970). Unfortunately because we have very few synonymous polymorphisms per gene, this method tends to underestimate the true value of . For most data sets this underestimation is small, but it can be large. We therefore estimated the extent of bias by simulating data under a range of parameter values, using the actual numbers of sites from the real data such that the expected numbers of polymorphisms and substitutions were equal to the mean values. For example, if we estimated to be 0.5 and σ_{μ} to be 0.1, we simulated data for σ_{μ}-values of 0.1, 0.2, and 0.3 and for -values between 0.4 and 1.0 in steps of 0.05. For each simulated data set we estimated and using linear regression we inferred the relationship between (estimated) and (true). Using this relationship we inferred the true value of from the value estimated from the real data (Figure S2 and Figure S3). To obtain a corrected SE we multiplied the observed standard error by the ratio of the corrected estimate of divided by the observed estimate of . This slightly underestimates the true SE since we have not taken into account the small amount of error associated with estimating the regression line. To test for heterogeneity in between species we assumed that the estimate of was normally distributed; under this assumption is χ^{2}-distributed with *k* − 1 d.f. for *k* species. was calculated as a weighted average, where the weights were inversely proportional to the variance of the estimate (Eyre-Walker 1996).

### Variation of efficiency of selection

We tested whether the strength of selection on nonsynonymous mutations was correlated to the effective population size, which can be seen as testing whether the fraction of deleterious mutations varies with *N*_{e}. This can be done by considering the correlation of *P*_{n}/*P*_{s} and θ_{s} or *P*_{n}/*P*_{s} and *N*_{e} (= θ_{s}/(4μ)), where *P*_{n} and *P*_{s} are the numbers of nonsynonymous and synonymous mutations, respectively, and *N*_{e} values are point estimates from the genetic diversity and mutation rates taken from the literature. However, *P*_{s} and θ_{s} are not independent. We overcome this problem by splitting *P*_{s} into two independent values by generating a random hypergeometric variable as follows:(9)(10)(Piganeau and Eyre-Walker 2009; Stoletzki and Eyre-Walker 2011). One of the *P*_{s} values is used to estimate *P*_{n}/*P*_{s} (see below) and the other one is used to estimate θ_{S}. There are two further problems to consider with this method: first, *P*_{n}/*P*_{s} can be an overestimate or an underestimate of the true value of *P*_{n}/*P*_{s} and, second, the ratio *P*_{n}/*P*_{s} is undefined if *P*_{s} = 0. Both of these problems can be overcome by considering the correlation between ψ and θ_{s}:(11)(Piganeau and Eyre-Walker 2009). Hence, using our method to split *P*_{s} into independent values, we have two independent pairs of θ_{S} and ψ; we present results from only one pair. Some of the data sets contain relatively little polymorphism, which results in substantial variance of ψ. To overcome this problem we sum data across loci. For this we ranked loci according to their neutral diversity obtained from θ_{s2} and binned them into groups of size *n* (*e.g.*, 2, 4, 8, and 16). For each group average θ_{s2} and corresponding *N*_{e2} values were calculated. Furthermore, for each group, the sums of *P*_{n} and *P*_{s1} were calculated to calculate ψ_{1}. Note that ψ_{2} can be obtained in a similar manner; however, results were qualitatively comparable and we therefore show only results for ψ_{1} *vs.* θ_{2} and *N*_{e2}. Also we show only results for group size 4 because results for group sizes >2 were similar. The correlations were performed by calculating Spearman’s rank correlation and probabilities were combined using the unweighted *Z* method (Whitlock 2005). For summary of symbols see File S1.

## Results

To investigate variation in the effective population size within genomes, we assembled protein-coding sequences from 10 species. The data sets are from 6 plant species, 3 animal species, and 1 fungus. The data sets range in size from 66 to 918 loci per species and from 8 to 40 sequences per gene (Table 1). In all analyses we assume that synonymous mutations are neutral.

### Variation of diversity and *N*_{e} within a genome using χ^{2}- and HKA tests

The level of genetic diversity appears to vary considerably within each genome (Figure 1); however, the number of polymorphisms per gene is generally quite low and hence this variation might be due to sampling error. To test whether the variation is significant we used two tests, which make different assumptions about the rate of recombination within loci: either free or no recombination. Both tests suggest that there is variation in the level of diversity in most species: all species are significant assuming free recombination and 6 of 10 are significant assuming no recombination (Table 2). This variation in diversity between loci could be due to variation in the effective population size or to variation in the mutation rate. To investigate whether variation in the mutation rate might be responsible, we estimated the number of synonymous substitutions for each locus (*D*_{S}), between the species of interest and an outgroup species. In many species there is a significant positive correlation between the numbers of synonymous substitutions *d*_{s} and polymorphisms *p*_{s} per site (Table 3), suggesting that part of the variation in diversity is due to variation in the mutation rate. However, if we test whether there is significant variation between loci, taking into account the mutation rate, as estimated from the divergence between species, using either a χ^{2}-test of independence or the more conservative HKA test, then we find significant evidence in the majority of species, whether or not we assume free or no recombination within loci: 9 of 10 loci for the free recombination test and 6 of 9 loci for the no recombination test (the HKA test could not be performed on *S. bicolor* due to the large number of genes in which the divergence was zero) (Table 2).

### Correlates of *N*_{e}

The variation in *N*_{e} across the genome is likely to be due to genetic hitchhiking and background selection. Both processes are expected to be stronger in regions of the genome with low rates of recombination and a high density of sites subject to natural selection. To investigate which or whether either of these factors is responsible for the variation in *N*_{e}, we investigated whether the variation in *N*_{e} was correlated to the rate of recombination and density of selected sites in four species for which these data were available: *D. melanogaster*, human, mouse, and *A. thaliana*. We measured the density of selected sites as either the number of nucleotides in annotated exons (genic density) or the number of nucleotides in conserved regions (conserved site density), as annotated in the UCSC conservation track, in windows of size 50 kb, 500 kb, and 5 Mb, where the window is centered on the gene from which the polymorphism data were taken (there is no conservation track for *A. thaliana*, so in this species we investigated just the density of genic sites). Results for the different window sizes were generally consistent, so we present the results from the 500-kb window size. We estimated *N*_{e} as the synonymous diversity divided by synonymous divergence.

In *D. melanogaster* we find, as others have done, that our estimate of *N*_{e} is positively correlated to recombination rate (Spearman’s correlation coefficient *r* = 0.45, *P* < 0.01). It is, however, also positively correlated to the density of conserved sites (*r* = 0.24, *P* < 0.01), which is unexpected, although not genic sites (*r* = 0.03, *P* = 0.65). The positive correlation with conserved site density might be due to the positive correlation that exists between the density of conserved sites and the rate of recombination (*r* = 0.56, *P* < 0.01), and indeed if we perform a multiple regression, we find that the correlation between *N*_{e} and the density of conserved sites disappears (*P* = 0.74), while the positive correlation between *N*_{e} and recombination rate remains (*P* < 0.01).

In humans we find, as others have done, that both diversity (*r* = 0.14, *P* = 0.02) and divergence (*r* = 0.18, *P* < 0.01) are positively correlated to the rate of recombination (Lercher and Hurst 2002; Hellmann *et al.* 2005), and there is, as a consequence, no correlation between estimates of *N*_{e} and the rate of recombination (*r* = 0.026, *P* = 0.69). *N*_{e} is significantly negatively correlated to the density of genic sites (*r* = −0.19, *P* < 0.01), but not conserved sites (*r* = −0.085, *P* = 0.17). Using multiple regression does not alter this picture: *N*_{e} is correlated only to the density of genic sites.

In mouse we see no significant correlations between estimates of *N*_{e} and the rate of recombination (*r* = 0.054, *P* = 0.72) and the density of genic (*r* = 0.089, *P* = 0.53) or conserved sites (*r* = 0.093, *P* = 0.51). This picture is unaffected by the use of multiple regression.

In *A. thaliana* we see a pattern like that in humans: both diversity (*r* = 0.10, *P* = 0.04) and to a lesser extent divergence (*r* = 0.064, *P* = 0.11) are positively correlated to recombination rate, and *N*_{e} is positively but not significantly correlated to recombination rate (*r* = 0.080, *P* = 0.11). *N*_{e} is significantly negatively correlated to genic density (*r* = −0.11, *P* = 0.02). Unfortunately there are no data on conserved sites in this species.

### Quantifying variation of *N*_{e}

Since we find evidence for variation in *N*_{e} in many of our species, we attempted to quantify the amount of variation using a hierarchical Bayesian model. We assume underlying distributions for *N*_{e} and μ (*e.g.*, log-normal distributions) and estimate the shape parameters and σ_{μ} and hence the variances of these distributions; the mean of each distribution is constrained to be equal to one (see *Materials and Methods*). We investigate two different models: in the first we assume free recombination and in the second we assume no recombination within loci, but free recombination between loci. These two models are likely to set the upper and lower bounds on the true level of variation in *N*_{e}. Under the free recombination model all the variation in diversity is attributed to variation in *N*_{e}, variation in the mutation rate, and sampling error. In the model with no recombination, variation in diversity may additionally be due to variation in the coalescent process. Hence, the free recombination model gives an upper estimate on the variation in *N*_{e} and the no recombination model gives a lower bound.

We applied our method to the polymorphism data from each of the 10 eukaryotic species to estimate the variation of *N*_{e} within each genome along with the variation in the mutation rate, σ_{μ} (Table 4). As expected, in all cases the estimate of is larger when free recombination is assumed, but the estimates from the two models are highly correlated (*r* = 0.95). The estimate of σ_{μ} is unaffected by the model of recombination assumed. We find evidence that the value of varies between species for both the free and the no recombination models (*P* = 2.5 × 10^{−9} and *P* = 4.2 × 10^{−8}, respectively). We find that the level of variation of *N*_{e} is the lowest for *M. musculus* and highest for *Capsella rubella* for both recombination models. The estimates of and σ_{μ} were of similar magnitude for each taxon, suggesting that overall variation in the mutation rate and variation in the effective population size contribute a similar amount to the variation in diversity.

The level of variation in *N*_{e} we estimate using our method is quite modest. For example, *C. rubella* has the highest estimate of , but under this distribution the genes in the 90th percentile have an *N*_{e} that is only 7.2-fold greater than that of those in the 10th percentile; *i.e.*, 80% of genes have an effective population size within 7.2-fold of each other. Four species have estimates of < 0.6, meaning that the difference between the 90th and the 10th percentile is <4-fold.

The estimated distribution appears to fit the data reasonably well (Figure 2). We would not expect the fit to be perfect, particularly at the lower end of the distribution, since this is where sampling error is a major issue; *e.g.*, many genes have no polymorphism because of sampling error, not because they have an effective population size of zero. It is possible that assuming a log-normal distribution places some unwanted constraints on the estimation procedure; in particular, the probability density tends to zero for low *N*_{e.} We therefore also fitted a gamma distribution to the data (Table S3); with this distribution the probability density does not necessarily decline to zero near the origin. However, the estimated distributions are very similar to those obtained assuming a log-normal distribution (Figure S4 and Figure S5). The species that show low variation in *N*_{e} are also those that tend to show little evidence of variation in *N*_{e}, as judged by the χ^{2}- and HKA tests. This implies that failure to detect variation in *N*_{e} is largely because there is limited variation in *N*_{e} rather issues with statistical power.

### Variation in the efficiency of selection

Although we estimate the variation in the effective population size to be modest, it is of interest to investigate whether this translates into significant differences in the efficiency of natural selection across the genome. To investigate this we tested whether there was a correlation between ψ = *P*_{n}/(*P*_{s} + 1) and either θ_{s} or *N*_{e} for each locus in a manner that controls for the obvious nonindependence of the two variables (see *Materials and Methods*). We remove the nonindependence by splitting *P*_{s} into two independent parts and we use ψ because it reduces the bias inherent in the estimation of *P*_{n}/*P*_{s}; furthermore it allows *P*_{n}/*P*_{s} to be calculated for all genes (Piganeau and Eyre-Walker 2009). This test is not very powerful since ψ has a large variance; furthermore, it is statistically biased in a manner that tends to generate a positive correlation between ψ and θ_{s} or *N*_{e}. We therefore follow the approach suggested by Piganeau and Eyre-Walker (2009) and grouped genes according to their θ or *N*_{e} value. The results are qualitatively similar for groupings of 4, 8, and 16 genes, so we present the results for groups of 4. There is a significant negative correlation between both θ_{s} and ψ and *N*_{e} and ψ in *A. lyrata* and *C. grandiflora*, and a marginally significant correlation between ψ and θ* _{s}* in

*D. melanogaster*, although only the correlations in

*Cryptostegia grandiflora*are significant after correction for multiple tests; otherwise the correlations are generally weak and nonsignificant. However, overall we find significant evidence for a negative correlation between ψ and θ

_{s}or

*N*

_{e}if we combine probabilities: between ψ and θ

_{s}

*P*= 0.043 and between ψ and

*N*

_{e}

*P*= 0.021.

The relationship between ψ and *N*_{e} can potentially yield information about the distribution of fitness effects (DFE) (Loewe and Charlesworth 2006; Loewe *et al.* 2006; Woolfit 2006; Elyashiv *et al.* 2010). If we assume that the DFE for nonsynonymous mutations is a gamma distribution and that synonymous mutations are neutral, then *P*_{n}/*P*_{s} is expected to be proportional to , where β is the shape parameter of the gamma distribution (Welch *et al.* 2008). Hence we can estimate β by considering the slope of the regression line between log(ψ) and log(*N*_{e}). Since the log of zero is undefined, we grouped genes in groups of size *n* such that no group had a zero estimate of ψ or *N*_{e}. We attempted to estimate β in the species that individually showed a significant correlation between ψ and *N*_{e}. However, we could not perform the analysis of *A. lyrata* because the diversity is so low that it was impossible to define groups that did not have zero values for both ψ and *N*_{e}. The estimates of β using this method are 0.41 (SE = 0.15) in *C. grandiflora* and 0.23 (0.15) in *D. melanogaster*; these are similar to those obtained using an independent method that uses the site frequency spectrum (Keightley and Eyre-Walker 2007): 0.27 (0.08) for *C. grandiflora* and 0.29 (0.07) for *D. melanogaster* (Table S2). This suggests that the gamma distribution is a reasonable approximation to the DFE, at least for mutations of weak effect.

## Discussion

The effective population size (*N*_{e}) is one of the most important parameters in population genetics and evolutionary biology. It has been shown that *N*_{e} varies across the genome of *D. melanogaster* and some plant species, and it is thought that it might vary across the human genome (Hellmann *et al.* 2005). Here we have shown that it varies in most species that we have considered. However, the variation in *N*_{e} is not consistently correlated to either the rate of recombination or the density of selected sites. This might in part be because the variation in *N*_{e} is quite limited: most genes in a genome have an *N*_{e} that is within a few fold of that of most other genes. Nevertheless the variation is sufficient to cause differences in the effectiveness of natural selection on segregating nonsynonymous polymorphisms.

There are a number of factors that might have led us to over- or underestimate the variation in *N*_{e}. First, we have assumed that there is either free recombination or no recombination within loci to estimate the variation in the effective population size. This is unsatisfactory since we know that recombination is one of the factors that generates variation in the effective population size, at least in species like *Drosophila*, in which there is a correlation between diversity and the rate of recombination. Unfortunately it is not easy to get around this problem. However, as we have noted earlier, the estimate assuming free recombination should give an upper estimate on the amount of variation, because under this method all variation in the diversity is assumed to arise from sampling error and variation in the mutation rate and *N*_{e}. In reality, some of the variation between genes is a consequence of variation in the length of the genealogy in genes with little or no recombination.

Second, we have used the divergence between species as an estimate of the mutation rate, but if the mutation rate at a locus changes through time, for which there is evidence (Aguileta *et al.* 2006; Hodgkinson and Eyre-Walker 2011), then we will tend to overestimate the variation in *N*_{e}: this is most easily seen by assuming there is variation in the mutation rate, but no variation in *N*_{e}; if the mutation rate has changed through time, then the divergence will not be a perfect measure of the recent mutation rate and there will appear to be variation in *N*_{e}.

Third, we have assumed that synonymous mutations are neutral, but there is evidence of selection in humans (Iida and Akashi 2000) and other species (Duret 2002; Pond and Muse 2005); although it is clear that selection has acted upon synonymous mutations in the past in *D. melanogaster*, the evidence of selection currently acting is contradictory (Akashi 1996; McVean and Vieira 2001; Zeng and Charlesworth 2010) and biased gene conversion may be acting (Galtier *et al.* 2006; Zeng and Charlesworth 2010). Most of the other species we have analyzed have not been investigated in any detail. We need to consider two models. In the first model, let us assume that there is no variation in *N*_{e} but that there is variation in the strength of selection on synonymous codons. Such a model would generate apparent variation in *N*_{e} with the genes subject to the strongest selection apparently having the highest *N*_{e}, because negative selection affects divergence to a greater extent than polymorphism (Kimura 1983). However, this would lead to the regions of the genome with the lowest diversity apparently having the highest effective population size. This is clearly not the case. If we split *P*_{s} into two independent samples, using a hypergeometric distribution, then we find a positive correlation between our estimate of *N*_{e} and *P*_{s} (Table 3). In the second model, let us imagine that there is variation in *N*_{e} and variation in the strength of selection on codon usage bias, but that they are uncorrelated to each other. In this case selection on codon usage bias will tend to generate an overestimate of the variation in *N*_{e}: as *N*_{e} increases, selection becomes more effective, but this reduces the divergence more than the level of polymorphism, yielding a higher apparent effective population size. So genes in regions of high *N*_{e} will tend to have an exaggerated *N*_{e}. There is also another effect that needs to be considered. We have estimated the level of synonymous divergence using the method of Goldman and Yang (1994; Yang and Nielsen 1998), which assumes that codon bias is due to mutation bias; however, this method will tend to overestimate the synonymous substitution rate if codon bias is due to selection, because it will incorrectly infer that genes with high bias have a small number of synonymous sites and hence a relatively large number of substitutions (Bierne and Eyre-Walker 2003; Yang 2006). As a consequence, the divergence in high-biased genes will be overestimated, but at the same time the mutation rate will tend to be underestimated because of the action of selection. These two factors may cancel each other out.

Fourth, we have applied our method only to protein-coding sequences, so we are estimating the variation in the effective population size that applies to the proteome. There might be further variation in *N*_{e} in regions that are relatively devoid of protein-coding sequences, such as heterochromatin. Whether this is important depends on whether there are functional sequences within these regions. We have also considered genes only on the autosomes and occasionally the homogametic sex chromosome (14 loci in *H. sapiens*). We have not considered genes on the heterogametic sex chromosome, which often appear to have much lower effective population sizes. However, the heterogametic sex chromosome usually has very few genes (Graves 2006).

Fifth, in estimating the variation in *N*_{e} we have assumed that there is either free recombination or no recombination and the population size has been stationary. Variation in population size can generate variation in diversity between loci, which may for example be mistaken for the signature of genetic hitchhiking (Tajima 1989; Pluzhnikov *et al.* 2002). In principle we could take this into account by estimating a demographic model from the polymorphism data while simultaneously estimating the variation in *N*_{e}. This is difficult and is beyond the scope of the current work.

Finally, we have not taken into account ancestral polymorphism within our method. Ignoring ancestral polymorphism will lead us to underestimate the variation in *N*_{e} because loci with large *N*_{e} will tend to have higher divergences than loci with small *N*_{e} and this will appear as though these loci have higher mutation rates; variation in *N*_{e} will therefore be underestimated because the mutation rate has been overestimated. In principle it is possible to include ancestral polymorphism within the method, but we observe a lack of convergence, probably because the number of polymorphisms for each gene was so low. However, we have chosen data sets in which divergence is generally considerably larger than diversity; for example, we chose macaque as the outgroup to humans because variation in *N*_{e} does appear to generate variation in the divergence between human and chimpanzee (McVicker *et al.* 2009).

Despite finding variation in *N*_{e} in many of the species we tested, we find no consistent evidence that *N*_{e} is correlated to either the rate of recombination or the density of selected sites, the two factors that we would have expected variation in *N*_{e} to depend upon. This is probably in part due to the fact that we are using synonymous diversity; as such, our estimates of diversity are subject to considerable error. The lack of a strong correlation between recombination rate and *N*_{e} may also be due to the fact that the genetic maps in *A. thaliana* and mouse are relatively crude. Furthermore, for our mouse species we are using an F_{2} genetic linkage map constructed from intercrosses between *M. m. domesticus* and *M. m. castaneus* to infer recombination rates for *M. m. castaneus*. In humans it has previously been shown that diversity over divergence is correlated positively to recombination rate (Hellmann *et al.* 2005) and that *d*_{n}/*d*_{s} is correlated to gene density (Bullaughey *et al.* 2008). In contrast to Hellmann *et al.* (2005), we do not find a significant correlation between *N*_{e} and recombination rate, but they used long noncoding sequences to investigate diversity over divergence; their estimates were therefore subject to much less error than ours. It is surprising that there is a correlation between genic density but not conserved site density in humans. This might be due to the fact that there is approximately twice as much variation in genic density as in conserved site density (coefficient of variation: 0.79 *vs.* 0.30). It might also be due to differences in the DFE between the two types of sites: background selection is most effective when the strength of selection acting upon deleterious mutations is similar in magnitude to the rate of recombination (Nordborg *et al.* 1996).

In contrast, genetic hitchhiking depends upon the rate of advantageous mutation and sequences undergoing considerable adaptive evolution may not appear as conserved. The correlation between *N*_{e} and the density of genic sites may therefore suggest that hitchhiking is more important in generating variation in *N*_{e} than background selection. The lack of a correlation between *N*_{e} and the density of selected sites in *Drosophila*, once correlations to the rate of recombination have been taken into account, may reflect the fact that the variation in *N*_{e} is generated by genetic hitchhiking and a lot of adaptive evolution goes on outside coding sequences (Andolfatto 2005).

Across species we find evidence that variation in *N*_{e} leads to variation in the effectiveness of natural selection on nonsynonymous mutations across the genome (Table 5). However, this is individually significant for just two genomes: *C. grandiflora* and *A. lyrata*. A lack of a correlation in other genomes may be due to the fact that we have little power to detect the correlation since (i) some of the data sets are quite small, (ii) there is limited variation in *N*_{e}, and (iii) in most of these species the DFE is very leptokurtic. The kurtosis of the DFE is such that changes in effective population size do not greatly change the proportion of mutations that are effectively neutral. It can be shown that under a gamma DFE the proportion of effectively neutral mutations is proportional to (Ohta 1977; Kimura 1979, 1983; Welch *et al.* 2008). Since β-values are typically between 0.1 and 0.3 in most species (Table S2), changes in *N*_{e} tend to cause small changes in the proportion of effectively neutral mutations; for example, a 10-fold increase in effective population size will reduce the proportion of effectively neutral mutations by only 37% if β = 0.2. We find no evidence of a significant negative correlation between ψ and either θ_{S} or *N*_{e} in humans, in agreement with the work of Bullaughey *et al.* (2008). They found no evidence that the ratio of the nonsynonymous (*d*_{N}) to the synonymous (*d*_{S}) substitution rate between human, chimpanzee, and macaque was correlated to the rate of recombination.

We find evidence that the amount of variation in *N*_{e} varies between species; however, there are no obvious correlates of this variation. Both plants and animals have species with high and low levels of variation. Surprisingly we find no obvious effect of self-fertilization as suggested by previous studies (Cutter and Payseur 2003; Roselius *et al.* 2005). *A. thaliana*, *C. rubella*, and *Boechera stricta* are all self-fertile with selfing rates of ∼0.95, 1, and 0.94, respectively (Charlesworth and Vekemans 2005; Song *et al.* 2006; Foxe *et al.* 2009), whereas the closely related species *A. lyrata* and *C. grandiflora* are obligate outcrossing species. However, the variation in *N*_{e} seems to be relatively low for *C. grandiflora* and *B. stricta* and similar for the two *Arabidopsis* species. It also should be noted that the confidence intervals on the estimate of *N*_{e} in *C. rubella* are very large and a substantial amount of variation is still shared between *C. grandiflora* and *C. rubella* so these estimates are not independent. Moreover, the lack of an effect for self-compatibility in our estimates of *N*_{e} for *Arabidopsis* may be not surprising as self-compatibility might evolved relatively recently in *Arabidopsis* (Bechsgaard *et al.* 2006; Tang *et al.* 2007). Furthermore, both *Arabidopsis* species have high sequence diversity in pericentromeric regions (Borevitz *et al.* 2007; Kawabe *et al.* 2008) that is not caused by varying mutation rates. Therefore this could be a major determinant of variation in *N*_{e} in those species and interfere with the effects of the breeding system.

Although the variation we observe in the effective population size appears to be modest, it does appear to influence both the level of neutral genetic diversity and the effectiveness of selection. This potentially has important implications. If slightly deleterious mutations contribute substantially to phenotypic traits, then variation in the effective population size may affect where the genetic variation underlying fitness and other traits is distributed. For example, Rockman *et al.* (2010) have recently shown that expression QTL (eQTL) tend to be present in regions of the *Caenorhabditis elegans* genome with the highest rates of recombination and lowest density of genes, where *N*_{e} is expected to be largest. However, population genetic theory also suggests that such weakly selected mutations are unlikely to contribute much to the overall genetic variance in fitness unless the proportion of mutations under such weak selection is large (Eyre-Walker 2010). Variation in the effective population might also affect the rate of adaptive evolution, as appears to be the case in *Drosophila* (Betancourt and Presgraves 2002). Advantageous mutations can potentially come from three sources. They can be generated *de novo*, in which case we expect regions of the genome with large *N*_{e} to adapt faster because the number of chromosomes an advantageous mutation can occur in is larger, and selection will be more effective on a greater proportion of the advantageous mutations. Advantageous mutations can also arise from standing genetic variation (Pritchard and Rienzo 2010; Pritchard *et al.* 2010). If these mutations were previously strongly deleterious, the genetic variation is not expected to depend upon *N*_{e}, unless the mutations are highly recessive. If, however, the advantageous mutations were previously neutral or weakly selected, regions of the genome with high *N*_{e} are expected to have more genetic variation and hence adapt more rapidly.

## Acknowledgments

We are grateful to several anonymous referees for comments. T.I.G. was financially supported by the John Maynard Smith studentship.

## Footnotes

*Communicating editor: M. A. Beaumont*

- Received July 12, 2011.
- Accepted September 23, 2011.

- Copyright © 2011 by the Genetics Society of America