ARTICLE https://doi.org/10.1038/s41467-020-17162-z OPEN Germline de novo mutation rates on exons versus introns in humans Miguel Rodriguez-Galindo 1, Sònia Casillas 2,3, Donate Weghorn 1,4 ✉ & Antonio Barbadilla 2,3 ✉ A main assumption of molecular population genetics is that genomic mutation rate does not depend on sequence function. Challenging this assumption, a recent study has found a 1234567890():,; reduction in the mutation rate in exons compared to introns in somatic cells, ascribed to an enhanced exonic mismatch repair system activity. If this reduction happens also in the germline, it can compromise studies of population genomics, including the detection of selection when using introns as proxies for neutrality. Here we compile and analyze published germline de novo mutation data to test if the exonic mutation rate is also reduced in germ cells. After controlling for sampling bias in datasets with diseased probands and extended nucleotide context dependency, we find no reduction in the mutation rate in exons compared to introns in the germline. Therefore, there is no evidence that enhanced exonic mismatch repair activity determines the mutation rate in germline cells. 1 Centre for Genomic Regulation (CRG), The Barcelona Institute of Science and Technology, Dr. Aiguader 88, 08003 Barcelona, Spain. 2 Institute of Biotechnology and Biomedicine, Universitat Autònoma de Barcelona, Bellaterra (Cerdanyola del Vallès), 08193 Barcelona, Spain. 3 Department of Genetics and Microbiology, Universitat Autònoma de Barcelona, Bellaterra (Cerdanyola del Vallès), 08193 Barcelona, Spain. 4 Universitat Pompeu Fabra (UPF), Barcelona, Spain. ✉email: dweghorn@crg.eu; antonio.barbadilla@uab.es NATURE COMMUNICATIONS | (2020)11:3304 | https://doi.org/10.1038/s41467-020-17162-z | www.nature.com/naturecommunications 1 ARTICLE NATURE COMMUNICATIONS | https://doi.org/10.1038/s41467-020-17162-z O ne of the most general and widely accepted predictions DNM densities on exons and adjacent introns, namely extended of the neutral theory of molecular evolution is that “the sequence context dependency and several chromatin features, more sequence conservation, the more functional including H3K36me3 epigenetic marks. Finally, we provide esti- (selective) constraint on the sequence”1. This principle explains mates of exonic and intronic DNM rates. why different functional regions in the genome have different levels of polymorphism and divergence, such as the lower var- Results iation at nonsynonymous vs synonymous sites in protein-coding No evidence of reduced exonic DNM rate compared to introns. genes or in exonic vs intronic sequences2. This relationship To study the distribution of germline mutations across exons and between constraint and variation constitutes one of the most introns, we collected a total of 679,547 SNV DNMs from seven powerful approaches in the current search for functional regions family-based WGS datasets, consolidating a high-density, high- in the genome and the detection of natural selection at the quality DNM map across the human genome (see “Methods”). molecular level. An integral part of estimating constraint, or The compiled datasets show highly similar mutation spectra and purifying selection, on functional genomic regions is the com- an enrichment with CpG > TpG transitions (Supplementary parison of the observed number of mutations to the expectation Fig. 1). We first analyzed whether exonic and intronic mutation under neutral evolution. In genes, this neutral expectation is densities differ among human DNMs after accounting for usually estimated from putatively nonfunctional regions or sites, sequence composition. For that purpose, we computed the including intronic sequence3,4. A main requirement for the observed total mutation burden at exonic and intronic sites by validation of this assumption is that mutation rate on exons and summing over 95,633 internal exon-centered sequences of size introns does not correlate with that sequence function. 2001 base pairs (bp), carrying a subset of 50,780 genic mutations. Mutation rate can vary strongly across the human genome, Since per-nucleotide mutation probability is influenced by the with regional differences up to threefold in the germline5 and at neighboring sequence context, we derived the expected mutation least up to fivefold in tumor cells6. It is influenced by several burden at each position of each 2001-bp internal exon-centered factors, including replication time, chromatin state, and expres- window from a context-dependent model (see “Methods”). We sion level6–8. A priori, none of these factors are expected to initially used a trinucleotide-context-dependent germline whole- correlate directly with genic sequence function (exonic vs intro- genome mutation signature model, in line with the analysis nic). Another important determinant of mutation rate is DNA presented in Frigola et al.17 for somatic mutations. sequence composition. Recent studies have addressed mutational Frigola et al.17 found that the mutation burden of POLE- processes and their associated sequence-dependent signatures, mutant tumors in positions dominated by exonic DNA is lower both in the soma and the germline9–14. Germline and many than expected (Fig. 1a). In contrast, in our study the observed cancer tumor signatures exhibit a higher relative rate of C > T germline exonic mutation burden was significantly increased by transitions for single nucleotide variants (SNV)6,9,15. Conse- 7.2% (s.d. 1.4%; P = 0.001, permutation-based test) compared to quently, due to their higher G/C content, exonic regions show a the expectation across introns and exons (Fig. 1b). This result is context-driven relative increase in mutation rate compared to robust to biases in mutation calling due to region mappability intronic regions16, which can be corrected for with the proper differences (Supplementary Table 1, Supplementary Fig. 2) and mutational model. effects of transcription-coupled repair on the mutational pattern A differential mutation rate between intronic and exonic DNA (Supplementary Fig. 3). This suggests that the hypothesized beyond the context dependence would require a mutational pro- mechanism of enhanced repair on exons compared to introns in cess that recognizes the difference between the two functional POLE-aberrant tumors is not determining mutation rate in the sequence categories. Surprisingly, Frigola et al.17 found in tumoral germline. DNA, primarily from skin melanomas and DNA-polymerase-ε (POLE)-mutant colorectal cancers, that mutation rates are lower in exons than in introns after accounting for the trinucleotide- Sampling bias explains exonic mutation density excess. Even in context-dependent mutational signature. This reduced mutation the absence of the supposed enhanced MMR effect on exonic rate in exons is similar both in synonymous and nonsynonymous DNA, we would expect a very slight deficit of exonic relative to sites, which rules out purifying selection as an explanation. The intronic DNMs, due to strong purifying selection on lethal de study suggests that the lower mutation rate in exons results from novo variants in the early stages of embryonic development23. an enhanced mismatch repair (MMR) activity in exons compared Therefore, we next investigated potential factors that could to introns. In turn, the increased repair activity is attributed to explain the observed increase in exonic relative to intronic different amounts of H3K36me3 epigenetic marks on exons and mutation burden in DNM datasets, such as technical differences introns17. in sequencing and calling (Supplementary Table 2) and In the germline, whether originating from replication errors or enrichment with diseased probands. As detailed in Table 1, the mediated by DNA damage, the dominant mutational processes analyzed DNM datasets are heterogeneous regarding their study are expected to produce mismatches18–22. Hence, any MMR- conditions, including disease cohorts. Diseased probands, e.g. related mechanism is expected to play an important role in those with autism spectrum disorder (ASD) or preterm birth, germline DNA damage repair. If the enhanced somatic exonic are more likely than average to carry mutations with functional MMR activity found by Frigola et al.17 could be extrapolated to impact24,25. This ascertainment bias in the data is expected to the germline, as the study suggests, then population and func- entail an enrichment with exonic nonsynonymous variants26,27. tional genomics studies would be compromised, and they should To test this, we classified the 4669 mutations in the internal include differential exonic and intronic mutation rates as an exons as synonymous and nonsynonymous changes, resulting integral part of their explanatory models. in 3488 nonsynonymous and 1170 synonymous DNMs (corre- Here, we investigate the relative mutation rates of exons and sponding to a ratio of 2.98:1). We then repeated the internal introns in the human germline using de novo mutation (DNM) exon-centered analysis for each of the two exonic mutation data. We show that DNM densities do not differ between exons categories. and introns after accounting for trinucleotide sequence compo- Figure 2a shows that the observed synonymous profile sition and an excess of nonsynonymous exonic variation arising matches the expected profile almost perfectly, with a slight from sampling bias. We further explore factors that can impact nonsignificant deficit (−1.1%, P = 0.353). However, exonic 2 NATURE COMMUNICATIONS | (2020)11:3304 | https://doi.org/10.1038/s41467-020-17162-z | www.nature.com/naturecommunications NATURE COMMUNICATIONS | https://doi.org/10.1038/s41467-020-17162-z ARTICLE a Colorectal POLE mutants (n = 6) b Germline cells (n = 11,237 trios) Observed muts Observed muts Expected muts 60 Expected muts 80 Number of mutations 60 40 40 20 20 –1000 –500 0 500 1000 –1000 –500 0 500 1000 100 100 Exon (%) 50 50 0 0 –1000 –500 0 500 1000 –1000 –500 0 500 1000 Distance from exon center (bp) Distance from exon center (bp) Exon Exon Intron Intron Intron Intron Fig. 1 Internal exon-centered analyses on somatic and germline de novo mutations. Exon-centered 2001-nt-wide observed and expected mutational profiles (top) and exon density (bottom) in a somatic and b germline cells. The light red line represents the observed number of mutations at each position, whereas the dark red and black lines represent smoothed numbers of observed and expected mutations, respectively, obtained from a polynomial fit. a Profile of mutations in six POLE-mutant colorectal tumors, reprinted by permission from Springer Nature Customer Service Centre GmbH: Springer Nature Genetics, Reduced mutation rate in exons due to differential mismatch repair, Frigola et al.17. b Profile of mutations in the germline of 11,237 trios. Table 1 Properties of analyzed DNM datasets including excess in exonic burden. Dataset Phenotypic condition Control/proband DNMs Exonic excess [%] Emp. p value Halldorsson36 Mixed (multiple diseases) 180,151 11.7 ± 2.8 0.001 Yuen38 Mostly ASD 127/117,612 8.0 ± 4.1 0.015 Goldmann35 Mixed (preterm birth) 35,793 13.0 ± 6.6 0.016 Sasani39 Random (mostly healthy) 27,454/0 12.3 ± 7.9 0.052 An37 ASD + healthy sibling 115,697/117,942 3.5 ± 2.2 0.070 GoNL16 Random (mostly healthy) 11,016/0 12.1 ± 14.5 0.181 Goldmann29 Healthy 73,755/0 2.5 ± 4.4 0.291 Autism probands ASD 0/235,554 8.4 ± 2.5 0.001 Representative sample Non-ASD + ASD 98,300/1,700 −2.5 ± 3.5 0.233 Healthy probands Healthy 189,579/0 −0.5 ± 2.5 0.417 The seven used datasets and their references are shown above the line, while below the results for the composed datasets (see “Methods”) are given. In each group, datasets are ordered from most to least significant exonic mutation excess. Errors of the exonic excess denote 1 s.d. from 1000 permutations (see “Methods”). nonsynonymous mutations show a large and statistically De novo variants show extended-context dependency. The significant excess compared to the expectation under the stratification into synonymous and nonsynonymous changes trinucleotide-context model (10.4%, P = 0.001, Fig. 2b). Note entailed a polarization of the exonic excess (Fig. 2), intensifying that this stratification by functional mutation category entails a the signal for nonsynonymous variants (10.4%) with a con- reduction of the number of both synonymous and nonsynon- comitant decrease for synonymous variants (−1.1%). This type of ymous mutations relative to flanking introns across the window polarization could be due to the incompleteness of our mutational of stacked sequences. This is why the number of mutations, model. Our mutational model for the expected number of when moving outwards from the center, converges to that of mutations was constructed using trinucleotide-context-dependent Fig. 1b. Overall, Fig. 2 suggests that disease ascertainment mutation probabilities. However, it has been shown that SNPs during data acquisition may be responsible for the overall segregating in the human population are affected by the extended excess of 7.2% of exonic variants. Therefore, we next repeated flanking sequence, with a heptameric context explaining a the analysis for all seven DNM datasets individually, as well as majority of the observed mutation rate variability13,14. Figure 3 for assembled samples with only healthy, only ASD, or a and Supplementary Figs. 4, 5 show that this is confirmed by representative mixture of probands (see “Methods”). We found DNMs based on the relative frequencies of all four nucleotides a significant exonic mutation excess only in cohorts with a high around mutations in our DNM dataset, although the effect of the fraction of diseased probands or those assembled purely from extended flanking sequence is small compared to the one diseased samples, while all cohorts with mostly or exclusively observed in POLE-mutated tumor genomes28. healthy probands show no signal (Table 1). Moreover, when we We assessed the impact of context dependency by expanding stratify exonic mutations by functional impact, we find a our mutational signature model to incorporate the pentameric and statistically significant excess only in nonsynonymous muta- heptameric mutational sequence context (based on exact computa- tions among ASD individuals (Supplementary Table 3). tion and a likelihood decomposition approach, respectively; see NATURE COMMUNICATIONS | (2020)11:3304 | https://doi.org/10.1038/s41467-020-17162-z | www.nature.com/naturecommunications 3 ARTICLE NATURE COMMUNICATIONS | https://doi.org/10.1038/s41467-020-17162-z a Synonymous germline b Nonsynonymous germline 40 40 Number of mutations Number of mutations 20 20 Observed muts Observed muts Expected muts Expected muts 0 0 –1000 –500 0 500 1000 –1000 –500 0 500 1000 Distance from exon center (bp) Distance from exon center (bp) Fig. 2 Internal exon-centered analysis for synonymous and nonsynonymous DNMs. The light red line represents the observed number of mutations at each nucleotide position, while the dark red and black lines represent averages in bins of size 25 positions for observed and expected mutations, respectively. a Synonymous DNM profile with observed exonic mutation difference of −1.1% (P = 0.353). b Nonsynonymous DNM profile, showing an exonic excess of 10.4% (P = 0.001). Due to the removal of nonsynonymous and synonymous mutations in panels (a) and (b), respectively, the total number of exonic mutations relative to flanking introns is reduced. The number of mutations, when moving away from the center, converges to that of Fig. 1b. Table 2 Extended sequence context dependency for complements the previously reported influence of nucleosome Goldmann et al.29. organization on human germline DNMs5,33. Model Exonic Emp. Log- # param AIC Estimation of exonic and intronic de novo mutation rate. We excess [%] p value likelihood estimated germline DNM rates for exons and introns separately. 1-mer 16.2 ± 5.4 0.001 −68,111 12 136,247 Using the largest dataset with only healthy probands29, we CpG 1.7 ± 4.3 0.348 −66,594 18 133,224 estimated 1.38 × 10−8 and 1.11 × 10−8 mutations per site per 3-mer 2.7 ± 4.3 0.281 −66,315 192 133,014 generation for exons and introns, respectively. This difference 5-mer 2.9 ± 4.6 0.258 −66,054 1344 134,797 reflects the higher mutability of exons with respect to introns 7-mer 2.8 ± 4.3 0.266 −65,981 2496 136,955 due to sequence differences, namely higher CpG and overall GC content16. These estimates are consistent with a previously Errors of the exonic excess denote 1 s.d. from 1000 permutations (see “Methods”). reported whole-genome based rate of 1.2 × 10−8 mutations per site per generation20. Conversely, estimates obtained from the “Methods”, Supplementary Fig. 6). We applied these extended- pooled dataset reflect the disease ascertainment bias, with a context models to the largest DNM dataset that had no diseased similar intronic mutation rate (1.15 × 10−8), but a much larger probands, Goldmann et al.29, as well as to this dataset and the exonic rate (1.52 × 10−8), in line with previous findings in dis- pooled dataset stratified by synonymous and nonsynonymous eased cohorts34. variants. We found that while the overall likelihood increases for increasing context size, penalization due to the additional Discussion parameters of the larger context models entails that the We compiled human DNM data to show that the rate of genera- trinucleotide-context-dependent model is found to be the best tion of new genetic variants, the mutation rate, does not sig- model for the current datasets based on the Akaike information nificantly vary between exons and adjacent introns when criterion (AIC) (Table 2, Supplementary Tables 4–6). accounting for sequence context. Moreover, we went beyond pre- vious analyses that used extreme rare variants as a proxy for H3K36me3 does not correlate with exonic mutation density. DNMs13,14 and described directly germline mutation patterns The enhanced exonic MMR activity compared to introns in based on a large aggregated DNM dataset. We corroborated earlier POLE-aberrant tumors was proposed to be mediated by the findings, in particular an extended-context dependence of germline H3K36me3 mark17. Using our dataset of mutations from healthy variants. At the same time, the internal exon-centered analysis, probands, we therefore investigated the relative exonic mutation with its relatively low number of mutations compared to the entire density as a function of H3K36me3 and nucleosome density. dataset, is still adequately described by a trinucleotide-context- These two features show differential coverage between (mainly dependent model. Beyond context dependence, the sampling bias internal) exons and introns (Supplementary Fig. 7), and had been introduced by enrichment with diseased probands is one of the previously described to contribute to the recognition of splice most important confounding factors of DNM analyses. We showed marks at internal exon–intron boundaries30,31. We observed no that its effects can lead to significant deviations from the null significant correlation (r = −0.03, P = 0.84) between the exonic model and, depending on the application, should be addressed mutation enrichment and the H3K36me3 mark (Fig. 4b). This through an informed choice of samples. contrasts with the recruitment mechanism described in somatic Our analysis shows that the results found in the soma cannot cells32, which is invoked as the mechanistic hypothesis behind the be directly extrapolated to the germline, and the MMR-dependent findings in Frigola et al.17 (Fig. 4a). Conversely, we find a nega- process that was proposed as an explanation for the decreased tive, nearly significant correlation (r = −0.28, P = 5.38 × 10−2) exonic mutation burden in somatic cells does not seem to with nucleosome coverage (Supplementary Fig. 8). This result determine germinal cell mutation rates. Last, this study provides a 4 NATURE COMMUNICATIONS | (2020)11:3304 | https://doi.org/10.1038/s41467-020-17162-z | www.nature.com/naturecommunications NATURE COMMUNICATIONS | https://doi.org/10.1038/s41467-020-17162-z ARTICLE a 100% 75% A % of bases 50% C G T 25% 0% –10 –9 –8 –7 –6 –5 –4 –3 –2 –1 0 1 2 3 4 5 6 7 8 9 10 Base position C > T (G > A) T > C (A > G) 100% 100% 75% 75% % of bases 50% 50% 25% 25% 0% 0% –10 –9 –8 –7 –6 –5 –4 –3 –2 –1 0 1 2 3 4 5 6 7 8 9 10 –10 –9 –8 –7 –6 –5 –4 –3 –2 –1 0 1 2 3 4 5 6 7 8 9 10 Base position Base position b –1 3-mer –2 –1 5-mer –3–2–1 7-mer C > T (G > A) Relative rate 1e–03 +1 +1 +1 +2 +2 +3 1e–04 1e–05 T > C (A > G) 1e–06 +1 +1 +1 +2 +2 +3 Fig. 3 Mutational classes with highest sequence context dependency around de novo mutations. a Frequencies of nucleotides neighboring our entire set of 679,547 DNMs and subsets belonging to C > T (G > A) and T > C (A > G) 1-mer classes in a window of size 21 bp. Black dashed lines represent the whole-genome background frequencies for the four nucleotides. The extended sequence context dependency varies across 1-mer mutation classes. b Heatmap of estimated relative mutation rates, corrected by abundance of reference k-mer, for C > T (G > A) (top) and T > C (A > G) (bottom) 1-mer classes, up to a 7-mer resolution. For each 1-mer class, each of the three grids delineates mutation contexts of different length, defined by the upstream sequence (y-axis) and downstream sequence (x-axis) from the central (mutated) nucleotide. Boxed regions indicate motifs previously identified as hypermutable (pink) or hypomutable (green). Supplementary Figs. 4 and 5 show the corresponding plots for the other 1-mer classes. NATURE COMMUNICATIONS | (2020)11:3304 | https://doi.org/10.1038/s41467-020-17162-z | www.nature.com/naturecommunications 5 ARTICLE NATURE COMMUNICATIONS | https://doi.org/10.1038/s41467-020-17162-z a Colorectal POLE mutants (n = 6) b Germline cells (n = 3195 trios) 0 r = –0.68 40 r = –0.03 Exonic mutations difference (%) Exonic mutations difference (%) –10 P = 6.71 × 10 –8 P = 8.41 × 10–1 –20 20 –30 0 –40 –50 –20 –60 –40 –70 1 2 3 4 5 1 2 3 4 5 H3K36me3 exon-to-intron ratio H3K36me3 exon-to-intron ratio Fig. 4 Deviation in the exonic mutation burden as a function of the H3K36me3 exon-to-intron ratio. Blue dots denote 50 groups of genes binned by their exon-to-intron ratio of H3K36me3 coverage (x-axis). The relative difference between the total observed and expected number of exonic mutations (computed using a 3-mer model) per group is shown on the y-axis. The trend line and its confidence interval were added using the seaborn package of Python, while the correlation coefficient and its significance were computed using the same iteratively re-weighted least-squares approach as used by Frigola et al.17 to ensure comparability. a POLE-mutant colorectal tumors. The H3K36me3 histone mark is derived from colonic mucosa (E075), reprinted by permission from Springer Nature Customer Service Centre GmbH: Springer Nature Genetics, Reduced mutation rate in exons due to differential mismatch repair, Frigola et al.17. b DNMs from healthy probands (a total of 3195 trios). The H3K36me3 histone mark is derived from H1 stem cells (E003). clear-cut answer to the challenge posed by Frigola et al.17: It Moreover, positions where mutation calling would be technically challenging validates the main assumption of molecular population genetics— because of mappability issues were removed, leaving a total of 10,237 genes for the gene by gene analyses. We also filtered out the internal exon-centered 2001-nt that genomic mutation rate does not depend on sequence func- windows that overlapped at least one nucleotide with regions with mappability tion—and demonstrates different mutational dynamics in issues (Supplementary Table 1) for the supplementary internal exon-centered somatic vs germinal cells. analysis restricted to highly mappable regions (Supplementary Fig. 2). Coordinates of unreliable regions42 were obtained from the UCSC Genome Browser, available at http://genome.ucsc.edu/cgi-bin/hgFileUi?db=hg19&g=wgEncodeMapability. Methods Meta-exon and meta-intron coordinates of genes with at least five meta-exons De novo mutation data. We aggregated DNMs from seven family-based WGS (not only internal) were extracted from GENCODE v19 for the analyses with datasets coming from multiple centers and projects: the Genomes of the Nether- chromatin features across genic regions (Supplementary Fig. 7). lands (GoNL) project16, the Inova Translational Medicine Institute Preterm Birth Study35, Inova Translational Medicine Institute’s Longitudinal Childhood Genome Sequence context model. We directly computed the probability of a mutation Study29, deCODE genetics36, Simons Simplex Collection (SCC) and Korean ASD into an alternative nucleotide, Ha where a ∈ {1, 2, 3}, given the reference nucleotide cohort37, the Autism Genetic Research Exchange (AGRE) repository38 and Centre Hr and its flanking sequence X ¼ ðX50 ; X30 Þ, where X50 and X30 are the 5′ and 3′ d’Etude du Polymorphisme Humain (CEPH)39. Mutation datasets were down- flanking sequences, respectively. Here, Hr ∈ M = {A, C, G, T}, where the latter loaded from the supplementary tables of the respective papers35–39 or by direct denote nucleotides adenine, cytosine, guanine, and thymine, and request to the authors29. Data from the GoNL were downloaded from http://www. H a 2 M0 ¼ fm 2 M ; m ≠ H r g. For example, for the 5-mer GCACG > GCTCG nlgenome.nl. Most of datasets were originally mapped to hg19, with exception of mutation Hr = A, Ha = T, X50 ¼ ð G; C Þ, X30 ¼ ð C; G Þ and X = (G, C, C, G). the data from An et al.37 and Halldorson et al.36, which were mapped to hg38. Therefore, the probability of each of the possible k-mer changes, normalized by the Subsequently, coordinates in these datasets were lifted over to hg19, the most abundance of each reference k-mer in the genome, was computed as follows: common reference genome in our data. To avoid possible biases arising from mutation calling on sexual chromosomes, only autosomal SNVs were used, leaving NðH a ; H r ; XÞ a total of 679,547 germline SNV DNMs coming from 11,237 trios. PðH a jH r ; XÞ ¼ ; ð1Þ GðH r ; XÞ Effect prediction of de novo mutations. The predicted consequence class of all where N(Ha, Hr, X) is the genome-wide number of observed mutations into DNMs was obtained using the Ensembl Variant Effect Predictor (VEP)40 for the alternate allele Ha with given reference allele Hr and flanking sequence X. G(Hr, X) GRCh37/hg19 assembly. Since some DNMs were reported to have more than one is the abundance of the reference k-mer with reference nucleotide Hr and flanking consequence, e.g. different transcripts or overlapping genes, only one predicted sequence X in the genome. We computed the relative abundance of each reference consequence for each DNM was retrieved (according to VEP criteria). Predictions k-mer in the autosomal genome using the pyFasta package. We also computed are classified from major to mild according to the Ensembl Variation hierarchy. strand-wise signatures restricted to mutations falling in genic regions (exons at the canonical CDS and the respective introns) for the supplementary analysis in Supplementary Fig. 3. To compute strand-wise signatures, we polarized mutations Whole-genome de novo mutation spectrum. All mutations were divided into according to the transcription strand on which the canonical CDS is annotated. nine classes, considering the fact that CpG sites are highly mutagenic. The number For some k-mer models, given limitations imposed by the amount of DNMs, we of mutations were corrected by the relative abundance of the context in the whole used a decomposition approach to compute the probability. For a k-mer model of genome, e.g. the total number of C > T (G > A) transitions occurring at CpG sites sequence length k, let Hr be a reference core h-mer of length h and Ha an alternate divided by the relative abundance of CpG sites in the genome. We performed the core h-mer, where 1 ≤ h < k. Let the tuple X = (x1, …, xg) with g = (k − h) mutational analysis across all used studies (Supplementary Fig. 1). For all sub- elements represent again the flanking sequence of the core h-mer, where xi ∈ sequent analyses, extended nucleotide-context-dependent mutational models M ∀i ∈ {1,…, g}. In other words, X ∈ Mg where Mg is the g-fold Cartesian product. were used. For example, with k = 7, h = 3 and the mutation ACTGACT > ACTCACT, then Hr = TGA, Ha = TCA and X = (x1 = A, x2 = C, x3 = C, x4 = T). We then approximate the probability P(Ha∣Hr, X) by: Genomic coordinates of internal exons and flanking introns. Coordinates for a total of 20,345 protein-coding genes were obtained from GENCODE v19 41. Genes Yg Pi ðH a jH r ; xi Þ without introns and overlapping genes were discarded, leaving a filtered set of PðH a jH r ; XÞ  PðH a jH r Þ  ; ð2Þ i¼1 PðH a jH r Þ 13,474 genes. Genes located on chromosomes X, Y and on the mitochondrial genome were removed from the analysis, leaving a total of 12,754 autosomal genes. with Finally, all transcripts per gene were merged into meta-exon and meta-intron coordinates, both 5′ and 3′ flanking exons were removed as well as UTRs. Only P PZ2Mg NðH a ; H r ; ZÞ internal exons (unfiltered by mappability issues, see below) were used for the main PðH a jH r Þ ¼ ; ð3Þ internal exon-centered mutational analyses. Z2Mg GðH r ; ZÞ 6 NATURE COMMUNICATIONS | (2020)11:3304 | https://doi.org/10.1038/s41467-020-17162-z | www.nature.com/naturecommunications NATURE COMMUNICATIONS | https://doi.org/10.1038/s41467-020-17162-z ARTICLE and healthy probands and 1700 classified as ASD and repeated the internal exon- P centered analysis. To generate the purely healthy and purely ASD cohorts, we used Z2Yxi NðH a ; H r ; ZÞ solely mutations from probands with the respective condition. Pi ðH a jH r ; xi Þ ¼ P ; ð4Þ Z2Yxi GðH r ; ZÞ where Nucleosome and H3K36me3 histone mark data. We downloaded narrow peak ( ( ) coordinates and genome-wide read-coverage of H3K36me3 from human if j ¼ i; yj ¼ xi embryonic stem cell H1-hESC (E003), as proxy for germline cells, from the Epi- Yxi ¼ ðy1 ; :::; yj ; :::; yg Þ 8j 2 f1; :::; gg : ð5Þ genome Roadmap consortium43 data portal (http://www.roadmapepigenomics.org/ otherwise ; yj 2 M data). The genome-wide nucleosome positioning density graph of ENCODE44 cell We implemented this framework using custom Python code. The composite line GM12878 (lymphoblastoid cell line) was obtained via the UCSC genome likelihood model was applied to the 7-mer analysis of the data pooled across all browser (https://hgdownload.soe.ucsc.edu/downloads.html). Nucleosome peak cohorts using k = 7 and h = 5. Also in the analysis of the largest single dataset regions were identified across the genome by using the bwtool program (with purely composed of healthy probands29 and the largest single dataset36, in each for parameters local-extrema -maxima -min-sep = 150). The window of 146 bp 5-mers with k = 5 and h = 3 and for 7-mers with k = 7 and h = 3. The rest of flanking the peak coordinate (73 bp per side) was considered the region covered by probabilities were computed using the direct approach. Supplementary Fig. 6 a nucleosome. shows the relationship between the exact computations of mutational probabilities and the composite likelihood model for the pooled dataset. The number of parameters for the direct approach increases exponentially as k Coverage of chromatin features across exons and introns. Exons and introns in increases following f(x) = 4x ⋅ 3, where x = k. For the decomposition approach they each gene were classified according to their position with respect to the tran- increase linearly as k − h increases from a fixed h following g(y, x) = f(x) ⋅ (1 + scription start site, where the ones that occupy different positions in different 3(y − x)), where y = k and x = h. transcripts were discarded. We also discarded exons and introns at the lower quartile of length to compute the coverage for a set of exons or introns of het- erogeneous lengths in a given position: the fraction of bases covered by H3K36me3 Comparison of sequence context models. We selected the sequence context and nucleosomes at the center of the stack corresponding to the window defined by dependency model that best explains mutations across our set of exonic and the shortest exon or intron remaining after the filtering. Finally, the difference intronic sequences by means of the AIC. We interrogated each of the 191,361,633 between the exonic and intronic coverage was computed via the two-tailed (~6.4% of the whole-genome length) exonic and intronic sites on the 95,633 2001- Mann–Whitney p value of the comparison of both distributions. nt windows for the state in the observed data: mutated (and type) or not mutated. We also computed the positions in the genome covered by H3K36me3 or For recurrent sites, we chose one observed mutation at random. For a given model, nucleosomes across 95,633 internal exon-centered 4001-nt windows. By stacking we computed the log-likelihood as the sum across sites of the logarithm of the sequences, we obtained middle exon-centered profiles of coverage across exons and estimated probability of the observed state at the site. Probabilities were estimated introns (Supplementary Fig. 7). with mutations from the entire dataset, through the direct or the decomposition approach as stated above. Nucleosome and H3K36me3 binned gene analysis. For each gene, we computed the readcount-based exonic enrichment of H3K36me3 or nucleosomes as the ratio Internal exon-centered mutational analysis. A total of 95,633 stacked 2001-nt between the exonic and intronic total number of bases covered by reads of the sequences centered on the middle position of internal meta-exons were used to chromatin feature. Genes with no exonic and intronic bases covered by reads were compare the observed and expected mutational profiles across exons and introns. removed from the analysis, as well as genes without any observed exonic or We computed the frequency of mutation at a site l with reference core sequence H lr intronic mutation. Thus, a total of 7215 and 6529 genes remained for the and flanking sequence Xl as H3K36me3 and nucleosome analysis, respectively. X 3 For a given gene, we computed the exonic expected number of mutations as fl ¼ PðH la jH lr ; Xl Þ; l 2 f1; :::; Lg ; ð6Þ follows: a¼1 ^e ¼ Pe  n ; n ð10Þ where L = 2001 denotes the total number of considered sites. Then, each frequency was normalized by the total frequency on the sequence: where n is the total number of mutations (both exonic and intronic) observed on the gene and Pe is the (binomial) probability of a mutation to fall on the exonic f region of the gene, which in turn is computed as: f lresc ¼ PL l : ð7Þ l0 ¼1 f l0 Le Finally, the total number of observed mutations ns on each of the 2001-nt Pe ¼ ; ð11Þ Le þ Li sequence s, of a total of S = 95,633 stacked sequences, was redistributed across both middle exonic and flanking intronic sites according to the normalized frequencies: with ^ls ¼ f lresc  ns ; n s 2 f1; :::; Sg ; ð8Þ X Le X 3 Le ¼ PðH lae jH lre ; Xle Þ ; ð12Þ thus yielding the expected number of mutations at site l of a given sequence s. By le ¼1 a¼1 adding up the values of all the stacked sequences, we obtain the cumulative number of expected mutations at site l, and X S Li X X 3 ^l ¼ n ^ls : n ð9Þ Li ¼ PðH lai jH lri ; Xli Þ : ð13Þ s¼1 li ¼1 a¼1 For the internal exon-centered analysis on synonymous or nonsynonymous mutations, we separated all possible exonic mutations in middle exon sequences Here, le ∈ {1,…, Le} and li ∈ {1,…, Li} denotes the set of all exonic and intronic into two groups: those with synonymous consequence and those with a positions of a given gene, respectively. Le and Li represent the exonic and intronic consequence ranking higher than synonymous in the Ensembl Variation hierarchy. target size, respectively, expressed as the sum of the probability PðH la jH lr ; Xl Þ of all Then we computed the expected numbers by only adding frequencies for either possible three mutations that can happen across all exonic or intronic sites of the synonymous or nonsynonymous mutations. gene. The probability was computed under a 3-mer model for each of the genes as explained above only with mutations from the composed dataset of healthy probands. Computation of effect size and statistical significance. We performed 1000 Afterwards, genes were grouped into 50 bins according to their exonic random permutations of the observed mutations in each stacked sequence based on enrichment of H3K36me3 or nucleosomes. Then, with the observed ne and the probability of each site to acquire a mutation. The effect size, defined as the expected n ^e exonic mutations over all genes in the bin, we computed the relative relative increase or decrease in observed exonic mutations with respect to the difference between the observed and expected number of exonic mutations per bin expected number, was computed based on the simulation mean expected value. as follows: The error of this estimate is given as one standard deviation derived from the 1000 permutations. Moreover, we computed an empirical one-sided p value as the ne  n^e Exonic mutations difference ½% ¼  100 : ð14Þ fraction of the simulations with more (or fewer) exonic mutations than the ^e n observed number of exonic mutations. Finally, we computed the correlation between the median exonic chromatin feature enrichment and the difference in exonic mutations across the bins. The Composed datasets. Mutations were only resampled from datasets with known trend line and its confidence intervals were added using the bootstrapping conditions of the probands, either from healthy probands or those with ASD. functions of the python seaborn package, which confers equivalent weights in the Given an ASD prevalence in humans of 1.7%, we created a random sample of regression to all points. In order to guarantee that the trend is not the result of a 100,000 whole-genome DNMs, taking 98,300 mutations classified as strictly from few outliers, the correlation coefficient and its significance were computed using an NATURE COMMUNICATIONS | (2020)11:3304 | https://doi.org/10.1038/s41467-020-17162-z | www.nature.com/naturecommunications 7 ARTICLE NATURE COMMUNICATIONS | https://doi.org/10.1038/s41467-020-17162-z iteratively re-weighted least-squares approach, letting the variance of exonic 15. Cooper, D. N. & Youssoufian, H. The CpG dinucleotide and human genetic chromatin feature enrichment of the bins influence the weight of each point. disease. Hum. Genet. 78, 151–155 (1988). 16. Francioli, L. C. et al. Genome-wide patterns and properties of de novo Estimation of absolute mutation rates. We estimated absolute mutation rate in mutations in humans. Nat. Genet. 47, 822–826 (2015). our set of 95,633 middle exons as a proxy of mean exonic mutation rate and 17. Frigola, J. et al. Reduced mutation rate in exons due to differential mismatch absolute mutation rate on the rest of the 2001-nt window as a proxy of mean repair. Nat. Genet. 49, 1684–1692 (2017). intronic mutation rate as follows: 18. Crow, J. F. The origins, patterns and implications of human spontaneous mutation. Nat. Rev. Genet. 1, 40–47 (2000). N obs μ¼ : ð15Þ 19. Ozturk, S. & Demir, N. DNA repair mechanisms in mammalian germ cells. N site  N gen Histol. Histopathol. 26, 505–517 (2011). Here, μ is the mutation rate per site and generation, Nobs is the number of observed 20. Milholland, B. et al. Differences between germline and somatic mutation rates mutations, Nsite is the number of sites (we used Nsite = 13,632,264 exonic sites and in humans and mice. Nat. Commun. 8, 1–8 (2017). Nsite = 177,729,369 flanking intronic sites) and Ngen is the number of generations. 21. Gao, Z. et al. Overlooked roles of DNA damage and maternal age in A total of Ngen = 2582 gametogeneses for the largest dataset with healthy pro- generating human germline mutations. Proc. Natl Acad. Sci. USA 116, bands29 and Ngen = 22,474 gametogeneses in the pooled dataset. 9491–9500 (2019). 22. García-Rodríguez, A., Gosálvez, J., Agarwal, A., Roy, R. & Johnston, S. DNA damage and repair in human reproductive cells. Int. J. Mol. Sci. 20, 1–22 (2019). Reporting summary. Further information on research design is available in 23. Lynch, M. Rate, molecular spectrum, and consequences of human mutation. the Nature Research Reporting Summary linked to this article. Proc. Natl Acad. Sci. USA 107, 961–968 (2010). 24. Sanders, S. J. et al. De novo mutations revealed by whole-exome sequencing Data availability are strongly associated with autism. Nature 485, 237 (2012). All the analyses in this study were based on published datasets. Mutation data from the 25. Li, J., Oehlert, J., Snyder, M., Stevenson, D. K. & Shaw, G. M. Fetal de novo Genomes of the Netherlands (GoNL) project16 was downloaded from (http://www. mutations and preterm birth. PLoS Genet. 13, e1006689 (2017). nlgenome.nl). The remaining mutation datasets were either by direct request to the 26. Iossifov, I. et al. The contribution of de novo coding mutations to autism authors29 or downloaded from the supplementary tables of their respective publications35– spectrum disorder. Nature 515, 216 (2014). 39. Coordinates of unreliable regions42 were obtained from the UCSC Genome Browser, 27. Kosmicki, J. A. et al. Refining the role of de novo protein-truncating variants available at http://genome.ucsc.edu/cgi-bin/hgFileUi?db=hg19&g=wgEncodeMapability. in neurodevelopmental disorders by using population reference samples. Nat. Narrow peak coordinates and genome-wide read-coverage of H3K36me3 from human Genet. 49, 504 (2017). embryonic stem cell H1-hESC (E003) were downloaded through the Epigenome Roadmap 28. Martincorena, I. et al. Universal patterns of selection in cancer and somatic consortium43 data portal (http://www.roadmapepigenomics.org/data). The genome-wide tissues. Cell 171, 1029–1041.e21 (2017). nucleosome positioning density graph of ENCODE44 cell line GM12878 (lymphoblastoid 29. Goldmann, J. M. et al. Germline de novo mutation clusters arise during oocyte cell line) was obtained via the UCSC genome browser (https://hgdownload.soe.ucsc.edu/ aging in genomic regions with high double-strand-break incidence. Nat. downloads.html). Genet. 50, 487–492 (2018). 30. Tilgner, H. et al. Nucleosome positioning as a determinant of exon recognition. Nat. Struct. Mol. Biol. 16, 996–1001 (2009). Code availability 31. Luco, R. F. et al. Regulation of alternative splicing by histone modifications. Custom scripts and associated files needed to reproduce all analyses described here are Science 327, 996–1000 (2010). provided together with the code at https://bitbucket.org/weghornlab/germline_intron_ 32. Li, F. et al. The histone mark H3K36me3 regulates human DNA mismatch exon_mutrate. Segments of the code are based on scripts shared by Frigola et al.17 to ensure repair through its interaction with MutSα. Cell 153, 590–600 (2013). comparability of the main analysis. 33. Li, C. & Luscombe, N. M. Nucleosome positioning stability is a modulator of germline mutation rate variation across the human genome. Nat. Commun. Received: 7 January 2020; Accepted: 2 June 2020; 11, 1–13 (2020). 34. Neale, B. M. et al. Patterns and rates of exonic de novo mutations in autism spectrum disorders. Nature 485, 242 (2012). 35. Goldmann, J. M. et al. Parent-of-origin-specific signatures of de novo mutations. Nat. Genet. 48, 935–939 (2016). 36. Halldorsson, B. V. et al. Characterizing mutagenic effects of recombination through a sequence-level genetic map. Science 363, eaau1043 (2019). References 37. An, J.-Y. et al. Genome-wide de novo risk score implicates promoter variation 1. Kimura, M. The Neutral Theory of Molecular Evolution (Cambridge in autism spectrum disorder. Science 362, eaat6576 (2018). University Press, 1983). 38. Yuen, R. K. et al. Whole genome sequencing resource identifies 18 new 2. Castle, J. C. SNPs occur in regions with less genomic sequence conservation. candidate genes for autism spectrum disorder. Nat. Neurosci. 20, 602–611 PLoS ONE 6, e20660 (2011). (2017). 3. Mikkelsen, T. S. et al. Initial sequence of the chimpanzee genome and 39. Sasani, T. A. et al. Large, three-generation human families reveal post-zygotic comparison with the human genome. Nature 437, 69 (2005). mosaicism and variability in germline mutation accumulation. eLife 8, 1–24 (2019). 4. Chamary, J. V., Parmley, J. L. & Hurst, L. D. Hearing silence: non-neutral 40. McLaren, W. et al. The Ensembl Variant Effect Predictor. Genome Biol. 17, evolution at synonymous sites in mammals. Nat. Rev. Genet. 7, 98–108 (2006). 122 (2016). 5. Smith, T. C., Arndt, P. F. & Eyre-Walker, A. Large scale variation in the rate of 41. Harrow, J. et al. GENCODE: the reference human genome annotation for the germ-line de novo mutation, base composition, divergence and diversity in ENCODE project. Genome Res. 22, 1760–1774 (2012). humans. PLoS Genet. 14, 1–30 (2018). 42. Derrien, T. et al. Fast computation and applications of genome mappability. 6. Lawrence, M. S. et al. Mutational heterogeneity in cancer and the search for PloS ONE. 7, e30377 (2012). new cancer-associated genes. Nature 499, 214–218 (2013). 43. Roadmap Epigenomics Consortium. Integrative analysis of 111 reference 7. Stamatoyannopoulos, J. A. et al. Human mutation rate associated with DNA human epigenomes. Nature 518, 317–329 (2015). replication timing. Nat. Genet. 41, 393 (2009). 44. ENCODE Project Consortium. An integrated encyclopedia of DNA elements 8. Schuster-Böckler, B. & Lehner, B. Chromatin organization is a major influence in the human genome. Nature 489, 57–74 (2012). on regional mutation rates in human cancer cells. Nature 488, 504 (2012). 9. Alexandrov, L. B. et al. Signatures of mutational processes in human cancer. Nature 500, 415–421 (2013). Acknowledgements 10. Alexandrov, L. B. et al. The repertoire of mutational signatures in human We thank especially Núria López-Bigas for her insight and suggestions during the cancer. Nature 578, 94–101 (2020). development of this work. We also thank her research group for sharing the software 11. Nik-Zainal, S. et al. Mutational processes molding the genomes of 21 breast code in all their contributions, especially their Frigola et al.17 paper. We would also like to cancers. Cell 149, 979–993 (2012). acknowledge Jakob Goldmann, Christian Gilissen, and Wendy Wong for sharing high- 12. Rahbari, R. et al. Timing, rates and spectra of human germline mutation. Nat. quality mutation data. We are grateful to David Castellano for helpful discussions and Genet. 48, 126–133 (2016). Vladimir Seplyarskiy for a critical reading of the manuscript. The project that gave rise to 13. Carlson, J. et al. Extremely rare variants reveal patterns of germline mutation these results received the support of a fellowship from “la Caixa” Foundation (ID rate heterogeneity in humans. Nat. Commun. 6, 3753 (2018). 100010434) with fellowship code LCF/BQ/DR19/11740019 (M.R.-G.); by the Ministerio 14. Aggarwala, V. & Voight, B. F. An expanded sequence context model broadly de Economía y Competitividad (Spain) (CGL2017-89160P to M.R.-G. and A.B.); and explains variability in polymorphism levels across the human genome. Nat. AGAUR (Generalitat de Catalunya) (2017SGR-1379 to A.B.). We acknowledge the Genet. 48, 349–355 (2016). support of the Spanish Ministry of Science and Innovation to the EMBL partnership, the 8 NATURE COMMUNICATIONS | (2020)11:3304 | https://doi.org/10.1038/s41467-020-17162-z | www.nature.com/naturecommunications NATURE COMMUNICATIONS | https://doi.org/10.1038/s41467-020-17162-z ARTICLE Centro de Excelencia Severo Ochoa and the CERCA Programme/Generalitat de Cata- Reprints and permission information is available at http://www.nature.com/reprints lunya (M.R.-G. and D.W.). Publisher’s note Springer Nature remains neutral with regard to jurisdictional claims in published maps and institutional affiliations. Author contributions M.R.-G. performed the analyses. M.R.-G. and D.W. designed the analyses. M.R.-G., D.W., and A.B. wrote the manuscript. D.W., S.C., and A.B. supervised the study. D.W. and A.B. conceived the study. All authors read and approved the paper. Open Access This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give Competing interests appropriate credit to the original author(s) and the source, provide a link to the Creative The authors declare no competing interests. Commons license, and indicate if changes were made. The images or other third party material in this article are included in the article’s Creative Commons license, unless Additional information indicated otherwise in a credit line to the material. If material is not included in the Supplementary information is available for this paper at https://doi.org/10.1038/s41467- article’s Creative Commons license and your intended use is not permitted by statutory 020-17162-z. regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this license, visit http://creativecommons.org/ Correspondence and requests for materials should be addressed to D.W. or A.B. licenses/by/4.0/. Peer review information Nature Communications thanks Edward Hollox and Wendy Wong for their contribution to the peer review of this work. Peer reviewer reports are © The Author(s) 2020 available. NATURE COMMUNICATIONS | (2020)11:3304 | https://doi.org/10.1038/s41467-020-17162-z | www.nature.com/naturecommunications 9