Rapid evolution of the human mutation spectrum

  1. Kelley Harris  Is a corresponding author
  2. Jonathan K Pritchard  Is a corresponding author
  1. Stanford University, United States
  2. Howard Hughes Medical Institute, Stanford University, United States

Abstract

DNA is a remarkably precise medium for copying and storing biological information. This high fidelity results from the action of hundreds of genes involved in replication, proofreading, and damage repair. Evolutionary theory suggests that in such a system, selection has limited ability to remove genetic variants that change mutation rates by small amounts or in specific sequence contexts. Consistent with this, using SNV variation as a proxy for mutational input, we report here that mutational spectra differ substantially among species, human continental groups and even some closely related populations. Close examination of one signal, an increased TCCTTC mutation rate in Europeans, indicates a burst of mutations from about 15,000 to 2000 years ago, perhaps due to the appearance, drift, and ultimate elimination of a genetic modifier of mutation rate. Our results suggest that mutation rates can evolve markedly over short evolutionary timescales and suggest the possibility of mapping mutational modifiers.

https://doi.org/10.7554/eLife.24284.001

eLife digest

DNA is a molecule that contains the information needed to build an organism. This information is stored as a code made up of four chemicals: adenine (A), guanine (G), cytosine (C), and thymine (T). Every time a cell divides and copies its DNA, it accidentally introduces ‘typos’ into the code, known as mutations. Most mutations are harmless, but some can cause damage. All cells have ways to proofread DNA, and the more resources are devoted to proofreading, the less mutations occur. Simple organisms such as bacteria use less energy to reduce mutations, because their genomes may tolerate more damage. More complex organisms, from yeast to humans, instead need to proofread their genomes more thoroughly.

Recent research has shown that humans have a lower mutation rate than chimpanzees and gorillas, their closest living relatives. Humans and other apes copy and proofread their DNA with basically the same biological machinery as yeast, which is about a billion years old. Yet, humans and apes have only existed for a small fraction of this time, a few million years. Why then do humans need to replicate and proofread their DNA differently from apes, and could it be that the way mutations arise is still evolving?

Previous research revealed that European people experience more mutations within certain DNA motifs (specifically, the DNA sequences ‘TCC’, ‘TCT’, ‘CCC’ and ‘ACC’) than Africans or East Asians do.

Now, Harris (who conducted the previous research) and Pritchard have compared how various human ethnic groups accumulate mutations and how these processes differ in different groups.

Statistical analysis of the genomes of thousands of people from all over the world did indeed show that the mutation rates of many different three-letter DNA motifs have changed during the past 20,000 years of human evolution. Harris and Pritchard report that when groups of humans left Africa and settled in isolated populations across different continents, each population quickly became better at avoiding mutations in some genomic contexts, but worse in others. This suggests that the risk of passing on harmful mutations to future generations is changing and evolving at an even faster rate than was originally suspected.

The results suggest that every human ethnic group carries specific variants of the genes which ensure that DNA replication and repair are accurate. These differences appear to influence which types of mutations are frequently passed down to future generations. An important next step will be to identify the genetic variants that could be controlling mutational patterns and how they affect human health.

https://doi.org/10.7554/eLife.24284.002

Introduction

Germline mutations not only provide the raw material for evolution but also generate genetic load and inherited disease. Indeed, the vast majority of mutations that affect fitness are deleterious, and hence biological systems have evolved elaborate mechanisms for accurate DNA replication and repair of diverse types of spontaneous damage. Due to the combined action of hundreds of genes, mutation rates are extremely low–in humans, about one point mutation per 100 MB or about 60 genome-wide per generation (Kong et al., 2012; Ségurel et al., 2014).

While the precise roles of most of the relevant genes have not been fully elucidated, research on somatic mutations in cancer has shown that defects in particular genes can lead to increased mutation rates within very specific sequence contexts (Alexandrov et al., 2013; Helleday et al., 2014). For example, mutations in the proofreading exonuclease domain of DNA polymerase ϵ cause TCTTAT and TCGTTG mutations on the leading DNA strand (Shinbrot et al., 2014). Mutational shifts of this kind have been referred to as ‘mutational signatures’. Specific signatures may also be caused by nongenetic factors such as chemical mutagens, UV damage, or guanine oxidation (Ohno et al., 2014).

Together, these observations imply a high degree of specialization of individual genes involved in DNA proofreading and repair. While the repair system has evolved to be extremely accurate overall, theory suggests that in such a system, natural selection may have limited ability to fine-tune the efficacy of individual genes (Lynch, 2011; Sung et al., 2012). If a variant in a repair gene increases or decreases the overall mutation rate by a small amount–for example, only in a very specific sequence context–then the net effect on fitness may fall below the threshold at which natural selection is effective. (Drift tends to dominate selection when the change in fitness is less than the inverse of effective population size). The limits of selection on mutation rate modifiers are especially acute in recombining organisms such as humans because a variant that increases the mutation rate can recombine away from deleterious mutations it generates elsewhere in the genome.

Given these theoretical predictions, we hypothesized that there may be substantial scope for modifiers of mutation rates to segregate within human populations, or between closely related species. Most triplet sequence contexts have mutation rates that vary across the evolutionary tree of mammals (Hwang and Green, 2004), but evolution of the mutation spectrum over short time scales has been less well described. Weak natural mutators have recently been observed in yeast (Bui et al., 2017) and inferred from human haplotype data (Seoighe and Scally, 2017); if such mutators affect specific pathways of proofreading or repair, then we may expect shifts in the abundance of mutations within particular sequence contexts. Indeed, one of us has recently identified a candidate signal of this type, namely an increase in TCCTTC transitions in Europeans, relative to other populations (Harris, 2015); this was recently replicated (Mathieson and Reich, 2016). Here, we show that mutation spectrum change is much more widespread than these initial studies suggested: although the TCCTTC rate increase in Europeans was unusually dramatic, smaller scale changes are so commonplace that almost every great ape species and human continental group has its own distinctive mutational spectrum.

Results

To investigate the mutational processes in different human populations, we classified each single nucleotide variants (SNV) in the 1000 Genomes Phase 3 data (Auton et al., 2015) in terms of its ancestral allele, derived allele, and 5’ and 3’ flanking nucleotides. We collapsed strand complements together to obtain 96 SNV categories. Since the detection of singletons may vary across samples, and because some singletons may result from cell-line or somatic mutations, we only considered variants seen in more than one copy. We further excluded variants in annotated repeats (since read mapping error rates may be higher in such regions) and in PhyloP conserved regions (to avoid selectively constrained regions) (Pollard et al., 2010). From the remaining sites, we calculated the distribution of derived SNVs carried by each Phase 3 individual. We used this as a proxy for the mutational input spectrum in the ancestors of each individual.

To explore global patterns of the mutation spectrum, we performed principal component analysis (PCA) in which each individual was characterized simply by the fraction of their derived alleles in each of the 96 SNV categories (Figure 1A). PCA is commonly applied to individual-level genotypes, in which case the PCs are usually highly correlated with geography (Novembre et al., 2008). Although the triplet mutation spectrum is an extremely compressed summary statistic compared to typical genotype arrays, we found that it contains sufficient information to reliably classify individuals by continent of origin. The first principal component separated Africans from non-Africans, and the second separated Europeans from East Asians, with South Asians and admixed native Americans (Figure 1—figure supplement 2) appearing intermediate between the two.

Figure 1 with 8 supplements
Global patterns of variation in SNV spectra.

(A) Principal component analysis of individuals according to the fraction of derived alleles that each individual carries in each of 96 mutational types. (B) Heatmaps showing, for pairs of …

https://doi.org/10.7554/eLife.24284.003
Figure 1—source data 1

This text file shows the number of SNPs in each of the 96 mutational categories that passed all filters in each 1000 Genomes continental group.

https://doi.org/10.7554/eLife.24284.004
Figure 1—figure supplement 8
Variation of the mutation spectrum with DNA replication timing.

We partitioned the genome into 10 equal replication timing quantiles using data obtained from (Woodfine et al., 2004), then computed mutation spectrum differences within each quantile. Although most …

https://doi.org/10.7554/eLife.24284.012
Figure 1—figure supplement 7
Mutation spectra of ChromHMM chromatin states (Part II of II).
https://doi.org/10.7554/eLife.24284.011
Figure 1—figure supplement 6
Mutation spectra of ChromHMM chromatin states (Part I of II).

To investigate whether any mutation spectrum shifts might be confined to particular chromatin states, we used chromHMM annotations of the human embryonic stem cell line HESC-H1 (Hoffman et al., 2013)…

https://doi.org/10.7554/eLife.24284.010
Figure 1—figure supplement 5
Mutation spectra of transcribed vs non-transcribed DNA.

Using the UCSC Genome Browser annotations of the human reference hg19, we determined whether each SNP occurs in a transcribed or non transcribed region. We further divided SNPs occurring in …

https://doi.org/10.7554/eLife.24284.009
Figure 1—figure supplement 4
The effects of biased gene conversion on mutation spectra.

When using segregating variation to study the mutation spectrum, one potential source of bias is that strong-to-weak mutations, where the ancestral allele is G or C and the derived allele is A or T, …

https://doi.org/10.7554/eLife.24284.008
Figure 1—figure supplement 3
Mutation spectrum comparison p-values.

Each left-hand plot shows all chi-squared p-values corresponding to the ratios from Figure 1A. In the absence of recent mutation spectrum evolution, only one out of 96 SNP categories is expected to …

https://doi.org/10.7554/eLife.24284.007
Figure 1—figure supplement 2
PCA of all 1000 Genomes continental groups.

All admixed North and South American individuals were omitted from Figure 1 in the main text to clarify the separation of other populations along an African vs non-African axis and an East vs West …

https://doi.org/10.7554/eLife.24284.006
Figure 1—figure supplement 1
Pairwise mutation spectrum comparisons among continental groups.

Each of these plots compares the mutation spectra of two populations P1 and P2. Letting fi denote the fraction of SNVs in population Pi that have a given triplet context, ancestral allele, and …

https://doi.org/10.7554/eLife.24284.005

Remarkably, we found that the mutation spectrum differences among continental groups are composed of small shifts in the abundance of many different mutation types (Figure 1B). For example, comparing Africans and Europeans, 43 of the 96 mutation types are significant at a p<105 threshold using a forward variable selection procedure. The previously described TCCTTC signature partially drives the difference between Europeans and the other groups, but most other shifts are smaller in magnitude and appear to be spread over more diffuse sets of related mutation types. East Asians have excess AT transversions in most sequence contexts, as well as about 10% more *AC*CC mutations than any other group. Compared to Africans, all Eurasians have proportionally fewer C* mutations relative to A* mutations.

Replication of mutation spectrum shifts

One possible concern is that batch effects or other sequencing artifacts might contribute to differences in mutational spectra. Therefore we replicated our analysis using 201 genomes from the Simons Genome Diversity Project (Mallick et al., 2016). The SGDP genomes were sequenced at high coverage, independently from 1000 Genomes, using an almost non-overlapping panel of samples. We found extremely strong agreement between the mutational shifts in the two data sets (Figure 2). For example, all of the 43 mutation types with a significant difference between Africa and Europe (at p<105) in 1000 Genomes also show a frequency difference in the same direction in SGDP (comparing Africa and West Eurasia). In both 1000 Genomes and SGDP, the enrichment of *AC*CC mutations in East Asia is larger in magnitude than any other signal aside from the previously described TCCTTC imbalance.

Figure 2 with 2 supplements
Concordance of mutational shifts in 1000 Genomes versus SGDP.

Each panel shows natural-log mutation spectrum ratios between a pair of continental groups, based on 1000 Genomes (x-axis) and SGDP (y-axis) data. Data points encoded by (+) symbols denote mutation …

https://doi.org/10.7554/eLife.24284.013
Figure 2—figure supplement 2
Regression of the SGDP heatmap coefficients versus the corresponding 1000 Genomes heatmap coefficients.
https://doi.org/10.7554/eLife.24284.015
Figure 2—figure supplement 1
Heatmap comparisons between continental groups in 1000 Genomes and the SGDP.

Here, each 1000 Genomes population is projected down to the sample size of the corresponding SGDP population in order to sample alleles with a similar distribution of ages and frequencies.

https://doi.org/10.7554/eLife.24284.014

The greatest discrepancies between 1000 Genomes and SGDP involve transversions at CpG sites, which are among the rarest mutational classes. These discrepancies might result from data processing differences or random sampling variation, but might also reflect differences in the fine-scale ethnic composition of the two panels.

Evidence for a pulse of TCCTTC mutations in Europe and South Asia

To investigate the timescale over which the mutation spectrum change occurred, we analyzed the allele frequency distribution of TCCTTC mutations, which are highly enriched in Europeans (Figure 3A; p<1×10300 for Europe vs. Africa) and to a lesser extent in South Asians. We calculated allele frequencies both in 1000 Genomes and in the larger UK10K genome panel (Walter et al., 2015). As expected for a signal that is primarily European, we found particular enrichment of these mutations at low frequencies. But surprisingly, the enrichment peaks around 0.6% frequency in UK10K, and there is practically no enrichment among the very lowest frequency variants (Figure 3B and Figure 3—figure supplement 1). CT mutations on other backgrounds, namely within TCT, CCC and ACC contexts, are also enriched in Europe and South Asia and show a similar enrichment around 0.6% frequency that declines among rarer variants (Figure 3C). This suggests that these four mutation types comprise the signature of a single mutational pulse that is no longer active. No other mutation types show such a pulse-like distribution in UK10K, although several types show evidence of monotonic rate change over time (Figure 3—figure supplements 3, 4 and 5).

Figure 3 with 6 supplements
Geographic distribution and age of the TCC mutation pulse.

(A) Observed frequencies of TCCTTC variants in 1000 Genomes populations. (B) Fraction of TCCTTC variants as a function of allele frequency in different samples indicates that these peak around …

https://doi.org/10.7554/eLife.24284.016
Figure 3—figure supplement 6
Mutation spectrum comparisons partitioned by allele frequency.

Each of these heatmaps shows a subset of the data used to construct Figure 1B, partitioned by allele frequency to show how rare variants are the most highly differentiated between populations. Black …

https://doi.org/10.7554/eLife.24284.022
Figure 3—figure supplement 5
Mutation type enrichment as a function of allele frequency in UK10K (Part III of III).

The 11 panels in this three-part figure show the full dependence of mutation spectrum on allele frequency in the UK10K data.

https://doi.org/10.7554/eLife.24284.021
Figure 3—figure supplement 4
Mutation type enrichment as a function of allele frequency in UK10K (Part II of III).

The eleven panels in this three-part figure show the full dependence of mutation spectrum on allele frequency in the UK10K data.

https://doi.org/10.7554/eLife.24284.020
Figure 3—figure supplement 3
Mutation type enrichment as a function of allele frequency in UK10K (Part I of III).

The eleven panels in Figure 3—figure supplements 2, 3 and 4 show the full dependence of mutation spectrum on allele frequency in the UK10K data. If we let F(f,m) denote the fraction of SNVs of frequency …

https://doi.org/10.7554/eLife.24284.019
Figure 3—figure supplement 2
Fraction of TCCTTC mutations as a function of allele frequency in ExAC.

Lek et al. compiled data from 60,706 exomes to create the Exome Aggregation Consortium dataset, which enables the analysis of ultra-rare human variation (Lek et al., 2016). The overall fraction of …

https://doi.org/10.7554/eLife.24284.018
Figure 3—figure supplement 1
TCCTTC mutation fraction as a function of allele frequency in all 1000 Genomes populations.

To enable better comparison with the 1000 Genomes data, the UK10K SNPs have been downsampled to 200 individuals. The age distribution of alleles of a given frequency varies as a function of the …

https://doi.org/10.7554/eLife.24284.017

We used the enrichment of TCCTTC mutations as a function of allele frequency to estimate when this mutation pulse was active. Assuming a simple piecewise-constant model, we infer that the rate of TCCTTC mutations increased dramatically 15,000 years ago and decreased again 2000 years ago. This time-range is consistent with results showing this signal in a pair of prehistoric European samples from 7000 and 8000 years ago, respectively (Mathieson and Reich, 2016). We hypothesize that this mutation pulse may have been caused by a mutator allele that drifted up in frequency starting 15,000 years ago, but that is now rare or absent from present day populations.

Although low frequency allele calls often contain a higher proportion of base calling errors than higher frequency allele calls do, it is not plausible that base-calling errors could be responsible for the pulse we have described. In the UK10K data, a minor allele present at 0.6% frequency corresponds to a derived allele that is present in 23 out of 3854 sampled haplotypes and supported by 80 short reads on average (assuming 7x coverage per individual). When independently generated datasets of different sizes are projected down to the same sample size, the TCCTTC pulse spans the same range of allele frequencies in both datasets (Figure 3—figure supplements 1 and 2), which would not be the case if the shape of the curve were a function of low-frequency errors.

Fine-scale mutation spectrum variation in other populations

Encouraged by these results, we sought to find other signatures of recent mutation pulses. We generated heatmaps and PCA plots of mutation spectrum variation within each continental group, looking for fine-scale differences between closely related populations (Figure 4 and Figure 4—figure supplement 1 through 6). In some cases, mutational spectra differ even between very closely related populations. For example, the *AC*CC mutations with elevated rates in East Asia appear to be distributed heterogeneously within that group, with most of the load carried by a subset of the Japanese individuals. These individuals also have elevated rates of ACAAAA and TATTTT mutations (Figure 4A and Figure 4—figure supplement 4). This signature appears to be present in only a handful of Chinese individuals, and no Kinh or Dai individuals. As seen for the European TCC mutation, the enrichment of these mutation types peaks at low frequencies, that is, 1%. Given the availability of only 200 Japanese individuals in 1000 Genomes, it is hard to say whether the true peak is at a frequency much lower than 1%.

Figure 4 with 6 supplements
Mutational variation among east Asian populations.

(A) PCA of east Asian samples from 1000 Genomes, based on the relative proportions of each of the 96 mutational types. See Figure 4—figure supplement 2 through 6 for other finescale population PCAs. …

https://doi.org/10.7554/eLife.24284.023
Figure 4—source data 1

This text file shows the number of SNPs in each of the 96 mutational categories that passed all filters in each finescale 1000 Genomes population.

https://doi.org/10.7554/eLife.24284.024
Figure 4—figure supplement 6
PCA of European populations.

Population abbreviations are: IBS = Iberian Population in Spain; TSI = Toscani in Italia; GBR = British in England and Scotland; CEU = Utah Residents (CEPH) with Northern and Western Ancestry; FIN = …

https://doi.org/10.7554/eLife.24284.030
Figure 4—figure supplement 5
PCA of South Asian populations.

Population abbreviations are: ITU = Indian Telugu from the UK; GIH = Gujarati Indian from Houston, Texas; PJL = Punjabi from Lahore, Pakistan; BEB = Bengali from Bangladesh; STU = Sri Lankan Tamil …

https://doi.org/10.7554/eLife.24284.029
Figure 4—figure supplement 4
PCA of East Asian populations.

Population abbreviations are: CDX = Chinese Dai in Xishuangbanna, China; JPT = Japanese in Tokyo, Japan; CHB = Han Chinese in Bejing, China; CHS = Southern Han Chinese; KHV = Kinh in Ho Chi Minh …

https://doi.org/10.7554/eLife.24284.028
Figure 4—figure supplement 3
PCA of African populations.

Population abbreviations are: MSL = Mende in Sierra Leone; LWK = Luhya in Webuye, Kenya; YRI = Yoruba in Ibadan, Nigeria; GWD = Gambian in Western Divisions; ESN = Esan in Nigeria.

https://doi.org/10.7554/eLife.24284.027
Figure 4—figure supplement 2
PCA of American populations.

Population abbreviations are: CLM = Colombians from Medellin, Colombia; MXL = Mexican Ancestry from Los Angeles, USA; PUR = Puerto Ricans from Puerto Rico; PEL = Peruvians from Lima, Peru; ACB = …

https://doi.org/10.7554/eLife.24284.026
Figure 4—figure supplement 1
Mutation spectrum differences within Africa, Europe, East Asia, and South Asia.

Figure 4B of the main text shows heat map comparisons between East Asian populations, which display fine-scale differences that are exceptionally well defined. For completeness, this figure shows …

https://doi.org/10.7554/eLife.24284.025

PCA reveals relatively little fine-scale structure within the mutational spectra of Europeans or South Asians (Figure 4—figure supplements 5 and 6). However, Africans exhibit some substructure (Figure 4—figure supplement 3), with the Luhya exhibiting the most distinctive mutational spectrum. Unexpectedly, a closer examination of PC loadings reveals that the Luhya outliers are enriched for the same mutational signature identified in the Japanese. Even in Europeans and South Asians, the first PC is heavily weighted toward *AC*CC, ACAAAA, and TATTTT, although this signature explains less of the mutation spectrum variance within these more homogeneous populations. The sharing of this signature may suggest either parallel increases of a shared mutation modifier, or a shared aspect of environment or life history that affects the mutation spectrum.

Mutation spectrum variation among the great apes

Finally, given our finding of extensive fine-scale variation in mutational spectra between human populations, we hypothesized that mutational variation between species is likely to be even greater. To compare the mutation spectra of the great apes in more detail, we obtained SNV data from the Great Ape Diversity Panel, which includes 78 whole genome sequences from six great ape species including human (Prado-Martinez et al., 2013). Overall, we find dramatic variation in mutational spectra among the great ape species (Figure 5 and Figure 5—figure supplement 1).

Figure 5 with 1 supplement
Mutational differences among the great apes.

(A) Relative abundance of SNV types in 5 ape species compared to Bornean Orangutan; data from (Prado-Martinez et al., 2013). Boxes indicate labels in (B). For additional comparisons see Figure …

https://doi.org/10.7554/eLife.24284.031
Figure 5—figure supplement 1
Mutation spectra of great apes.

These heatmap comparisons demonstrate that closely related great apes such as Chimpanzees and Bonobos have more similar mutation spectra than more distantly related apes do.

https://doi.org/10.7554/eLife.24284.032

As noted previously (Moorjani et al., 2016a), one major trend is a higher proportion of CpG mutations among the species closest to human, possibly reflecting lengthening generation time along the human lineage, consistent with previous indications that species closely related to humans have lower mutation rates than more distant species (Goodman, 1961; Li and Tanimura, 1987; Scally and Durbin, 2012). However, most other differences are not obviously related to known processes such as biased gene conversion and generation time change. The AT mutation rate appears to have sped up in the common ancestor of humans, chimpanzees, and bonobos, a change that appears consistent with a mutator variant that was fixed instead of lost. It is unclear whether this ancient AT speedup is related to the AT speedup in East Asians. Other mutational signatures appear on only a single branch of the great ape tree, such as a slowdown of AC mutations in gorillas.

Discussion

The widespread differences captured in Figures 1 and 2 may be footprints of allele frequency shifts affecting different mutator alleles. But in principle, other genetic and non-genetic processes may also impact the observed mutational spectrum. First, biased gene conversion (BGC) tends to favor C/G alleles over A/T, and BGC is potentially more efficient in populations of large effective size compared to populations of smaller effective size (Galtier et al., 2001). However, despite the bottlenecks that are known to have affected Eurasian diversity, there is no clear trend of an increased fraction of C/GA/T relative to A/TC/G in non-Africans vs Africans, or with distance from Africa (Figure 1—figure supplement 7), and previous studies have also found little evidence for a strong genome-wide effect of BGC on the mutational spectrum in humans and great apes (Do et al., 2015; Moorjani et al., 2016a). For these reasons, we think that evolution of the mutational process is a better explanation than BGC or selection for differences that have been observed between the spectra of ultra-rare singleton variants and older human genetic variation (Carlson et al., 2017);

It is also known that shifts in generation time or other life-history traits may affect mutational spectra, particularly for CpG transitions (Martin and Palumbi, 1993; Amster and Sella, 2016). Most CpG transitions result from spontaneous methyl-cytosine deamination as opposed to errors in DNA replication. Hence the rate of CpG transitions is less affected by generation time than other mutations (Hwang and Green, 2004; Moorjani et al., 2016b; Gao et al., 2016). We observe that Europeans have a lower fraction of CpG variants compared to Africans, East Asians and South Asians (Figure 1B), consistent with a recent report of a lower rate of de novo CCGCTG mutations in European individuals compared to Pakistanis (Narasimhan et al., 2016). Such a pattern may be consistent with a shorter average generation time in Europeans (Moorjani et al., 2016b), although it is unclear that a plausible shift in generation-time could produce such a large effect. Apart from this, the other patterns evident in Figure 1 do not seem explainable by known processes.

In summary, we report here that, mutational spectra differ significantly among closely related human populations, and that they differ greatly among the great ape species. Our work shows that subtle, concerted shifts in the frequencies of many different mutation types are more widespread than dramatic jumps in the rate of single mutation types, although the existence of the European TCCTTC pulse shows that both modes of evolution do occur (Harris, 2015; Moorjani et al., 2016b; Mathieson and Reich, 2016).

At this time, we cannot exclude a role for nongenetic factors such as changes in life history or mutagen exposure in driving these signals. However, given the sheer diversity of the effects reported here, it seems parsimonious to us to propose that most of this variation is driven by the appearance and drift of genetic modifiers of mutation rate. This situation is perhaps reminiscent of the earlier observation that genome-wide recombination patterns are variable among individuals (Coop et al., 2008), and ultimate discovery of PRDM9 (Baudat et al., 2010); although in this case it is unlikely that a single gene is responsible for all signals seen here. As large datasets of de novo mutations become available, it should be possible to map mutator loci genome-wide. In summary, our results suggest the likelihood that mutational modifiers are an important part of the landscape of human genetic variation.

Materials and methods

Data availability

Request a detailed protocol

All datasets analyzed here are publicly available at the following websites:

Human mutation spectrum processing

Request a detailed protocol

Mutation spectra were computed using 1000 Genomes Phase 3 SNPs (Auton et al., 2015) that are biallelic, pass all 1000 Genomes quality filters, and are not adjacent to any N’s in the hg19 reference sequence. Ancestral states were assigned using the UCSC Genome Browser alignment of hg19 to the PanTro2 chimpanzee reference genome; SNPs were discarded if neither the reference nor alternate allele matched the chimpanzee reference. To minimize the potential impact of ancestral misidentification errors, SNPs with derived allele frequency higher than 0.98 were discarded. We also filtered out regions annotated as ‘conserved’ based on the 100-way PhyloP conservation score (Pollard et al., 2010), download from http://hgdownload.cse.ucsc.edu/goldenPath/hg19/phastCons100way/, as well as regions annotated as repeats by RepeatMasker (Smit et al., 2013), downloaded from http://hgdownload.cse.ucsc.edu/goldenpath/hg19/database/nestedRepeats.txt.gz. To be counted as part of the mutation spectrum of population P (which can be either a continental group or a finer-scale population from one city), a SNP should not be a singleton within population P–at least two copies of the ancestral and derived alleles must be present within that population.

An identical approach was used to extract the mutation spectrum of the UK10K ALSPAC panel (Walter et al., 2015), which is not subdivided into smaller populations. The data were filtered as described in Field et al. (2016). The filtering procedure performed by Field et al. (2016) reduces the ALSPAC sample size to 1927 individuals.

We also computed mutation spectra of the Simons Genome Diversity Panel (SGDP) populations (Mallick et al., 2016). Four of the SGDP populations, West Eurasia, East Asia, South Asia, and Africa, were compared to their direct counterparts in the 1000 Genomes data. Three additional SGDP populations, Central Asia and Siberia, Oceania, and America, had no close 1000 Genomes counterparts and were not analyzed here (although each project contained a panel of people from the Americans, the composition of the American panels was extremely different, with the 1000 Genomes populations being much more admixed with Europeans and Africans). SGDP sites with more than 20% missing data were not utilized. All other data processing was done the same way described for the 1000 Genomes data.

The following table gives the same size of each population panel, as well as the total number of SNPs segregating in the panel that are used to compute mutation type ratios:

DatasetPopulationNumber of individualsNumber of SNPs
1 kgAfrica50416,870,400
1 kgEurope5038,508,040
1 kgEast Asia5047,895,925
1 kgSouth Asia4899,552,781
SGDPAfrica456,569,658
SGDPWest Eurasia694,201,571
SGDPEast Asia493,312,645
SGDPSouth Asia383,449,624

Great ape diversity panel data processing

Request a detailed protocol

Biallelic great ape SNPs were extracted from the Great Ape Diversity Panel VCF (Prado-Martinez et al., 2013), which is aligned to the hg18 human reference sequence. Ancestral states were assigned using the Great Ape Genetic Diversity project annotation, which used the Felsenstein pruning algorithm to assign allelic states to internal nodes in the great ape tree. In the Great Ape Diversity Panel, the most recent common ancestor (MRCA) of the human species is labeled as node 18; the MRCAs of chimpanzees, bonobos, gorillas, and orangutans, respectively, are labeled as node 16, node 17, node 19, and node 15. We extracted the state of each MRCA at each SNP in the alignment and used it to polarize the ancestral and derived allele at that site; a SNP was discarded whenever the ancestral node was assigned an uncertain or polymorphic ancestral state. As with the human data, SNPs with derived allele frequency higher than 0.98 were not used, and both repeats and PhyloP-annotated conserved regions were filtered away.

Visual representation of mutation spectra

Request a detailed protocol

The mutation type of an SNP is defined in terms of its ancestral allele, its derived allele, and its two immediate 5’ and 3’ neighbors. Two mutation types are considered equivalent if they are strand-complementary to each other (e.g. ACGATG is equivalent to CGTCAT). This scheme classifies SNPs into 96 different mutation types, each that can be represented with an A or C ancestral allele.

To compute the frequency fP(m) of SNP m in population P, we count up all SNPs of type m where the derived allele is present in at least one representative of population P (which can be either a specific population such as YRI or a broader continental group such as AFR). After obtaining this count CP(m), we define fP(m) to be the ratio CP(m)/mCP(m), where the sum in the denominator ranges over all 96 mutation types m. The enrichment of mutation type m in population P1 relative to population P2 is defined to be fP1(m)/fP2(m); these enrichments are visualized as heat maps in Figures 1B, 3B and 4A.

To track changes in the mutational spectrum over time, we compute fP(m) in bins of restricted allele frequency. This involves counting the number of SNPs of type m that are present at frequency ϕ in population P to obtain counts CP(m,ϕ) and frequencies fP(m,ϕ)=CP(m,ϕ)/mCP(mϕ). Deviation of the ratio fP(m,ϕ)/fP(m) from one indicates that the rate of m has fluctuated recently in the history of population P. To make the sampling noise approximately uniform across alleles of different frequencies, alleles of derived count greater than five were grouped into approximately log-spaced bins that each contained similar numbers of UK10K SNPs. More precisely, we defined a set of bin endpoints b1,b2, such that the total number of SNPs ranging in derived allele count between bi and bi+11 is greater than or equal to the number of 5-ton SNPs, while the total number of SNPs ranging in derived allele count from bi to bi+1-2 is less than the number of 5-ton SNPs.

In some cases, for example Figure 2, Figure 2—figure supplement 1B and Figure 3—figure supplement 1, site frequency spectra were projected down to a smaller sample size before counting SNPs in order to more accurately compare datasets of different sample sizes. A binomial sampling approach was used to project a sample of N haplotypes does to a smaller sample size n. Letting CP(N)(m,ϕ) denote the SNP counts in the large sample of N haplotypes, effective SNP counts CP(n)(m,ϕ) in a sample of n haplotypes are computed as follows:

CP(n)(m,k/n)=(nk)=1N-1(/N)k(1-/N)n-kCP(N)(m,/N)

Significance testing

Request a detailed protocol

One central goal of this paper is to test whether many mutation types differ in rate between human populations or whether mutation spectrum shifts have been rare events affecting only a small proportion of mutation types. A simple statistical method for answering this question would be to perform 96 separate chi-square tests, one for each triplet-context-dependent mutation type, as follows:

Let Si denote the total number of SNPs segregating in population Pi, and let Si(m) denote the number of SNPs of mutation type m. If mutation type m is more prevalent in population P1 than in population P2, a chi-square test provides a natural way of assessing the significance of this difference. As described in Harris (2015), this test is performed on the following two-by-two contingency table:

S1(m)P1-S1(m)
S2(m)P2-S2(m)

It would be appealing to conclude that every mutation type ‘passing’ this chi-square test is a mutation type that has changed in rate during recent human history. However, if we were to perform the full set of 96 tests, they would not be independent. A sufficiently large increase in the rate of one mutation type m1 in population P1 after divergence from P2 could cause another mutation type m2, whose rate has remained constant, to comprise significantly different fractions of the SNPs from P1 and P2. To minimize this effect, we formulate the following iterative procedure of conditionally independent tests: first, compute a chi-square significance value punordered(m) for each mutation type m using the two-by-two chi-square table above. We then use these values to order the SNPs from lowest p value to highest and compute a set of ordered p values pordered(m). For the mutation type m0 with the lowest unordered p value, punordered(m0)=pordered(m0). For mutation type mi, which has the ith lowest unordered p value and i<96, pordered(mi) is computed from the following contingency table:

S1(mi)j=i+196S1(mj)
S2(mi)j=i+196S2(mj)

For mutation type m96, which has the highest unordered p value, the ordered p value is computed from the contingency table

S1(m96)S1(m95)
S2(m96)S2(m95)

This procedure is guaranteed to find fewer mutation types to differ significantly in rate between populations compared to separate chi-square tests.

Principal component analysis

Request a detailed protocol

The python package matplotlib.mlab.PCA was used to perform PCA on the complete set of 1000 Genomes diploid genomes. First, the triplet mutational spectrum of each haplotype h was computed as a 96-element vector encoding the mutation frequencies (fh(m))m of the non-singleton derived alleles present on that haplotype. The mutational spectrum of each diploid genome was then computed by averaging together the spectra of its two constituent haplotypes. In the same way, a separate PCA was performed on each of the five continental groups to reveal finescale components of mutation spectrum variation.

Dating of the TCCT mutation pulse

Request a detailed protocol

We estimated the duration and intensity of TCCT rate acceleration in Europe by fitting a simple piecewise-constant rate model to the UK10K frequency data. To specify the parameters of the model, we divide time into discrete log-spaced intervals bounded by time points t1,,td, assigning each interval a TCCT mutation rate r0,rd. In units of generations before the present, the time discretization points were chosen to be: 20, 40, 200, 400, 800, 1200, 1600, 2000, 2400, 2800, 3200, 3600, 4000, 8000, 12,000, 16,000, 20,000, 24,000, 28,000, 32,000, 36,000, 40,000. We assume that the total rate r of mutations other than TCCT stays constant over time (a first-order approximation).

In terms of these rate variables, we can calculate the expected shape of the TCCT pulse shown in Figure 2B of the main text. The shape of this curve depends on both the mutation rate parameters ri and the demographic history of the European population, which determines the joint distribution of allele frequency and allele age. To account for the effects of demography, we use Hudson’s ms program to simulate 10,000 random coalescent trees under a realistic European demographic history inferred from allele frequency data (Tennessen et al., 2012) and condition our inference upon this collection of trees as follows:

Let A(m,t) be the function for which titi+1A(m,t)dt equals the coalescent tree branch length, averaged over the sample of simulated trees, that is ancestral to exactly m lineages and falls between time ti and ti+1. Given this function, which can be empirically estimated from a sample of simulated trees, the expected frequency spectrum entry k/n is

E(k/n)=i=1dti-1tiA(k,t)dtj=1ni=1dti-1tiA(j,t)dt

and the expected fraction of TCCT mutations in allele frequency bin k/n is

E(fTCCT(k/n))=i=1driti-1tiA(k,t)dtri=1dti-1tiA(k,t)dt.

The expected value of the TCCT enrichment ratio being plotted in Figure 2B is

E(rTCCT(k/n))=i=1driti-1tiA(k,t)dtj=1ni=1dti-1tiA(j,t)dti=1dti-1tiA(k,t)dtj=1ni=1driti-1tiA(j,t)dt

In Figure 2B, enrichment ratios are not computed for every allele frequency in isolation, but for allele frequency bins that each contain similar numbers of SNPs. Given integers 1km<km+1n, the expected TCCT enrichment ratio averaged over all SNPs with allele frequency between km/n and km+1/n is:

E(rTCCT(km/n))=i=1driti-1tik=kmkm+1A(k,t)dtj=1ni=1dti-1tiA(j,t)dti=1dti-1tik=kmkm+1A(k,t)dtj=1ni=1driti-1tiA(j,t)dt

We optimize the mutation rates r1,,rd using a log-spaced quantization of allele frequencies k1/n,,km/n defined such that all bins contain similar numbers of SNPs. The chosen allele count endpoints k1,,km are: 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 20, 30, 40, 50, 60, 70, 80, 90, 100, 200, 300, 400, 500, 600, 700, 800, 900, 1000, 2000, 3000, 4000. Given this quantization of allele frequencies, we optimize r1,,rd by using the BFGS algorithm to minimize the least squares distance D(r0,,rd) between E(rTCCT(km/n)) and the empirical ratio rTCCT(km/n) computed from the UK10K data. This optimization is subject to a regularization penalty that minimizes the jumps between adjacent mutation rates ri and ri+1:

D(r0,,rd)=m=1d(E(rTCCT(km/n))-rTCCT(km/n))2+0.25i=1d(ri-1-ri)2

Although the underlying model of mutation rate change assumed here is very simple, it still represents an advance over the method used in (Harris, 2015) to estimate of the timing of the TCCTTC mutation rate increase. That method relied upon explicit estimates of allele age from a dataset of less than 100 individuals, which are much noisier than integration of a joint distribution of allele age and frequency across a sample of thousands of haplotypes.

Data availability

The following previously published data sets were used
    1. Swapan Mallick
    2. David Reich
    3. et al
    (2016) Simons Genome Diversity Project
    Publicly available from the Simons Foundation. Directions for downloading available here: http://simonsfoundation.s3.amazonaws.com/share/SCDA/datasets/2014_11_12/StepstodownloadtheSGDPdataset_v4.docx.
    1. Prado-Martinez
    2. Tomas Marques-Bonet
    3. et al
    (2013) Whole genome sequences for a set of 79 great ape individuals. Genome sequencing
    Publicly available at NCBI BioProject, submitted as part of the Great Ape Genome Diversity Project (accession no: PRJNA189439).
    1. Monkol Lek
    2. Daniel MacArthur
    3. et al
    (2016) Exome Aggregation Consortium
    Summary data publicly available for download at http://exac.broadinstitute.org/downloads.
    1. Prado-Martinez
    2. Tomas Marques-Bonet
    3. et al
    (2013) Great Ape Genome Diversity Project
    Publicly available at NCBI Sequence Read Archive (accession no: SRP018689).

References

    1. Galtier N
    2. Piganeau G
    3. Mouchiroud D
    4. Duret L
    (2001)
    GC-content evolution in mammalian genomes: the biased gene conversion hypothesis
    Genetics 159:907–911.
    1. Goodman M
    (1961)
    The role of immunochemical differences in the phyletic development of human behavior
    Human Biology 33:131–162.

Article and author information

Author details

  1. Kelley Harris

    Department of Genetics, Stanford University, Stanford, United States
    Contribution
    KH, Conceptualization, Data curation, Software, Formal analysis, Funding acquisition, Validation, Investigation, Visualization, Methodology, Writing—original draft, Writing—review and editing
    For correspondence
    kelleyh@stanford.edu
    Competing interests
    The authors declare that no competing interests exist.
    ORCID icon "This ORCID iD identifies the author of this article:" 0000-0003-0302-2523
  2. Jonathan K Pritchard

    1. Department of Genetics, Stanford University, Stanford, United States
    2. Department of Biology, Stanford University, Stanford, United States
    3. Howard Hughes Medical Institute, Stanford University, Stanford, United States
    Contribution
    JKP, Conceptualization, Funding acquisition, Visualization, Writing—original draft, Writing—review and editing
    For correspondence
    pritch@stanford.edu
    Competing interests
    The authors declare that no competing interests exist.

Funding

National Institutes of Health (NRSA-F32 Grant GM116381)

  • Kelley Harris

Howard Hughes Medical Institute (Investigator Grant)

  • Jonathan K Pritchard

National Institutes of Health (R01 Grant HG008140)

  • Jonathan K Pritchard

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

Acknowledgements

This work was funded by NIH grants GM116381 and HG008140, and by the Howard Hughes Medical Institute. We thank Jeffrey Spence and Yun S Song for technical assistance. We also thank Ziyue Gao, Arbel Harpak, Molly Przeworski, Joshua Schraiber, and Aylwyn Scally for comments and discussion, as well as two anonymous reviewers.

Copyright

© 2017, Harris et al.

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.

Metrics

  • 20,433
    views
  • 1,787
    downloads
  • 159
    citations

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

Download links

Share this article

https://doi.org/10.7554/eLife.24284

Further reading

    1. Evolutionary Biology
    Mattias Siljestam, Claus Rueffler
    Research Article Updated

    The majority of highly polymorphic genes are related to immune functions and with over 100 alleles within a population, genes of the major histocompatibility complex (MHC) are the most polymorphic loci in vertebrates. How such extraordinary polymorphism arose and is maintained is controversial. One possibility is heterozygote advantage (HA), which can in principle maintain any number of alleles, but biologically explicit models based on this mechanism have so far failed to reliably predict the coexistence of significantly more than 10 alleles. We here present an eco-evolutionary model showing that evolution can result in the emergence and maintenance of more than 100 alleles under HA if the following two assumptions are fulfilled: first, pathogens are lethal in the absence of an appropriate immune defence; second, the effect of pathogens depends on host condition, with hosts in poorer condition being affected more strongly. Thus, our results show that HA can be a more potent force in explaining the extraordinary polymorphism found at MHC loci than currently recognised.

    1. Ecology
    2. Evolutionary Biology
    Rebecca D Tarvin, Jeffrey L Coleman ... Richard W Fitch
    Research Article

    Understanding the origins of novel, complex phenotypes is a major goal in evolutionary biology. Poison frogs of the family Dendrobatidae have evolved the novel ability to acquire alkaloids from their diet for chemical defense at least three times. However, taxon sampling for alkaloids has been biased towards colorful species, without similar attention paid to inconspicuous ones that are often assumed to be undefended. As a result, our understanding of how chemical defense evolved in this group is incomplete. Here, we provide new data showing that, in contrast to previous studies, species from each undefended poison frog clade have measurable yet low amounts of alkaloids. We confirm that undefended dendrobatids regularly consume mites and ants, which are known sources of alkaloids. Thus, our data suggest that diet is insufficient to explain the defended phenotype. Our data support the existence of a phenotypic intermediate between toxin consumption and sequestration — passive accumulation — that differs from sequestration in that it involves no derived forms of transport and storage mechanisms yet results in low levels of toxin accumulation. We discuss the concept of passive accumulation and its potential role in the origin of chemical defenses in poison frogs and other toxin-sequestering organisms. In light of ideas from pharmacokinetics, we incorporate new and old data from poison frogs into an evolutionary model that could help explain the origins of acquired chemical defenses in animals and provide insight into the molecular processes that govern the fate of ingested toxins.