Mutation rate evolution in cetacean genomes: understanding drivers of variability and observation bias in estimates of genomic mutation rate
Amy M. Van Cise1, Sophie Garrote1, Phillip A. Morin2*
- School of Aquatic and Fishery Sciences, University of Washington, Seattle, WA
- Southwest Fisheries Science Center, National Marine Fisheries Service, National Oceanic and Atmospheric Administration, La Jolla, CA
*Corresponding author email: phillip.morin@noaa.gov
Running page head: Genomic mutation rate in cetaceans
Abstract
Genomic mutation is the fundamental process driving evolutionary processes that result in contemporary species, subspecies, and locally adapted populations. Thus, genomic studies of divergence, speciation, and evolutionary demographics rely on accurate estimates of mutation rate. Historically, mutation rates have been estimated using a small portion of the genome and/or averaged among multiple species within a family or order, but growing evidence indicates potentially high variability in mutation rate both across the genome and among species. Here, we examine variability in mutation rate within cetaceans by mapping reference-quality genomes using a range of pipeline parameters in order to determine which factors have the greatest effect on variability in substitution rate or its estimation. Our results indicate that mutation rate has varied over evolutionary time and across the genome in cetaceans. Notably, we show that generational mutation rate increases with maximum lifespan despite decreases in annual mutation rate with maximum lifespan. Mutation rate also decreases with chromosome size and is higher in coding regions than non-coding regions. Similarly, we found significant effects of pipeline parameter choice on mutation rate, which is most affected by the choice of reference genome used in the alignment pipeline. These results indicate the importance of thoughtful selection of mutation rate used in future cetacean genomic studies, in order to ensure accurate downstream analyses of evolutionary processes.
Introduction
As the fundamental process underlying evolution, mutation drives processes such as evolutionary adaptation and divergence among species. The rate at which mutations accrue will ultimately determine the rate at which these processes occur; conversely, reconstruction of evolutionary processes as well as genomic estimates of historical population demography have the potential to be biased by inaccurate estimates of mutation ratee.g. 1. As reliance on genomics increases for the conservation and management of threatened populations2, it is increasingly important that we comprehensively understand how mutation rate varies within taxa and with observation bias.
Broad scale trends in mutation rate variability across taxa have been correlated with phylogenetic inertia as well as various demographic factorse.g. 3,4–7. Three hypotheses have been proposed as potential drivers of mutation rate variability among species. (1) The longevity hypothesis posits that mutation rate decreases with increased lifespan, allowing individuals to live for longer periods of time without adverse effects of the accumulation of harmful mutations over their lifetimes8. (2) The generation time hypothesis posits that mutation rate decreases with respect to generation time likely due to mutation accumulation in gamete DNA9,10. (3) The metabolic rate hypothesis posits that mutation rate increases with increasing metabolic rate (often proxied by body length) due to the oxidative effects of metabolism on the cell11.
Mutation rate has been shown to vary over evolutionary timescales4. At a finer scale, mutation rate varies across the genome (e.g. among chromosomes and coding vs. non-coding regions), and site-specific variability is only partially context-dependent12,13. In addition to these factors, there are multiple decision points along the bioinformatic pipeline from sequence data to aligned genome that may affect the accurate tally of genomic mutations and downstream estimation of mutation rate. However, there is limited research on how the selection of evolutionary timescale, genomic region, or alignment pipeline parameters affect estimates of mutation rate.
Mutation rate estimates form the basis of many conservation genomic and evolutionary genomic studies; for example, as a parameter included in estimates of divergence among species, historical population size, or patterns of selective divergence as they relate to neutral variation. Biases in this estimate may be reflected in many downstream analyses that shape our understanding of evolutionary biology; they may also affect conservation priorities and management actions. In one example, a recent analysis of baleen whale mutation rates based on wild pedigrees resulted in higher estimates for these species than when estimated using phylogenetic divergence; using the new mitochondrial mutation rate reduces estimates of pre-exploitation North Atlantic humpback whale abundance by 86%14 and considerably shifts our understanding of the recovery of that population. Possibly owing in part to the diversity of factors affecting mutation rate, some studies suggest a lack of direct relationship between mtDNA genetic diversity and population size3.
The increasing accessibility of genomic sequencing technologies and high-quality reference genomes for marine mammals16,17 has precipitated a rise in the use of genomic techniques to address previously unanswered questions about the ecology and conservation of these speciese.g. 18,19–21. Because these species are rare, remote, elusive, and protected, it is often difficult to obtain the sample size needed to accurately detect and describe ecological processes and patterns. In cases such as this, whole genome analyses increase sensitivity and power to describe genetic patterns using a smaller sample set, thus making it a potentially powerful tool for studies of evolutionary ecology in marine mammals. These genomic studies in marine mammals, e.g. descriptions of evolutionary phylogenies or historical demographics of species, currently rely on estimates of mutation rate that are averaged across a family or even suborder of species and often based on only a small portion of the genomee.g. 22,23.
Here, we explore variability in phylogeny-based estimation of genomic mutation rate within cetaceans, and describe the major drivers affecting that variability with the aim of providing insight and context for the estimation of an accurate and appropriate mutation rate across cetacean genomics studies. Note that some harmful mutations will be filtered out of a population over time and therefore not accounted for in phylogeny-based estimates of mutation rate. We therefore estimate the substitution rate; i.e. the rate at which new mutations accumulate over time, which is considered equal to mutation rate in cases of neutral evolution and may be lower than mutation rate in cases of positive or stabilizing selection. In the context of phylogeny-based rates, we use mutation rate and substitution rate interchangeably.
To understand how evolutionary processes and observation bias affect contemporary estimates of substitution rate within cetaceans, we examine several factors related to each. First, we examine process driven variability by estimating taxon- and chromosome-specific substitution rates, and comparing species-specific substitution rates with key demographic traits considered to affect substitution rate (i.e., longevity, generation time, and body length). We also examine potential drivers of observation error by varying key parameters in the pipeline used to align genomes and estimate a substitution rate. These input parameters include the choice of reference genome used for alignment, maximum read depth cutoff for genotyping a locus, inclusion of coding or non-coding regions in the mutation rate estimate, divergence date estimate used in the mutation rate calculation, and whether or not ancestral heterozygosity is included in the mutation rate calculation.
Methods
Sample and reference genome preparation
We use publicly available whole genome sequence (WGS) data from 14 cetacean species representing the breadth of cetacean phylogenetic diversity. Shotgun WGS or Arima Hi-C (Arima Genomics, Inc) chromatin-linked data were obtained from NCBI Genbank using SRAToolkit (v. 3.0.7, https://github.com/ncbi/sra-tools), or from the European Nucleotide Archive (Table 1). Sequence files were quality filtered using BBduk in BBtools (https://jgi.doe.gov/data-and-tools/software-tools/bbtools/REF)prior to read mapping.
In order to test for the effect of reference genome on the resulting substitution rate, we selected four reference genomes to represent three divergence times: a deeply diverged taxon, a recently diverged taxon, and an intermediate taxon (Table 2). The deeply diverged taxon is represented by the hippopotamus (Hippopotamus amphibius, divergence from all cetacea = approx. 53 MYA), and the intermediate taxon is represented by the pygmy sperm whale (Kogia breviceps, divergence from other cetacea approx. 34-37 MYA). We selected two species to represent the recently diverged taxon: the killer whale (Orcinus orca, divergence from other odontocetes approx. 11-20 MYA) is used for all odontocetes, and the North Atlantic right whale (Eubalaena glacialis, divergence from other mysticetes approx. 23-27 MYA) is used for all mysticetes. Exact divergence dates for each species-reference combination are listed in Supplemental Table S1.
We downloaded all reference genomes as fastq files from NCBI Genbank, and created index files using bwa-mem224 and samtools25. We masked repeats with RepeatMasker (version 4.1.4) and bedtools26. All steps described here, making up the pre-alignment pipeline, are outlined at: https://github.com/UW-WADE-lab/Whole-Genome-Alignment
| Species | Abbreviation | Data.source | Accession | HiC.SRA | Recent reference | Average Generation Time | Maximum Length (m) | Maximum Lifespan |
|---|---|---|---|---|---|---|---|---|
| minke whale | Bacu | ENA | ERR11040176 | HiC | Egla | 22.1 | 8.65 | 50 |
| common dolphin | Ddel | ENA | ERR11040185 | HiC | Oorc | 14.8 | 1.97 | 30 |
| Risso’s dolphin | Ggri | ENA | ERR13132929 | HiC | Oorc | 19.6 | 3.2 | 42.5 |
| long-finned pilot whale | Gmel | ENA | ERR11837528 | HiC | Oorc | 24 | 3.65 | 41.08 |
| northern bottlenose whale | Hamp | ENA | ERR10908646 | HiC | Oorc | 17.8 | 9.25 | 37 |
| Pacific white-sided dolphin | Lalb | ENA | ERR11042964 | HiC | Oorc | 18.1 | 1.74 | 30 |
| false killer whale | Pcra | NA | NA | SRA | Oorc | 25 | 3.4 | 62.5 |
| striped dolphin | Scoe | ENA | ERR11042965 | HiC | Oorc | 22.5 | 1.8 | 57.5 |
| blue whale | Bmus | NCBI | SRR5665644 | SRA | Egla | 30.8 | 25.75 | 131.7 |
| Rice’s whale | Bric | NCBI | SRR10331559 | SRA | Egla | 18.4 | 12 | 71.56 |
| North Atlantic right whale | Egla | NCBI | SRR11097130 | SRA | NA | 35.7 | 15.5 | 73 |
| humpback whale | Mnov | NCBI | SRR17854487 | SRA | Egla | 21.5 | 15.5 | 93 |
| gray whale | Erob | NCBI | SRR12437599 | SRA | Egla | 22.9 | 12.5 | 77 |
| killer whale | Oorc | NCBI | SRR574970, SRR574978, SRR574980, SRR574981 | SRA | NA | 25.7 | 5.55 | 90 |
| Amazon River dolphin | Igeo | NCBI | SRR11430609 | SRA | Oorc | 10.2 | 1.998 | 37.93 |
| pygmy sperm whale | Kbre | NCBI | SRR11430611 | SRA | Oorc | 12.1 | 3.125 | 32.43 |
| Blainville’s beaked whale | Mden | NCBI | SRR13167975 | SRA | Oorc | 23.4 | 4.11 | 60 |
| melon-headed whale | Pele | NCBI | SRR11097149 | SRA | Oorc | 19.6 | 2.75 | 47 |
| vaquita | Psin | NCBI | SRR15435906 | SRA | Oorc | 11.4 | 1.37 | 21 |
| Species | Abbreviation | Database | Accession | Contig.N50 | Scaffold.N50 | Coverage |
|---|---|---|---|---|---|---|
| pygmy sperm whale | Kbre | NCBI | GCA_026419965.1 | 42.8 Mb | 122.4 Mb | 23.0x |
| common hippopotamus | Hamb | NCBI | GCA_030028045.1 | 50.3 Mb | 149.2 Mb | 33.0x |
| killer whale | Oorc | NCBI | GCA_937001465.1 | 45.6 Mb | 114.2 Mb | 34.0x |
| North Atlantic right whale | Egla | NCBI | GCA_028564815.2 | 37.5 Mb | 129.6 Mb | 32.0x |
Sample alignment
Genomic datasets from 19 cetacean species were first aligned to the the pygmy sperm whale, which we chose to represent an “intermediate” divergence time from both mysticetes and odontocetes. We then repeated sequence alignment for each species using the hippopotamus as a deeply diverged reference, and either the killer whale (odontocetes) or the North Atlantic right whale (mysticetes) as a recently diverged reference (Table 1). This resulted in a total of three sequence alignments for each of the 19 species included in the study. The pygmy sperm whale, North Atlantic right whale, and killer whale were each excluded from self-alignment, so that each of these three species have two alignments rather than three.
The alignment pipeline, regardless of sequencing approach, comprised six steps: 1) align sample sequence data to reference, 2) sort aligned genome, 3) remove duplicate sequences, 4) mask repeat regions, 5) index genome, 6) call genotypes and filter variants by read depth. All alignment scripts can be found at https://github.com/UW-WADE-lab/Whole-Genome-Alignment.
All sequence data were aligned using the Burrows-Wheeler Aligner implemented by bwa-mem224 in Step 1. Forward and reverse Illumina short-read WGS sequence data were aligned simultaneously. Arima Hi-C genomes, which comprise longer reads that contain chimeric reads from chromatin cross-linking prior to DNA extraction27, require an additional filtering step post-alignment. Therefore, we aligned Arima Hi-C forward and reverse reads separately, then filtered forward and reverse read files using custom perl scripts generated by Arima to identify and remove chimeras, before combining forward and reverse read alignments (https://github.com/ArimaGenomics/mapping_pipeline).
Steps 2-6 were identical for both types of sample genome. We sorted genome alignments using samtools25 and removed duplicates using Picard (http://broadinstitute.github.io/picard). We removed previously identified repeat regions with bedtools26. We indexed the sample genome using samtools, called genotypes using bcftools28, and filtered variants to remove those with a read depth less than 10 or greater than 500, also using bcftools. Note that this read depth filter is less conservative than is normally used in an alignment pipeline29,30; this allowed us to directly examine the effect of read depth on mutation rate later in the study.
Data analysis: calculating substitution rate
Substitution rate (\(\mu\)) can be estimated as genomic divergence (i.e. genomic heterozygosity, \(d_{xy}\)) over divergence time (\(t\)). Following Robinson et al. (2022), we use the following equation:
\[\begin{equation} \mu=d_{xy}/2t \end{equation}\]
Genomic divergence (\(d_{xy}\)), in turn, is estimated as:
\[\begin{equation} d_{xy}=(h_1+2h_2)/2g \end{equation}\]
where \(h_1\) is the number of sites in which the sample genome is heterozygous, \(h_2\) is the number of sites in which the sample genome is homozygous for the alternate allele, and \(g\) is the number of sites in which both the sample and reference genomes had called genotypes.
This calculation may result in an overestimate of substitution rate because ancestral (pre-divergence) heterozygosity is not taken into account, which we can account for in the following:
\[\begin{equation} \mu=(d_{xy}-\pi_{anc})/2t \end{equation}\]
The above equations produce a per-year substitution rate estimate. Substitution rates are often reported per generation in the literature, so to consider the effect of including generation time in substitution rate estimates we converted our per-year substitution rates to an estimate of substitution rate per generation. For most species, generation time was available from Taylor et al.31. All generation times are reported in Table 1; the published source of generation time for all species is listed in Supplemental Table S1.
Data analysis
Computation of substitution rate, and all downstream data analyses, were completed in R32 using Rstudio. We subset our substitution rate estimates according to each test described below. Unless otherwise noted below, for consistency among analyses we ensured that all datasets included genomes from all 19 species considered in the study, using estimates only from reads aligned to the pygmy sperm whale reference genome, divergence date estimates from McGowen et al.33 (6-partition AR model), and without correcting for ancestral heterozygosity.
Data analysis: phylogenetic or demographic effects on substitution rate
To examine the effect of phylogenetic inertia on substitution rate, we test for significant differentiation among taxonomic families in whole-genome substitution rate per year (substitution/site/year) as well as whole-genome substitution rate per generation substitution/site/generation) using an ANOVA. We compared observed patterns using both annual and generational estimates of substitution rate because they are both commonly reported in marine mammal genomic literature14,e.g. 22.
We also consider potential effects of lifespan and body size on annual whole-genome substitution rate using data on species-specific maximum lifespan and maximum body length. Most of these data were acquired from Groot et al.34, and missing data were filled in from published literature (Table 1, Supplemental Table S1). We first tested for correlation between maximum lifespan and maximum body length, and found them to be correlated in our dataset (Spearman’s \(\rho\) = 0.83). Next, we fit a generalized linear regression using both maximum lifespan and maximum body length as predictor variables. Because phylogenetic inertia may be a significant factor affecting variability in mutation rate among species, we used a phylogenetic generalized least squares (PGLS) modeling approach implemented using the phytools package in R35, and chose the best model and significant parameters based on reported AIC and p values, respectively.
Data analysis: pipeline effects on substitution rate
We first consider the effect of the following sources of potential variability in the genome alignment pipeline on downstream estimates of substitution rate:
- Choice of reference species for genomic alignment
- Chromosome
- Inclusion of coding vs. non-coding regions in estimate
- Maximum read depth cutoff for called genotypes
We also consider the effect of variability in two estimated parameters that are used in the calculation of substitution rate:
- Choice of divergence date estimates
- Choice of whether and how to correct for ancestral heterozygosity
Choice of reference species: Using all alignments from deep, intermediate, and recently diverged reference species, we estimated a whole genome substitution rate (\(\mu\)) following the approach outlined above. We used a Shapiro-Wilk normality test36,37 to determine that this dataset is normally distributed (p = 0.01). We then grouped whole genome substitution rates in two stratification schemes: four strata (by reference species) and three strata (by divergence time: deep divergence, intermediate divergence, and recent divergence), and used an ANOVA38 analysis to determine whether substitution rate varies significantly among groups in either stratification scheme.
Chromosome: We indexed the genomic position of each chromosome using a bed file generated from the annotated pygmy sperm whale reference genome, extracted SNP data on a per-chromosome basis, and estimated a substitution rate (\(\mu\)) for each chromosome of each genome included in the study as described above. We used a Shapiro-Wilk normality test to determine that this substitution rate dataset is not normally distributed (p = 1.6^{-6}). We therefore used a non-parametric Kruskal-Wallis analysis39 to determine whether substitution rate varies significantly among chromosomes across 18 species included in the study (excluding pygmy sperm whale), and a post hoc Dunn test to identify chromosomes driving any signal of differentiation. Finally, we tested for a relationship between chromosome-level substitution rates across species and chromosome length using a linear regression.
Inclusion of coding vs non-coding regions: We indexed the genomic position of coding and non-coding regions using a bed file generated from the annotated pygmy sperm whale reference genome, extracted SNP data in each region, and estimated substitution rates (\(\mu\)) in coding and non-coding regions of each genome included in the study. We used a Shapiro-Wilk normality test to determine that this substitution rate dataset is not normally distributed (p = 0.04). We then used a Kruskal-Wallis test to compare substitution rate estimates between these two regions across 18 species, excluding the pygmy sperm whale.
Maximum read-depth cutoff: We first grouped all substitutions into bins between read depths of 11-500 reads (i.e. 11-20, 21-30, … 491-500) and calculated the cumulative number of substitutions at each depth bin. Following this, we calculated substitution rate (\(\mu\)) for each species at maximum read depth cutoffs from 10-500. Finally, we use a generalized additive model to test whether read depth is a significant driver of variability in substitution rate estimates and a Kruskal-Wallis test to estimate the size of the effect of read depth on substitution rate estimates.
Divergence date estimates: We generated estimates of substitution rate (\(\mu\)) using fossil-dated divergence times (\(t\)) from McGowen et al.33 (6-partition AR model) and Lloyd and Slater40 (SAFE metatree analysis), as well as an average across all published divergence times from TimeTree41 (timetree.org, Supplemental Table S1). We then re-estimated whole genome substitution rate using each of the species-specific divergence dates. We used a Shapiro-Wilk normality test to determine that this substitution rate dataset is normally distributed (p = 0.15) and used an ANOVA to compare substitution rate differentiation across the three strata.
Ancestral heterozygosity: Because the ancestral heterozygosity (\(\pi_{anc}\)) of a species is unknown and unknowable using currently available technology, we generated substitution rate estimates under five different states of \(\pi_{anc}\). First, we estimated substitution rate (\(\mu\)) without correcting for \(\pi_{anc}\), and by correcting for \(\pi_{anc}\) under the assumption that the contemporary genomic heterozygosity of each species is equivalent to its ancestral heterozygosity, i.e., using a whole-genome estimate of heterozygosity from each species from Morin et al. (2025) to correct for \(\pi_{anc}\). We further tested the sensitivity of our estimates of \(\mu\) to a range of values of \(\pi_{anc}\) by re-estimating \(\mu\) for every species using the maximum and minimum reported values of \(\pi\) within cetaceans from Morin et al. (2025). We used a Shapiro-Wilk normality test to determine that this substitution rate dataset is not normally distributed (p = 3.29^{-5}) and used a non-parametric Kruskal-Wallis test with an ad hoc Dunn test and Bonferroni correction to compare substitution rate differentiation across the five different \(\pi_{anc}\) estimates.
Results
For each combination of species and reference genome, we generated 18 estimates of whole-genome substitution rate by varying estimates of divergence date, ancestral heterozygosity, and annual vs. generational substitution rate estimates. Using the pygmy sperm whale reference genome, divergence date from McGowen et al.33, and without correcting for ancestral heterozygosity, we also generated an estimate of substitution rate at each chromosome (n = 20 autosomes plus X chromosome), in both coding and non-coding regions (n = 2 per species), and at maximum read depths of 100-500 (n = 50 per species).
Whole genome substitution rate varies among individual species included in the study; rates largely fall between previously published average estimates for mysticetes and odontocetes (Figure 1).
Figure 1. Left: Whole genome annual substitution rate for each cetacean species included in the study, colored by species. Warm colors = mysticetes. Cool colors = odontocetes. Reference distance shown by shape. Species abbreviation codes as in Table 1. Average annual mysticete and odontocete substitution rate estimates from Dornburg et al. 2012 (UCLN model) shown as dashed and dot-dash lines, respectively. Images from phylopic.org. Right: Difference in substitution rates between odontocetes (blue) and mysticetes (yellow) across three divergence time strata.
Phylogenetic or demographic effects on substitution rate
We first examined the effect of phylogenetic inertia on substitution rate, and found that annual substitution rate was significantly different between the two cetacean suborders (ANOVA p = 0.077, F = 3.57, Figure 2) and accounts for approximately 18% of annual substitution rate variability, with substitution rates that were significantly lower in mysticetes than in odontocetes.
We then re-analyzed substitution rate between suborders using generational substitution rates, and found no significant differentiation (ANOVA p = 0.15, Figure 2). Further, the observed pattern in substitution rate differentiation between suborders is flipped, with a higher mean generational substitution rate in the mysticetes than in the odontocetes.
When examining the effect of maximum lifespan, asymptotic length, and generation time on annual substitution rate, while controlling for phylogenetic inertia, the model with the lowest AIC included only lifespan but it was not a significant parameter (p = 0.97). However, when considering generational substitution rate using maximum lifespan and asymptotic length as the predictor variables - again controlling for phylogenetic inertia - we found that the AIC-selected model included a significant positive effect of maximum lifespan (p = 1.4^{-5}, Figure 2), such that generational substitution rate increases with increasing maximum lifespan. Models including asymptotic length had higher AIC values than the model with lifespan as the only parameter; it follows that asymptotic length was not a significant predictor of generational substitution rate when considered as the only model parameter or in combination with maximum lifespan.
Figure 2. Top: Variability in annual (left) and generational (right) substitution rate estimated in two suborders. Bottom: Phylogeny-controlled linear trend in generational substitution rate by species’ asymptotic lifespan.
Pipeline effects on substitution rate
Reference species
The reference genome used had the greatest overall magnitude of effect among all factors tested in this study (Figure 1, Table 3). Reference genomes for this analysis were chosen to reflect three periods of evolutionary divergence: deep divergence, intermediate divergence, and recent divergence. Substitution rate per year was highest in most species when compared to the deep divergence reference genome. Alignment to the intermediate divergence reference genome resulted in the lowest substitution rate estimates across all species in the study, and alignment to the recent divergence reference genome resulted in substitution rate estimates that were between the deep and intermediate reference genomes in most cases. Yearly substitution rate estimates did not differ significantly across mysticetes and odontocetes when aligned to the deep reference genomes, but did vary significantly when aligned to the intermediate and recently diverged genomes, with higher annual substitution rates in odontocetes than in mysticetes (Figure 1).
Chromosome
Annual substitution rates varied significantly among chromosomes (Kruskal-Wallis p = 0). The range in substitution rates across chromosomes within each species (mean range = 1.5^{-10}) was larger than the range in whole genome substitution rates across all species included in the study (6.2^{-11}. Chromosome-level substitution rates were significantly higher in at least half pairwise comparisons in chromosomes 15, 16, and 19. The substitution rate was significantly lower in chromosome X than all other chromosomes (Figure 3). Broadly, relative patterns in substitution rate were consistent across chromosomes, with odontocete rates higher than mysticete substitution rates in all chromosomes, while the Amazon river dolphin and long-finned pilot whale substitution rates remained lower than all species across all chromosomes. Across all species, variability in substitution rate correlated significantly with chromosome size (linear regression p < 0.0001), with lower substitution rates in larger chromosomes (Figure 3).
Coding vs. non-coding regions
Annual substitution rate estimates were significantly different in coding vs. noncoding regions (Kruskal-Wallis p=9.4^{-7}, Figure 3). Rates were higher in coding regions than non-coding regions across all species included in the study, increasing by an average of 5.6^{-11} substitution/site/year (min = 4.8^{-11}, max = 6.4^{-11}). Similar to the relative patterns observed among species’ chromosome-level substitution rates, relative substitution rates were similar in coding and non-coding regions, with most odontocetes having higher substitution rates than most mysticetes and the Amazon river dolphin and long-finned pilot whale having the lowest substitution rates in both regions.
Figure 3. Top: Annual substitution rate across chromosomes aligned to the pygmy sperm whale reference genome. Red stars indicate chromosomes that were significantly different from at least 10 chromosomes. Bottom: Annual substitution rate in coding vs. non-coding regions of the genome. Points are colored by species in both subplots; warmer colors are mysticets and cooler colors are odontocetes.
Genotyping read depth
Across all species, annual substitution rate estimates increase rapidly with increasing maximum read depth cutoff to an asymptotic maximum substitution rate estimate (Figure 4). Most species reached their asymptotic substitution rate at a maximum read depth cutoff between 50 and 150 reads. GAM modeling indicated that read depth was a significant predictor of substitution rate estimate (p < 0.0001), and a Kruskal-Wallis test confirmed that substitution rates varied significantly among cumulative read depth bins (e.g. read depth of 10-20, 10-30, 10-40 … 10-500). This effect was driven primarily by significant pairwise differences between read depth bins with a cutoff < 100 and read depth bins with a cutoff > 100, where mutation rate estimates stabilized. Similar to previous tests, relative patterns among species remained the same with increasing read depth.
Divergence date estimates
Varying divergence date estimates also had a significant effect on annual substitution rate estimates, although the magnitude of the effect is smaller than the effect of reference genome or coding vs. non-coding regions (ANOVA p=0.01, F = 5.53, Figure 4). Divergence dates from Lloyd and Slater40 resulted in the highest within-group variability in substitution rate estimates among species, and relative substitution rates among species differed slightly among all three strata included in this test.
Ancestral heterozygosity
Varying ancestral heterozygosity also had a significant effect on annual substitution rate estimates (Kruskal-Wallis p=10^{-9}, Figure 4). This was driven by the group of substitution rates that were corrected using the maximum reported cetacean heterozygosity in Morin et al. (2025), which were significantly different from all other correction methods. We did not find any other significant differences across all other pairwise comparisons of correction methods.
Figure 4. Left: Annual substitution rate (substitutions/site/year) estimated with varying upper read depth cutoff for genotyping variants along the genome. Right top: Annual substitution rate (substitutions/site/year) estimated using divergence times from Lloyd and Slater 2021, McGowen et al. 2020, and Timetree.org. Right bottom: Annual substitution rate (substitutions/site/year) estimated using five approaches to accounting for ancestral heterozygosity: 1) no \(\pi_{anc}\), 2) \(\pi_{anc}\) equal to contemporary heterozygosity, 3) \(\pi_{anc}\) equal to maximum cetacean heterozygosity, 4) \(\pi_{anc}\) equal to mean cetacean heterozygosity, and 5) \(\pi_{anc}\) equal to minimum cetacean heterozygosity. Red stars indicate chromosomes that were significantly different in > 5 post hoc pairwise comparisons.
Comparison among all pipeline effects
We compared effect sizes (\(\eta^2\)) among strata for all pipeline tests (Table 3) to determine which parameters have the greatest effect on substitution rate estimates. Variability in reference species had the greatest effect on substitution rate estimates, and variability in phylogenetic distance to reference made up the majority of that effect. This was followed by variability among genomic locations (chromosomes) and coding vs. non-coding regions. Varying divergence time, maximum read depth, and the method for handling ancestral heterozygosity each had a comparatively modest effect on substitution rates, except when using the maximum value for heterozygosity reported by Morin et al. (2025).
| Pipeline effect | \(\eta^2\) |
|---|---|
| reference distance | 0.69 |
| reference species | 0.81 |
| chromosome position | 0.61 |
| coding/noncoding | 0.68 |
| read depth | 0.27 |
| divergence date | 0.18 |
| \(\pi_{anc}\) | 0.51 |
| \(\pi_{anc}\) (ex max \(\pi_{anc}\)) | 0.08 |
Discussion
Evolutionary patterns in mutation rate
Species-specific mutation rate has varied considerably over the evolutionary history of cetaceans, and at any point in time varies considerably across chromosomes and between coding and non-coding regions. In addition to ecological and evolutionary variability in mutation rate, we found that the way mutation rate is estimated will yield significantly different results, both in specific mutation rate estimates as well as in observed relational patterns among species.
Perhaps most notably, we observed phylogenetic patterns in annual substitution rate supporting previous assertions that mysticetes have lower annual mutation rates than odontocetes22,23. While the same phylogenetic patterns were not observed using generational mutation rates, we instead found support for higher generational substitution rates in longer-lived species, supporting recent pedigree-based estimates of high generational mutation rates in mysticetes14. While these results may seem contradictory, they highlight the broad range of lifespans within cetaceans and suggest that mutation rate has evolved in response to the evolution of longer lifespans. There may be evolutionary selection for lower annual mutation rates within the long-lived Mysticeti, but that phylogenetic trend does not fully counteract the effect of longer lifespans in these species, resulting in higher generational mutation rates in mysticetes even with lower annual mutation rates.
Another notable pattern is the large variability in substitution rate estimates driven by the phylogenetic distance to the reference species. Our results indicate that substitution rate has not varied linearly over evolutionary time; rather, we found that the lowest substitution rates were estimated in reference to intermediate divergence times, with higher rates in both recent and deep divergence times (Figure 1). This could be caused by an excess of mapping errors in the hippopotamus reference genome. However, scaffold N50 indicates similar quality in all reference genomes used in the study; further, if this observed pattern were purely an effect of decreased mapping accuracy to references of increasing divergence, we would expect that estimated substitution rate increases linearly with increasing divergence. Based on this evidence we hypothesize a period of high mutation rates deep within the evolutionary history of cetaceans and hippos, followed by a period of low mutation rates. Although cetacean mutation rate seems to have increased in recent evolutionary history, it remains lower than rates based on deeper divergence within the cetacean phylogeny. Further, substitution rate estimates integrated over deep evolutionary time were not significantly different between the two cetacean suborders, but a pattern of differentiation emerges when considering intermediate divergence time and increases in magnitude in recent divergence time (Figure 1).
Logistical considerations
In addition to the evolutionary and demographic effects we observed on substitution rate, we found that many pipeline parameters can affect substitution rate estimates. In general, altering pipeline parameters tended to shift all estimates of substitution rate in the same direction, with minimal observed shifts in the relative position of substitution rates among species.
Among pipeline effects on substitution rate, the choice of reference genome had the greatest effect on species-specific estimates of whole genome substitution rate, and the time since divergence from the reference genome made up a large portion of this variability. To our knowledge, this is the first report of a difference in substitution rate estimates based on the choice of reference sequence. Reference-driven differences in substitution rate observed in this study may be caused by several factors. Our results suggest that mutation rate accelerates and slows over time, and that cetacean mutation rates were higher in deep evolutionary time, slowed in mid-evolutionary time, and are marginally higher in shallow or contemporary evolutionary time. By comparing this pattern in mysticetes and odontocetes, we observed that the annual substitution rates of these two evolutionary lineages was not significantly different over deep evolutionary time, but differed significantly over intermediate and recent evolutionary time (Figure 1), suggesting that the annual substitution rates of each of these suborders may have embarked on differing evolutionary trajectories once they split from each other, resulting in the significantly lower annual substitution rates observed in contemporary mysticetes than in contemporary odontocetes.
The degree of variability driven by choice of reference genome - up to an order of magnitude in some cetacean species - brings with it the question of how to choose the correct reference genome for alignment of new genomic data. Because substitution rate estimates represent an average over the evolutionary time period since divergence from the reference species it may be best to choose a reference genome that reflects the time period of the question. For example, population genomics or other studies concerned with recent evolutionary time might choose a reference genome from a recently-diverged species, and estimate substitution rate based on the reference species for use in downstream analyses, e.g. PSMC or BEAST-generated divergence dates.
Compared to whole-genome substitution rate estimates among reference genomes, estimates of substitution rates among chromosomes varied less, but still significantly. This variability in substitution rate among chromosomes was driven at least in part by chromosome size, with a pattern of lower substitution rates in larger chromosomes that has been observed in other taxa as well42. The significantly lower substitution rate observed in Chromosome X is likely driven at least in part by its shorter residence time in males, which contribute more mutations to their offspring than females43,44. Interestingly, the opposite pattern is observed in birds, where the Z chromosome that has a shorter residence time in females nonetheless has a higher mutation rate than autosomal chromosomes45. Finally, we observed higher substitution rates in coding regions than non-coding regions. These results, taken together, highlight the importance of choosing appropriate genomic regions when estimating mutation rates for use in conservation or evolutionary genomic studies.
Variability caused by other pipeline parameters (i.e. read depth, divergence date, and ancestral heterozygosity) were among the smallest effect sizes in the study. The effect of maximum read depth cutoff was primarily limited to maximum read depths < 150. The plateau in substitution rates at maximum read depths > 150 indicates that false polymorphisms are not likely to be introduced by genotyping error in regions with anomalously high read depths (up to 500). Conversely, capping read depth at a cutoff less than ~100-150 may negatively bias mutation rate estimates. Variation in divergence date estimates also had comparatively minimal effect on substitution rate variability, which indicates a high level of precision and convergence among our current estimates of divergence dates among cetacean species. However, divergence date was the only parameter to affect the relative position of species-specific mutation rate estimates, caused primarily by disagreement between published models and data regarding the divergence of Amazon river dolphin and long-finned pilot whale from the pygmy sperm whale.
Correcting for species-specific ancestral heterozygosity has a linearly predictable effect on substitution, because \(\pi_{anc}\) is simply subtracted from estimated genomic divergence before dividing by divergence time. We expect that subtracting any value from genomic divergence will cause a significant effect on the estimated substitution rate; the magnitude of the effect is determined by the size of \(\pi_{anc}\). Our dataset confirmed this expectation. Further, we found that correcting for species-specific genomic heterozygosity was comparable with mean and minimum contemporary cetacean genomic heterozygosity; correcting for heterozygosity in this way has a minimal, but significant, effect on substitution rate estimates. However, our estimates of contemporary cetacean genomic heterozygosity ranged by over an order of magnitude (min = 1.1^{-4}, max = 0.005037), and correcting for ancestral heterozygosity using larger values had a much larger effect on substitution rate estimates. Because species-specific ancestral heterozygosity is currently unknowable, it is important to carefully consider whether and how to account for ancestral heterozygosity when estimating species-specific substitution rates, and to clearly describe this correction in the methods of genomic studies.
Although not directly targeted by this study, we found that substitution rates for the same species, estimated using the same method, differ based on whether they are estimated annually vs. generationally, with potentially profound effects on downstream interpretation of evolutionary relationships or conservation applications. Recent disagreement about whether mutation rates in mysticetes are larger or smaller than odontocetes14,e.g., 22 may be at least partially explained by differences in the use of annual vs. generational mutation rate, and future studies will need to carefully consider the implications of each of these for the interpretation of their results.
Subsitution rate vs. pedigree-based mutation rate estimates
Recently, the development of a method to estimate mutation rate by quantifying the number of de novo germline mutations passed down through multi-generational family pedigrees has produced estimates that are considered more accurate and robust than estimates based on phylogenetic relationshipse.g. 4,14. While this method remains subject to some of the same sources of error (e.g. sampling bias and postmortem DNA damage, genotyping error), pedigree-based mutation rate estimates allow for direct quantification of inherited mutations without the need for the estimation of parameters such as divergence date or generation time. Although the pedigree approach is a more direct way of estimating mutation rate, for most taxa - particularly elusive marine species - parent and offspring genomes are not available for direct comparison. Further, while using this method with a large sample size can produce an accurate estimate of generational mutation rate, translating that to an annual mutation rate estimate would require knowledge of the age of parents, which is often difficult or impossible to obtain. Pedigree-based estimates will reflect contemporary population demographics and the age of the parents, and may differ significantly from historical mutation rates; however, recent studies comparing the two methods indicate general agreement when comparisons are drawn between generational mutation rates only46.
Our study includes two species - the blue and humpback whales - for which pedigree-based mutation rate has also been estimated14. While sample size was limited in this study (N = 1 blue whale pedigree and N = 5 humpback whale pedigrees), the study provides good context for comparison among pedigree-based mutation rate estimates and phylogeny-based substitution rate estimates.
Interestingly, the phylogeny-based estimate of generational substitution rate most similar to pedigree-based estimates was obtained when blue whale sequence data was were aligned to the deeply-diverged hippopotamus (\(\mu_{ped}\) = 1.24 x 10-8, CI: 0.85 x 10-8 to 1.63 x 10-8; mean hippo \(\mu_{phy}\) = 1.4220435^{-8}). All phylogeny-based estimates of generational whole genome subsitution rate in blue whales aligned to the North Atlantic right whale, regardless of divergence date or ancestral heterozygosity, were within the 95% confidence interval of the pedigree-based estimate from Suarez-Menendez et al.14. When aligned to the pygmy sperm whale, phylogeny-based generational substitution rate estimates (mean \(\mu_{phy}\) = 6.7435825^{-9})) were just below the 95% confidence interval from Suarez-Menendez et al.14. Given that the estimates in Suarez-Menendez et al.14 were made using only one pedigree, it is possible that adding additional pedigrees will increase variability in pedigree-based estimates of mutation rate.
Similarly, across five humpback whale pedigrees, Suarez-Menendez et al.14 estimated a pedigree-based generational mutation rate of 1.12 × 10−8 (CI: 0.94 × 10−8 to 1.30 × 10−8). Again, the mean phylogeny-based generational substitution rate estimate for the humpback whales genome to the hippopotamus reference (1.17 x 10-8) was in close alignment with the pedigree-based mutation rate estimate. Mean phylogeny-based estimates of substitution rate were smaller by half when the humpback whale was aligned to the pygmy sperm whale (5.03 x 10-9) or North Atlantic right whale (6.9 x 10-9).
Conclusions and Recommendations
Species-specific mutation rates evolve over time; cetacean mutation rates are affected by phylogeny and maximum lifespan, which are highly correlated within the order. Mysticetes, which trend toward longer-lived species, have lower annual mutation rates but similar generational mutation rates to odontocetes, highlighting the importance of using the same type of mutation rates when drawing comparison across species. Because mutation rate varies over time within an evolutionary lineage, the choice of reference species will have a large affect on estimates of heterozygosity and mutation rate; therefore, reference species should be carefully selected to reflect the timespan of the research question. While divergence date and ancestral heterozygosity estimates had smaller effects on mutation rate, we found evidence that the mutation rate varies along the genome; therefore the genomic regions included in a study - e.g. coding vs. non-coding regions or choice of chromosomes - will also affect our estimates of heterozygosity and mutation rate. In most population genetic studies, it is recommendable to include a large number of sites across the genome to better represent the genomic average. In targeted studies, it may be more effective to estimate heterozygosity or mutation rate across a carefully selected subset of the genome. Given the phylogenetic, demographic, genomic, and pipeline effects on variability in mutation rate demonstrated in this study, it is vital that genomic studies carefully consider how to process sequence data within the context of the targeted research question, and give consideration to the potential effects of pipeline parameters on the interpretation of their results.
Because substitution rate variability among species within each suborder was minimal, we report mean substitution rates estimated using the pygmy sperm whale reference genome, divergence from McGowen et al. (2020), and no correction for ancestral heterozygosity for both odontocetes (annual = 2.49 x 10-10, generational = 4.8 x 10-9) and mysticetes (annual = 2.3 x 10-10, generational = 5.8 x 10-9). These estimates may be appropriate for use in studies of population genomics or “recent” evolutionary genomics. However, we recommend that any studies involving species not included in this study estimate a mutation rate using the method employed here; scripts for counting substitutions and estimating substitution rate are available at https://github.com/UW-WADE-lab/Mutation-rate-variability-in-cetaceans. Specific rates for each species and reference genome included in this study are reported in Supplemental Table S2.
ACKNOWLEDGEMENTS
We appreciate the expert input and review of early drafts of this manuscript by multiple colleagues, including Dr. Annabel Beichman and Dr. Michael McGowen. High-quality genomic sequence data used in this study were made publicly available by the Cetacean Genomes Project in collaboration with the Vertebrate Genomes Project.
DATA AVAILABILITY
All genomic sequence data and reference genomes were accessed via NCBI or ENA using the accession numbers listed in this manuscript. All genome alignment scripts are available at https://github.com/UW-WADE-lab/Whole-Genome-Alignment, and all scripts used in mutation rate estimation, data analysis, and manuscript generation are available at https://github.com/UW-WADE-lab/Mutation-rate-variability-in-cetaceans.