Discovery of an unusual high number of de novo mutations in sperm of older men using duplex sequencing

De novo mutations (DNMs) are an important player in heritable diseases and evolution. Of particular interest are highly recurrent DNMs associated with congenital disorders that have been described as selfish mutations expanding in the male germline, thus becoming more frequent with age. Here, we have adapted duplex sequencing (DS), an ultra-deep sequencing method that renders sequence information on both DNA strands; thus, one mutation can be reliably called in millions of sequenced bases. With DS, we examined ∼4.5 kb of the FGFR3 coding region in sperm DNA from older and younger donors. We identified sites with variant frequencies of 10−4 to 10−5, with an overall mutation frequency of the region of ∼6×10−7. Some of the substitutions were re-current and were found at a higher variant frequency in older donors than in younger ones, or exclusively, in older donors. Also, older donors harbored more mutations associated with congenital disorders. Other mutations were present in both age groups suggesting that these might result from a different mechanism (e.g., post-zygotic mosaicism). We also observed that independent of age, the frequency and deleteriousness of the mutational spectra was more similar to COSMIC than to gnomAD variants. Our approach is an important strategy to identify mutations that could be associated with a gain-of-function of the receptor tyrosine kinase activity, with unexplored consequences in a society with delayed fatherhood.

There are certain de novo mutations (DNMs) that are highly recurrent with mutation rates orders of 2 magnitude higher than genome average. These mutations have been discovered because of their 3 association with congenital disorders. Moreover, these mutations have several other associated 4 characteristics (reviewed in ( properties (activating mutations); they occur exclusively in the male germline; and older men have a 7 higher probability of having an affected child than younger males, known as the paternal age-effect 8 (PAE) described already decades ago (Risch et al. 1987;Crow 2000;Crow 2012). In the literature, 9 these mutations have been termed as PAE-mutations, selfish, or RAMP (recurrent, autosomal 10 dominant, male biased, and paternal age effect) mutations, as reviewed in (Arnheim and Calabrese 11 2009; Goriely and Wilkie 2012; Arnheim and Calabrese 2016). 12 Over the last decades, it was shown that the high incidence levels and steep increase with age of 13 these PAE-mutations are not solely due to errors occurring during the continuous replicative process 14 of spermatogenesis (Risch et al. 1987;Crow 2000;Crow 2012). Instead, it was suggested that these 15 behave like driver mutations, well known in cancer, which are mutations that promote their own 16 clonal expansion. In particular, PAE-mutations have been described in genes such as FGFR2, FGFR3, The clonal expansion of some PAE-mutations in the testis with age also results in more mutant 4 sperm in older individuals, and an increased incidence of the associated congenital disorder with 5 paternal age, as observed for example for achondroplasia (ACH), a growth disorder with a defect in 6 FGFR3 signaling (Risch et al. 1987 In spite of the importance of PAE-mutations in the male germline which is highlighted by their high 11 incidence, increased frequency with paternal age, and phenotypes potentially leading to early or 12 late-onset disorders in children of older men, we know very little about this mutagenic mechanism. could harbor as yet unknown mutations that could be expanding with paternal age but have gone 22 undetected so far, though, with important health consequences in children of older fathers. To 23 overcome this, a more recent study used a sequencing strategy to discover prospective driver 24 mutations in genes of the RTK-RAS pathway by enriching for mutant SSCs identified in testes biopsies 25 by a high RTK activity that identified 61 mutations at a frequency of 0.06% or larger (Maher et al. 1 2018) . 2 In this work, we further developed the discovery of prospective RTK mutations in the male germline 3 to include ultra-low frequency variants (with a variant frequency of <2 x 10 -5 ). For this purpose, we 4 used sperm since this biological material is easy to sample and provides information on whether 5 mutations are passed on throughout the different maturation and differentiated stages of the male 6 germline. Accurately identifying ultra-low frequency mutations by sequencing is still technically very 7 challenging, especially if no enrichment of mutants is possible (e.g., high RTK signaling). To date, only 8 a few methods can achieve the required sensitivity and accuracy. One of these methods is duplex 9 sequencing (DS), a strategy that organizes sequence reads derived from the same DNA molecule into 10 families with information on the forward and reverse strands (Schmitt et al. 2012). Given the 11 excessive sequencing depth required in DS, in this work we focused on the coding region of has been categorized as an oncogene with many missense mutations and high oncogenic score 16 (Tomasetti and Vogelstein 2015). 17 We sequenced the coding regions of the FGFR3 gene at an average coverage depth of ~17,000x (up 18 to ~38,000x) in young and an old donor pools. We identified 75 unique variants at 72 different 19 genomic positions in the exonic FGFR3 sequence (exon 3-15 without exon 11); of those, 20 have 20 been reported in other databases to be associated with congenital disorders (Human Gene Mutation 21 database (HGMD)) and/or tumors (Catalogue Of Somatic Mutations In Cancer (COSMIC)); 34 DNMs 22 have never been reported before (Table S1). Some DNMs were found exclusively or at a higher 23 frequency in the older donor set. Interestingly, we also identified many higher frequency variants 24 independent of the age of the donors, suggesting that these might not be testis-exclusive mutations,

6
In order to confirm that variants called by our pipeline represent true mutations, six of the 7 mutations identified with DS were measured with other ultra-sensitive detection methods, namely 8 bead emulsion amplification (BEA) as described in (Shinde et al. 2013) or droplet digital PCR (ddPCR); 9 see Supplemental Table S2. BEA is an in-house digital PCR method based on the amplification of 10 single molecules on magnetic beads within an emulsion (Boulanger et al. 2012;Shinde et al. 2013). 11 The ddPCR technology (Bio-Rad Laboratories) is based on the partitioning of PCR reactions into 12 thousands of individual reactions within water-in-oil droplets that are read out in an automated 13 droplet flow-cytometer after the PCR reaction is concluded (Hindson et al. 2011); see Methods for 14 further information. Our DS measurements were similar among methods; albeit, given the low 15 number of events captured in DS the confidence intervals are larger for DS than for ddPCR. 16 Moreover, in BEA and ddPCR more starting molecules were used (~300,000 genomes) making these 17 two methods more accurate for the quantification of very low mutation frequencies compared to 18 DS. A Fischer's exact test for all possible comparisons did not reveal any significant difference 19 between the measured mutation frequencies. 20 As a further demonstration of the high sensitivity and accuracy of DS, we carried out a mixing 21 experiment. For this purpose, different genomic DNA obtained from Coriell Cell Repositories 22 (Camden, NJ) each carrying a point mutation (c.742C>T, c.746C>G, c.749C>G, and c.1620C>A) were 23 spiked into DNA from a younger donor at ratios from 1/10 to 1/10,000 in steps of one order of 24 magnitude ( Figure 2). With this orthogonal assay, we demonstrated that the variant frequency 25 detected with DS was equivalent to the expected input of the 10-fold dilutions. Note that the 26 starting input amounts varied due to pipetting errors and the occasional inaccuracy of 27 spectrophotometric measurements used to establish the initial DNA concentration. Fourteen out of 28 16 variants (four variants per library) were identified. Two replicates of the 1:10,000 dilutions were 1 not detected likely due to sampling drop-out at the Poisson distribution level. The Pearson's 2 correlation coefficient between observed and expected input dilutions of the 14 measurements was 3 very high (R 2 = 0.96, p-value = 9.8 x 10 -10 ) demonstrating the power of DS to measure ultra-rare 4 variants ( Figure 2).

13
Substitutions detected in the coding region of FGFR3 in sperm DNA 14 In total, we sequenced 12 libraries, each prepared from sperm DNA pooled from five healthy 15 individuals, mainly of European ancestry. Six libraries featured younger donors with ages ranging 16 from 26 to 30 years, and the remaining six featured older donors with ages ranging from 49 to 59 17 years; see Supplemental Table S3 and Supplemental Table S4. The distribution of donor ages in the 18 young and old pool was significantly different (Supplemental Figure S2). Further, in spite of 19 differences in the library preparation protocols, we did not observe differences in the overall 20  The mutants detected in this study had a VAF (estimated as the ratio of the alternate allele divided 9 by the reference allele or mutation count per depth of coverage at a given genomic position) 10 between 1.5 x 10 -5 to 7.9 x 10 -4 ( Figure 3 and Supplemental Table S1). Of the 75 discovered variants, 11 34 of them have not been reported before in any public databases, 13 were unique substitutions 12 associated with tumors (mainly bladder, skin, and multiple myeloma, which are the most common 13 cancer associations within FGFR3), and 32 were substitutions reported in the gnomAD database 14 (likely viable mutants). We also calculated the overall mutation frequency, which is an estimate of the total number of 1 variants per number of sequenced base pairs, and can be used to compare overall mutation 2 frequencies with other methods, or differences between donor groups or domains ( Figure 4). Note 3 that this mutation frequency or substitution/base pair represents de novo mutation (DNM) events, 4 and was estimated as described in Methods. As shown in Figure 4A, the mutation frequency of the 5 exonic region of FGFR3 is ~4.58 x 10 -7 , which is one order of magnitude higher than measured  Figure S4), suggesting a unique mutagenic 9 mechanism in this gene explored in more detail in the next sections. Interestingly, the number of 10 missense mutations was significantly higher compared to synonymous, stop-gained or splice region 11 variants with approximately two-thirds of the missense variants being deleterious based on the SIFT 12 score (Supplemental Figure S5). Also, deleterious mutations or mutations with unknown effect based 13 on the SIFT score were overall more frequent than tolerated/benign variants ( Figure 4A). 14 A closer look into the younger and older sperm pools did not show differences in mutation 15 frequencies between these two pools ( Figure 4B); although, the mutation frequency in the younger 16 donor pool was slightly higher than in the older donor pool (5.0 x 10 -7 vs 4.1 x 10 -7 , respectively). 17 However, some interesting trends can be observed in the older donor pools such as a higher 18 mutation frequency in the Immunoglobulin-like domains I, II and III (IgI-III), but a lower frequency in 19 the TK domain, as well as a higher frequency of substitutions with a reported germline disorder. 20 These trends suggest that older donors might harbor more substitutions associated with a 21 phenotype that follow a paternal age effect. 22 across the three main domains (IgI-III n=38, TM_Inter-domain n=21, and Tyrosine Kinase n=33); yet 1 normalized by mutation frequency, the tyrosine kinase (TK) domain had the lowest mutation 2 frequency. Interestingly, this domain had the largest proportion of deleterious mutations (Figure 4C  3 and Supplemental Table S1). The lower frequencies detected in the TK domain might indicate that 4 the high deleteriousness in this domain is not tolerated in the male germline. Similarly, in a study 5 analyzing all possible substitutions in codon p.K650 found that highly activating substitutions 6 associated with cancer (e.g., seborrheic keratoses) were underrepresented in sperm compared to 7 less "oncogenic" substitutions (e.g., p.K650E and p.K650T resulting in thanatophoric dysplasia type II 8 (TDII) and familial acanthosis nigricans, respectively) (Goriely et al. 2009). 9 We captured in the transmembrane (TM) and surrounding regions the highest mutation frequency 10 of DNM associated with germline disorders ( Figure 4C). Also, the proportion of missense 11 substitutions compared to synonymous substitution was highest in this domain. Interestingly, in this 12 domain we also observed substitutions with the highest VAF (e.g., p.Y373C, p.G380R, and p.E447K); 13 almost 50% higher in magnitude than the median of the other variants (~4 x 10 -5 ) measured in 14

Mutational mechanism 9
Aiming to explore the possible mutational mechanisms behind our variants, we tested if the 10 observed variants have a selective advantage linked to the clonal expansion of driver mutations. In 11 evolutionary genetics, the measure of dN/dS (ratio of non-synonymous versus synonymous 12 substitutions per non-synonymous or synonymous sites, respectively) is indicative of positive, 13 neutral or negative selection. The interpretation of this ratio is that a number close to the reference 14 value (usually one) means no selection; a ratio below the reference means negative selection, and a 15 ratio above the reference implies positive selection (revised in (Nielsen 2005)). The reference value 16 is calculated by the distribution of dN/dS occurring at random with respect to nonsynonymous vs. We further analyzed the mutational spectra, transcriptional bias, and mutational signatures in our 26 data (Supplemental Figure S6A). Based on the cosine similarity of the FGFR3 mutation spectra, there 27 is a greater similarity with the COSMIC variants than with variants reported in gnomAD (cosine 28 similarity of 0.93 vs 0.87, respectively) shown in Supplemental Table S7 and Supplemental Figure reported in tumors for FGFR3 (COSMIC) than compared to random DNMs captured in the general 1 population (gnomAD). Also, note that we do not observe differences in the mutational spectra 2 between younger and older donor groups (Supplemental Figure S8) with both groups showing a high 3 cosine similarity score (0.92). The lack of a strong strand bias in the substitution types also indicates 4 that transcription cannot explain our data (Supplemental Figure S6B). Finally, the mutational 5 signature that compares the variants also in the context of the 5ʹ and 3ʹ base adjacent to the variant 6 shows that the main substitutions were strong to weak transitions in the context of CpG sites (VCG). 7 Further, the comparison of our observed mutational signature with COSMIC indicate an overlap of 8 ~51% of the variants explained by pattern SBS24 (Aflatoxin exposure), 30% by SBS1 (spontaneous or 9 enzymatic deamination of 5-methylcytosine to thymine) and 18.6% by SBS5 (unknown) (see 10 Supplemental Figure S9). However, it is likely that driver genes have each their own mutational 11 signature and needs to explored further. 12 Overall, the higher overall DNM mutation frequencies, the elevated VAFs at individual positions, and 13 the trends of dN/dS for positive selection indicate that substitutions captured in FGFR3 of sperm 14 cannot be explained by a mutagenic mechanism similar to genome-wide DNM, but instead might be 15 the result of deleterious/activating mutations that promote the clonal expansion of mutant cells. 16 However, a larger dataset of variants and more sequencing data might show the underlying trends 17 with better power. 18

Testis specific clonal expansion or post-zygotic mosaic? 19
It is possible that the captured variants are expanding clonally with age only in the testis. It is possible that variants occurring only in the older donor pool might be testis-exclusive mosaics 9 growing larger over time (see Table 1 whereas, in younger donors libraries, they were found only once (VAF of 1.6 x 10 -5 ). This amino acid 13 substitution is known to be associated with hypochondroplasia, which also suggests that this variant 14 could expand in the male germline with age. A more in-depth analysis of more sperm samples or the 15 distribution within testis with a less expensive approach will provide more details in this regard. 16

Younger Donors
Older

6
Albeit challenging, PZM (labelled as PZM-candidates in Table 1) might be distinguished from testis-7 specific mosaics since they occur at a relatively high VAF already in young sperm donors or testis;   Figure 6A, PAE candidates 2 had higher CADD scores similar to COSMIC variants implying that some mutations expanding in 3 sperm can be as deleterious as tumor-associated variants. Interestingly, PZM candidates and 4 gnomAD variants had the lowest CADD scores. This suggests that at least part of these PAE candidate 5 substitutions might have a stronger functional impact than most of the variants currently 6 segregating in the population. Whether these are associated with a paternal-age effect and might 7 have a deleterious effect like a congenital disorder, a late-onset disease (e.g., cancer), or instead a 8 normal phenotype, is not known. 9 Using the same grouping of mutations, we investigated the association with germline disorders, 10 tumors or gnomAD ( Figure 5B). PAE-candidates contained the largest proportion of variants 11 associated with germline disorders and tumors in comparison to PZM-candidates highlighting further 12 differences between these two groups. Variants found in the general population (gnomAD) were 13 associated more often with PZM-candidates. 14 15 16   targets. Further enrichment of targeted regions is performed by two rounds of hybridization capture 21 followed by PCR, as previously described (Schmitt et al. 2015). This approach rendered an 22 unprecedented DCS coverage obtained for every library. Supplemental Figure S3 and S10 shows 23 variation across different libraries, which can be explained by the different parameters used (e.g., input DNA, PCR cycles, sequencing read number). We also observed variable DCS coverage values 1 between different targeted regions within each library that can be attributed to differences in 2 performance of each restriction enzyme, amplification efficiency of each amplicon, and capture oligo 3 efficiency Supplemental Figure S10. The lower DCS coverage obtained in the central parts of the 4 restriction fragments is originated by trimming low quality read ends. However, with our DS 5 approach we achieved a median DCS coverage of ~11,000 allowing to call mutations at a frequency 6 of at least 1/10,000. 7

Mutation accumulation in FGFR3 8
The higher incidence of mutations in FGFR3 increasing with age in the paternal germline was first In this study, we were able to detect 7 variants (6 amino acid substitutions) that we classified as PAE-20 candidates and are likely expanding in the male germline given the higher variant frequencies in 21 older donor pools (Table 1). Three out of six of these amino acid substitutions have been associated 22 with PAE disorders and are suspected to increase with paternal age. One of these PAE-candidates is 23 the c.1118A>G (p.Y373C) measured at a VAF of ~3 x 10 -4 ; one of the highest frequencies measured in 24 our dataset. This variant is causal for the autosomal dominant congenital disorder TDI, also described as a PAE-disorder (Rousseau et al. 1996b;Wilcox et al. 1998). Similarly, the c. respectively. These two variants were also identified in a study measuring prospective mutations 12 enriched in testis biopsies with a high RTK activity (Maher et al. 2018). The p.N540K amino acid 13 substitution is related to the skeletal dysplasia associated with a paternal age effect, the autosomal 14 dominant congenital disorder hypochondroplasia (Bellus et al. 1995;Rousseau et al. 1996a). The 15 higher frequency detected in older donors coupled with an association to a PAE-disorder suggests 16 that the p.N540K substitution might be also a testis-specific mosaic expanding with age. 17 Apart from the few characterized PAE-candidates, our study identified further variants that could be 18 potentially expanding with paternal age and have not been described before. These occur mainly 19 within or downstream of the TM domain (e.g., p.T394T, p.R399H, p. R421Q, p.A429E) and are 20 associated with a high deleteriousness and tumors ( Figure 6). Moreover, they are viable since they 21 have been detected in the general population. Thus, our approach is an important tool to identify 22 prospective mutations expanding with paternal age, but more evidence is required to demonstrate a 23 selective advantage of these variants in the male germline and whether these form testis-specific 24 mosaics with age. 25

Testis-exclusive mosaics versus post-zygotic mosaicism 1
Driver mutations and their association with congenital disorders gain a special importance in 2 industrial societies in which parenthood is delayed. The full understanding of driver mutations 3 accumulating in the male germline is key to assess the risks and the impact on future generations. 4 DNMs can occur in adult spermatogonia that are positively selected and are therefore restricted to 5 the seminiferous tubules forming testis-exclusive mosaics. These are often associated with PAE-6 disorders; the numbers of mutant sperm increase with paternal age and they can be transmitted and We captured several substitutions with VAF over 1 x 10 -5 . These are viable variants found in the 10 human population (gnomAD) and some of them have also been associated with tumors (COSMIC). 11 Since they were not detected in multiple libraries and/or their frequency was not higher in older 12 donors, we considered that these might come from a different source: DNMs occurring post- Moreover, we noticed a difference in the association with germline disorders and tumors between 18 PZM-candidates and PAE-candidates; although, our mutation deleteriousness analysis did not show 19 any significant differences. 20 The expansion of PZM variants with age is not yet fully understood. In theory, the fitness of these 21 variants conferred to the host cell should dictate their development over time. Therefore, some 22 mosaic variants might expand clonally, and others might remain at low frequencies. But similar to 23 gonadal-exclusive mosaics, these would also be associated with a high recurrence risk in the 24 offspring. Depending on how early in development they occur, the range and expansion might vary occurred before the separation of the somatic and gonadal lineages, then these would be equally 3 likely to be found at high levels in younger sperm donor pools or in more pieces of testis biopsies 4 (see such an example in (Maher et al. 2018)). Yet, in blood, the higher number of aberrant clonal 5 expansions observed in older individuals is not fully explained by this theory. In the germline, driver 6 mutations expand clonally with age and its frequency is higher in older individuals, but in the case of 7 post-zygotic mutations occurring early in development, this is currently unknown. As reviewed in 8 development. An expansion over time is observed in some tissues, but the opposite can occur in 10 others due to a deleterious effect and subsequent decrease in relative fitness in that particular 11 tissue. 12

Effect of mutations on the RTK 13
To date the majority of known PAE mutations affect the RTK or components of its downstream RAS- Thus, the observed PAE with consequences in neurological, heart, or cancer development, might 5 also be partially the result of DNMs in driver genes that lead to abnormalities in the signaling 6 pathway (e.g., RTK-RAS). However, in order to better understand these effects, we need to further 7 investigate DNMs in driver genes in the male germline, how these DNMs affect the self-propagation 8 of the mutation, and the nature of their phenotypic consequences. 9 with substantial modifications including an initial restriction enzyme digest to pre-select the target 16 regions. We used in the different libraries some modifications to our main protocol either starting 17 with distinct input DNA amounts, different adapters or amplification strategies. Supplemental Table  18 S3 summarizes the differences between each of our library preparation protocols. Genomic DNA 19 (500-5000 ng, Supplemental Table S3) was subject to an overnight restriction enzyme digest that 20 targeted 5 regions of the FGFR3 gene (Supplemental Table S8). Resulting fragments are expected to 21 be between ~300 to ~600 bp. A double size selection was performed using SPRIselect beads 22 (Beckman Coulter) in order to exclude fragments outside of this size range. We used 3 different 23 adapter synthesis protocols to produce Adapter 1, Adapter 2 and Adapter 3. The details of the 1 protocols used to produce these adapters can be found in SM. Supplemental Figure S11 shows the 2 structural differences between adapters: Adapter 1 is identical to the original open-hairpin structure, 3 whereas adapters 2 and 3 have a closed loop hairpin structure that is opened up after the ligation 4 step by the cleavage of the uracil with the USER enzyme mix (New England Biolabs). It is suggested 5 that closed loops adapters reduce adapter dimers during ligation (New England Biolabs). Structurally, 6 adapter 3 is similar to adapter 2, but is 9 bp longer and contains a phosphorothioate bond before the 7 T-overhang on the 3' end such that removal of this overhang is reduced. 8 Size selected genomic DNA was end-repaired and A-tailed using the NEBNext® Ultra™ II End 9 Repair/dA-Tailing Module (New England Biolabs) according to the manufacturer's instructions 10 described in SM. DNA fragments with A-3' end overhangs were then ligated to the DS adaptors with 11 T-3' overhangs using the NEBNext® Ultra™ II Ligation Module (New England Biolabs) following the 12 manufacturer's instructions and then purified by 0.8 volumes of AMPure XP beads (Beckman 13 Coulter). 14 Amplification of ligated DNA was executed using KAPA HiFi HotStart ReadyMix (KAPA Biosystems). 15 Except for 4 libraries, a strategy of 12 cycles of single primer extensions was adopted before the PCR 16 reaction in order to get multiple copies of the initial genomic template to improve the 17 representation of both forward and reverse SSCS and prevent the exponential amplification of 18 templates resulting in very large family sizes. Reaction volumes and PCR conditions are described in 19 SM and Supplemental Table S9. Primer sequences are shown in Supplemental Table S10. After the 20 extension/PCR step, PCR products were purified with AMPure XP beads followed by 2-3 rounds of 21 targeted capture steps to enrich further for templates of interest. A third targeted capture was done 22 for two of the libraries using the same procedure as the second one. The number of cycles of both 23 post-capture PCRs varied across libraries (Supplemental Table S3). Pools were verified using Libraries were diluted and pooled for sequencing according to the concentration of dsDNA measured 1 with DeNovix® dsDNA High Sensitivity. Finally, the libraries were subject to further quantification 2 with the KAPA Library Quantification kit (KAPA Biosystems). Sequencing reactions were performed 3 on the MiSeq Illumina platform using the kit MiSeq Reagent v3 600 cycles (Illumina) at the Center for 4 Medical Research of the Johannes Kepler University and at the VBCF NGS Unit. 5

Data processing and variant filtering 6
FastQ files were analyzed in Galaxy (in both public (usegalaxy.org) and private (zusie.jku.at) servers)) 7 according to a DS specific pipeline that includes an error correction tool (Stoler et al. 2020 approach (two-stage set-up method of Benjamini, Krieger and Yekutieli). Germline association was 3 investigated using ClinVar data extracted from wANNOVAR output. Tumor association was 4 investigated by consulting the presence or absence of variants in the extracted COSMIC data. 5

Mutation analysis 6
Mutation frequencies vs VAF 7 The VAF measures the mutant count on a single position and is estimated as the alternate allele 8 count divided by the coverage (reference allele) on that position. The mutation frequency considers 9 all positions of a specific region (region size) and represents the ratio of the number of different 10 variants (variant count) per number of sequenced or analyzed nucleotides (given by the formula 11 below). The latter is estimated as the sum of the mean coverage of a given library times the size of 12 the sequenced/analyzed region (e.g., the region size was adjusted per domain or upstream-

Mutational spectra 18
After counting each substitution type (e.g., A>C) we either represented the mutational spectra by 19 the measure of relative count (variant divided by the number of all variants) or as a mutation 20 frequency. The latter is calculated as described before, except that here we normalize for the 21 sequenced nucleotides of the reference allele, which is estimated by the frequency of the reference 22 allele in the targeted regions multiplied by all sequenced nucleotides. 23 For the mutational spectra of the databases COSMIC v94 and gnomAD v3.1.1, all variants except for 3 indels of the targeted regions were downloaded. Next, we compared them to the spectra of the 4 duplex sequencing data using the measure of relative counts. For comparisons within duplex 5 sequencing data (e.g., Younger and Older group), we used the measure of mutation frequencies to 6 account also for differences in the coverage. We also categorized the mutational spectra after the 7 (un)transcribed strand and the mutational signature within the trinucleotide context (includes the 5ʹ 8 and 3ʹ adjacent nucleotide to the variant) with the tools SigProfilerMatrixGenerator and 9 SigProfilerPlotting

Cosine similarities 13
Comparing two mutational spectra is frequently done by estimating the cosine similarity. The 14 expected cosine similarity for a reference spectrum is calculated by bootstrapping (Abascal et al. 15 2021). In short, we randomly sampled n mutations, where n is the number of mutations in the query 16 spectra, from the reference spectra and calculated the cosine similarity between the sample and the 17 original reference spectra (1,000 iterations). The mean and 95% confidence intervals were 18 estimated from the bootstrapped cosine similarities. 19 The rate of nonsynonymous vs. synonymous mutations 20 To test for selection at protein-coding sequences, we calculated the dN/dS ratio analogous to what 21 was described in (Nielsen 2005) and (Li et al. 2015). The value of dN represents the number of 22 nonsynonymous mutations per sites in the sequence where a mutation would produce a 23 nonsynonymous difference. The value of dS is the number of synonymous mutations per sites in the 1 sequence where a mutation would produce a synonymous difference. The numbers of 2 nonsynonymous and synonymous sites was calculated using the Nei-Gojobori method (Nei and 3 Gojobori 1986). Under neutrality, the dN/dS ratio is expected to be equal to 1, a ratio <1 is 4 suggestive of purifying selection, a ratio >1 is suggestive of positive selection. In order to assess the 5 probability of a random hN/hS ratios, we generated a distribution of dN/dS ratios that would be 6 expected if the observed spectrum of mutational changes were occurring at random with respect to 7 nonsynonymous vs. synonymous sites (1,000 iterations). 8

9
Raw sequencing data is available at the NCBI Sequence Read Archive (BioProject ID: PRJNA684907). 10 No Competing interests 11 The authors declare no competing interests. 12 13

14
Open access funding was provided by the Austrian Science Fund (FWF). This work was also supported 15 by the "Austrian Science Fund" (FWFP30867000), by the FWF Doctoral College "NanoCell" 16 (FWFW1250) and the European Regional Development Fund (REGGEN ATCZ207). Further funding was 17 provided by the Linz Institute of Technology (LIT213201001). We would like to thank Nick Stoler, Anton