The distribution of fitness effects (DFE) describes the selection coefficients of newly arising mutations and fundamentally influences population genetic processes. However, the extent and mechanisms of differences in the DFE for non-synonymous mutations have not been systematically investigated across species with divergent phylogenetic histories and ecologies. Here, we inferred the DFE in natural populations of 11 animal (sub)species, including humans, mice, fin whales, vaquitas, wolves, collared flycatchers, pied flycatchers, halictid bees, Drosophila, and mosquitoes. We found that mammals have a higher proportion of strongly deleterious mutations (defined as ; 22% to 47% in mammals; 0.0% to 5.4% in insects and birds) and a lower proportion of weakly deleterious mutations than insects and birds. Further, the DFE co-varies with phylogeny, such that the mean mutation effects are more similar in closely related species (Pagel’s = 0.84, P = 0.01). Next, we investigated whether various summary statistics of the DFE were related to variation in life-history traits across these organisms. We found some support for genome size, body mass, and long-term effective population size being correlated with the DFE. Overall, our findings are consistent with predictions derived independently from the Fisher’s Geometric Model (FGM), which defines organismal complexity as the number of phenotypes under selection. FGM predicts that mutations are more deleterious in complex organisms, while strongly deleterious mutations occur more frequently in smaller populations. Our study demonstrates strong phylogenetic signal in the evolution of a fundamental population genetics parameter, and proposes that, through mechanisms of epistasis, long-term population size and organismal complexity could be underlying variation in the DFE across animals.

Population genetic inference for 11 animal (sub)species

We retrieved high-quality population-level polymorphism datasets for 11 animal (sub)species. These data include three insect species (mosquitoes [33], Drosophila [7], and halictid bees [34]), two bird species (pied flycatchers and collared flycatchers [35]), and six mammal (sub)species from five species (arctic wolves [36,37], gray wolves [38], vaquitas [28], fin whales [39], mice [22], and humans [7]; arctic wolves and gray wolves are subspecies of wolves; Table 1, Text A in S1 Appendix). To further understand the mechanisms of DFE evolution, we also compiled several life history traits and genomic attributes for each species from literature, including mutation rate (Text B in S1 Appendix), long-term population size, body mass, generation time, age at maturity, maximum longevity, genome size and recombination rate (sources and citations in S1 Data). The 11 (sub)species we evaluated are distinctive in phylogenetic position (Fig 1A) and the compiled biological features that could affect the DFE. For example, the literature-derived long-term population sizes () varied from 4,485 (vaquitas) to (Drosophila) individuals. All datasets were filtered with the same stringency to retain only high-confidence genotype calls and at least eight diploid high-quality individuals in coding regions. We tallied the number of variants at different minor allele frequencies and generated folded site frequency spectra (SFS) for synonymous and nonsynonymous/missense variants (SYN-SFS and MIS-SFS, respectively) for each species (Fig A in S1 Appendix, S2 Data). To implement a robust but flexible workflow to compare the DFE across species, we developed varDFE (this study), a Python package based on the and Fit packages (Fig B in S1 Appendix; [20,40]).

We inferred demographic parameters from the SYN-SFS (Fig C in S1 Appendix, S3 Data). Although the demographic models we infer from synonymous variants could be biased by selection at linked sites, these demographic models nevertheless control for the distortion in the MIS-SFS caused by selection at linked sites and population history when inferring the DFE parameters [20]. Demographic models were selected using a stepwise likelihood-based procedure (Fig D and Text A in S1 Appendix). A demographic model with a single change in population size (i.e., a two-epoch model) provided a good fit to SYN-SFS and reasonable parameter estimates for most species except for the mosquitoes, halictid bees, gray wolves, and collared flycatchers. For the gray wolf and collared flycatcher datasets, the three-epoch demographic model (i.e., two changes in population size) performed better with more realistic parameter estimates or drastically improved fit to the SYN-SFS. For the mosquitoes and halictid bees datasets, no demographic model using the full SYN-SFS fit the data well, and using the singleton-masked SYN-SFS substantially improved model fit (S4 Data, Text A in S1 Appendix). All the inferred demographic models are qualitatively consistent with previous estimates (S5 Data). There is a strong correlation between literature-based estimates of and estimates of long-term population size from synonymous variants ( = 0.94, Fig E in S1 Appendix), suggesting that selection at linked sites has a limited effect on demographic inferences. We utilized the best-fit demographic parameters from the chosen models (Table 1) for downstream DFE inference.

Mutations are more deleterious in mammals compared with insects and birds

Conditional on the inferred demography, we estimated the DFE for new nonsynonymous mutations for each species, assuming that the DFE follows a gamma distribution [20], mutations are neutral or deleterious (), and mutations act in an additive manner (h = 0.5). The gamma distribution yielded a robust fit to the observed MIS-SFS across all datasets (Fig F in S1 Appendix, S6 Data).

Overall, mutations are 4.8 to 4,923 times more deleterious in mammals (the mean selection coefficient, |s|_mean, is 0.010 in mice to 0.67 in arctic wolves) compared to insects (|s|_mean = 0.00014 in Drosophila to 0.00056 in mosquitoes) and birds (|s|_mean = 0.00052 in pied flycatchers and |s|_mean = 0.0021 in collared flycatchers; Fig 1B). Noticeably, the scale () parameter for the arctic wolves reached the upper boundary of the parameter space during inference. To rule out potential data artifacts, we analyzed an independently generated, whole-genome dataset from Russian Karelia gray wolves [41], which also yielded a high value, confirming |s|_mean is very deleterious in wolves. Since |s|_mean may be heavily influenced by the strongly deleterious tail of the DFE, which we have less ability to accurately estimate, we also examined the median |s| (|s|_median) values across species. Consistent with patterns in |s|_mean, median selective effects are more deleterious in mammals (|s|_median = 0.00034 in vaquitas to 0.0062 in arctic wolves) compared to insects and birds (|s|_median = in Drosophila to 0.00026 in collared flycatchers; Fig G in S1 Appendix). The shape () parameter of the gamma distribution differs between species groups as well. On average, mammals have lower estimates for the shape parameter ( = 0.11 in arctic wolves to = 0.21 in mice) compared to insects and birds ( = 0.22 in collared flycatchers to = 0.36 in Drosophila; Fig 1C). The extent of variation in that we detected is in line with previous findings [29].

We also observed that mammals have higher proportions of strongly deleterious mutations () than insects and birds, consistent with these species having more deleterious |s|_mean estimates (Fig 1D). The maximum-likelihood gamma distribution for each species revealed that 22% (vaquitas) to 47% (arctic wolves) of mutations in mammals are strongly deleterious (), compared to 0.0% (Drosophila) to 5.4% (collared flycatchers) in insects and birds. On the other hand, insects and birds have a larger proportion of weakly deleterious mutations (, 41% in collared flycatchers to 68% in Drosophila), compared to mammals (16% in arctic wolves to 30% in mice). The proportions of neutral to nearly neutral mutations () are similar in all 11 (sub)species. To directly estimate the proportion of (nearly) neutral mutations, we considered a different functional form, where DFE followed a mixture of gamma distribution with an additional point mass at neutrality (neugamma DFE [7], eq 1, S7 Data). This neugamma DFE improved model fit for most datasets (lower AIC scores) compared to the gamma DFE (Figs H and I in S1 Appendix), except for mosquitoes, gray wolves, and vaquitas. Similar to the patterns obtained assuming a gamma DFE (Fig 1D), the inferred proportion of neutral mutations (p_neu) varied across lineages (Figs 1E, Fig J in S1 Appendix). p_neu varied from 0.081% (mosquitoes) to 16% (pied flycatchers) in insects and birds, and from 0.089% (vaquitas) to 28% (humans) in mammals. We also attempted to directly estimate the proportion of lethal mutations () [20,42]. However, our approach likely lacks sufficient power to confidently infer the proportion of lethal mutations, due to poor model fit and unrealistic estimates of the proportion of lethal mutations (Figs I–K and Text C in S1 Appendix, S8 Data).

To evaluate the robustness of our analysis, we tested whether several technical treatments affected the DFE inference. First, we examined different probability distributions for the DFE, including the gamma, neugamma, and lognormal distributions (Figs L and M in S1 Appendix, S9 Data), and compared the median selective effects (|s|_median) across distributions. When assuming a lognormal DFE, we focused on comparing median (|s|_median) rather than mean selective effects (|s|_mean), as the inferred distribution includes a long tail of strongly deleterious mutations. In this case, the median provides a more robust and interpretable summary. When assuming the neugamma DFE, |s|_median remained more deleterious in mammals compared with insects and birds (Fig G in S1 Appendix). When assuming the lognormal DFE, this pattern was largely maintained, although the |s|_median estimate in vaquitas (|s|_median = 0.00029) became lower than those in mosquitoes (0.00034) and collared flycatchers (0.00037; Fig G in S1 Appendix). However, the AIC for the lognormal DFE is the lowest only in four datasets (mosquitoes, gray wolves, arctic wolves, and vaquitas; Fig N in S1 Appendix). Furthermore, the lognormal distribution exhibits a long tail of strongly deleterious mutations (Figs L and O in S1 Appendix), where there is the least confidence of -based DFE inference [42]. Both of these factors suggest that a lognormal-distributed DFE is less supported by the data. The gamma distribution was confirmed as an adequate functional form for DFE comparisons, given its parsimony, and its second-lowest AIC across most datasets (Fig N in S1 Appendix), and good visual fit to the observed SFS (Fig F in S1 Appendix). Next, we tested the potential impacts of demographic model misspecification, singleton-masking treatments, mutation rate variation, differences in nonsynonymous to synonymous mutation rate ratio (), and GC-biased gene conversion, on the inference of the DFE. The average mutation effects and proportion of each mutation category remain unchanged qualitatively across these different technical treatments (Figs P-U and Text D in S1 Appendix, S10 Data - S14 Data).

To confirm our DFE inference can confidently detect strongly deleterious mutations () even in very large populations, we performed a SLiM simulation using demographic parameters inferred from Drosophila and DFE parameters inferred from humans. The long-term population sizes () inferred from synonymous and non-coding mutations in the simulation are both lower than the true Drosophila population size used, as expected due to selection at linked sites (Fig V in S1 Appendix). Using a sample size of 100, was underestimated by 3.8%–5.0% when using the SFS from synonymous variants and by 2.8%–3.5% when using the SFS from non-coding variants. With a sample size of 10, underestimation was 3.1%–5.2% and 1.7%–2.5%, respectively. The DFE inferred from the simulation closely matches the true human DFE used in the simulation, where the assumed human |s|_mean is 0.014 (Table 1). The average inferred |s|_mean is 0.014 (95% CI: 0.012–0.017) when using an SFS sample size of 100, and 0.015 (95% CI: 0.010–0.023) with a sample size of 10 (Fig W in S1 Appendix). This simulation indicates that our methods can confidently detect strongly deleterious mutations () and correctly infer |s|_mean in large populations, even when estimates of long-term population size from synonymous variants are biased due to selection at linked sites. To assess whether accurate inference is robust to misspecification of the assumed functional form of the DFE, we also fit a lognormal distribution to the simulated data, and compared its median to the true human DFE used to generate the data. The average inferred |s|_median assuming a lognormal-distributed DFE is 0.014 (95% CI: 0.012–0.017) when using an SFS sample size of 100, and 0.015 (95% CI: 0.010–0.023) with a sample size of 10. These estimates are larger than the true gamma-distributed human |s|_median of 0.0012. However, the direction of the bias is away from the more neutral Drosophila-like DFE, suggesting that estimates of |s|_median are not under-estimated in extremely large populations (Fig X in S1 Appendix). Thus, the lack of strongly deleterious mutations inferred in species with large population sizes is likely not due to biases in our inference procedure, and instead reflects true differences in the DFE.

In summary, the distribution of fitness effects is more similar in closely related species, with mammals harboring more strongly deleterious and less weakly deleterious variations compared with insects and birds.

The DFE co-varies with phylogeny in animals

The apparent co-variation of the DFE with phylogeny motivated us to test for a phylogenetic signal more formally. Pagel’s quantifies how trait evolution depends on the underlying phylogenetic tree across species, where = 0 suggests no phylogenetic signal and = 1 indicates strong phylogenetic dependence [43]. Assuming a gamma-distributed DFE, we calculated Pagel’s for log10 transformed mean selection coefficients [(|s|_mean)], log10 transformed median selection coefficients [(|s|_median)], and the shape parameter (, Fig 1B – 1C). For (|s|_mean), Pagel’s was 0.84 (P = 0.01, likelihood ratio test [LRT]). A strong phylogenetic signal was also observed for (|s|_median) ( = 0.76, P = 0.03) and the shape parameter ( = 0.86, P = 0.004). Since the DFE parameters for arctic wolves reached upper boundaries, Pagel’s was recalculated excluding the arctic wolves dataset. Strong phylogenetic signal was consistently observed for all metrics ( = 0.86, P = 0.008 for (|s|_mean); = 0.81, P = 0.02 for (|s|_median); = 0.81, P = 0.01 for ). To account for the long phylogenetic branches between mammals and insects, we repeated Pagel’s analysis only within mammals. Signals of phylogenetic dependence were reduced ( = 0.45, P = 0.49 for (|s|_mean); = 0.15, P = 0.88 for (|s|_median); = 0.90, P = 0.08 for ), which could either be due to less of a signal within mammals or reduced statistical power due to the small number of species analyzed. Overall, these analyses suggest that the evolution of DFE is influenced by the shared evolutionary history determined by the phylogenetic structure, but phylogenetic dependence is less evident within mammals (S15 Data).

To formally test whether the parameters of a gamma-distributed DFE vary across species, we examined a null model where the shape and scale parameters of the gamma distribution were constrained to be the same across the 11 (sub)species datasets, while still accounting for their individual demographic histories. Using a grid-search approach, we calculated the Poisson likelihood of the expected MIS-SFS given each parameter combination in each species. The maximum likelihood estimate (MLE) for this null model is found by summing the individual likelihood surfaces of each species’ DFE (Fig 2A). The alternative model allows each species to have its own shape and scale parameters, as described in previous sections (Fig 2B – 2C). This alternative model fits the data (MIS-SFS) significantly better than the null model (LRT, = 18,203, df = 20, P < 10⁻¹⁶; Figs 2A – 2C, Fig Y in S1 Appendix). Therefore, the variation of DFE across taxa is statistically supported by our analyses.

To further investigate the relationship between phylogenetic proximity and DFE similarity, we next used the likelihood ratio test framework to compare pairs of species. Here, the null model is relaxed to constrain DFE parameters to be identical only for the two species compared, rather than across all species. The resulting LRT statistic () serves as a measure of DFE similarity between species pairs, with lower values indicating higher similarity (eq 2). To visualize pairwise comparisons, we conducted hierarchical clustering of log10-scaled pairwise LRT statistics []. The insects and mammals clustered into two main lineages, with birds being split between the two: the pied flycatchers fell into the insect lineage while the collared flycatchers clustered with vaquitas (Fig 2D). In addition, is significantly lower (P = 0.02, Wilcoxon test) in pairwise comparisons within the same taxonomic class (average = ), such as mammal-mammal pairs, compared to inter-class pairwise comparisons (average = ), such as mammal-insect pairs (Fig 2E). In summary, species that are more closely related phylogenetically, tend to have more similar DFEs compared to species that are more distantly related to each other.

The relationship between life history traits and the DFE

We next explored the relationship between the DFE and candidate life history traits related to long-term population size and organismal complexity, using the phylogenetic generalized least squares model ( [44]). estimates the amount of phylogenetic signal () in the regression residuals, and incorporates in regression parameter estimates, to control for interspecific autocorrelation due to phylogeny. Using this framework, we examined the relationship between the gamma-distributed DFE parameters (|s|_mean, |s|_median and the shape parameter) and nine different life history traits. First, we explored the relationship between the DFE and long-term effective population size. While it would seem natural to test for a correlation between |s|_mean and , was used to convert to s, which could induce a negative correlation between |s|_mean and (Text E in S1 Appendix). To partially mitigate this effect, we instead used published independently-derived estimates of long-term population size made from noncoding variants (Fig E in S1 Appendix). We observed a marginally significant negative correlation (Fig 3A) between the literature-based long-term population sizes and log10 transformed mean selection effects , (, slope (b) = −0.90, = 0, P = 0.02), though this correlation was no longer significant after Bonferroni correction for nine traits (P_adj = 0.2). On the other hand, statistically stronger positive correlations were observed for body mass and genome size with |s|_mean (body mass: , b = 0.26, = 0, P = 0.002 [0.02]; genome size: , b = , = 0, P = 0.002 [0.02]; Fig 3B – 3C). Other parameters tested, including the generation time, age at maturity, maximum longevity, recombination rates, and mutation rates, were not significantly correlated with |s|_mean after considering phylogenetic dependence (Fig Z in S1 Appendix). The two wolf (sub)species’ DFE deviated the most from the regression, likely due to their very deleterious estimated DFEs.

As |s|_mean is sensitive to the value of to unscale the scale parameter in the gamma distribution, we next focused on the shape () parameter of the gamma distribution. The shape parameter is not influenced by (Text E in S1 Appendix), and often estimated more confidently than the scale parameter in DFE inference [20,29]. We performed the analysis between and life history traits (Fig AA in S1 Appendix). Literature-based long-term population size was positively correlated with (P = 0.003 [0.02 after bonferroni-adjustment], Fig 3D), whereas body mass and genome size were negatively correlated with (P = [0.001], P = 0.001 [0.01], respectively; Fig 3E – 3F). Although not significantly correlated with |s|_mean, age at maturity is negatively correlated with (P = 0.004 [0.04], Fig AA in S1 Appendix). We performed the analysis between and life history traits as well (Fig AB in S1 Appendix). The correlation with is less pronounced, with synonymous-site-based, but not literature-based, long-term population size weakly correlating with (P = 0.04 [0.4]). Body mass and genome size remain correlated with (P = 0.006 [0.06], P = 0.001 [0.01], respectively). We also repeated the analyses excluding the arctic wolves outlier. Literature-based long-term population size, body mass, and genome size remain correlated at various degrees with |s|_mean, |s|_median and (S16 Data). To account for the long phylogenetic branches between mammals and insects, we repeated the analyses only within mammals. None of the traits tested were significantly correlated with |s|_mean, |s|_median or after Bonferroni-adjustment. Only a marginally significant correlation between the shape parameter and genome size was observed (P = 0.01 [0.1]; Figs AC–AE in S1 Appendix).

The marginally significant negative correlation of population size and mean mutation effects (|s|_mean) suggests potential constancy of population-scaled mutation effects across species () as previously hypothesized [8]. Therefore, we next sought to determine if population size by itself can explain the variation in the DFE observed across species by testing whether the population-scaled DFE varies across species. Assuming the population-scaled DFE follows a gamma distribution, the mean mutation effects () still vary across species (Fig 3G–3H). While |s|_mean co-varied with phylogeny, this same phylogenetic signal was not observed in estimates of ( = 0, P = 1 for ()). The mean population-scaled mutation effects are the highest in two wolves populations ( = in arctic wolves and 8,120 in gray wolves), and the lowest in the vaquitas (161) and humans (195; Fig 3G). The proportion of mutations with different population-scaled selection coefficients varied substantially within Classes as well (Fig 3H). For example, humans were estimated to have only 5.12% mutations with , while mice have 41.8% of mutations in this category. Adopting the grid-search analysis for population-scaled DFE comparisons, we found that the alternative model where each species has its own population-scaled shape and scale parameters fit the data significantly better than the null model assuming a constant in different species (LRT, = 99,508, df = 20, P < 10⁻¹⁶; Fig AF in S1 Appendix). When comparing the population-scaled DFE between species pairs, for remained significantly lower (P = 0.01, Wilcoxon test) in pairwise comparisons within the same taxonomic class ( = ) compared to inter-class comparisons (; Fig AG in S1 Appendix), as observed in the statistics for s (Fig 2E). However, we repeated the analyses for , and none of the nine traits tested was significantly correlated with (Fig AH in S1 Appendix).

Fisher’s geometric model enables concurrent examination of the influence of population size and organismal complexity on the DFE

Previous work has suggested that Fisher’s geometric model (FGM) can explain patterns of DFE variation across divergent species. FGM models the fitness landscape with a minimal set of parameters to study the emergent properties of mutational effects [14]. In FGM, the fitness of an organism is characterized by a multi-dimensional Gaussian decay function with a local optimum in n-dimension phenotype space, where the phenotype dimensionality n is termed “organismal complexity". Mutations create new phenotypes, modeled as random vectors with effect scale that act within a subset of m phenotype dimensions (, m is defined as mutation pleiotropy) [10]. Both m and reflect intrinsic features of an organism’s genotype–phenotype-fitness architecture, and serve as proxies for organismal complexity from different perspectives [7,16]. When the population is at the fitness optimum and mutations are universally pleiotropic (m = n), the derived DFE is a gamma distribution with only neutral or deleterious mutations () [16]. When the population is under mutation-selection-drift equilibrium, the increase in fitness from beneficial mutations should counteract the drift load, i.e., the decrease in fitness caused by fixed deleterious mutations, and the population has an equilibrium phenotypic distance () to the fitness optimum. The formula for this FGM-derived DFE, which models the full DFE (), is a function of three key parameters in FGM (eq 3): long-term population size (), mutation pleiotropy (m), and scale of mutation effects ().

We inferred the parameters for this FGM-derived DFE for 11 animal (sub)species (eq 3, derived in [10]). The FGM-derived DFE provided an equally good or better fit to the MIS-SFS compared to the gamma DFE (Figs AI and AJ in S1 Appendix, S17 Data). In addition, the inferred long-term population sizes from the FGM-derived DFE () largely agree with the estimates of derived from demographic inference using genetic variation data ( = 0.36, linear regression; P = 0.03, ; Fig 4A), providing additional validation of this model. However, estimates are often higher than , which could be due to the different biological interpretations in these two estimates. is inferred from demographic modeling using putatively neutral polymorphisms, whereas is inferred from deleterious variation at nonsynonymous sites. reflects the equilibrium distance from the fitness optimum, and likely captures a much longer time period than the neutral coalescent time that represents [45]. The vaquita population’s estimate deviated the most between the two methods, likely reflective of its long-term small population size and recent bottleneck [28].

We found some support that the mutation effect scale () is reflective of a species’ phylogenetic position. Mammals have larger (0.077 to 9.5) compared with insects and birds (0.0037 to 0.052; Fig 4B). The mutation pleiotropy (m) is negatively correlated with the inferred ( = 0.76, linear regression, P < 0.001, ; Fig AK in S1 Appendix), suggesting potential auto-correlations during inference. Although less conserved within lineages, mammals tend to have smaller m (0.17 to 0.60) compared with insects and birds (0.49 to 4.5). However, Pagel’s did not suggest phylogenetic signal for either ( = 0, P = 1) or m ( = 0.51, P = 0.67; S15 Data).

Next we examined the proportion of mutations with various ranges of s predicted from the FGM-derived DFE (Fig 4D). Similar to the findings assuming the gamma DFE, strongly deleterious mutations are scarce in insect and bird populations (0.0% to 2.0% in collared flycatchers), but much more abundant in mammals (9.1% in mice to 55% in arctic wolves). Additionally, we found that the proportion of beneficial mutations varied among populations. FGM predicts that as the population size decreases, the proportion of beneficial mutations increases to counteract the increased drift load [10]. Indeed, assuming FGM, we found that the proportion of beneficial mutations is negatively correlated with the long-term effective population size taken from the literature ( = 0.73, linear regression; P < 0.001, ; Fig 4C), providing further support for this model. Note that here we are relating the proportion of beneficial mutations to the effective population sizes estimated from the literature, and as such, this correlation is not directly induced by the FGM being fit to the data.

Data sets and site frequency spectra

Polymorphism datasets from 11 animal (sub)species were compiled from various whole genome resequencing projects using a standardized curation workflow designed to retain only high-quality variants and samples for DFE inference (Text A in S1 Appendix). Most datasets had pre-annotated variant effects. If not, we annotated the variant call files (VCF) using default options in snpEff (v.5.1 [75]). To retain only high-quality variants, we selected genomic sites that: 1. passed the filters set in the original dataset; 2. were in canonical coding region (CDS) but not in soft-masked, unknown or CpG islands of the genome; 3. were monomorphic or biallelic single-nucleotide polymorphisms (SNP) that were annotated as “synonymous" or “nonsynonymous/missense". To retain only high-quality samples from as close to a panmictic population as possible, we selected at least eight diploid samples for each species that: 1. had at least 8x average sequencing depth; 2. belonged to one natural population with no population substructure or admixture identified by principal component analysis (PCA); 3. were not related to other samples identified by kinship analysis.

To summarize polymorphism data, the obtained VCFs were further filtered to exclude sites with more than 20% missing or 75% heterozygous genotypes. We projected down the sample size and computed folded synonymous and nonsynonymous/missense site frequency spectra (SYN-SFS and MIS-SFS, respectively) for each species. Computing a folded SFS mitigates uncertainties in ancestral state classifications. The projected sample size was calculated using the hypergeometric probability distribution to maximize the number of SNPs available and account for the sporadic missing genotypes (S1 Data). The total synonymous and nonsynonymous sequence lengths (L_SYN and L_MIS) were obtained using previous estimates of nonsynonymous to synonymous mutation rate ratio [7].

We also compiled several life history traits and genomic attributes for each species from the literature (sources and citations in S1 Data).

Demographic inference

To control for how population size changes influence the SFS, we first inferred demographic parameters from the SYN-SFS using the Demog1D_sizechangeFIM module in varDFE (v.0.1.0, Fig B in S1 Appendix, based on v.2.1.1 [20,40] and Python v.3.10.2). For each species, we fitted four models, including the standard neutral, two-epoch, three-epoch, and four-epoch models. To account for the uncertainties in genotype calls, we repeated the demographic inference by masking the singleton entries in the SYN-SFS. In total, eight demographic inference runs (four demographic models with or without masking singletons) for each species were conducted. All settings were the same across datasets except for the starting parameter positions, which were set based on knowledge from prior demographic inference (S1 Data, S3 Data). For each run, the best-fit parameters with the maximum multinomial log-likelihood in 100 replicates were chosen, and their uncertainties were estimated through Fisher’s Information Matrix (FIM).

To choose the best-fit demographic model and singleton-masking treatment out of the eight runs per (sub)species for downstream DFE inference, we performed a stepwise likelihood-based model selection procedure (Fig D and Text A in S1 Appendix, S4 Data). At the start of model selection, the two-epoch demographic model based on full SYN-SFS run was chosen for all datasets. However, if it resulted in poor model fit, we chose a more complex demographic model or masked the singletons in the SYN-SFS. We calculated the population-scaled synonymous mutation rate , using the optimal_sfs_scaling function in . is only equivalent to Watterson’s under the standard neutral model, and changes with demographic model choice [40]. As exon mutation rates are hard to estimate in non-model organisms, and existing estimates can vary [64], we compiled a list of literature-reported mutation rates, and selected the most widely used, phylogenetic-based and whole-genome-based mutation rate estimates in each species (Text B in S1 Appendix). For each dataset, synonymous-site-based ancestral population size, or long-term population size, is defined as the reference population size in , and is estimated using in the chosen demographic run for downstream analyses (Table 1).

DFE inference

Conditional on the inferred demographic models, we estimated the DFE for new nonsynonymous mutations using methods from Fit [20]. Briefly, Fit takes advantage of the pattern that deleterious mutations are less likely to segregate in the sample of individuals and those that are segregating are more likely to be at lower frequency compared to neutral mutations. Assuming the MIS-SFS is under selection while the SYN-SFS is neutral, the best DFE parameters should fit the differences between MIS-SFS and SYN-SFS in each minor allele frequency bin. We first computed and stored the expected MIS-SFS for a range of population-scaled selection coefficients (−10,000 −10⁻⁵ and 10⁻⁵ 100) in each species, given their best fit demographic scenarios, using the DFE1D_refspectra module in varDFE. Note that varDFE assumes semi-dominance of the mutations involved (h = 0.5), and s refers to the selection coefficient against heterozygotes.

To infer the deleterious-only DFE, we assumed that the DFE follows a gamma distribution, given its strong theoretical and empirical support [7,16,19]. We parameterized the DFE for each species as , where (; shape) and (; scale) are parameters to be inferred. For each species, 100 replicates were run from a permuted starting parameter of = 0.2 and = 4,000 using the DFE1D_inferenceFIM module in varDFE [7]. The best-fit parameters with the maximum Poisson log-likelihood were chosen, and the corresponding Akaike information criterion (AIC) was calculated. To quantify uncertainty of the parameter estimates, we Poisson resampled the MIS-SFS 200 times, and repeated the DFE inference to obtain the 95% confidence interval [19,20]. Recall that we optimized the gamma distribution parameters for the population-scaled selection coefficient . To obtain the distribution of s, we unscaled the gamma distribution by , therefore, , where and had been inferred from demography estimation. The mean selection coefficient |s|_mean for each species was calculated given the gamma distribution’s property that . Similarly, . The median selection coefficient |s|_median was found as the value of s where the cumulative distribution function equals 0.5 using the pgamma function in R (v.4.4.2 [76]). To compute the proportion of mutations with different values of s or , we found the cumulative probability in the given range using the pgamma function.

To estimate the proportion of effectively neutral mutations, we added a point mass at neutrality () to the original gamma distribution. The DFE function can be written as .

(1)

The distribution of s is , where . The mean selection coefficient |s|_mean for each species is . The median selection coefficient |s|_median was found where the cumulative distribution function equals 0.5.

To examine the robustness of DFE inference, we conducted additional analyses assessing the effects of DFE functional forms, demographic model misspecification, singleton-masking treatments, mutation rate variation, differences in nonsynonymous to synonymous mutation rate ratio (), and GC-biased gene conversion (detailed in Text D in S1 Appendix). To account for uncertainty in the functional form of the DFE, we performed DFE inference assuming three different functional forms (gamma, neugamma, and lognormal) under the best-fit demographic models and singleton masking treatments. For the lognormal distribution, we parameterized the DFE for each species as . The distribution of s was unscaled as , where . The median selective effects |s|_median for each species was calculated using the lognormal distribution’s property that .

Phylogenetic signal detection

To evaluate whether the DFE co-varies with phylogeny, we calculated Pagel’s for the log10 transformed mean selection coefficient (), the log10 transformed median selection coefficient (), and the shape parameter () in the gamma-distributed DFE, life history traits and genomic attributes using the phytools package (v.2.4.4 [77]). The underlying phylogenetic structure and divergence times for 11 (sub)species were compiled from the TimeTree database (http://timetree.org/). To test if the candidate life history traits and genomic attributes (X) are correlated with the observed DFE parameters, and correct for phylogenetic dependency when appropriate, we utilized the method in phytools. The regression was specified as , , or to normalize the data, except when X was genome size, in which case X was not log-transformed. The estimated in does not equate to the Pagel’s estimated individually for X or Y, but represents whether the correlation of X and Y can be accounted for by phylogenetic dependence in residual error. Candidate explanatory variables (X) include synonymous-site-based and literature-based long-term population size ( and ) in diploids, body mass in grams, generation time in years per generation, age at maturity in days, maximum longevity in years, genome size in Mb, recombination rate in recombination events per bp per generation, and mutation rate in mutations per bp per generation. P-values were adjusted using a Bonferroni correction across the nine explanatory variables considered. To control for outliers or account for the long branches between mammals and insects, we repeated the Pagel’s and analysis excluding the arctic wolves data point, or only within mammals (S15 Data, S16 Data).

To formally test for variation in the DFE across species, we used a grid-search approach (DFE1D_gridsearch in varDFE) to compare each species’ DFE estimates to a null model where the DFE was constrained to be the same across species by likelihood-ratio tests (LRT). 2000 grid points were evenly spaced in biologically meaningful ranges for shape ( = 0.1 to 0.5) and scale ( = 0.0001 to 0.5) parameters. Using integrate in   we obtained the expected SFS for each pair and calculated the Poisson log-likelihood relative to the empirical MIS-SFS in each species. For the null model, we summed log-likelihoods across species to find the MLE. To test whether the shape () and scale () are different in any x number of species, the LRT was constructed as , where LL₀ is the log-likelihood for the null model, with two DFE parameters inferred (; params = 2) given the inferred demographic parameters for each species (), and LL₁ is the log-likelihood for the alternate model, where each species is allowed to have its own DFE parameters (; params = 2x). Across the 11 (sub)species tested, asymptotically, should follow a distribution with . However, as is a composite likelihood, its empirical distribution under the null model of constrained DFE parameters is broader than the distribution [7]. Although small P-values (e.g., P < 10⁻¹⁶) should not be interpreted literally, simulations showed that did not exceed 34 under the null model in a two-species comparison, suggesting that values of like those seen here would be unlikely under the null hypothesis [7].

We also calculated the pairwise LRT statistics in all possible pairs of species (x = 2). For example, the LRT statistics comparing humans (population 1) and Drosophila (population 2) DFE can be written as follows:

(2)

Hierarchical clustering for was performed by the hclust function in R. To test the significance of the variation of across species, we repeated the LRT by evenly spacing 2000 grid points in shape ( = 0.1 to 0.5) and population-scaled scale ( = 100 to ) parameters. We also calculated Pagel’s for the log10 transformed population-level mean selection coefficient (). We repeated analysis with as the response variable, following methods described above.

Fitting the FGM-derived DFE

To investigate the genetic mechanisms underlying DFE variation, we implemented a functional form of the DFE directly derived from the Fisher’s geometric model [10]. The FGM proposes that fitness (w) can be described as a Gaussian function of the distance from the optimum (z), , in an n-dimension phenotypic space. Random mutations affect a subset of m phenotypes (), with fitness effect size r following a zero-mean Gaussian distribution with scale . The DFE for a population at mutation-selection-drift equilibrium can be described using the mutation pleiotropy (m), scale of mutation effects () and long-term population size (). It was derived originally by Lourenco and colleagues [10], therefore we refer to it as . Here is the gamma function and K(.) is the modified Bessel function of the second kind [10].

(3)

While complex in its expression, this FGM-derived distribution is functionally equivalent to the previously described gamma and lognormal distributions for generating mutational effects. Therefore, we estimated its parameters (m, , and ) following methods described in the DFE inference section, through maximum-likelihood optimization by fitting the observed differences between MIS-SFS and SYN-SFS across minor allele frequency bins. To compute the proportion of mutations with different values of s, we found the cumulative probability for s in given ranges using the scipy.integrate.quad function (v.1.8.0) in Python. The proportion of beneficial mutations was obtained by calculating the cumulative probability for s > 0. We performed two regressions using methods described above: first, inferred in this FGM-derived DFE and the obtained in demographic inference; second, the proportion of beneficial mutations and . To test for phylogenetic signal in the inferred m and parameters, we calculated Pagel’s (S15 Data).

Simulating a human DFE with the Drosophila demography

To investigate whether Fit can accurately detect strongly deleterious mutations in species with large population sizes, we performed a simulation study using Drosophila demography and the human DFE. We used the stdpopsim library to obtain coding sequence (CDS) annotations for Drosophila melanogaster chromosomes 2L, 2R, 3L, and 3R from the FlyBase_BDGP6.32.51_CDS [78]. CDS regions were labeled as “exon", and non-coding regions as “bkgd". To simulate realistic genomic regions, we randomly sampled 108 1-Mb segments from these chromosomes with probabilities proportional to chromosome lengths. Annotations for each segment were subset and converted to SLiM’s coordinate format. Next, we implemented a two-epoch Drosophila demographic model in which the ancestral population size was set to 2,766,461, and the population expanded to 7,482,090 over 507,689 generations. We used a uniform recombination rate of which is an average of recombination rates for Drosophila chromosomes from the stdpopsim catalog. The mutation rate was set to per bp per generation. Deleterious mutations followed a human-like gamma DFE (mean = –0.0276, shape = 0.189) with a nonsynonymous-to-synonymous ratio of 2.85. Background (noncoding) regions were assigned neutral mutations only. All mutations were additive (semi-dominant). From each simulation run on the randomly sampled 1-Mb Drosophila chromosome segments, we output SFS using a sample size of 100 or a sample size of 10. To reduce computational cost, we applied a scaling factor Q = 100, adjusting population size, mutation rate, recombination rate, and selection coefficients accordingly. One out of 108 simulations did not complete due to time constraints on the cluster and was discarded. To match the total length of chromosomes 2L, 2R, 3L, and 3R in Drosophila, we sampled the SFS outputs from 107 successful simulations with replacement, and summed 107 SFS outputs to generate one bootstrapped SFS that mimics a resampled version of the simulated Drosophila genome ( 107 Mb). To assess inference uncertainty, we repeated the sampling 50 times, producing 50 bootstrapped SFSs for each sample size (2n = 100 or 2n = 10). We then inferred the DFE on the 50 bootstrapped SFSs as done for the empirical data. We assumed either a gamma-distributed or a lognormal-distributed DFE during inference. To evaluate the influence of selection at linked sites on synonymous mutations, we also inferred demographic parameters from the “bkgd" (non-coding) regions to compare with the demographic parameters inferred from synonymous mutations.