Somatic mutation rates scale with lifespan across mammals

Data reporting

No statistical methods were used to predetermine sample size. The experiments were not randomized and the investigators were not blinded to allocation during experiments and outcome assessment.

Ethics statement

All animal samples were obtained with the approval of the local ethical review committee (AWERB) at the Wellcome Sanger Institute and those at the holding institutions.

Sample collection

We obtained colorectal epithelium and skin samples from a range of sources (Supplementary Table 1). For comparability across species an approximately 1-cm biopsy of the colorectal epithelium was taken from the terminal colon during necropsy. All necropsies occurred as soon as possible post-mortem to minimize tissue and DNA degradation. Tissue samples taken later than 24 h post-mortem typically showed extensive degradation of the colorectal epithelium, making the identification of colorectal crypts challenging. These samples were also associated with poor DNA yields and so were not included in the study. Sampled tissue was fixed in PAXgene FIX (PreAnalytiX, Switzerland), a commercially available fixative, during the necropsy. After 24 h in the fixative at room temperature, samples were transferred into the PAXgene STABILIZER and stored at –20 °C until further processing.

Sample processing

Samples were processed using a workflow designed for detection of somatic mutations in solid tissues by laser-capture microdissection (LCM) using low-input DNA sequencing. For a more detailed description see the paraffin workflow described in another study29. In brief, PAXgene-fixed tissue samples of the colorectal epithelium were paraffin-embedded using a Sakura Tissue-Tek VIP tissue processor. Sections of 16 µm were cut using a microtome, mounted on PEN-membrane slides and stained with Gill’s haematoxylin and eosin by sequential immersion in the following: xylene (two minutes, twice), ethanol (100%, 1 min, twice), deionized water (1 min, once), Gill’s haematoxylin (10 s, once), tap water (20 s, twice), eosin (5 s, once), tap water (20 s, once), ethanol (70%, 20 s, twice) and xylene or Neo-Clear, a xylene substitute (20 s, twice).

High-resolution scans were obtained from representative sections of each species. Example images are shown in Fig. 1a, Extended Data Fig. 2. Individual colorectal crypts were isolated from sections on polyethylene naphthalate (PEN) membrane slides by LCM with a Leica LMD7 microscope. Haematoxylin and eosin histology images were reviewed by a veterinary pathologist. For some samples we also cut a section of muscle tissue from below the colorectal epithelium of the section to use as a germline control for variant calling (Supplementary Table 2). Pre- and post-microdissection images of the tissue were recorded for each crypt and muscle sample taken. Each microdissection was collected in a separate well of a 96-well plate.

Crypts were lysed using the Arcturus PicoPure Kit (Applied Biosystems) as previously described8,29. Each crypt then underwent DNA library preparation, without a quantification step to avoid loss of DNA, following a protocol described previouslyl29. For some animals, a PAXgene fixed bulk skin biopsy was used as the germline control. For these skin samples, DNA was extracted using the DNeasy Blood & Tissue Kit (Qiagen).

Library preparation and sequencing

Libraries from microdissected samples were prepared using enzymatic fragmentation, adapter ligation and whole-genome sequencing following a method described previously29. Libraries from skin samples were prepared using standard Illumina whole-genome library preparation. Samples were multiplexed and sequenced using Illumina XTEN and Novaseq 6000 machines to generate 150 base pair (bp) paired-end reads. Samples were sequenced to around 30× depth (Supplementary Table 2).

Sequence read alignment

For each species, sequences were aligned to a reference assembly (Supplementary Table 2) using the BWA-MEM algorithm59 as implemented in BWA v.0.7.17-r1188, with options ‘-T 30 -Y -p -t 8’. The aligned reads were sorted using the bamsort tool from the biobambam2 package (v.2.0.86; gitlab.com/german.tischler/biobambam2), with options ‘fixmates=1 level=1 calmdnm=1 calmdnmrecompindetonly=1 calmdnmreference=<reference_fasta> outputthreads=7 sortthreads=7’. Duplicate reads were marked using the bammarkduplicates2 tool from biobambam2, with option ‘level = 0’.

Variant calling

Identification of somatic substitutions and short indels was divided into two steps: variant calling, and variant filtering to remove spurious calls (see ‘Variant filtering’ below). For human colorectal crypts, we obtained previously sequenced and mapped reads from a study in which colorectal crypts were isolated by LCM8, and processed them using the sample variant calling and filtering process that was applied to the non-human samples.

Substitutions were identified using the cancer variants through expectation maximization (CaVEMan) algorithm60 (v.1.13.15). CaVEMan uses a naive Bayesian classifier to perform a comparative analysis of the sequence data from a target and control sample from the same individual to derive a probabilistic estimate for putative somatic substitutions at each site. The copy number options were set to ‘major copy number = 5’ and ‘minor copy number = 2’, as in our experience this maximizes the sensitivity to detect substitutions in normal tissues. CaVEMan identifies and excludes germline variants shared in the target (colorectal crypt) and matched normal (skin or muscle tissue) samples, and produces a list of putative somatic mutations that are present only in the target sample. CaVEMan was run separately for each colorectal crypt, using either bulk skin or muscle microdissected from the sample colorectal biopsy as the matched normal control (Supplementary Table 2). For two human donors for whom an alternative tissue was not available, a colonic crypt not included as a target sample was used as the matched normal control.

Indels were identified using the Pindel algorithm61 (v.3.3.0), using a second sample from the same individual as a matched control. The indel calls produced by Pindel were subsequently re-genotyped using the vafCorrect tool (https://github.com/cancerit/vafCorrect), which performs a local sequence assembly to address alignment errors for indels located at the end of sequence reads, and produces corrected counts of sequence reads supporting the indel and corrected estimates of variant allele fraction (VAF; the fraction of reads supporting the alternate allele at the variant site).

Variant filtering

A number of post-processing filters were applied to the variant calls to remove false positives (Supplementary Fig. 1a, b).

Quality flag filter

CaVEMan and Pindel annotate variant calls using a series of quality flags, with the ‘PASS’ flag denoting no quality issues affecting the call60,61. Variant calls presenting any flag other than ‘PASS’ were discarded.

Alignment quality filter

Variants were excluded if more than half of the supporting reads were clipped. The library preparation methods create short insert size libraries that can result in reads overlapping. To avoid the risk of double counting mutant reads we used fragment-based statistics. Variants without at least four high-quality fragments (alignment score ≥ 40 and base Phred quality score ≥ 30) were excluded. Variants were excluded if reads supporting the variant had a secondary alignment score that was greater than the primary alignment score. This filter was not applied to indel calls.

Hairpin filter

To remove variants introduced by erroneous processing of cruciform DNA during the enzymatic digestion, we applied a custom filter to remove variants in inverted repeats29. This filter was not applied to indel calls.

Chromosome and contig filter

For species with chromosome-level assemblies, we discarded variants located in non-chromosomal contigs, including the mitochondrial genome (calling of mitochondrial variants is described in the section ‘Mitochondrial variant calling and filtering’). For males, variants on the Y chromosome were excluded for species in which the Y chromosome was annotated in the assembly.

N-tract and contig-end filter

To reduce artefactual calls due to read misalignment, we discarded variants located within 1 kb of a tract of 50 or more consecutive N bases in the reference assembly, as well as variants within 1 kb of the start or end of a contig (this implies discarding all variants in contigs shorter than 2 kb).

Sequencing coverage filter

A sample-specific read depth filter was designed to exclude sites with coverage above the 99th coverage percentile in the sample or its matched normal control, as well as sites with coverage of less than 10× in the sample or its matched normal control.

Allelic strand bias filter

We discarded variants without any supporting reads on either the forward or the reverse strand.

Indel proximity filter

We discarded variants for which the total number of reads supporting the presence of an indel within 10 bp of the variant was more than three times larger than the number of reads supporting the presence of the variant. This filter was not applied to indel calls.

Spatial clustering filter

Visual assessment of variant calls and mutational spectra showed spatially clustered variants to be highly enriched for artefacts. Therefore, we discarded groups of two or more variants located within 1 kb of each other.

Beta-binomial filter

For each crypt, an artefact filter based on the beta-binomial distribution was applied, which exploits read count information in other crypts from the same individual. More specifically, for each sample, we fitted a beta-binomial distribution to the variant allele counts and sequencing depths of somatic variants across samples from the same individual. The beta-binomial distribution was used to determine whether read support for a mutation varies across samples from an individual, as expected for genuine somatic mutations but not for artefacts. Artefacts tend to be randomly distributed across samples and can be modelled as drawn from a binomial or a lowly overdispersed beta-binomial distribution. True somatic variants will be present at a high VAF in some samples, but absent in others, and are hence best captured by a highly overdispersed beta-binomial. For each variant site, the maximum likelihood estimate of the overdispersion factor (ρ) was calculated using a grid-based method, with values ranging between 10−6 and 10−0.05. Variants with ρ > 0.3 were considered to be artefactual and discarded. The code for this filter is based on the Shearwater variant caller62. We found this to be one of the most effective filters against spurious calls (Supplementary Fig. 1b).

Minimum VAF filter

For each sample, we discarded variants for which the VAF was less than half the median VAF of variants passing the beta-binomial filter (see above) in that sample.

Maximum indel VAF filter

For each sample, we discarded indels that presented a VAF of greater than 0.9, as such indels were found to be highly enriched in spurious calls in some species. This filter was not applied to substitution calls.

To validate our variant calling strategy, we used LCM to microdissect two sections from the same mouse colorectal crypt. We expected to detect a high fraction of shared somatic variants in these two sections, as their cells should be derived from the same ancestral epithelial stem cell. Both sections were submitted for independent library preparation, genome sequencing, variant calling and filtering using our pipeline. The majority of substitution variant calls (2,742 of 2,933, 93.5%) were shared between both sections (Supplementary Fig. 1c). By contrast, when comparing five separate crypts from a mouse, a maximum of two variants were shared between two crypts, and no variants were shared by three or more crypts (Supplementary Fig. 1d).

Sample filtering

Our method for estimation of mutation rates assumes monoclonality of colorectal crypt samples. This assumption can be violated owing to several causes, including contamination from other colorectal crypts during microdissection or library preparation, contamination with non-epithelial cells located in or near the crypt, insufficient time for a stem cell to drift to clonality within the crypt, or the possibility that in some species, unlike in humans8, polyclonal crypts are the norm. Therefore, a truncated binomial mixture model was applied so as to remove crypts that showed evidence of polyclonality, or for which the possibility of polyclonality could not be excluded. An expectation–maximization (EM) algorithm was used to determine the optimal number of VAF clusters within each crypt sample, as well as each cluster’s location and relative contribution to the overall VAF distribution. The algorithm considered a range of numbers of clusters (1–5), with the optimal number being that which minimized the Bayesian information criterion (BIC). As the minimum number of supporting reads to call a variant was four, the binomial probability distribution was truncated to incorporate this minimum requirement for the number of successes, and subsequently re-normalized. The EM algorithm returned the inferred optimal number of clusters, the mean VAF (location) and mixing proportion (contribution) of each clone, and an assignment of each input variant to the most likely cluster. After applying this model to the somatic substitutions identified in each sample, sample filtering was performed on the basis of the following three criteria.

Low mutation burden

We discarded samples that presented fewer than 50 somatic variants, which was indicative of low DNA quality or sequencing issues.

High mutation burden

We discarded samples with a number of somatic variants greater than 3 times the median burden of samples from the same individual (excluding samples with fewer than 50 variants). This served to exclude a small minority of samples that presented evident sequencing quality problems (such as low sequencing coverage), but which did not fulfil the low-VAF criterion for exclusion (see below).

Low VAF

We discarded samples in which less than 70% of the somatic variants were assigned to clusters with VAF ≥ 0.3. However, this rule was not applied to those cases in which all the samples from the same individual had primary clusters with mean VAF < 0.3; this was done to prevent the removal of samples from individuals presenting high fractions of non-epithelial cells, but whose crypts were nonetheless dominated by a single clone.

These criteria led to the exclusion of 41 out of 249 samples. On the basis of visual assessment of sequencing coverage and VAF distributions, we decided to preserve three samples (ND0003c_lo0004, ND0003c_lo0011, TIGRD0001b_lo0010) that we considered to be clonal, but which would have been discarded on the basis of the criteria above.

Mitochondrial variant calling and filtering

For six species whose reference genome assemblies did not include the mitochondrial sequence, mitochondrial reference sequences were obtained from the GenBank database (Supplementary Table 5). For each species, alignment to the reference genome was performed using BWA (v.0.7.17-r1188), as described above (see ‘Sequence read alignment’). Pileup files were generated for mtDNA genomes using the ‘bam2R’ function in the deepSNV (v.1.32.0) R package62,63. The mapping quality cut-off was set to 0, taking advantage of the fact that the mitochondrial genome coverage for most samples was more than 100-fold higher than the nuclear genome coverage, and hence most reads with poor mapping scores should be of mitochondrial origin. Mitochondrial variants were called using the Shearwater algorithm62 (deepSNV package v.1.32.0). Multiple rounds of filtering were applied to identify and remove false positives. The first set of filters removed germline polymorphisms, applied a maximum false discovery rate (FDR) threshold of q > 0.01, required that mismatches should be supported by at least one read on both the forward and reverse strands, and merged consecutive indel calls. Further filtering steps were as follows.

Minimum VAF filter

Only variants with VAF > 0.01 were considered for analysis, based on the quality of the mutational spectra.

Sequencing coverage filter

Owing to species-specific mtDNA regions of poor mappability, we discarded sites with a read coverage of less than 500×.

D-loop filter

Analysis of the distribution of mutations along the mitochondrial genome revealed clusters of mutations within the hypervariable region of mtDNA known as the D-loop. To obtain estimates of the mutation burden in mtDNA unaffected by hypermutation of the D-loop, mutations in the D-loop region (coordinates MT:1–576 and MT:16,024–16,569 in human) were excluded from this analysis.

High mutation burden

We discarded samples that had a number of somatic mtDNA variants greater than four times the mean mtDNA burden across all samples. This served to exclude a small minority of samples that were suspected of enrichment in false positive calls. Visual inspection of these samples in a genome browser confirmed the presence of high numbers of variants found on sequence reads with identical start positions and/or multiple base mismatches, suggestive of library preparation or sequencing artefacts.

We examined the mutational spectra of somatic mtDNA substitutions on a trinucleotide sequence context (Extended Data Fig. 14b). The specificity of the filtered variant calls was supported by the observation that the mutational spectra across species were broadly consistent with those previously observed in studies of human tissues46, with a dominance of C>T and T>C transversions and a strong replication strand bias.

Mitochondrial copy number analysis

Sequence reads from each sample were separately mapped to the species-specific mtDNA reference sequence to estimate average mtDNA sequencing coverage. Excluding nuclear reference sequences from the alignment enabled even coverage to be obtained across the mitochondrial genome by preventing the mismapping of sequence reads to inherited nuclear insertions of mitochondrial DNA (known as NuMTs). Next, coverage information from individual mtDNA and whole-genome alignment (BAM) files was obtained using the genomecov tool in the bedtools suite (v.2.17.0)64. Mitochondrial copy number was calculated according to the formula

$${{rm{depth}}}_{{rm{mtDNA}}}times {rm{ploidy}}/{{rm{depth}}}_{{rm{gDNA}}},$$

where depthmtDNA and depthgDNA are the mean coverage values for mtDNA and the nuclear genome, respectively, and ploidy = 2 (assuming normal somatic cells to be diploid). For simplicity, the sex chromosomes were excluded from the calculation of the mean nuclear genome coverage.

Calculation of analysable genome size

To estimate the somatic mutation rate, it was first necessary to establish the size of the analysable nuclear genome (that is, the portion of the genome in which variant calling could be performed reliably) for each sample (Supplementary Table 4). For both substitutions and indels, the analysable genome of a sample was defined as the complement of the union of the following genomic regions: regions reported as ‘not analysed’ by the CaVEMan variant caller; regions failing the ‘chromosome and contig’ filter; regions failing the ‘N-tract and contig-end’ filter; and regions failing the ‘sequencing coverage’ filter (see ‘Variant filtering’). For the analysis of mitochondrial variants, the analysable genome of a sample was defined as the portion of mtDNA that satisfied the ‘sequencing coverage’ filter (see ‘Mitochondrial variant calling and filtering’), after subtracting the hypervariable region (D-loop).

Life-history data

Obtaining accurate lifespan estimates is challenging; although point estimates of maximum lifespan are available for many species, their veracity is often difficult to assess and estimates can vary widely for the same species (Supplementary Table 6). There can be many causes for this variation, including errors in recording and real variation in longevity between populations (that is, captive versus wild). As we were interested in whether the somatic mutation burden has an association with lifespan in the absence of extrinsic mortality, we sought to obtain estimates of longevity from individuals under human care, to minimize the effect of external factors such as predation or infection.

Mortality records for 14 species were obtained from the Species360 database, authorized by Species360 research data use agreement no. 60633 (Species360 Zoological Information Management System (ZIMS) (2020), https://zims.species360.org). This database contains lifespan data of zoo animals from international zoo records. Using records from 1980 to the present, we excluded animals for which the date of birth or death was unknown or uncertain. To avoid infant mortality influencing the longevity estimates for each species, we removed animals that died before the age of female sexual maturity, as defined by the AnAge database65. This resulted in a mean of 2,681 animal lifespan records per species for the species in the study (minimum 309, maximum 8,403; Supplementary Table 6). For the domestic dog, we combined records for domestic dogs (Canis lupus familiaris) and wolves (Canis lupus), because of the paucity of records for domestic dogs in Species360. Although the data are curated, they are still vulnerable to the presence of inaccurate records, which can bias the lifespan estimates. To reduce the effect of these outliers, for each species lifespan was estimated as the age at which 80% of the adults from that species had died66 (Supplementary Table 6).

Human longevity estimates were obtained using census birth and death record data from Denmark, (1900–2020), Finland (1900–2019) and France (1900–2018), retrieved from the Human Mortality Database (University of California, Berkeley (USA), and Max Planck Institute for Demographic Research (Germany); https://www.mortality.org, https://www.humanmortality.de). We selected these countries because they had census records going back at least 100 years. To remove the effect of infant mortality, we excluded individuals who died before the age of 13. For each country, we selected the cohort born in 1900 and calculated the age at which 80% of the individuals had died (Denmark, 87 years; Finland, 83 years; France, 81 years). We then used the mean of the three countries as our estimate of the human 80% lifespan (83.7 years) (Supplementary Table 6).

To test the effects of different estimates of lifespan on our results, we also obtained maximum longevity estimates for each species from a range of databases67 and a survey of the literature (Supplementary Table 6). Other life-history metrics were obtained from the AnAge database65 (Supplementary Table 6).

Mutational signature analysis

Mutational signatures of substitutions on a trinucleotide sequence context were inferred from sets of somatic mutation counts using the sigfit (v.2.1.0) R package31. Initially, signature extraction was performed de novo for a range of numbers of signatures (N = 2,…,10), using counts of mutations grouped per sample, per individual and per species. To account for differences in sequence composition across samples, and especially across species, mutational opportunities per sample, per individual and per species were calculated from the reference trinucleotide frequencies across the analysable genome of each sample (see ‘Calculation of analysable genome size’), and supplied to the ‘extract_signatures’ function in sigfit. The ‘convert_signatures’ function in sigfit was subsequently used to transform the extracted signatures to a human-relative representation (Fig. 2b), by scaling the mutation probability values using the corresponding human genome trinucleotide frequencies. The best-supported number of signatures, on the basis of overall goodness-of-fit31 and consistency with known COSMIC signatures (https://cancer.sanger.ac.uk/signatures/), was found to be N = 3. The cleanest deconvolution of the three signatures was achieved when using the mutation counts grouped by species, rather than by sample or individual. The three extracted signatures (labelled SBSA, SBSB and SBSC) were found to be highly similar to COSMIC signatures SBS1 (cosine similarity 0.96), SBS5 (0.89) and SBS18 (0.91), respectively. These signatures were independently validated using the MutationalPatterns (v.1.12.0) R package68, which produced comparable signatures (respective cosine similarities 0.999, 0.98 and 0.89).

This de novo signature extraction approach, however, failed to deconvolute signatures SBSA and SBSB entirely from each other, resulting in a general overestimation of the exposure to SBSA (Extended Data Fig. 15). To obtain more accurate estimates of signature exposure, the deconvolution was repeated using an alternative approach that combines signature fitting and extraction in a single inference process31. More specifically, the ‘fit_extract_signatures’ function in sigfit was used to fit COSMIC signature SBS1 (retrieved from the COSMIC v,3.0 signature catalogue; https://cancer.sanger.ac.uk/signatures/) to the mutation counts grouped by species (with species-specific mutational opportunities), while simultaneously extracting two additional signatures de novo (SBSB and SBSC). Before this operation, COSMIC SBS1 was transformed from its human-relative representation to a genome-independent representation using the ‘convert_signatures’ function in sigfit. By completely deconvoluting SBS1 and SBSB, this approach yielded a version of SBSB that was more similar to COSMIC SBS5 (cosine similarity 0.93); the similarity of SBSC to COSMIC SBS18 was the same under both approaches (0.91).

Finally, the inferred signatures were re-fitted to the mutational spectra of mutations in each sample (using the ‘fit_signatures’ function in sigfit with sample-specific mutational opportunities) to estimate the exposure of each sample to each signature. The fitting of the three signatures yielded spectrum reconstruction similarity values (measured as the cosine similarity between the observed mutational spectrum and a spectrum reconstructed from the inferred signatures and exposures) with median 0.98 and interquartile range 0.96–0.99. Although the purely de novo extraction approach and the ‘fitting and extraction’ approach yielded comparable versions of signatures SBSB and SBSC, the fixing of COSMIC SBS1 in the latter approach resulted in lower SBS1 exposures and higher SBSB exposures in most samples, owing to the cleaner deconvolution of these two signatures (Fig. 2, Extended Data Fig. 15).

To examine potential variation in the spectrum of signature SBS5 across species, the following procedure was conducted for each species: individual-specific mutation counts and mutational opportunities were calculated for each individual in the species, and the ‘fit_extract_signatures’ function was used to fit COSMIC signatures SBS1, SBS18 and SBS34 (transformed to a genome-independent representation using the ‘convert_signatures’ function) to the mutational spectra of each individual, while simultaneously inferring one additional signature (corresponding to signature SBS5 as manifested in that species; Extended Data Fig. 6).

To assess the presence in non-human colorectal crypts of mutational signatures caused by APOBEC or colibactin, which have been previously observed in human crypts8, we used an expectation–maximization algorithm for signature fitting, in combination with likelihood ratio tests (LRTs). More specifically, for each non-human sample, we tested for exposure to colibactin (signature SBS88, COSMIC v.3.2) by comparing the log-likelihoods of (i) a model fitting COSMIC signatures SBS1, SBS5, SBS18, SBS34 and SBS88, and (ii) a reduced model fitting only the first four signatures. Benjamini–Hochberg multiple-testing correction was applied to the P values that resulted from the LRTs, and colibactin exposure was considered significant in a sample if the corresponding corrected q-value was less than 0.05. We followed the same approach to assess exposure to APOBEC (SBS2 and SBS13), using two separate sets of LRTs for models including either SBS2 or SBS13, in addition to SBS1, SBS5, SBS18 and SBS34. APOBEC exposure was considered significant in a sample if its q-values for the models including SBS2 and SBS13 were both less than 0.05. This analysis identified 1/180 crypts with significant exposure to each of colibactin and APOBEC (although the evidence for the presence of the relevant signatures in these two crypts was not conclusive). To test for depletion of colibactin or APOBEC exposure in non-human crypts relative to human crypts, we first applied the LRT-based method described above to a published set of 445 human colorectal crypts8, identifying 92 colibactin-positive and 9 APOBEC-positive crypts. We then compared the numbers of colibactin- and APOBEC-positive crypts in the human and non-human sets using two separate Fisher’s exact tests (‘fisher.test’ function in R). This revealed the difference in colibactin exposure to be highly significant (P = 7 × 10–14), unlike the difference in APOBEC exposure (P = 0.30).

Mutational spectra of somatic indels identified in each species were generated using the ‘indel.spectrum’ function in the Indelwald tool for R (24/09/2021 version; https://github.com/MaximilianStammnitz/Indelwald).

Selection analysis

Evidence of selection was assessed using the ratio of nonsynonymous to synonymous substitution rates (dN/dS) in the somatic mutations called in each species. The dNdScv (v.0.0.1.0) R package38 was used to estimate dN/dS ratios for missense and truncating substitutions in each species separately. Reference CDS databases for the dNdScv package were built for those species with available genome annotation in Ensembl (https://www.ensembl.org; Supplementary Table 2), using the ‘buildref’ function in dNdScv. For each species, the ‘dndscv’ function was applied to the list of somatic substitutions called in samples of that species, after de-duplicating any substitutions that were shared between samples from the same individual to avoid counting shared somatic mutations multiple times. In addition, the analysis was restricted to genes that were fully contained in the analysable genomes of all samples from the species (a condition satisfied by the vast majority of protein-coding genes). Genome-wide and gene-specific dN/dS ratios were obtained for missense and truncating substitutions in each species; no genes with statistically significant dN/dS ≠ 1 were observed.

Copy number analysis

For species with chromosome-level assemblies (cat, cow, dog, horse, human, mouse, rabbit and rat), the total and the allele-specific copy number (CN) was assessed in each sample, adapting a likelihood model that was previously applied to the detection of subclonal CN changes in healthy human skin19. This method exploits two sources of evidence: relative sequencing coverage and B-allele fraction (BAF; the fraction of reads covering a heterozygous single-nucleotide polymorphism (SNP) that support one of the alleles). Human samples PD36813x15 and PD36813x16 were excluded from this analysis owing to the poor quality of their SNP data.

For each sample, sequencing coverage was measured in non-overlapping 100-kb bins along the reference genome of the species, using the coverageBed tool in the bedtools suite (v.2.17.0)64. For each bin, the coverage per base pair was calculated by dividing the number of reads mapping to the bin by the bin length, and multiplying the result by the read length (150 bp). A normalized sample–normal coverage ratio was then calculated for each bin by dividing the bin coverage in the sample by the corresponding coverage in its matched normal control (see ‘Sample processing’). Heterozygous SNPs were isolated for each sample by selecting germline SNPs with a BAF between 0.4 and 0.6 in the matched normal sample, and a coverage of at least 15 reads in both the target sample and its matched normal sample. After assigning each SNP to its corresponding 100-kb genome bin, the bins in each sample were divided into two sets: (i) bins with coverage ≥ 10 in both the target sample and its matched normal, and at least one heterozygous SNP; and (ii) bins with coverage ≥ 10 in both the target sample and its matched normal, and no heterozygous SNPs. For the first set, estimates of total and allele-specific CN were inferred by maximizing the joint likelihood of a beta-binomial model for BAF and a negative binomial model for relative coverage, as previously described19. The most likely combination of allele CN values was obtained for each bin by conducting an exhaustive search of CN values between 0 and 4, and selecting the combination maximizing the joint likelihood (calculated on the basis of expected BAF and relative coverage values). A penalty matrix was used to penalize more complex solutions over simpler ones, as previously described19. For the second set of bins (bins without SNPs), only estimates of total CN were inferred, by maximizing the likelihood of a negative binomial model for relative coverage. The most substantial differences between these methods and the one previously published are: (i) SNPs were obtained from the variant calling output, instead of from a public database; (ii) relative coverage was calculated per 100-kb bin, rather than per SNP; (iii) SNPs were not phased within each gene, but within each bin; (iv) no reference bias was assumed (that is, the underlying BAF of heterozygous SNPs was assumed to be 0.5); (v) the minimum sample purity was raised to 0.85; (vi) putative CN changes were not subjected to significance testing, but selected according to their likelihood, and subsequently filtered by means of a segmentation algorithm (see below).

Estimates of total and allele-specific CN per bin were merged into CN segments, which were defined as contiguous segments composed of five or more bins with identical CN states. Segmentation was performed separately for total and allele-specific CN estimates in each sample. After this process, any pair of adjacent segments with the same CN assignment, and separated by a distance shorter than five bins, was merged into a single segment. Finally, within each species, segments presenting CN values other than 2 (or 1/1 for allele-specific CN), and being either shorter than 10 bins (1 Mb), or shared among two or more samples, were discarded, resulting in the removal of nearly all spurious CN changes.

Estimation of mutation rate

For each sample, the somatic mutation density (mutations per bp) was calculated by dividing the somatic mutation burden (total number of mutations called) by the analysable genome size for the sample (see ‘Calculation of analysable genome size’). The adjusted somatic mutation burden (number of mutations per whole genome) was then calculated by multiplying the mutation density by the total genome size of the species (see below). The somatic mutation rate per year (mutations per genome per year) was obtained by dividing this adjusted mutation burden by the age of the individual, expressed in years (Supplementary Table 2). The expected ELB for each sample was calculated by multiplying the somatic mutation rate by the estimated lifespan of the species (see ‘Life-history data’).

The total genome size of a species was estimated as the total size of its reference genome assembly. Across species, the mean genome size was 2.67 Gb, ranging between 2.41 Gb and 3.15 Gb and with a standard deviation of 221 Mb (Supplementary Table 4). This suggests that inter-species variation in genome size should not have a substantial influence on the somatic mutation rate estimates. For an assessment of alternative measures of mutation rate, see ‘Association of mutation rate and end-of-lifespan burden with lifespan’.

Association of mutation rate with life-history traits

The association of the somatic mutation rate with different life-history traits was assessed using LME models. In particular, associations with the following traits were examined: lifespan (in years), adult mass (or adult weight, in grams), BMR (in watts), and litter size (see ‘Life-history data’). Associations for lifespan, adult mass and BMR were assessed using the following transformed variables: 1/lifespan, log10(adult mass) and log10(BMR). To account for the potentially confounding effect of the correlation between metabolic rate and body mass, the residuals of a fitted allometric regression model of BMR on adult mass (equivalent to a simple linear regression of log10(BMR) on log10(adult mass)) were used as a mass-adjusted measure of metabolic rate, referred to as ‘BMR residuals’.

For each variable, an LME model was implemented for the regression of somatic mutation rates per sample on the variable of interest, using the ‘lme’ function in the nlme R package (v.3.1-137; https://cran.r-project.org/web/packages/nlme). To account for non-independence of the samples, both at the individual level and at the species level, the model included fixed effects (intercept and slope parameters) for the variable of interest, and random effects (slope parameters) at the individual and species levels. In addition, to account for the heteroscedasticity of mutation rate estimates across species, the usual assumption of constant response variance was replaced with explicit species-specific variances, to be estimated within the model.

To determine the fraction of inter-species variance in mutation rate explained by each life-history variable individually, the LME model described above was used to produce predictions of the mean mutation rate per species; only fixed effects were used when obtaining these predictions, random effects being ignored. The variance of these predictions was then compared to the variance in observed mean mutation rates; the latter were calculated for each species as the mean of the observed mean rates per individual, to avoid individuals with larger numbers of samples exerting a stronger influence on the species mean. The fraction of inter-species variance explained by the model was calculated using the standard formula for the coefficient of determination,

$${R}^{2}={rm{ESS}}/({rm{ESS}}+{rm{RSS}}),$$

where ESS is the explained sum of squares, and RSS is the residual sum of squares:

$${rm{ESS}}={sum }_{i}{({hat{y}}_{i}-overline{y})}^{2},{rm{RSS}}={sum }_{i}{({y}_{i}-{hat{y}}_{i})}^{2}.$$

In this formulation, ({y}_{i}) and ({hat{y}}_{i}) denote the observed and predicted mutation rates for species i, respectively, and (bar{y}) is the overall mean rate. This definition of R2 coincides with the fraction of variance explained (FVE), defined as 1 minus the fraction of variance unexplained (FVU):

$${rm{FVE}}=1{rm{mbox{–}}}{rm{FVU}}=1{rm{mbox{–}}}[{rm{RSS}}/({rm{ESS}}+{rm{RSS}})]={rm{ESS}}/({rm{ESS}}+{rm{RSS}})={R}^{2}.$$

As the predicted and observed values correspond to mean mutation rates per species, rather than mutation rates per sample, FVE provides a measure of the fraction of inter-species variance explained by the fixed effects of the LME model. Among the variables considered, 1/lifespan was found to have the greatest explanatory power (FVE = 0.84, using a free-intercept model).

To compare the explanatory power of variables other than 1/lifespan when considered either individually or in combination with 1/lifespan, the method described above was also applied to two-variable combinations of 1/lifespan and each of the remaining variables, using an LME model with fixed effects for both variables and random effects for 1/lifespan only. The R2 formula above was used to measure the fraction of inter-species variance explained by each model. In addition, to test whether the inclusion of a second explanatory variable was justified by the increase in model fit, LRTs between each two-variable LME model and a reduced LME model including only 1/lifespan were performed using the ‘anova’ function in the nlme R package.

To further assess the potential effects of body mass and lifespan on each other’s association with the somatic mutation rate, allometric regression models (equivalent to simple linear models under logarithmic transformation of both variables) were fitted to the mean somatic mutation rate per species, using either adult mass or lifespan as the explanatory variable. In addition, the ‘allometric residuals’ of mutation rate, adult mass and lifespan (that is, the residuals of pairwise allometric regressions among these three variables) were used to examine the associations between somatic mutation rate and either body mass or lifespan, after accounting for the effect of the third variable (partial correlation analysis). For example, to account for the potential influence of body mass on the relationship between somatic mutation rate and lifespan, the residuals of an allometric regression between mutation rate and adult mass, and the residuals of an allometric regression between lifespan and adult mass, were analysed using simple linear regression. This analysis supported a strong association between somatic mutation rate and lifespan (independently of the effect of mass; FVE = 0.82, P = 3.2 × 10–6; Fig. 3c), and a non-significant association between somatic mutation rate and body mass (independently of the effect of lifespan). Therefore, the relationship between somatic mutation rate and lifespan does not appear to be mediated by the effect of body mass on both variables. Of note, this result remains after excluding naked mole-rat: after removing this species, partial correlation analysis still reveals a strong association between somatic mutation rate and lifespan (FVE = 0.77, P = 4.1 × 10–5), and a non-significant association between somatic mutation rate and body mass (P = 0.84). This demonstrates that the observed relationships are not dependent on the presence of naked mole-rat in the study.

To assess the robustness of the LME regression analyses described above, we performed bootstrap analysis on each LME model, at the level of both individuals and species. More specifically, for each level we used each of the LME models to perform regression on 10,000 bootstrap replicates, produced by resampling either species or individuals with replacement. We then assessed the distributions of FVE across bootstrap replicates (Extended Data Fig. 13c). In addition, we performed a similar bootstrap analysis using a collection of maximum longevity estimates obtained from the literature (see ‘Life-history data’). We applied the zero-intercept LME model described above (for regressing mutation rate on inverse lifespan) on a set of 5,000 bootstrap replicates, each of which used a set of species lifespan estimates randomly sampled from the collection of literature-derived estimates (Extended Data Fig. 12).

The results obtained with the LME models were additionally validated using an independent hierarchical Bayesian model, in which the mean somatic mutation burden of each individual was modelled as following a normal distribution with mean defined as a linear predictor containing a species-specific slope parameter and a multiplicative offset (corresponding to the individual’s age; inclusion of this offset minimizes the heteroscedasticity of the mutation rate across species, which results from dividing mutation burdens by age). Species-specific slope parameters were in turn modelled as normally distributed around a global slope parameter, equivalent to the fixed-effect slope estimated by the LME model. This hierarchical model produced very similar results to those of the LME model for all life-history variables (Extended Data Fig. 13a).

We note that samples CATD0002b_lo0003 and MD6267ab_lo0003 were excluded from all regression analyses, owing to the fact that each shared the most of its somatic variants with another sample from the same individual (indicating, in each case, that both samples were closely related), hence violating the assumption of independence among samples. The inclusion of these two samples, however, had no effect on the outcome of the analyses.

Association of mutation rate and end-of-lifespan burden with lifespan

The relationship between somatic mutation rate and species lifespan was further explored by adapting the LME model described in the previous section to perform constrained (zero-intercept) regression of the adjusted mutation rate per year on the inverse of lifespan, 1/lifespan (see ‘Life-history data’, ‘Estimation of mutation rate’ and ‘Association of mutation rate with life-history traits’). The use of zero-intercept regression was motivated by the prediction that, if somatic mutation is a determinant of maximum lifespan, then it would be expected for all species to end their lifespans with a similar somatic mutation burden. Indeed, this was confirmed by simple linear regression of the species mean end-of-lifespan mutation burden against species lifespan (slope P = 0.39). Thus, if m is the mutation rate per year, and L is the species’ lifespan, the expected relationship is of the form.

where k is a constant representing the typical end-of-lifespan mutation burden across species. According to this relationship, the mutation rate per year is linearly related to the inverse of lifespan,

Therefore, the cross-species average end-of-lifespan burden (k), can be estimated as the slope parameter of a zero-intercept linear regression model with the mutation rate per year (m) as the dependent variable, and the inverse of lifespan (1/L) as the explanatory variable. To this purpose, the LME model described in the previous section was altered by removing the fixed-effect intercept parameter, thus considering only fixed- and random-effect slope parameters for 1/Lifespan.

The zero-intercept LME model estimated a value of k = 3,210.52 (95% confidence interval 2,686.89–3,734.15). The fraction of inter-species variance explained by the zero-intercept model (FVE) was 0.82, whereas the LME model described in the previous section (which estimated k = 2,869.98, and an intercept of 14.76) achieved FVE = 0.84 (see ‘Association of mutation rate with life-history traits’). To test whether the increase in model fit justifies the inclusion of an intercept, both models were compared using a LRT (as implemented by the ‘anova’ function in the nlme R package (v.3.1-137)). This yielded P = 0.23, indicating that the free-intercept model does not achieve a significantly better fit than the zero-intercept model. Similarly, the zero-intercept model yielded lower values for both the Bayesian information criterion (BIC) and the Akaike information criterion (AIC). Notably, equivalent analyses using somatic mutation rates per megabase and per protein-coding exome (instead of per whole genome) yielded comparable results (Extended Data Fig. 11).

To investigate the possibility of phylogenetic relationships between the species sampled confounding the analysis, a phylogenetic generalized linear model was used to regress the mean mutation rate of each species on the inverse of its lifespan (1/L), while accounting for the phylogenetic relationships among species. A phylogenetic tree of the 15 species examined was obtained from the TimeTree resource69, and the phylogenetic linear model was fitted using the ‘pgls’ function in the caper R package (v.1.0.1; https://cran.r-project.org/web/packages/caper). The estimates produced by zero-intercept regression of mean mutation rates per species on 1/lifespan were compared between this phylogenetic generalized linear model and a simple linear model (‘lm’ function in R). The use of this simple model, as well as the use of mean mutation rates per species (rather than mutation rates per sample), was necessary owing to the impossibility of replicating the heteroscedastic mixed-effects structure of the LME model used for the main association analyses (see ‘Association of mutation rate with life-history traits’) within the phylogenetic linear model. Both the phylogenetic linear model and the simple linear model produced similar estimates (Extended Data Fig. 13b), suggesting that phylogenetic non-independence of the samples does not have a substantial effect on the association analyses.

Cell division analysis

To investigate the extent to which differences in cell division rate could explain differences in mutation rate and burden across species, we obtained estimates of intestinal crypt cell division rates from mouse70, rat71 and human72,73 (Supplementary Table 7). Using these cell division rates, our lifespan estimates and the observed substitution rates, we calculated the number of cell divisions at the end of lifespan and the corresponding number of mutations per cell division expected under a simple model assuming that all mutations occur during cell division (Supplementary Table 7).

We investigated whether differences in the number of cell divisions among species could explain the observed differences in mutation burden. Although colorectal cell division rate estimates are lacking for most species, existing estimates from mouse, rat and human indicate that the total number of stem cell divisions per crypt in a lifetime varies greatly across species—for example, there are around 6- to 31-fold more divisions per intestinal stem cell in a human than in a rat over their respective lifetimes, depending on the estimate of cell division rate used (Supplementary Table 7). Mouse intestinal stem cells are estimated to divide once every 24 h (ref. 70), whereas estimates of the human intestinal stem cell division rate vary from once every 48 h (ref. 72) to once every 264 h (ref. 73). Thus, mouse intestinal stem cells divide 2–11 times faster than human intestinal stem cells. By the end of lifespan, an intestinal stem cell is predicted to have divided around 1,351 times in a mouse, around 486 times in a rat and 2,774–15,257 times in a human (depending on the estimate of cell division rate used). Applying our somatic mutation burden and lifespan data, this implies that the somatic mutation rate per cell division in a mouse is around 1.5- to 8.4-fold higher than in a human. However, the observed fold difference in somatic mutation rate between these two species is 16.9 (Table 1). Therefore, differences in cell division rate appear unable to fully account for the observed differences in mutation rate across species. Nevertheless, we note that accurate cell division rate estimates for basal intestinal stem cells are lacking for most species.

Reporting summary

Further information on research design is available in the Nature Research Reporting Summary linked to this paper.

Leave a Comment