Human genomes carry hundreds of mutations that are predicted to be deleterious in some environments, potentially affecting the health or fitness of an individual. We characterize the distribution of deleterious mutations among diverse human populations, modeled under different selection coefficients and dominance parameters. Using a new dataset of diverse human genomes from seven different populations, we use spatially explicit simulations to reveal that classes of deleterious alleles have very different patterns across populations, reflecting the interaction between genetic drift and purifying selection. We show that there is a strong signal of purifying selection at conserved genomic positions within African populations, but most predicted deleterious mutations have evolved as if they were neutral during the expansion out of Africa.

The Out-of-Africa (OOA) dispersal ∼50,000 y ago is characterized by a series of founder events as modern humans expanded into multiple continents. Population genetics theory predicts an increase of mutational load in populations undergoing serial founder effects during range expansions. To test this hypothesis, we have sequenced full genomes and high-coverage exomes from seven geographically divergent human populations from Namibia, Congo, Algeria, Pakistan, Cambodia, Siberia, and Mexico. We find that individual genomes vary modestly in the overall number of predicted deleterious alleles. We show via spatially explicit simulations that the observed distribution of deleterious allele frequencies is consistent with the OOA dispersal, particularly under a model where deleterious mutations are recessive. We conclude that there is a strong signal of purifying selection at conserved genomic positions within Africa, but that many predicted deleterious mutations have evolved as if they were neutral during the expansion out of Africa. Under a model where selection is inversely related to dominance, we show that OOA populations are likely to have a higher mutation load due to increased allele frequencies of nearly neutral variants that are recessive or partially recessive.

It has long been recognized that a human genome may carry many strongly deleterious mutations; Morton et al. (1) estimated that each human carries on average four or five mutations that would have a “conspicuous effect on fitness” if expressed in a homozygous state. Empirically estimating the deleterious mutation burden is now feasible through next-generation sequencing (NGS) technology, which can assay the complete breadth of variants in a human genome. For example, recent sequencing of over 6,000 exomes revealed that nearly half of all surveyed individuals carried a likely pathogenic allele in a known Mendelian disease gene (i.e., from a disease panel used for newborn screening) (2). Although there is some variation across individuals in the number of deleterious alleles per genome, we still do not know whether there are significant differences in deleterious variation among populations. Human populations vary dramatically in their levels of neutral genetic diversity, which suggests variation in the effective population size, N e . Theory suggests that the efficacy of natural selection is reduced in populations with lower N e because they experience greater genetic drift (3, 4). In an idealized population of constant size, the efficacy of purifying selection depends on the relationship between N e and the selection coefficient s against deleterious mutations. If 4N e s << 1, deleterious alleles evolve as if they were neutral and can, thus, reach appreciable frequencies. This theory raises the question of whether human populations carry differential burdens of deleterious alleles due to differences in demographic history.

Several recent papers have tested for differences in the burden of deleterious alleles among populations; these papers have focused on primarily comparing populations of western European and western African ancestry. Despite similar genomic datasets, these papers have reached a variety of contradictory conclusions (4⇓⇓⇓⇓–9). Initially, Lohmueller et al. (10) found that a panel of European Americans carried proportionally more derived, deleterious alleles than a panel of African Americans, potentially as the result of the Out-of-Africa (OOA) bottleneck. More recently, analyses using NGS exome datasets from samples of analogous continental ancestry found small or no differences in the average number of deleterious alleles per genome between African Americans and European Americans—depending on which prediction algorithm was used (11⇓–13). Simulations by Fu et al. (11) found strong bottlenecks with recovery could recapitulate patterns of differences in the number of deleterious alleles between African and non-African populations, supporting Lohmueller et al. (10), but in contrast to work by Simons et al. (12).

It is important to note two facts about these contradictory observations. First, these papers tend to use different statistics, which differ in power to detect changes across populations, as well as the impact of recent demographic history (6, 11). Lohmueller et al. (10) compared the relative number of nonsynonymous to synonymous (or “probably damaging” to “benign”) SNPs per population in a sample of n chromosomes, whereas Simons et al. (12) examined the special case of n = 2 chromosomes, namely, the average number of predicted deleterious alleles per genome (i.e., heterozygous + 2 * homozygous derived variants per genome). One way to think about these statistics is that the total number of variants, S, gives equal weight, w = 1, to an SNP regardless of its frequency, p. The average number of deleterious variants statistic gives weights proportional to the expected heterozygous and homozygous frequencies or w = 2p(1 − p) + p2 = 2p − p2. The average number of deleterious alleles per genome is fairly insensitive to differences in demographic history because heterozygosity is biased toward common variants. In contrast, the proportion of deleterious alleles has greater power to detect the impact of recent demographic history for large n across the populations because it is sensitive to rare variants that tend to be more numerous, younger, and enriched for functionally important mutations (14⇓–16). Second, empirical comparisons between two populations have focused primarily on an additive model for deleterious mutations, even though there is evidence for pathogenic mutations exhibiting a recessive or dominant effect (17, 18), and possibly an inverse relationship between the strength of selection s and the dominance parameter h (19).

There remains substantial conceptual and empirical uncertainty surrounding the processes that shape the distribution of deleterious variation across human populations. We aim here to clarify three aspects underlying this controversy: (i) Are there empirical differences in the total number of deleterious alleles among multiple human populations? (ii) Which model of dominance is appropriate for deleterious alleles (i.e., should zygosity be considered in load calculations)? (iii) Are the observed patterns consistent with predictions from models of range expansions accompanied by founder effects? We address these questions with a new genomic dataset of seven globally distributed human populations.

Results

Population History and Global Patterns of Genetic Diversity. We obtained moderate coverage whole-genome sequence (median depth 7×) and high coverage exome sequence data (median depth 78×) from individuals from seven populations from the Human Genome Diversity Panel (HGDP) (20). Unrelated individuals (no relationship closer than first cousin) were selected from seven populations chosen to represent the spectrum of human genetic variation from throughout Africa and the OOA expansion, including individuals from the Namibian San, Mbuti Pygmy (Democratic Republic of Congo), Algerian Mozabite, Pakistani Pathan, Cambodian, Siberian Yakut, and Mexican Mayan populations (Fig. 1A). The 2.48-Gb full genome callset consisted of 14,776,723 single nucleotide autosomal variants, for which we could orient 97% to ancestral/derived allele status (SI Appendix). Fig. 1. Decrease in heterozygosity and estimated N e with distance from southern Africa. (A) Locations of HGDP populations sampled for genome and exome sequencing are indicated on the map. Putative migration paths after the origin of modern humans are indicated with arrows (adapted from ref. 46). (B) PSMC curves for individual genomes, corrected for differences in coverage. Whereas populations experiencing an OOA bottleneck have substantially reduced N e , African populations also display a reduction in N e between ∼100 kya and 30 kya (see SI Appendix for simulations of population history with resulting PSMC curves). (C) For each individual’s exome, the number of putatively deleterious variants (equivalent to number of heterozygotes + twice the number of derived homozygotes) is shown by population. Heterozygosity among the seven populations decreases with distance from southern Africa, consistent with an expansion of humans from that region (21). The Namibian San population carried the highest number of derived heterozygotes, ∼2.39 million per sample, followed closely by the Mbuti Pygmies (SI Appendix, Table S1 and Fig. S5). The North African Mozabites carry more heterozygotes than the OOA populations in our dataset (2 million) but substantially fewer than the sub-Saharan samples, likely reflecting a complex history of an OOA migration, followed by reentry into North Africa and subsequent recent gene flow with neighboring African populations (22). The Maya have the lowest median number of heterozygotes in our sample, ∼1.5 million, which may be inflated due to recent European admixture (23). Two Mayan individuals displayed substantial recent European admixture (>20%) as assessed with local ancestry assignment (24) (SI Appendix, Fig. S6); these individuals were removed from analyses of deleterious variants. When we recalculated heterozygosity in the Maya, it was reduced by 3.5%. The decline in heterozygosity in OOA populations with distance from Africa strongly supports earlier results based on SNP array and microsatellite data for a serial founder effect model for the OOA dispersal (25, 26). We analyzed population history for individuals having sufficient coverage from five of the studied populations using the pairwise sequential Markovian coalescent software (PSMC) to estimate changes in N e (11, 12, 27). Because dating demographic events with PSMC is dependent on both the assumed mutation rate and the precision with which a given event can be inferred, we compare relative bottleneck magnitudes and timing among the seven HGDP populations. Consistent with previous analyses (27), the OOA populations show a sharp reduction in N e , with virtually identical population histories (Fig. 1B and SI Appendix). Simulations indicate that the magnitude of the 12-fold bottleneck is accurately estimated (SI Appendix, Fig. S7), even if the time of the presumed bottleneck is difficult to estimate precisely using PSMC. Interestingly, both the Mbuti and the Namibian San show a moderate reduction in N e relative to the ancestral maximum, with the San experiencing an almost twofold reduction in N e and the Mbuti displaying a reduction intermediate between the San and OOA populations (see also refs. 20, 28, and 29). These patterns are consistent with multiple population histories (e.g., both short and long bottlenecks) and multiple demographic events, including a reduction in substructure from the ancestral human population rather than a bottleneck per se (27).

Differences in Deleterious Alleles per Individual Genomes. Owing to differences in coverage among the whole genome sequences, our subsequent analyses focus on the high-coverage exome dataset (78× median coverage) to minimize any bias in comparing populations (Materials and Methods). We classified all mutations discovered in the exome dataset into categories based on Genomic Evolutionary Rate Profiling (GERP) Rejected Substitution (RS) scores. These conservation scores reflect various levels of constraint within a mammalian phylogeny (Materials and Methods) and are used to categorize mutations by their predicted deleterious effect (30, 31). Importantly, the allele present in the human reference genome was not used in the GERP RS calculation, avoiding the reference-bias effect previously observed in other algorithms (11, 12) (SI Appendix, Fig. S8A). Variants were sorted into four groups reflecting the likely severity of mutational effects: “neutral” (−2 < GERP < 2), “moderate” (2 ≤ GERP < 4), “large” (4 ≤ GERP < 6), and “extreme” (GERP ≥ 6) (SI Appendix, Fig. S9). GERP categories were concordant with ANNOVAR functional annotations (SI Appendix, Table S2 and Fig. S8B). When considering the total number of derived alleles per individual, defined here as A I = (1 × HET) + (2 × HOM der ), we observe an increase of predicted deleterious alleles with distance from Africa (Fig. 1C). The number of predicted deleterious alleles per individual increases along the range expansion axis (from San to Maya), consistent with theoretical predictions for expansion load (32). The maximal difference in the number of deleterious alleles between African and OOA individuals is ∼150 alleles. This result is consistent with theoretical predictions; the rate at which deleterious mutations accumulate in wave-front populations is limited by the total number of mutations occurring during the expansion (32). Assuming an exomic mutation rate of u = 0.5 per haploid exome and an expansion that lasted for t = 1,000 generations, a very conservative upper limit for the excess of deleterious alleles in OOA individuals would be 2*u*t = 1,000. The cline in A I is most pronounced for large-effect alleles (4 ≤ GERP < 6, Fig. 2E), whereby the San individuals carry A I = 4,450 large-effect alleles on average, increasing gradually to 4,550 in Yakut. The Mayans carry slightly fewer large-effect mutations per individual than the Yakut, which may be influenced by the residual European ancestry (between 5–20%) in our sample. For extreme alleles (GERP ≥ 6), each individual in the dataset carries on average 110–120 predicted highly deleterious alleles with no significant differences among populations (Fig. 2F). The average additive GERP score—obtained by counting the GERP scores at homozygous sites twice—for all predicted deleterious variants per individual is lowest in the San (∼3.3) and highest in the Maya (∼3.8). Fig. 2. Individual counts of deleterious variants. (A–C) For each individual’s exome, the number of derived homozygotes is plotted by population for moderate-, large-, and extreme-effect GERP categories. (D–F) For each individual’s exome, the number of derived variants (equivalent to number of heterozygotes + twice the number of homozygotes) is plotted by population for moderate-, large-, and extreme-effect GERP categories. Similar patterns are found when we consider the number of derived homozygous sites per individual. We find that individuals from OOA populations exhibit significantly more homozygotes for moderate, large, and extreme variants than African populations (Fig. 2 A–C). In addition, we observe a clear increase in the number of derived homozygotes with distance from Africa for moderate (2 ≤ GERP < 4) and large (4 ≤ GERP < 6) mutation effects categories, whereas the number of derived “extreme” homozygotes (GERP ≥ 6) is similar among OOA populations: All OOA genomes possess 30–40 extremely deleterious alleles in homozygous state (Fig. 2C). These patterns are in excellent agreement with theoretical predictions for the evolution of genetic variation during range expansions (7). The average GERP score per individual for derived homozygous variants is less differentiated than the additive model (above), varying between 2.43–2.49. It is important to note that A I is strongly influenced by common variants. Goode et al. (33) observed that as much as 90% of deleterious alleles in a single genome have a derived allele frequency greater than 5%, suggesting that the bulk of mutational burden using this metric will come from common variants. To explore this idea, we randomly chose an individual in each population and calculated the proportion of deleterious variants that are rare (<10%, i.e., a singleton within our population samples) and common (>10%), for each GERP category (Fig. 3A). Common deleterious alleles contribute to more than 90% of an individual’s A I , and the proportion of common deleterious variants increases with distance from Africa, as can be seen by the decrease of rare deleterious variants. This includes common large-effect variants, which make up proportionally more of A I for an OOA individual than for an African individual. For example, in a Mayan individual, 93% of large-effect variants are common compared with a San individual, where only 85% of large-effect variants are common (SI Appendix, Fig. S12). Given the small number of chromosomes in each population (n = 14–16), estimates of allele frequencies are subject to sampling effects. We recently performed the same analysis on exome data from the 1000 Genome Phase 1 Project (34). We find a similar pattern as in our HGDP data: On a per-genome basis, common variants represent a majority of the alleles predicted to be deleterious (5). Fig. 3. Differences in the proportion of deleterious alleles by frequency class. (A) The proportion of rare versus common deleterious variants per individual. For a given individual, deleterious variants were divided into common (>10%, solid colors) and rare (<10%, white space). The contribution of common deleterious variants to an individual’s burden is much greater than rare variants. (B) For each population, we calculated the proportional site frequency spectrum by plotting the proportion of deleterious large-effect alleles in each frequency class (translucent coloring) along with the proportion of neutral alleles for each frequency class (opaque coloring). African populations have proportionally fewer rare deleterious alleles than expected from neutrality. Populations with OOA ancestry have proportionally more fixed deleterious mutations.

Differences in Deleterious Alleles at the Population Level. To further elucidate the relationship between predicted mutation effect and allele frequencies, we compared the site frequency spectrum (SFS) for neutral and large- (4 ≤ GERP < 6) effect variants (Fig. 3B; see SI Appendix, Fig. S14 for a comparison between neutral and extreme variants). For all populations, singletons are enriched for deleterious variants (compared with neutral variants), consistent with the effect of purifying selection against deleterious variants (15, 35). However, the SFSs of OOA and African populations show marked differences. The neutral and deleterious SFSs of OOA populations show a global shift toward higher frequencies, consistent with the effects of serial bottlenecks/founder effects. It follows that OOA populations have fewer rare deleterious variants than Africans, as well as a larger proportion of fixed deleterious alleles; almost 7.9% of large-effect variants are fixed in the Maya, whereas the San have only 1.8% of deleterious variants fixed (Fig. 3B).

Simulations of Purifying Selection Under a Range Expansion. We sought to interpret the population-specific patterns of genetic diversity for each GERP category under a model including serial founder effects across geographic space and purifying selection. We simulated the evolution of both neutral and deleterious mutations under a simple model of range expansion in a 2D habitat (SI Appendix, Fig. S21). At selected loci, the ancestral allele was assumed selectively neutral and mutants reduced an individual’s fitness by a factor 1 − s only if it was present in homozygous state, that is, deleterious mutations were assumed to be completely recessive. Three thousand generations (corresponding to about 75 kya) after the onset of the range expansion, we computed the average expected heterozygosity for all populations. Computational limitations of individual-based simulations prohibit a complete exploration of the parameter space for this model, but, by varying migration rates and selection coefficients, we identified parameter values that fit the observed clines in heterozygosity reasonably well (Fig. 4B). Specifically, we first identified selection coefficients that yield the same relative differences between observed neutral and selected heterozygosities (Fig. 4A). Then, the migration rate was adjusted to fit the observed clines in heterozygosities, assuming that the distance between two demes is 250 km (Fig. 4B). The fit selection coefficients were 0, 1.25 × 10−4, 1 × 10−3, and 2 × 10−3 for neutral, moderate, large, and extreme GERP scores categories, respectively; the GERP ≥ 6 category showed the worst fit and observed counts indicate that even stronger selection coefficients should be considered for these extreme mutations (16). We performed the same analysis using a model in which mutations are codominant and, as expected, we found that the fit selection coefficients are smaller than those obtained a recessive model. These coefficients are estimated as s = 0, 0.5 × 10−4, 1.2 × 10−4, and 2 × 10−4, respectively (SI Appendix, Fig. S16) (16). Fig. 4. Heterozygosity under range expansion simulations with different selection coefficients. (A) Observed and simulated patterns of the reduction of heterozygosity (RH). Selection coefficients used in the simulations are s = 0 (black), s = −0.000125 (lavender), s = −0.001 (red), and s = −0.002 (orange). (B) Colored circles show average expected heterozygosity for populations with ancestry from the OOA bottleneck. Solid lines show the regression lines obtained from simulations and dashed lines indicate 95% confidence intervals for the regression. The boxplots and colored circles on the left show the simulated heterozygosities in ancestral (i.e., African) populations, and the observed heterozygosity in our African dataset (San/Mbuti), respectively. (C) Comparison of the distribution of RH between African and non-African individuals for different GERP categories, tested with a two-tailed Student t test (SI Appendix, Fig. S15).

Evolutionary Forces Acting on Heterozygosity. To better understand which evolutionary forces have acted in different populations to shape their levels of genetic diversity, we define a new statistic, RH. RH measures the reduction in heterozygosity at conserved sites relative to neutral heterozygosity, RH = (H neu − H del )/H neu , where H neu indicates heterozygosity at neutral sites and H del at GERP score categories >2. RH can be seen as a way to quantify changes of functional diversity across populations relative to neutral expectations. For instance, a constant RH value across populations would suggests that average functional diversity is determined by the same evolutionary force(s) as neutral diversity, that is, genetic drift and migration. In contrast, if RH changes across populations, it suggests that different evolutionary forces have shaped neutral and functional diversity, that is, selection has changed functional allele frequencies. In our dataset, RH is significantly larger in sub-Saharan Africans than in OOA populations across all functional GERP categories (Fig. 4C), indicating that selection has acted differently relative to drift between the two groups. The correlation between RH value and predicted mutation effect observed in Africa (Fig. 4A) confirms that purifying selection has kept strongly deleterious alleles at lower frequencies than in OOA populations. We then asked whether there were significant differences across OOA population, as oriented by their distance from eastern Africa. Interestingly, we see that the OOA RH values do not depend on their distance from Africa for predicted moderate-effect alleles (P = 0.82; SI Appendix, Fig. S15), suggesting that the frequencies of moderate mutations have evolved mainly according to neutral demographic processes during the range expansion out of Africa. In contrast, for strongly deleterious variants (large and extreme GERP categories) we see a significant cline in RH (P = 0.01 and P = 1.12 × 10–6, respectively; SI Appendix, Fig. S15), which implies that purifying selection has also contributed to their evolution relative to demographic processes.

Models of Dominance. We next considered whether there is empirical evidence for nonadditive effects for deleterious variants. Prior studies generally calculated “mutation load” by assuming an additive model, summing the number of deleterious alleles per individual, without factoring in whether an SNP occurs in a homozygous or heterozygous state. Determining an individual’s mutation load is, however, highly dependent on the underlying model of dominance (36) (a formal definition of mutation load is given below). For humans, Mendelian diseases tend to be overrepresented in endogamous populations or consanguineous pairings, indicating that many of these mutations are recessive (37); Gao et al. (38) estimate 0.58 lethal recessive mutations per diploid genome in the Hutterite population. Gene conversion can also lead to differential burden of derived, recessive diseases alleles among populations (39). Even height, a largely quantitative trait, seems to be affected by the architecture of recessive homozygous alleles in different populations (40). To further clarify the impact of dominance, we compared the distribution of deleterious variants across genes associated with dominant or recessive disease as reported in Online Mendelian Inheritance in Man (OMIM) (41). We expect to see a lower proportion of large- and extreme-effect variants in genes with dominant OMIM mutation annotations, compared with genes with recessive OMIM mutation annotations. We tested this hypothesis with the HGDP as well as the much larger 1000 Genomes Phase 1 dataset (SI Appendix, Fig. S18B). We averaged the proportion of variants within each effect category and performed a Wilcoxon test to determine whether the distribution of the proportion of large-effect variants was different between dominant and recessive genes. In the HGDP dataset, we observed P = 0.06, and for the larger 1000 Genomes dataset, P = 0.03. Our results indeed show a significantly higher proportion of large-effect variants in genes with recessive annotations, compared with genes with dominant annotations, suggesting that deleterious variants in the genome may tend to be recessive. However, we caution that OMIM genes are here annotated as dominant or recessive, whereas dominance is a property of specific mutations, and therefore all deleterious variants in a gene will not necessarily have the same dominance coefficient. Nonetheless, our results are consistent with an interpretation that genes may have certain properties, for example negative selection against dominant mutations in crucial housekeeping or developmental genes, that influence the tolerable distribution of dominance among variants. We consider the effect of dominance (summarized by h, which measures the effect of selected mutations in heterozygotes relative to homozygotes) on mutation load in the HGDP population samples given the observed differences in heterozygosity.