ARTICLE https://doi.org/10.1038/s41467-020-16296-4 OPEN Extreme differences between human germline and tumor mutation densities are driven by ancestral human-specific deviations José María Heredia-Genestar 1, Tomàs Marquès-Bonet 1,2,3,4, David Juan 1✉ & Arcadi Navarro 1,2,3,5 ✉ 1234567890():,; Mutations do not accumulate uniformly across the genome. Human germline and tumor mutation density correlate poorly, and each is associated with different genomic features. Here, we use non-human great ape (NHGA) germlines to determine human germline- and tumor-specific deviations from an ancestral-like great ape genome-wide mutational land- scape. Strikingly, we find that the distribution of mutation densities in tumors presents a stronger correlation with NHGA than with human germlines. This effect is driven by human- specific differences in the distribution of mutations at non-CpG sites. We propose that ancestral human demographic events, together with the human-specific mutation slowdown, disrupted the human genome-wide distribution of mutation densities. Tumors partially recover this distribution by accumulating preneoplastic-like somatic mutations. Our results highlight the potential utility of using NHGA population data, rather than human controls, to establish the expected mutational background of healthy somatic cells. 1 Institute of Evolutionary Biology (CSIC-UPF), Department of Experimental and Health Sciences, Universitat Pompeu Fabra, 08003 Barcelona, Spain. 2 CRG- CNAG, Centre for Genomic Regulation (CRG), Barcelona Institute of Science and Technology (BIST), 08003 Barcelona, Spain. 3 Institució Catalana de Recerca i Estudis Avançats (ICREA), 08010 Barcelona, Spain. 4 Institut Català de Paleontologia Miquel Crusafont, Universitat Autònoma de Barcelona, 08193 Cerdanyola del Vallès, Barcelona, Spain. 5 Barcelonaβeta Brain Research Center (BBRC), Pasqual Maragall Foundation, 08005 Barcelona, Spain. ✉email: david.juan@upf.edu; arcadi.navarro@upf.edu NATURE COMMUNICATIONS | (2020)11:2512 | https://doi.org/10.1038/s41467-020-16296-4 | www.nature.com/naturecommunications 1 ARTICLE NATURE COMMUNICATIONS | https://doi.org/10.1038/s41467-020-16296-4 M utation density, at different scales, has been shown to Results correlate with different genomic features, such as Mutation density distribution. We compared the mutation regional GC-content or recombination rate1–5. In can- density distribution in human (1KGP29, SGDP50 (ref. 30)), non- cer, mutation density shows a different behavior than in human human great apes (NHGA: chimpanzee31,32, gorilla31,33), and germline and has been linked to chromatin states6, with higher human cancer34 data sets. We focused on high-quality orthologous mutation accumulation in closed chromatin. It has been sug- regions shared between human, chimpanzee and gorilla genomes, gested that the tumor mutation density distribution is mediated measuring the number of variants per 1 Mbp independently of the by the different replication timing of euchromatin and hetero- frequency of each variant (see Methods, Supplementary Note 1). chromatin regions. Closed-chromatin, late-replicating regions In agreement with previous reports1,3,4,6, we observe a variable have poorer accessibility or recruitment of the mismatch repair distribution of the mutation density across the genome in all data machinery, decreasing the efficiency of DNA damage repair in sets (Fig. 1a) and a weak correlation of mutation densities them7,8. The correlation between tumor mutation density and of the human germline and tumors1,6 (Table 1). Strikingly, the chromatin state is tissue-dependent, allowing the identification of NHGA–tumor correlations are much stronger than the the tissue-of-origin of metastatic tumor samples9,10. Recent stu- human–tumor correlation and are similar to the human–NHGA dies have detected mutation patterns in human healthy somatic germline correlations (Table 1 and Supplementary Tables 1, 2). tissues resembling those seen in tumors of that same tissue11,12. We compared the distribution of mutation density between These patterns suggest the presence of inherent differences pairs of data sets (Supplementary Fig. 1). Interestingly, we between the mutation landscapes of the human germline and observed that mutation density in tumors is higher in windows soma. where NHGAs have higher mutation density than humans At a smaller scale, sequence context is a good predictor of the (Fig. 1b, c). To control for differences in the shapes of mutation rate13, even beyond hypermutable CpG sites14–17. distributions, we ranked each set of windows according to their Sequence context has been widely used in cancer analyses to mutation density (Fig. 1d, e). These ranked distributions show a detect signatures of mutation associated with mutagens such as clear pattern: tumor mutation densities are higher in windows UV-light, tobacco smoke, or APOBEC activity18,19. These effects with higher ranks in NHGAs than in human (two-sided have also been detected in healthy somatic tissues11,12. Interest- Mann–Whitney U test p value human–chimpanzee = 3.7e−216; ingly, it has been shown that the human mutation spectrum can human–gorilla = 2.8e−161). This human-NHGA-tumor diago- be recapitulated by a combination of two signatures: SBS1, nal behavior is exclusive to human-NHGA comparisons, as it associated with CpG>T transitions and SBS5, associated with a cannot be observed when comparing chimpanzee to gorilla slight increase in many non-CpG mutations20,21. These two sig- (p value chimpanzee–gorilla = 0.28) (Supplementary Fig. 1), and natures are detected in virtually all cancer samples19. The sig- can be replicated under different conditions, in orangutan, and natures’ mutation load is correlated with age of diagnosis and when using other human data sets (Supplementary Notes 2, 3, their prevalence in most tumors suggests that many of the cancer Supplementary Tables 3–12, and Supplementary Fig. 1). mutations occur in preneoplastic stages19,21. De novo mutations are also affected by sequence context20,22,23. Subspecies’ diversity. One possible explanation for this obser- The rates of some particular mutation types have changed recently vation is the higher genetic diversity of chimpanzees and gorillas across ancestries24–27, and seem to have been under selection in relative to human populations. We explored this by analyzing the human lineage24. Sequence context studies have shown dif- individual NHGA subspecies with varying levels of diversity. ferences in the relative proportion of certain mutation types High-diversity NHGA subspecies have stronger correlations with between great ape species26. Furthermore, studies of de novo both human and tumor than the low-diversity subspecies (Sup- mutations in great ape samples revealed a slowdown of the overall plementary Note 3, Supplementary Table 5). In fact, high- mutation rate in humans relative to chimpanzees and gorillas28. diversity subspecies present NHGA-tumor correlations as strong These previous works suggest that the mutation density land- as when analyzing the whole species. This suggests that genetic scape in the human germline has been influenced by specific diversity is one of the drivers of the NHGA-tumor correlations. phenomena and that part of the mutations accumulating in Furthermore, the diagonal pattern is only characteristic of com- tumors might be driven by the same processes that occur in parisons between the germlines of humans and high-diversity healthy cells. However, little is known about the forces driving NHGA subspecies. A comparison of high and low-diversity mutation accumulation in the human germline and their impact chimpanzee and gorilla subspecies showed a clear horizontal split on the pronounced differences in mutation densities between (Supplementary Fig. 2). Mutation density in tumors co-localizes tumors and the human germline. Here, we use NHGA population with the most diverse NHGA subspecies, regardless of the data to understand the origin of these differences by analyzing mutation density in the least diverse. However, the diagonal their evolution through the human lineage. We analyzed muta- pattern could never be reconstructed by comparing two NHGA tion distribution (at the 1 Mbp scale) in human and NHGAs (our subspecies. In other words; while a lack of diversity distorts the closest relatives), and compared them with the distribution of distribution of the genome-wide mutation densities, the diagonal somatic mutations in tumors. This evolutionary analysis allows us pattern is caused by effects intrinsic to the human lineage. to study the conservation status of mutation density, as well as to Consistent with this idea, we observed a weak intermediate detect species-specific mutation hotspots and coldspots, thus pattern when comparing NHGA to three archaic hominid shedding light on the processes governing differential mutagen- genomes (Supplementary Note 3, Supplementary Fig. 1, Supple- esis across the genome. mentary Table 13). Although the effects of small sample sizes and We observe that the mutation distribution in NHGAs and in of ancient DNA damage should not be disregarded, these results tumors present striking similitudes. The differences with human suggest that at least part of the differentiation process in the mutation distribution are driven by changes in human’s past distribution of mutation densities was already established before demography affecting the distribution of CpG>T transitions and the human–Neanderthal split. mutations at non-CpG sites. We do not detect any mutational process driving the NHGA-tumor similitudes, suggesting that these similitudes are driven by the normal accumulation of Chromatin state and genomic features. Previous studies have mutations in healthy human cells. shown that tumor mutation distribution is associated with the 2 NATURE COMMUNICATIONS | (2020)11:2512 | https://doi.org/10.1038/s41467-020-16296-4 | www.nature.com/naturecommunications NATURE COMMUNICATIONS | https://doi.org/10.1038/s41467-020-16296-4 ARTICLE a b c 0.6 2.5 2 Tumor 1 kGP − 1.40x Chimpanzee 3 0.4 Density SGDP50 − 1.69x Gorilla 0.0 2 Chimpanzee − 1.55x 0 1 Gorilla − 2.11x 0 0.2 Tumor − 3.50x −2.5 −2 −1 0.0 −5 0 5 10 0 5 10 0 5 10 Standardized mutation density 1KGP 1KGP d e Most Most mutated mutated Tumor Most mutated Chimpanzee Gorilla Least mutated Least Least mutated mutated Least Most Least Most mutated 1KGP mutated mutated 1KGP mutated Fig. 1 Distribution of mutation density across data sets. a Distribution of the standardized mutation density in 1 Mbp windows in human, NHGA, and tumor data sets. Color denotes data set. The numbers next to the legend are the fold-enrichment between the 95th and 5th quantiles. b Distribution of the standardized mutation density in humans, chimpanzee, and tumor. Each point represents a 1 Mbp window. The x axis represents the human mutation density, the y axis the chimpanzee mutation density, and the point color, the tumor mutation density. The black line corresponds to the diagonal where the mutation density is equal in human and chimpanzee. c Same as (b) but comparing human and gorilla. d Distribution of the ranked mutation density in humans, chimpanzee, and tumor. Each point represents a 1 Mbp window. The x and y axis represent the ranking in mutation density in human and chimpanzee, respectively. Color represents the ranked mutation density in the tumor data set. The solid black line corresponds to the diagonal where the ranked mutation density is equal in human and chimpanzee. The dashed lines represent a 25% difference in ranking in both species. e Same as (d), comparing human and gorilla. Note 4, Supplementary Fig. 3). However, and strikingly, NHGAs Table 1 Correlation between distributions of mutation density. show strong positive correlations with epigenomic features associated with closed chromatin (derived from lymphoblastoid cell 1KGP Chimpanzee Gorilla lines), just as tumors do, but humans do not (Fig. 2a, Supplementary Chimpanzee 0.65 – – Table 16). We also observe consistent associations with human Gorilla 0.53 0.84 – chromatin states35 (Fig. 2b, c). GC-content, CpG-content, and Tumor 0.16 0.55 0.58 H3K36me1 (specifically recruited in the gene bodies of genes Pairwise Pearson’s correlation R of the standardized mutation density of 5040 1 Mbp windows regulated by CpG islands)36 show a clear positive correlation with between data sets. human but negative with NHGAs and tumors, suggesting that they might be contributing to the diagonal pattern (Figs. 2d, e, 3a, b). CpG sites. Intrigued by the connection of several CpG-related tumor’s tissue-of-origin chromatin conformation9. When analyz- features with the human-NHGA-tumor diagonal pattern that ing individual tumor types separately, we detect the human- implies a stronger correlation between mutation densities in NHGA-tumor differential pattern across all tumor types. tumors and NHGAs than with the human germline (Fig. 3a, b), Although the mutation density of tumors across different tissues is we analyzed separately CpG>T transitions and mutations at non- highly correlated, some high-sample size tumors (skin melanoma CpG sites. CpG>T transitions present very strong correlations and breast adenocarcinoma) present slightly lower correlations between all germline data sets and very poor correlations with (two-sided Mann–Whitney U test p value skin melanoma = tumor, driven by their strong association with CpG-content 1.8e−118; breast adenocarcinoma = 5e−86) (Supplementary (Fig. 3c, d). The relationship between CpG-content and mutation Note 3, Supplementary Tables 14, 15). density at non-CpG sites in NHGAs is more similar to tumors Interestingly, correlations between a variety of genomic features than to the human germline. Moreover, the correlations of non- and tumor mutation density are consistently more similar to the CpG mutation densities between human, NHGA, and tumors are correlations with NHGAs than with humans (Fig. 2a). Mutation similar to those observed using all sites (Fig. 3e, f, Supplementary densities in NHGAs have, like in humans, strong correlations with Tables 17–20). The correlation of genomic features with CpG>T sequence conservation and recombination rate (Supplementary transitions are similar between germline data sets, while the NATURE COMMUNICATIONS | (2020)11:2512 | https://doi.org/10.1038/s41467-020-16296-4 | www.nature.com/naturecommunications 3 ARTICLE NATURE COMMUNICATIONS | https://doi.org/10.1038/s41467-020-16296-4 a 0.8 Correlation with genomic features Tumor 1KGP SGDP50 0.4 Pearson's correlation R Chimpanzee Gorilla 0.0 −0.4 −0.8 H3K9me3 H3K9me2 H4K20me3 mappability H3K27me2 conservation H3K27me3 recombination rate genes replication time H2AZ HiC H3K14ac H3R2me1 nucleosomes H3K36me1 CpG GC−content H4K20me1 H3K27me1 PolII H3K79me3 H3K79me2 CTCF H2AK9ac H3K79me1 H3K9me1 H3K4me2 H4K16ac H4K91ac H3K4me1 H3K36me3 H4K12ac H3K23ac H3K4me3 H3K9ac H4K8ac H4K5ac H2BK5ac H3K27ac H2BK12ac H3K36ac H2AK5ac H3K18ac H2BK120ac H2BK20ac H3K4ac b Heterochromatin c Promoters d GC−content e H3K36me1 Most Most Most Most mutated mutated mutated mutated 1.00 0.20 1.2 Chimpanzee Chimpanzee Chimpanzee Chimpanzee 0.75 0.15 0.50 0.9 0.45 0.50 0.10 0.40 0.6 0.25 0.05 0.35 0.00 0.00 0.3 Least Least Least Least mutated mutated mutated mutated Least Most Least Most Least Most Least Most mutated 1KGP mutated mutated 1KGP mutated mutated 1KGP mutated mutated 1KGP mutated Fig. 2 Genomic Features. a Pearson’s correlation R of different data sets with human genomic features (data sources in Supplementary Table 16). Color denotes data set. b Human and chimpanzee ranked mutation density distribution. Color represents the overlap of heterochromatin in human lymphoblastoid cell lines (LCLs) measured by chromHMM states35, c the aggregate chromHMM states associated with the presence of promoters, d the window’s GC-content, and e the density of H3K36me1 histone mark ChIP-seq reads57. correlations of non-CpG sites are weaker in humans than in species correlation of CpG>T transitions and mutations at non- NHGA and tumors (Supplementary Note 5, Supplementary CpG sites is weaker in humans than in NHGAs and in tumors Fig. 2). Consequently, while tumor-specific differences in corre- (Supplementary Fig. 3, Supplementary Table 24). Therefore, the lations of mutation densities with genomic features are associated combination of human-specific dynamics in both, non-CpG and with CpG>T mutations, human germline-specific differences are CpG>T mutations, causes the diagonal pattern observed when related to the different behavior of human non-CpG mutations. comparing all SNVs. Despite the overall mutation slowdown previously described in human28,37, all data sets present similar proportions of recurrent Mutational signatures. Studies of cancer mutations have asso- CpG>T transitions, discarding a possible effect of CpG>T ciated certain mutational signatures with the action of specific saturation on these human-exclusive differences (Supplementary mutagens or cellular processes18. We benefited from these asso- Note 5, Supplementary Tables 21 and 22). Interestingly, we ciations to explore the contribution of different mutation observed that even a small number of de novo mutations38 can mechanisms to the observed differences in mutation densities by recover much higher correlations with all germline data sets for analyzing the trinucleotide context of SNV in all germline CpG>T than for non-CpG mutation densities (Supplementary data sets. Table 23). This observation suggests that CpG>T mutations At the whole-genome level, the triplet mutation spectra of require much less diversity than non-CpG mutations to recover human, chimpanzee, and gorilla are very similar (Supplementary its genome-wide distribution, probably due to their strong Fig. 4, Supplementary Table 25), although different from tumors. dependence of the skewed genome-wide distribution of CpG- It has been shown that the human mutation spectrum can be content. In contrast, the distribution of human non-CpG de novo recapitulated by a combination of the cancer signatures SBS1 and mutations replicates the general trend of all germline data sets SBS5 (refs. 20,21). We were able to replicate this association in showing intermediate but higher correlation with NHGAs than NHGA and in another primate species (Vervet monkey) with tumors (Supplementary Note 5, Supplementary Table 23), (Supplementary Note 6, Supplementary Tables 26 and 27). however, the small size of this data set makes difficult to discern Previous studies have shown that both, germline and somatic, whether their correlation with tumors is lower than expected mutation spectra of mice differ from the human mutation for NHGAs. spectra39. Our results show that mutation spectra and their When comparing the distribution of non-CpG mutations in specific association with signatures SBS1 and SBS5 are conserved human, NHGA, and tumors, we detect a horizontal pattern across the primate lineage. (Supplementary Fig. 3) similar to those observed in comparisons A subset of trinucleotides show significant species-specific of high- and low-diversity subspecies. This suggests that the enrichments (Chi-Squared test p value < 10−5). We detected no effects of loss of diversity in a given species can be detected mostly association between these trinucleotides and known mutation in non-CpG sites but not in CpG>T mutations. In fact, within- mechanisms (Fig. 4a, see Methods, Supplementary Note 6, 4 NATURE COMMUNICATIONS | (2020)11:2512 | https://doi.org/10.1038/s41467-020-16296-4 | www.nature.com/naturecommunications NATURE COMMUNICATIONS | https://doi.org/10.1038/s41467-020-16296-4 ARTICLE Most a mutated c CpG>T 0.03 CpG−content 0.02 0.01 Least Ranked mut. density Most CpG mutated mutated content ee nz Chimpanzee 0.05 P pa rill a d 1K G Ch im Go 0.04 Pearson’s R 0.03 Tumor −0.13 −0.14 −0.01 1 0.02 0.01 Gorilla 0.93 0.94 0 Chimpanzee 0.98 −1 e Non−CpG sites Least 0.025 CpG−content mutated Least Most 0.020 mutated 1KGP mutated 0.015 b All sites 0.010 Least Ranked mut. density Most mutated mutated 0.025 ee nz GP pa rill a f 1K Ch im Go CpG−content 1 kGP 0.020 SGDP50 Pearson’s R Tumor 0.38 0.66 0.66 1 Chimpanzee 0.015 Gorilla Tumor Gorilla 0.59 0.86 0 0.010 Least Most Chimpanzee 0.7 mutated mutated −1 Ranked mutation density Fig. 3 CpG-content. a Ranked mutation density in human and chimpanzee. Color represents the distribution of CpG-content (fraction of the window overlapping CpG sites) in the human reference hg19. b Loess smoothers of ranked mutation density and CpG-content for the different data sets. Color denotes data set. Shaded area represents the 95% CI. c CpG>T transitions corrected by the whole window size; loess smoothers same as in (b); d Pearson’s correlation R of the standardized mutation density of CpG>T transitions in different species; e same as in (b, c), but using only mutations at non-CpG sites; f Pearson’s correlation R of the standardized mutation density of mutations at non-CpG sites in different species. Supplementary Fig. 4, Supplementary Table 28). However, linear sample and with preneoplastic mutations in tumors19,40. This regression models show a positive and significant (p value < 10−4) suggests that the strong correlation between NHGA and tumor effect of the triplet’s GC-content and its fold-enrichment in the mutation densities is driven by the same mechanisms in healthy human–chimpanzee comparison (Supplementary Fig. 4). human cells and in the great ape lineage, while the genome-wide Although no specific known mutational mechanism could be distribution of mutations has been altered in the human germline. detected having a whole-genome effect in the mutation spectrum, other processes might have an effect on the mutation distribution. To test this, we analyzed the differential distribution of Discussion trinucleotide mutations in 1 Mbp windows across the genome in We analyzed the mutation density distribution at the 1 Mbp scale humans, NHGA, and tumors. Only trinucleotides with similar in the human and NHGA germlines, as well as in human tumors. enrichment between species (non-CpG, mainly C>G and T>C) We observed a moderate similitude between human and NHGA show differences in their distribution across the genome between germlines, and surprisingly a higher resemblance between human human, NHGA, and tumor (trinucleotide-difference test p value < tumors with the germlines of NHGAs than with those of humans. 10−5, see Methods, Supplementary Note 6, Fig. 4a). These discrepancies in mutation density in the human and We further analyzed this trinucleotide distribution across the NHGA germlines are differently associated with genomic and genome in individual tumor types. We compared the association epigenomic features. In fact, human germlines show a lower of the number of mutations caused by each cancer signature19 in correlation than NHGA germlines and tumors with most of these each individual tumor type to the human-NHGA-tumor pattern features, suggesting that they present a distorted mutational dis- (Supplementary Data 1). Signatures SBS5 and SBS40 show a tribution. Regions more densely mutated in humans than in significant association (signature-difference test p value < 10−4, NHGAs tend to be GC-rich, exon-rich, promoter and enhancer- see Methods) of the pattern with the tumor’s signature mutation rich, open chromatin and early replicating. Particularly, CpG- load (Fig. 4b). Both SBS5 and SBS40 are flat signatures whose related features show a positive correlation with human and a mutation load is associated with the age of diagnosis of the negative correlation with NHGA and tumor mutation densities. NATURE COMMUNICATIONS | (2020)11:2512 | https://doi.org/10.1038/s41467-020-16296-4 | www.nature.com/naturecommunications 5 ARTICLE NATURE COMMUNICATIONS | https://doi.org/10.1038/s41467-020-16296-4 a ATA>C b SBS1 Deamination of 5−methylcytosine 0.12 SBS2 APOBEC activity Tumor types TTA>C ATT>C with signature SBS3 Defective homologous recomb.–based repair ATG>C Mutation present in SBS4 Tobacco smoking 2+ samples TTT>C SBS5 Age of sample CCT>G C>A CCA>G 10 SBS6 Defective DNA mismatch repair difference between standard deviations TTG>C CTG>C CTA>C C>G CCC>G SBS9 Polymerase eta activity 20 0.08 ACA>T ACA>G CTT>C C>T SBS12 Unknown TCG>T Window−distribution test 30 GTG>C ACT>G ACG>T SBS13 APOBEC activity TCA>G GCA>G T>A SBS14 Polymerase epsilon mutation TTG>G TCC>G GCT>G Signature TTC>G T>C SBS16 Unknown #mutations TCC>T TCT>G ATC>C SBS17a Unknown GTT>C T>G 10M 0.04 SBS17b Unknown SBS18 Damage by reactive oxygen species 1M SBS19 Unknown TTA>A Significance window−distribution SBS22 Aristolochic acid exposure 100K CCG>G SBS26 test Defective DNA mismatch repair Chimpanzee ACC>A SBS28 Unknown 10K GTT>A Signif. 0.00 SBS29 Tobacco chewing 1 KGP GCG>G ACG>A Non.Sign. 1K SBS37 Unknown CCG>A GTC>A SBS40 Age of sample 0 Chimpanzee 1KGP GCG>A SBS41 Unknown ATC>A SBS44 Defective DNA mismatch repair −0.2 −0.1 0.0 0.1 0.2 0 1.3 3 5 log10 Fold enrichment 1 KGP vs chimpanzee −log10 p−value Fig. 4 Trinucleotide analysis. a Contribution to the higher chimpanzee-tumor mutation distribution similarity vs. genome-wide enrichment in human compared with chimpanzee. x axis: log10 of the enrichment of trinucleotides comparing human and chimpanzee. Left: enriched in chimpanzee; right: enriched in human. y axis: effect size (difference between the standard deviations of human–tumor and chimpanzee–tumor) of the trinucleotide-difference test (see Methods). Positive values: tumor distribution more similar to chimpanzee; negative values: tumor distribution more similar to human. Color represents the central nucleotide mutation type. Filled dots denote mutation types significant (p value < 10−5) in the trinucleotide-difference test (Kolmogorov–Smirnov test; multiple testing significance threshold <5 × 10−4). b Chi-squared -log10 p values of the association of each cancer signature mutation load to the trinucleotide-difference test (signature-difference test; see Methods). The dashed line represents the Bonferroni multiple correction significance threshold. Color represents the number of mutations associated with each signature in the whole data set. Dot size represents the number of tumor types with two or more samples showing the signature. Only non-artifact mutation signatures present in two or more tumor types are shown. The possible functional implications in human evolution require recovered as a whole. In contrast, the hypermutability of CpG further study. sites and its concentration in specific regions caused CpG>T These observations are driven by the different dynamics of transitions to have already recovered diversity levels similar to mutation density at CpG>T transitions (very similar in all those of high-diversity NHGAs. However, low-diversity NHGA germlines and very different in tumors) and at non-CpG sites subspecies show lower CpG>T and higher non-CpG degrees of (more similar in NHGAs and human tumors than in the human recovery of their mutation densities than the human germline. germline). This effect is exclusive of the human germline and, This suggests that some human-specific phenomena are influ- thus, must have been caused by human-specific conditions. encing the different behavior of human mutational landscapes. We observed that human and other primates present a very We pose that this human-specific behavior can be the result of the similar global triplet mutation spectrum. We detected an combination of at least two very different contributions. enrichment of certain trinucleotide mutations in humans and First, the recent human-exclusive population expansions31,42 NHGAs consistent with previous results (non-CpG, GC-rich are expected to cause an increase of clock-like CpG>T mutations mutations are enriched in humans)26. The enriched trinucleotides in the population37,43, leaving signatures akin to positive selec- are not associated with mutation signatures with known causes, tion, as it has been described in Native Americans25. The results nor do they contribute significantly to the higher similitude of presented in our study suggest that these signatures of apparent human tumors to NHGA germlines. This suggests the absence of positive selection may not be population exclusive and might be strong mechanistic changes biasing the accumulation of muta- pervasive through human evolution, especially in CpG-rich tions in any of the studied germlines. exons. Moreover, the highly skewed genome-wide distribution As previously described for human20,21, we observed that of CpG-content can contribute to the easy recuperation of mutation rates of three non-human primates are explained by mutation densities of the fast mutating CpG>T transitions (as mutation signatures SBS1 (mostly CpG>T transitions) and SBS5 observed in a reduced data set of de novo mutations)38. The (associated with normal accumulation of mutation in healthy decoupling of the CpG>T/non-CpG mutation rates within the somatic and germline cells)18,20,21,41. Moreover, the lower same region is stronger in humans than in NHGAs and tumors human-tumor than NHGA-tumor mutation distribution corre- and could imply that these population expansions introduce lation is driven by the accumulation of mutations associated with additional distortions on the non-CpG densities associated to signatures SBS5 and SBS40 (similar to SBS5 and recently dis- specific biases inherited from previously contracted populations. covered)19. These results suggest that the poor human-tumor We cannot disregard an additional contribution of human- correlation is caused by the fact that human (but not NHGAs) specific shifts in CpG>T transitions mutation rates, although they germline (and de novo mutations) do not currently reflect the have been suggested to be similar across all great apes37. expected mutation densities of healthy (and preneoplastic-like) Second, the recent slowdown in mutation rates reported in human somatic cells. humans28,37 seems to have affected differently CpG>T and non- We observed that the moderate human-NHGA and the low CpG sites. A relative deceleration of non-CpG mutations could human-tumor correlations of mutation densities at non-CpG also contribute to the lesser capability of non-CpG mutations to sites could be caused by losses in population diversity (as recover their expected genome-wide distribution. In addition, it observed in low-diversity NHGA subspecies). We propose that could contribute to the CpG>T/non-CpG differentiation, which is successive bottlenecks during human evolution removed a sub- observed more strongly in the human germline. Interestingly, our stantial part of nucleotide variation that still remains to be analyses with a data set of de novo mutations38 show an increase 6 NATURE COMMUNICATIONS | (2020)11:2512 | https://doi.org/10.1038/s41467-020-16296-4 | www.nature.com/naturecommunications NATURE COMMUNICATIONS | https://doi.org/10.1038/s41467-020-16296-4 ARTICLE in the contribution of CpG>T mutations compared with all the we filtered out any region of the human reference genome hg19 failing one of the other data sets. This suggests that the recent slowdown in humans following criteria: poor mappability of the human reference split into 35 bp k-mers, poor callability in ≥25% of the chimpanzee or gorilla samples, or, matching a might have been less pronounced in CpG>T mutations. However, known Copy Number Variable region in NHGA samples55 (Supplementary the small size of the de novo data sets makes it difficult to Note 1). 2052 Mbp of autosomal sequence passed this filtering (76.54% of the non- interpret the relation of their genome-wide distributions with the N base pairs in the human reference autosomes). We divided the autosomes into other data sets. Further studies of de novo mutations in both 1 Mbp overlapping (500 kb) windows and kept all windows where ≥50% of its bases humans and NHGAs are needed to elucidate possible differences passed our filtering. This left 5040 1 Mbp windows to analyze (Supplementary Fig. 1, Supplementary Table 3). in the evolution of mutation rates. We propose that the combi- These filters were applied to all data sets used, including both our callings and nation of population bottlenecks and expansions, together with external data sets used as released. All SNV counts, trinucleotide counts, and the specific nature of the different mutation types, drives the genomic features measurements through this study used only regions passing this differences observed in the distributions of human mutation filtering. For the analysis of archaic samples, we combined this filtering with the intersection of the callability mask of all three archaic samples. This specific densities. Future additional work will be necessary to establish the filtering was applied to all data sets when compared with the archaic samples. relevance of these and other phenomena in the observed human- specific behavior of the genome-wide distribution of mutation Mutation density. We used BASH scripts to measure mutation density of each densities. window in each data set by counting either the number of non-fixed segregating sites Our results imply that accumulated mutations in human (in the human, chimpanzee, and gorilla data sets) or the number of somatic mutations populations are a poor proxy of the expected mutational back- (in the tumor and human de novo data sets, accounting repeated mutations as independent mutational events). We divided this count of single nucleotide variants ground in healthy somatic cells. In fact, contrary to the common (SNV) by the fraction of the window passing our filtering. This results in a measure of assumption that extreme differences in the mutational landscape mutations per megabase pair (Mbp) of sequence for each window. We standardized of human germline and tumors are due to the abnormal behavior the resulting distribution within each data set deeming it as the mutation density. We of tumors, our results show that they are driven by human- ranked all windows within a data set by their distribution of mutation density to control for the different shapes of the data sets distributions. specific particularities in the accumulation of mutations in the germline. In fact, accumulated mutations in NHGAs (at least at Correlations between distributions. All correlations used in this analysis are non-CpG sites) happen to be more informative about the normal Pearson’s correlation coefficient (using the R function cor.test) between the stan- occurrence of mutations in healthy somatic cells than those dardized mutation densities (unranked) of the two data sets unless otherwise accumulated by humans. specified. Partial correlations, when used, were calculated using the pcor function Our results do not affect the well-established association from the ppcor R package. between the mutation’s observed population frequencies and their relative functional impact44. However, they suggest that higher Significance of the diagonal split. To measure the significance of the diagonal precision in the definition of damaging mutations could be split pattern observed when comparing the human and NHGA data sets, we divided all windows into two groups depending on if the ranked mutation density achieved by incorporating data from less biased populations. is higher in human than NHGAs or vice-versa. We calculated the two-sided Interestingly, a recent work exploring missense mutations in Mann–Whitney U test on the variable of interest (usually, the tumor mutation population data from NHGAs for the prediction of the impact of density) on both groups using the R function wilcox.test. human mutations has already given a first major step in that direction45. Genomic features. The genomic features used were filtered using the same These results have also paramount implications for the mappability, callability, and copy-number filters used for the mutation density data. The features used were either the overlap of the feature’s genomic coordinates understanding of tumors as often normally mutating primate cell with the fraction of the 1 Mbp window passing our filtering (e.g., GC-content, populations. This implies that the definition of positively selected CpG-content) or the average value or intensity of the feature in the passing fraction shifts in regional mutational burdens could be much refined by of the window (e.g., histone marks), depending on the original data (Supplemen- correcting mutation expectations by making use of primate tary Table 16). populations or even of tumoral cell populations. Moreover, dis- entangling the contribution of normal mutation processes in Trinucleotides. We classified each SNV into the 96 possible combinations of trinucleotides (12 different mutation types, by 16 combinations of the adjacent tumors (either if they correspond to preneoplastic cells or not) nucleotides, divided by two when folding them). We determined the adjacent opens the door to establishing them as models to understand reference sequence of each SNV using the getfasta option of bedtools56. We filtered somatic mutations in healthy cells, with particular interest in the out any variant where the liftOver tool46 could not map them to the chimpanzee experimental research of aging and somatic mosaicism. panTro5 or the gorilla gorGor5 reference genomes, or the trinucleotide sequence differed in one of the three reference genomes (accounting for strand). This filter was applied to all windows and we used for our analysis only windows where ≥50% Methods of it passed both the original high-quality orthologous regions and this three- Data sets used. For the human data sets we used the release variant calling of 2504 reference filter, leaving 4920 windows to use. We applied additional filters humans from the 1000 Genomes Project29 (1KGP), our own calling of 50 addi- requiring the trinucleotide to be species exclusive and to not overlap variants in tional human samples from the Simons Genome Diversity Panel30 (SGDP50), and other species, the borders of our orthologous regions filter, multinucleotide var- de novo mutations from 1548 trios38 that were mapped to the human reference iants, InDels within the same species, or InDels with frequency ≥50% in another hg19 using the liftOver tool46. We used our own mapping and calling of 69 species. This resulted in a high-confidence set of species-exclusive trinucleotides chimpanzees and bonobos (59 chimpanzees and 10 bonobos, referred to as where the ancestral and derived alleles could be reliably inferred. This filtering chimpanzees in short)31,32 and 43 gorillas31,33. We used the release variant calling affected more CpG>T than non-CpG sites, due to the recurrent nature of CpG>T of three archaic samples: Altai and Vindija 33.19 Neanderthals47,48, and Deni- transitions (Supplementary Note 6, Supplementary Table 25). sova49. Finally, for the tumor data set, we used the release variant calling of 2583 human tumors from the Pan-Cancer Analysis of Whole Genomes Consortium34. Mutation spectra. We calculated each species’ mutation spectra as the fraction of all trinucleotides in a data set belonging to one of the 96 trinucleotides. We cal- culated correlations between data sets using Pearson’s correlation (cor.test function Definition of high-quality orthologous regions. We mapped and called chim- in R). We measured the correlation of the mutation spectrum of each species and panzee and bonobo, gorilla, and human (SGDP50) samples to the human reference the combined effect of cancer mutation signatures SBS1 and SBS5 (refs. 19,40) by hg19 using BWA MEM50 and GATK51 following the best practices protocols52,53. the formula: 0.1 × SBS1 + 0.9 × SBS5, as CpG>T transitions are the main compo- In addition, we removed variants with ≥20% missings, any variant (SNV or InDel) nents of signature SBS1 and they represent ~10% of the trinucleotides in both the within 5 bp of an InDel (suggesting potential misalignments), and variants where human and NHGA data sets. heterozygous samples represented more than 80% of the calls (suggesting potential mismappings or duplicated regions). To avoid mismappings to the human reference and erroneous estimates of Whole-genome enrichment of trinucleotides. We calculated the enrichment of mutation density in the NHGA samples (too low density caused by lack of trinucleotides and its significance in each germline data set pair (human–chimpanzee, mapping reads or deletions or too high density caused by collapsed duplications)54 human–gorilla, and chimpanzee–gorilla). We calculated the enrichment of NATURE COMMUNICATIONS | (2020)11:2512 | https://doi.org/10.1038/s41467-020-16296-4 | www.nature.com/naturecommunications 7 ARTICLE NATURE COMMUNICATIONS | https://doi.org/10.1038/s41467-020-16296-4 trinucleotide T between species A and B by dividing fraction of T in species A/fraction References of T in species B. We calculated a chi-squared test using a contingency table with the 1. Hodgkinson, A., Chen, Y. & Eyre-Walker, A. The large-scale distribution of trinucleotide count in species A, in species B and the count of the rest of trinucleotides somatic mutations in cancer genomes. Hum. Mutat. 33, 136–143 (2012). in species A, and in species B. As the counts of trinucleotides are not independent of 2. Tyekucheva, S. et al. Human-macaque comparisons illuminate variation in each other, we sorted all trinucleotides from most to least significant, and rerun the neutral substitution rates. Genome Biol. 9, R76 (2008). test by decreasing significance order, while removing the previously used trinucleo- 3. Ananda, G., Chiaromonte, F. & Makova, K. D. A genome-wide view of tides from the count of total trinucleotides26. mutation rate co-variation using multivariate analyses. Genome Biol. 12, R27 CpG>T transitions are highly affected by the sample size of the data sets. We (2011). ran all the tests using both 1KGP and SGDP50 as the human data set. We detected 4. Hodgkinson, A. & Eyre-Walker, A. Variation in the mutation rate across incoherences on the significance and direction of the results in two CpG>T mammalian genomes. Nat. Rev. Genet. 12, 756–766 (2011). trinucleotides. We report the results using 1KGP where tests using both 1KGP and 5. Makova, K. D. & Hardison, R. C. The effects of chromatin organization on SGDP50 are coherent in both significance and direction of the enrichment. variation in mutation rates in the genome. Nat. Rev. Genet. 16, 213–223 The top 10% most enriched trinucleotides in each species pairwise comparison (2015). were compared with cancer mutation signatures40 and reported when the 6. Schuster-Böckler, B. & Lehner, B. Chromatin organization is a major influence trinucleotide represented ≥5% of the mutations within a signature. on regional mutation rates in human cancer cells. Nature 488, 504–507 (2012). Trinucleotide-difference test. We developed a method to determine which tri- 7. Supek, F. & Lehner, B. Differential DNA mismatch repair underlies mutation nucleotides contribute significantly to the difference between NHGAs-tumors and rate variation across the human genome. Nature 521, 81–84 (2015). human-tumors mutation density correlations: 8. Supek, F. & Lehner, B. Scales and mechanisms of somatic mutation rate For each trinucleotide T and each pair of species (human–chimpanzee, variation across the human genome. DNA Repair 81, 102647 (2019). human–gorilla, and chimpanzee–gorilla), we subtract the ranked mutation density 9. Polak, P. et al. Cell-of-origin chromatin organization shapes the mutational of T in species A minus the ranking in tumor, and in species B minus tumor. We landscape of cancer. Nature 518, 360–364 (2015). calculate the two-sided Kolmogorov–Smirnov test (using the R function ks.test) of 10. Kübler, K. et al. Tumor mutational landscape is a record of the pre-malignant the two resulting distributions. We use the p value of the ks test as the significance state. Preprint at bioRxiv https://doi.org/10.1101/517565 (2019). of the test and the difference between the standard deviation of both distributions 11. Martincorena, I. et al. High burden and pervasive positive selection of somatic (as both have a mean of 0) as the test’s effect size. The results when using 1KGP or mutations in normal human skin. Science 348, 880–886 (2015). SGDP50 as the human data sets are concordant in the direction of the association, 12. Martincorena, I. et al. Somatic mutant clones colonize the human esophagus but we discarded the SGDP50 results because the smaller number of SNVs (and of with age. Science 362, 911–917 (2018). each trinucleotide type) results in lower power when using SGDP50. 13. Aggarwala, V. & Voight, B. F. An expanded sequence context model broadly explains variability in polymorphism levels across the human genome. Nat. Association of GC-content in the trinucleotide sequence. We counted the Genet. 48, 349–355 (2016). number of cytosine and guanine bases in each trinucleotide and built a linear 14. Hwang, D. G. & Green, P. Bayesian Markov chain Monte Carlo sequence regression (using the R function glm). The GC-content of the triplet acted as a analysis reveals varying neutral substitution patterns in mammalian evolution. predictor of the result of the test (the log10 fold-enrichment in the whole-genome Proc. Natl Acad. Sci. USA 101, 13994–14001 (2004). enrichment analysis or the difference between the standard deviation of both 15. Keightley, P. D., Eöry, L., Halligan, D. L. & Kirkpatrick, M. Inference of distributions in the trinucleotide-difference test). mutation parameters and selective constraint in mammalian coding sequences by approximate Bayesian computation. Genetics 187, 1153–1161 (2011). Signature-difference test. In order to determine the contribution of each 16. Siepel, A. & Haussler, D. Phylogenetic estimation of context-dependent mutation signature to the difference between NHGAs-tumors and human-tumors substitution rates by maximum likelihood. Mol. Biol. Evol. 21, 468–488 (2003). mutation density correlations, we rerun the trinucleotide-difference test using the 17. Nachman, M. W. & Crowell, S. L. Estimate of the mutation rate per nucleotide 1KGP and chimpanzee data sets, while using the different individual tumor types in humans. Genetics 156, 297–304 (2000). (Supplementary Data 1). For each trinucleotide, tumor type and mutation sig- 18. Alexandrov, L. B. et al. Signatures of mutational processes in human cancer. nature, we built a linear regression (using R’s glm function) where the mutation Nature 500, 415–421 (2013). load of that signature in that tumor type19 predicted the effect size in the 19. Alexandrov, L. B. et al. The repertoire of mutational signatures in human trinucleotide-difference test for that tumor type (Supplementary Note 6). For each cancer. Nature 578, 94–101 (2020). signature, we built a contingency table where all 96 trinucleotides where classified 20. Rahbari, R. et al. Timing, rates and spectra of human germline mutation. Nat. by whether being significant or not (p value < 0.05) in the trinucleotide-difference Genet. 48, 126–133 (2016). test and the significance of the mutation load in the linear regression model. We 21. Alexandrov, L. B. et al. Clock-like mutational processes in human somatic ran a chi-squared test on that contingency table and obtained its significance. cells. Nat. Genet. 47, 1402–1407 (2015). 22. Smith, T. C. A., Arndt, P. F. & Eyre-Walker, A. Large scale variation in the Reporting summary. Further information on research design is available in rate of germ-line de novo mutation, base composition, divergence and the Nature Research Reporting Summary linked to this article. diversity in humans. PLOS Genet. 14, e1007254 (2018). 23. Carlson, J. et al. Extremely rare variants reveal patterns of germline mutation rate heterogeneity in humans. Nat. Commun. 9, 3753 (2018). Data availability 24. Harris, K. Evidence for recent, population-specific evolution of the human All the analyses were performed using publicly available data obtained from their original mutation rate. Proc. Natl Acad. Sci. USA 112, 3439–3444 (2015). publications. Human data sets: 1000 Genomes Project (ftp://ftp.1000genomes.ebi.ac.uk/ 25. Mathieson, I. & Reich, D. Differences in the rare variant spectrum among vol1/ftp/release/20130502/), Simons Genome Diversity Project (EBI European human populations. PLOS Genet. 13, e1006581 (2017). Nucleotide Archive accession numbers PRJEB9586 [https://www.ebi.ac.uk/ena/data/ 26. Harris, K. & Pritchard, J. K. Rapid evolution of the human mutation spectrum. view/PRJEB9586] and ERP010710 [https://www.ebi.ac.uk/ena/data/view/PRJEB9586]), eLife 6, e24284 (2017). International Cancer Genome Consortium (https://dcc.icgc.org/releases/PCAWG/ 27. Narasimhan, V. M. et al. Estimating the human mutation rate from germline_variations), and Human de novo mutations (European Variant Archive autozygous segments reveals population differences in human mutational accession number PRJEB15197 [https://www.ebi.ac.uk/ena/data/view/PRJEB15197]). processes. Nat. Commun. 8, 303 (2017). Non-human great ape data sets: the Great Apes Genome Project (Sequence Read Archive 28. Besenbacher, S., Hvilsom, C., Marques-Bonet, T., Mailund, T. & Schierup, M. (SRA) accession numbers PRJNA189439 [https://www.ncbi.nlm.nih.gov/bioproject/? H. Direct estimation of mutations in great apes reconciles phylogenetic dating. term=PRJNA189439] and SRP018689 [https://www.ncbi.nlm.nih.gov/sra/? Nat. Ecol. Evol. 3, 286–292 (2019). term=SRP018689]), European Nucleotide Archive (ENA) accession numbers 29. The 1000 Genomes Project Consortium. A global reference for human genetic PRJEB15086 [https://www.ncbi.nlm.nih.gov/bioproject/?term=PRJEB15086] variation. Nature 526, 68–74 (2015). (Chimpanzee), PRJEB3220 [https://www.ebi.ac.uk/ena/data/view/PRJEB3220] (Gorilla), 30. Mallick, S. et al. The Simons Genome Diversity Project: 300 genomes from and PRJEB19688 [https://www.ncbi.nlm.nih.gov/bioproject/?term=PRJEB19688] 142 diverse populations. Nature 538, 201–206 (2016). (Orangutan). Tumor data sets: the Pan-Cancer Analysis of Whole Genomes (https://dcc. 31. Prado-Martinez, J. et al. Great ape genetic diversity and population history. icgc.org/pcawg/). Archaic samples (http://cdna.eva.mpg.de/neandertal/Vindija/). The Nature 499, 471–475 (2013). remaining data are available within the Article, Supplementary Information or available 32. de Manuel, M. et al. Chimpanzee genomic diversity reveals ancient admixture from the authors upon reasonable request. with bonobos. Science 354, 477–481 (2016). 33. Xue, Y. et al. Mountain gorilla genomes reveal the impact of long-term Received: 9 September 2019; Accepted: 22 April 2020; population decline and inbreeding. Science 348, 242–245 (2015). 34. The ICGC/TCGA Pan-Cancer Analysis of Whole Genomes Consortium. Pan- cancer analysis of whole genomes. Nature 578, 82–93 (2020). 8 NATURE COMMUNICATIONS | (2020)11:2512 | https://doi.org/10.1038/s41467-020-16296-4 | www.nature.com/naturecommunications NATURE COMMUNICATIONS | https://doi.org/10.1038/s41467-020-16296-4 ARTICLE 35. Ernst, J. et al. Mapping and analysis of chromatin state dynamics in nine Acknowledgements human cell types. Nature 473, 43–49 (2011). We thank B. Lehner, D. Weghorn, and I. Lobón for their insights discussing the analyses, 36. Vavouri, T. & Lehner, B. Human genes with CpG island promoters have a and C. Warembourg for her statistical analysis help. J.M.H.G is supported by MDM- distinct transcription-associated chromatin organization. Genome Biol. 13, 2014-0370. T.M.B. is supported by BFU2017-86471-P (MINECO/FEDER, UE), Howard R110 (2012). Hughes International Early Career and Obra Social “La Caixa”. D.J. is supported by Juan 37. Moorjani, P., Amorim, C. E. G., Arndt, P. F. & Przeworski, M. Variation in the de la Cierva fellowship (FJCI-2016-29558) from MICINN. A.N. is supported by AEI- molecular clock of primates. Proc. Natl Acad. Sci. USA 113, 10607–10612 PGC2018-101927-BI00(FEDER/UE), MINECO-BFU2015-68649-P (MINECO/FEDER, (2016). UE), the Spanish National Institute of Bioinformatics of the Instituto de Salud Carlos III 38. Jónsson, H. et al. Parental influence on human germline de novo mutations in (PT17/0009/0020), and by FEDER (Fondo Europeo de Desarrollo Regional)/FSE (Fondo 1,548 trios from Iceland. Nature 549, 519–522 (2017). Social Europeo). This work was supported by Unidad de Excelencia María de Maexto, 39. Milholland, B. et al. Differences between germline and somatic mutation rates funded by the AEI (CEXS2028-000792-M) and Secretaria d’Universitats i Recerca and in humans and mice. Nat. Commun. 8, 15183 (2017). CERCA Programme del Departament d’Economia i Coneixement de la Generalitat de 40. Tate, J. G. et al. COSMIC: the catalogue of somatic mutations in cancer. Catalunya (GRC 2017 SGR 880). Nucleic Acids Res. 47, D941–D947 (2019). 41. Blokzijl, F. et al. Tissue-specific mutation accumulation in human adult stem cells during life. Nature 538, 260–264 (2016). Author contributions J.M.H.G. performed all the analyses. J.M.H.G and D.J wrote the paper. T.M.B., D.J., and A.N. 42. Li, H. & Durbin, R. Inference of human population history from individual conceived and supervised this work. All the authors read and approved the final paper. whole-genome sequences. Nature 475, 493–496 (2011). 43. Harpak, A., Bhaskar, A. & Pritchard, J. K. Mutation rate variation is a primary determinant of the distribution of allele frequencies in humans. PLoS Genet. Competing interests 12, e1006489 (2016). The authors declare no competing interests. 44. Charlesworth, B. & Charlesworth, D. Elements of Evolutionary Genetics (Roberts, 2010). 45. Sundaram, L. et al. Predicting the clinical impact of human mutation with Additional information deep neural networks. Nat. Genet. 50, 1161–1170 (2018). Supplementary information is available for this paper at https://doi.org/10.1038/s41467- 46. Hinrichs, A. S. et al. The UCSC Genome Browser Database: update 2006. 020-16296-4. Nucleic Acids Res. 34, D590–D598 (2006). 47. Prüfer, K. et al. The complete genome sequence of a Neanderthal from the Correspondence and requests for materials should be addressed to D.J. or A.N. Altai Mountains. Nature 505, 43–49 (2014). 48. Prüfer, K. et al. A high-coverage Neandertal genome from Vindija Cave in Peer review information Nature Communications thanks Amy Boddy, Elaine Guevara, Croatia. Science 358, 655–658 (2017). and the other, anonymous, reviewer(s) for their contribution to the peer review of this 49. Meyer, M. et al. A high-coverage genome sequence from an archaic Denisovan work. Peer review reports are available. individual. Science 338, 222–226 (2012). 50. Li, H. Aligning sequence reads, clone sequences and assembly contigs with Reprints and permission information is available at http://www.nature.com/reprints BWA-MEM. Preprint at arXiv https://arxiv.org/abs/1303.3997 (2013). 51. McKenna, A. et al. The Genome Analysis Toolkit: a MapReduce framework Publisher’s note Springer Nature remains neutral with regard to jurisdictional claims in for analyzing next-generation DNA sequencing data. Genome Res. 20, published maps and institutional affiliations. 1297–1303 (2010). 52. DePristo, M. A. et al. A framework for variation discovery and genotyping using next-generation DNA sequencing data. Nat. Genet. 43, 491–498 Open Access This article is licensed under a Creative Commons (2011). Attribution 4.0 International License, which permits use, sharing, 53. Van der Auwera, G. A. et al. From FastQ data to high confidence variant calls: adaptation, distribution and reproduction in any medium or format, as long as you give the Genome Analysis Toolkit best practices pipeline. Curr. Protoc. Bioinform. appropriate credit to the original author(s) and the source, provide a link to the Creative 43, 11.10.1–11.10.33 (2013). Commons license, and indicate if changes were made. The images or other third party 54. Hartasánchez, D. A., Brasó-Vives, M., Heredia-Genestar, J. M., Pybus, M. & material in this article are included in the article’s Creative Commons license, unless Navarro, A. Effect of collapsed duplications on diversity estimates: what to indicated otherwise in a credit line to the material. If material is not included in the expect. Genome Biol. Evol. 10, 2899–2905 (2018). article’s Creative Commons license and your intended use is not permitted by statutory 55. Sudmant, P. H. et al. Evolution and diversity of copy number variation in the regulation or exceeds the permitted use, you will need to obtain permission directly from great ape lineage. Genome Res. 23, 1373–1382 (2013). the copyright holder. To view a copy of this license, visit http://creativecommons.org/ 56. Quinlan, A. R. & Hall, I. M. BEDTools: a flexible suite of utilities for licenses/by/4.0/. comparing genomic features. Bioinformatics 26, 841–842 (2010). 57. Barski, A. et al. High-resolution profiling of histone methylations in the human genome. Cell 129, 823–837 (2007). © The Author(s) 2020 NATURE COMMUNICATIONS | (2020)11:2512 | https://doi.org/10.1038/s41467-020-16296-4 | www.nature.com/naturecommunications 9