A population’s ability to adapt is determined by its levels of additive genetic variance (V_A), and while it is agreed that most organisms have genetic variation for most traits, the extent to which it varies between species is poorly characterized. Here, we investigate this question by compiling 3,209 and 1,852 estimates of heritability and evolvability (the additive genetic variance divided by the square of the mean), respectively, for a variety of traits from 220 and 172 multicellular eukaryotic species. Using phylogenetic generalized linear mixed models, we find substantial and highly significant interspecific variation in evolvability. Much of the variation is explained by phylogenetic relatedness, with plants in our data having substantially higher evolvability than animals. While heritability also varies between species, the differences are more subtle, and plants are not exceptional. We investigate whether the variation in evolvability and heritability between species is due to variation in the mutation rate, effective population size, genome size, ploidy, and recombination rate, but find little evidence of any factor being important. However, the confidence intervals are large suggesting that we have little power to detect any associations between these factors and our estimates of V_A.

Literature search and QGV compilation

Heritability and evolvability estimates were compiled from the peer-reviewed literature indexed in the Web of Science for the period between 1992 and 2022. Searches were conducted using “Heritability” or “Evolvability” as topic keywords, with the aim of capturing estimates from a diverse range of species and from natural populations. Therefore, journals focusing on domesticated, crop species, or human populations were excluded. The search was further filtered by only selecting journals with over 200 papers meeting the search criteria and then scanning the first 50 articles by titles when sorted by relevance. This resulted in the choice of four journals: “Journal of Evolutionary Biology,” “Evolution,” “Heredity,” and “Proceedings of the Royal Society B.” When multiple estimates were available for the same trait within a population, we selected a single value based on a standardized evaluation of reliability. Preference was given to estimates derived from larger sample sizes and model structures minimizing potential biases. When estimates were otherwise comparable in these respects, we applied a hierarchy based on the relatedness structure used to estimate additive genetic variation. Animal models were preferred, as they utilize information from all known relationships and are generally the least biased and most precise [28]. When animal model estimates were unavailable, we selected those from parent-offspring regressions, which are less affected by dominance or common environmental effects than sibling-based designs ([29]).

The resulting dataset was then merged with the datasets of Mittell and colleagues [8] and Young and Postma [9]. We additionally conducted a targeted search for estimates from species that would increase overlap with our predictor variables. We estimate that we have captured approximately 19% of all relevant heritability and evolvability estimates in our data (S1 Appendix). All recorded estimates of evolvability and heritability were numerically checked where possible and recalculated if sufficient raw data were available. Each estimate was also evaluated for suitability in terms of the trait scale and transformation, as evolvability is only meaningful for traits measured on ratio or log-interval scales, where proportional differences are interpretable [ 4,6]. Heritability estimates were likewise checked for compatibility with the underlying trait scale, as h² is interpretable for interval and ratio-scale traits, but not for ordinal or nominal traits unless modeled and interpreted appropriately ([2]; [29]). If an evolvability estimate was not given, this was calculated if the necessary data were provided. We removed estimates where artificial selection had been performed on the focal trait, along with estimates from inbred lines, since inbreeding inflates V_A. However, estimates were retained from species that are naturally clonal or selfing. We log₁₀-transformed estimates of I_A because the estimates are not normally distributed. Consequently, 146 negative and zero evolvability values were removed, as these are undefined under log transformation. Removing these undefined values could, in principle, bias the results by truncating the lower tail of the I_A distribution, so we also ran a censored model that retained these. The results were qualitatively similar, but the censored model fit the data substantially worse, so we report only the log-transformed analysis

Because I_A was log-transformed while h² was modeled on its natural [0,1] scale, the two response variables are not directly comparable across analyses. Results are therefore interpreted separately. The bounded nature of h² additionally reduces total variance available to be explained, likely reducing statistical power relative to the I_A analyses.

For cases in which both I_A and h² have been estimated, we can calculate the residual variance, the sum of the non-additive and environmental variances, which we normalized by dividing it by the square of the mean: I_R = I_A(1 − h²)/h² [4].

Interspecific variation models

To partition the variation in recorded estimates of V_A, a phylogenetic generalized linear mixed model (PGLMM) was fitted using Markov chain Monte Carlo (MCMC) implemented in MCMCglmm [30] in R [31]. Weakly informative priors were used throughout. Random-effect variances were modeled using scaled F_1,1 priors (scale √1,000), while the residual variance followed an inverse-gamma prior with shape and scale set to 0.001. Fixed effects were assigned normal priors with a mean 0 and variance 10⁸. The MCMC chains were run for 780,000 iterations, with a burn-in of 180,000 and thinning every 250 iterations. Model convergence was assessed through visual inspection of chain mixing. Additionally, we used posterior predictive checks to visually assess model fit. We included a number of fixed and random effects in our model which aim to capture the likely major sources of variation in V_A (Table 1). Critical to our analysis are two species-level random effects that capture the variation in V_A among species, one associated with phylogenetic inertia (“Phylogeny” in Table 1) and one which is independent of this (“Non-Phylogenetic” in Table 1). The phylogenetic term takes into account that closely related species may covary in their levels of V_A; the tree was estimated using Time Tree [32].

The model was fitted to evolvability and heritability estimates separately. Model predictions for the log evolvability and heritability of each species was calculated by extracting the intercept, phylogeny, and species estimates for the ith species at the jth sampling of the chain

(1)

where X is either log (I_A) or h². These were then averaged across samples of the chain to yield the posterior mean evolvability and heritability for a linear morphological trait estimated with an animal model in a wild population (these are the reference categories in our analysis), along with the 95% credible intervals.

To test whether there is significant variation between specific phylogenetic groups represented in our dataset, we averaged the species estimates within each taxonomic subgroup for each sampling of the MCMC. We then calculated pairwise differences between these group means (e.g., plants versus animals, arthropods versus mammals, etc.). From these contrasts, we derived the posterior distribution of the difference in means, quantifying the variation among the phylogenetic groups sampled in our data.

To further investigate whether phylogenetic groups differ in their average evolvability or heritability, we estimated the evolvabilities and heritabilities for the internal nodes (n) of our phylogeny, providing posterior estimates for ancestral nodes.

(2)

These node-level estimates were then used to calculate pairwise differences among major clades, following the same procedure as above. Finally, we extended our original interspecific variation models by including taxonomic groups as fixed effects, allowing us to formally test whether differences in log (I_A) or h² between groups deviate more than expected under the null model of Brownian evolution (the phylogenetic term) with added independent and identically distributed species effects (the non-phylogenetic term).

To quantify the amount of variation between our species in terms of evolvability and heritability we took the estimated variance associated with the phylogenetic (V_phylogenetic) and non-phylogenetic (V_{non-phylogenetic}) species terms and then assuming the terms are normally distributed, we derived the quartiles of the distribution of the combined distribution. We quantify the variation in terms of the ratio of the upper and lower quartiles. This allows us to state, for example, that the top 25% of species have an evolvability that is at least x-fold greater than the bottom 25% of species. The 95% credible intervals on these ratios were derived from the samples of the chain.

Correlation analysis

We investigated whether our estimates of I_A and h² could be predicted by a number of factors: the point mutation rate, the point mutation rate multiplied by the number of protein coding sites in the genome, genome size, number of genes, ploidy levels, the effective population size, the rate of recombination per nucleotide, the mode of reproduction, body mass, body length, longevity, and generation time. Body length was defined as the longest linear dimension; we restricted this analysis to animals since body length is hard to assess in plants because of the root system. Data on the predictors were taken from various databases, including Amniote [33], Pantheria [34], Animal Diversity Web [35] as well as some meta-analytic studies such as Wang and Obbard [36], Lynch and colleagues [15], Stapley and colleagues [37], Lewin and Eyre-Walker [38], and the broader literature (see S1 Table for details). All predictors were log₁₀-transformed prior to analysis.

To assess the predictive power of each of our predictors, we used the same fixed and random effect structure, prior, and iteration settings as for the interspecific variation models. Since we did not have estimates of each predictor for every species, we ran the analysis for each predictor separately with heritability and evolvability, respectively, as the response variable. The amount of between-species variance explained by each predictor was quantified using a signed R² [8].

(3)

where V(x) is the variance in the focal predictor and β is the regression coefficient. The product β²V(x) represents the between-species variation in the respective measure of V_A that is attributable to variation in the predictor, with the sign of R² taken from the sign of β. The denominator quantifies the total between-species variation, including phylogenetic and species-level components.

pMCMC values were used to assess the degree to which coefficients overlapped zero [39], and for multi-level factors, we used the logic of the Wald test to perform an omnibus test of the factor’s overall effect. The Wald statistic was calculated from the posterior mean vector and covariance matrix [8].

Evolvability versus heritability

To investigate the relationship between heritability and evolvability, we fitted a bivariate PGLMM with the phylogenetic and non-phylogenetic terms as random effects (Table 1) in which h² and I_A were modeled jointly. The phylogenetic effects, non-phylogenetic effects, and residuals were each allowed to covary across response variables and priors were placed on the resulting covariance matrices such that the marginal prior for the variances was identical to that in the univariate models. All MCMC settings were identical to those of the univariate models.

There is an expectation for V_A to be positively correlated with the non-additive genetic variance and the environmental variance [4]. The sum of these two components is the residual variance (V_R), which we scale by dividing by the square of the mean to yield I_R. To test whether I_A and h² are correlated to I_R we ran bivariate PGLMMs under the same conditions as the h² and I_A model.

Fig 1. The predicted evolvability from our model for a linear, morphological trait estimated with an animal model from a wild population for each species, plotted across the phylogeny with their 95% credible intervals.

Organism silhouettes were obtained from PhyloPic (https://www.phylopic.org) and are in the public domain (CC0 1.0 Universal Public Domain Dedication or Public Domain Mark 1.0). The data underlying this figure can be found in the data repository accompanying this paper under the DOI: https://doi.org/10.5281/zenodo.20021814.

https://doi.org/10.1371/journal.pbio.3003819.g001

Fig 2. The difference in mean evolvability between groups of animals and plants.

Each point represents the posterior mean difference between group means, back transformed to the original data scale. Horizontal bars indicate the 95% credible intervals. The dashed line at 1 indicates no difference between groups. Differences greater than 1 indicate that the first node (first-mentioned group in the y-axis) has a x-fold higher predicted evolvability than the second node. All comparisons are plotted so that the means are positive. The data underlying this figure can be found in the data repository accompanying this paper under the DOI: https://doi.org/10.5281/zenodo.20021814.

https://doi.org/10.1371/journal.pbio.3003819.g002

Evolvability (I_A)

To investigate whether I_A varies between species, we fitted a phylogenetic GLMM with two species terms—one capturing differences due to phylogenetic relatedness and another capturing species-specific deviations independent of phylogeny; the total effect due to species is the sum of these terms. The fit of the model to the data was adequate, although the model predicted a deficit of small values and an excess of large values compared to that observed (Fig A in S2 Appendix). We find that I_A varies significantly between species (Fig 1), with 39% [95% credible intervals: 16%–58%] of the total variance in the random effects (including the residual variance) in the model attributable to interspecific differences (Table 2). The magnitude of this variation is substantial: the top 25% of species exhibit 9.3-fold [3.5, 23] greater evolvability than the bottom 25% (note that we have transformed evolvability back to the natural, rather than log scales, when discussing fold-differences in evolvability). Most of the interspecific variation reflects phylogenetic structure, with phylogeny accounting for 38% [15%, 58%] (Table 2) of the random effect variance.

A substantial part of this phylogenetic effect appears to be a difference between plants and animals, with plants having 4.7-fold [2.7, 8.3] higher I_A than animals. Plants also have significantly higher evolvability than each of the individual animal subgroups (Fig 2). This difference between plants and animals, and individual animal groups is also evident in the estimated ancestral values for each group; the ancestor to our sample of plants is estimated to have had 4.3-fold [1.1, 18] higher evolvability than the ancestor of the sample of animals (Fig B in S2 Appendix). However, the difference between plants and animals is consistent with the phylogenetic model. When group identity (plant versus animal) is included as a fixed effect alongside the phylogenetic and non-phylogenetic species effects in the model, the fixed effect coefficient was not significantly different from zero (mean effect = 0.73 [−1.46, 2.52]; Table B in S2 Appendix); i.e., plants are no more divergent from animals than we might expect given the time that has elapsed since the groups diverged and the rate at which evolvability changes.

However, the variation in evolvability between species is not simply due to the difference between plants and animals. If we restrict our interspecific variation model to animals only, we still find significant variation between species, with 24% [7%, 40%] of the variance attributable to species differences (Table 2). We also find significant variation in I_A between species within plants, mammals, arthropods, and fish individually (Table 2). In each of these cases, however, the phylogenetic and non-phylogenetic effects are not individually significant. This indicates that while there is a significant species effect, the data lack sufficient power to partition the species effect between the phylogenetic and non-phylogenetic effects.

Given the complexity of the full interspecific variation model, which includes many fixed and random effects, and an assumption that log (I_A) is linearly related to the predictors in our model, we performed a simplified analysis using only linear morphological traits estimated with an animal model. This results in 493 estimates from 63 species. Results were qualitatively unchanged: interspecific variation accounted for 58% [32%, 81%], and the phylogeny for 56% [28%, 80%] of the random effect variance. Moreover, if we simply plot the species’ mean I_A values for this data the phylogenetic effect remains evident (Fig C in S2 Appendix).

Although, our principal focus is whether evolvability varies between species, a number of other fixed effects are significant in our model. We observe substantial differences in evolvability among different trait types (Omnibus test, P(>x²) < 0.001), with behavioral traits exhibiting substantially greater I_A than other traits. We also find a significant effect of dimension (P(>x²) < 0.001; Table C in S2 Appendix), where traits that scale quadratically (e.g., area) and cubically (e.g., mass) show 1.7 [0.83, 3.5] and 3.5 [2.2, 5.5] times greater evolvability than those that scale linearly (Fig D in S2 Appendix). Counts exhibit 3.4 [1.9–5.9] times greater evolvability than linear traits. The method of estimation also has a significant effect (P(>x²) = 0.012; Table C in S2 Appendix), with full-sib designs yielding higher estimates that are 4.0 [1.6, 9.8] times greater on average than animal models, and midparent-offspring regression yielding the lowest values that are 0.40 [0.17, 0.92] of the values from the animal model (Fig D in S2 Appendix). Estimates of I_A also depend on the environment (P(>x²) = 0.003) in which the estimation was made. However, no individual level differed significantly from the reference level, and direct comparisons between laboratory and wild estimates were not significant (Table C in S2 Appendix; Fig D in S2 Appendix). This indicated that the overall significance arises from moderate differences across multiple environment categories rather than a single strong pairwise contrast.

Although these fixed effects account for some of the observed variation in I_A, they primarily serve to control for differences in methodology and trait choice across publications, and do not explain the interspecific variation in I_A that we detect. To investigate what might be causing this variation, we conducted PGLMMs including a suite of predictors: nucleotide mutation rate, estimated total point mutation rate of the proteome, genome size, number of genes, ploidy level, recombination rate per site, mode of reproduction, and estimates of effective population size. In addition, we assessed associations between evolvability and life-history traits such as body mass, body length, longevity, and generation time.

Among all variables tested, only body length (which we have only assessed for animals) exhibited a statistically significant association with evolvability after applying a Bonferroni correction for multiple testing. Body length accounted for approximately 46% [4.9%, 94%] of the between-species variance in evolvability, and the model predicts that doubling body size increases I_A by 1.3-fold [1.1, 1.5]. No other predictors were significant predictors of evolvability after correction for multiple tests; however, the credible intervals are broad, indicating we have little power to detect associations (Table 3).

Heritability

When we run our model on heritability, we find that the model fits the data reasonably well except that there is a large excess of h² estimates at zero compared to what the model predicts (Fig E in S2 Appendix). This excess of zero estimates is almost certainly methodological, arising because many of the analyses return zero estimates of V_A when relatives exhibit negatively correlated trait values.

As with evolvability, we find evidence of significant variation between species—the proportion of the random effects variance explained by between-species differences is 24% [4.9%, 52%], though these proportions are not directly comparable because we have log-transformed I_A and the two statistics have different bounds. However, this variation cannot be clearly attributed to either phylogenetic or non-phylogenetic components when considered separately (Table 4). The lack of a clear phylogenetic signal is evident when the predicted h² from our model for each species are plotted on a phylogeny (Fig 3). We also find significant variation in h² within plants and animals, respectively; however, the smaller taxonomic subgroups as mammals and fish only harbor marginally significant interspecific variation, while the remaining groups tested do not show significant differences (Table 4). Overall, we find that h² varies quite substantially between species; the ratio of the upper and lower quartiles is 1.8-fold [1.3, 2.6].

All the fixed effects (Table 1) in our interspecific variation model (i.e., the model without any predictors such as the nucleotide mutation rate) are significant, with the exception of the number of fixed effects (p_MCMC = 0.118; Table D in S2 Appendix). Consistent with previous research, h² varies significantly between trait types (Omnibus test, P(>x²) < 0.001) with morphological traits having the highest h² and fitness traits the lowest (the average difference in heritability between morphological and fitness traits is 0.15 [0.013, 0.28] (Table D in S2 Appendix; Fig F in S2 Appendix). We also find that the method of estimation has a significant effect (p_MCMC < 0.001, Table D in S2 Appendix; Fig F in S2 Appendix) with full-sib designs yielding estimates of heritability that are 0.21 [0.11, 0.30] higher than animal models, which yield similar estimates to other methods. We also find a weak but significant effect of dimension (P(>x²) = 0.012; Table D in S2 Appendix; Fig F in S2 Appendix) and the environment (P(>x²) = 0.002, Table D in S2 Appendix; Fig F in S2 Appendix). The effect of the environment seems to largely reflect higher heritability estimates from laboratory studies compared to those from natural populations (p_MCMC = 0.003, Table D in S2 Appendix; Fig F in S2 Appendix).

Although we observe significant variation in narrow-sense heritability (h²) across species, most predictors tested do not exhibit statistically significant associations after Bonferroni correction (Table 5). The exception is ploidy level (P(>x²) = 0.004), with polyploids showing lower heritability than diploids.

Evolvability versus heritability

Although estimates of I_A and h² are significantly correlated across individual estimates, the correlation is weak (r = 0.24 [0.20, 0.29]; Fig G in S2 Appendix), and there is no significant correlation between I_A and h² across species (r = 0.13 [−0.13, 0.38]; Fig G in S2 Appendix). There are several reasons why I_A and h² might be weakly correlated. In particular, it has previously been shown that the residual variance—the sum of the non-additive genetic and environmental variances—is positively correlated to V_A across individual estimates [4]. Consistent with this, we find a strong positive correlation between residual variance (mean-scaled, I_R) and I_A across individual estimates (r = 0.73 [0.70, 0.75]; Fig H in S2 Appendix) and species (r = 0.69 [0.47, 0.89]; Fig H in S2 Appendix). Similarly, we find a significantly negative relationship between h² and I_R across individual estimates (r = −0.33 [−0.36, −0.29]; Fig I in S2 Appendix) and species (r = −0.51 [−0.68, −0.28]; Fig I in S2 Appendix).

In addition to the correlations, we also analyzed the interspecific variation in I_R itself using an interspecific variation model and find that the variation in I_R follows a similar pattern to I_A; there is a strong phylogenetic component (interspecific variation = 29% [8.3%, 53%]) (Table E in S2 Appendix) with plants having substantially more residual variance than other phylogenetic groups (Fig J in S2 Appendix).

Why does evolvability vary between species?

Although we find substantial variation in evolvability, we do not find that it is predicted by any variable that we have tested except body length. Particularly surprising is the lack of a correlation between evolvability and the mutation rate since most models of the maintenance of the additive genetic variance predict that it should depend on the mutational variance [13]. There are several reasons why we might not have detected an effect. First, available proxies for mutational input are relatively crude. We assessed whether evolvability correlates with the nucleotide mutation rate alone, as well as with the mutation rate scaled by proteome size; for these variables, we have few estimates. We also considered the number of protein-coding genes (without multiplying this by the mutation rate, which is only known for a limited number of species). However, the number of genes likely represents a limited proxy for the genomic regions in which mutations can generate phenotypic effects. It does not account for the amount of regulatory DNA, which is likely to vary among species and indeed, most GWAS hits fall in non-coding regions [43,44], meaning that proxies based on proteome size or gene number likely miss much of the mutational target relevant for quantitative traits. Moreover, under most models of quantitative genetic variation maintenance, mutational variance is influenced not only by mutation rate but also by the distribution of effect sizes, which is also known to vary across species [45,46]. Unfortunately, we have very few estimates for the distribution of effect sizes. However, possibly the biggest omission in our estimate of the mutation rate are structural variants (SVs) including transposable element (TE) insertions. There is evidence from Drosophila [47,48] and yeast [49] that SVs can contribute substantially to additive genetic variation. Unfortunately, estimates for the rate at which SVs occur have been estimated for very few species. It is also possible that I_A is not correlated to our estimates of the mutational variance because it does not differ substantially between species, in which case it could not account for the variation in I_A we observe.

We also find no correlation between evolvability and the effective population size. This might be due to a lack of power; both our estimates of N_e and species mean I_A are subject to reasonable levels of error. However, in the light of current theory, a strong I_A ~ N_e relationship is not generally expected. Classical models show that V_A scales with N_e under neutrality [50] and that under mutation-selection-drift balance, V_A becomes sensitive to N_e only when drift is strong relative to selection [19,51]. Under stabilizing selection, most mutational input is deleterious and removed by selection before drift can act, making V_A largely insensitive to N_e except at very small N_e [20,52]. Moreover, the effective neutral mutational input (N_eu) does not necessarily increase with N_e: as N_e increases, the fraction of mutations that are effectively neutral declines, reducing the effective mutational heritability [13]. Recent work reframing these results in terms of nucleotide diversity similarly emphasizes that strong V_A ~ N_e scaling is expected only when traits behave neutrally or when N_e is sufficiently small that drift dominates selection [18]. The absence of a detectable I_A ~ N_e correlation in our analysis is therefore consistent with theoretical expectations.

The one significant correlation we detect is a positive relationship between evolvability and body size in animals (measured as the maximum linear dimension), which remains statistically significant after correcting for multiple tests (p = 0.003, multiple‐test–corrected significance threshold = 0.004). This might arise because larger-bodied species tend to have longer generation times, and species with long generation times have higher mutation rates per generation [36,53] although we tested for a correlation between evolvability and the mutation rate, we had few estimates in our analysis and hence low power. There is also some evidence that larger-bodied species have larger ranges [54,55] and might therefore be exposed to different environments, leading to maintenance of locally adaptive genetic variation. Finally, the observed positive correlation may also be influenced by trait scaling. If the additive genetic variance does not increase linearly with the square of the trait mean, I_A will tend to increase with larger trait values. While this could in principle be tested by excluding size traits, the majority of our dataset consists of such traits, and removing them would result in a severe loss of power, preventing reliable inference.

Although this study explored a range of factors that could explain interspecific variation in quantitative genetic parameters, there are additional important forces we have not explored. One such factor is selection, which is expected to affect levels of additive genetic variance (though not always under directional selection—see [56]). However, the extent to which selection varies among species is largely unknown, because meta-analyses of selection gradients rarely account for species identity or phylogenetic relatedness [57]. Finally, population structure—such as sub-division and migration—could influence V_A, but we have little information on these for the species we have analyzed.

Why is evolvability higher in plants?

We find that plants have significantly higher evolvability than animals. Although our analysis suggests that the difference between plants and animals is no more than we might expect under a simple phylogenetic model of trait evolution, the difference is substantial, with plants estimated to have ~4.7-fold [2.7, 8.3] more evolvability than animals.

Unfortunately, understanding why plants and animals differ is challenging because there are many fundamental differences between plants and animals, and both groups are monophyletic; as a result, we have limited statistical power to test hypotheses. For example, plants generally have many more genes than animals (in our dataset, the mean number of genes is 2.3-fold greater than in animals), but because there is relatively little variation in gene number within plants and animals, we have few independent contrasts to detect an association. We find no evidence that this pattern might be attributable to differences in the mutation rate, effective population size, mating system, genome size, or ploidy. However, there are two conspicuous differences between plants and animals that might contribute to the difference in evolvability between the groups. First, plants have an additional genome, in the chloroplast. While both plastid and mitochondrial genomes are small relative to the nuclear genome, studies have shown that genetic variation in the mitochondrial genome can significantly contribute to phenotypic variation in animals [58,59]. A meta-analysis of studies in which mtDNA was introgressed between nuclear backgrounds in animals suggested the effects can be substantial, with nearly 50% of cases changing the mean up or down by >10% [60] although this is likely to be an overestimate due to sampling error [61]. Similarly, the limited number of studies that have considered the effect of the organellar DNA on phenotypic variation in plants have found substantial effects [62–65].

The second major difference between plants and animals is the potential for transgenerational epigenetic inheritance (TEI) [66]. There are examples of TEI in both plants and animals [67], but no comparative study of how much epigenetic variation is inherited in the two groups of organisms. However, there are some reasons for believing that inherited epigenetic variation might be greater in plants than animals. There are several pathways by which TEI can be propagated [67], but the simplest is through DNA methylation, and most methylation marks are not removed during gametogenesis in plants, whereas extensive remodeling occurs in animals [66,67]. Experiments using epigenetic recombinant inbred lines (epiRILs) of Arabidopsis thaliana suggest that heritable differences in methylation can generate substantial phenotypic variation [68,69] comparable in magnitude to differences observed among natural accessions. Broad sense heritability amongst epiRILs can be as high as 90% for some traits, such as flowering time, but essentially absent for others [64,69]. However, an epiRIL experiment gives an upper estimate for the variation induced by heritable methylation since it involves crossing methylated and largely unmethylated parental strains, hence generating more diversity in methylation than would be seen in a natural population. The fact that largely unmethylated strains are viable suggests that DNA methylation might not play an essential or widespread role in generating large amounts of additive genetic variance, even if it can influence particular phenotypes. Although natural populations do contain epigenetic variation in methylation [70,71], much of this appears to be genetically determined [71]. It therefore seems unlikely that epigenetic variation is why plants have higher evolvability than animals.

Evolvability—Other effects

Hansen and colleagues [4] reported that the evolvability of cubic traits was roughly nine times that of linear traits, implying strong genetic correlations among trait dimensions (e.g., between height and width). In contrast, our analysis finds that the ratio of linear, quadratic, and cubic traits is closer to 1:2:3, consistent with variation in different dimensions being largely independent. The difference between our study and that of Hansen and colleagues [4] might be due to our larger sample size and the fact that we have modeled many effects they did not, including the method of estimation, trait type and the relatedness of the species in our sample; if we simply take the median estimates of the observed I_A values we find a pattern closer to that originally reported by Hansen and colleagues [4]. An additional difference is that our analyses are based on log-transformed I_A values, whereas Hansen and colleagues [23] used raw I_A, which might also contribute to the differences we observe. Biologically, one might have expected an intermediate scaling between the different dimensions, given that our dataset is extremely heterogeneous spanning diverse species and traits that likely vary in their degree of dimensional integration. However, a more robust test of these scaling relationships will require a dataset with a greater representation of higher-dimensional traits.

We also show, for the first time, that count traits exhibit significantly higher evolvability (I_A) and lower heritability (h²) compared to linear traits. While Hansen and colleagues [4] did consider counts they did so only for size traits. The relatively high I_A of count traits may reflect a combination of biological and statistical factors. Biologically, count traits often arise from modular or discrete developmental processes, where variation among semi-independent units (e.g., number of leaves or offspring) can accumulate additively. Methodologically, the scaling of variance by the mean may contribute to elevated I_A values when counts are low or distributions are skewed, as we might expect under a Poisson distribution, overdispersion or zero-inflated distributions, although not all count traits conform to such distributions.

Interestingly, when we examine evolvability across different trait types, we find no significant difference between life-history and morphological traits, contrary to previous analyses [4,9]. Our analysis, however, varies from previous studies in several ways. Firstly, we separate fitness and life-history, yet neither is significantly different from morphological traits in I_A. We also include a random effect grouping traits that, although labeled differently across publications, measure the same underlying property (e.g., body size). Lastly, we explicitly account for trait dimensions across all trait types. If we remove trait dimension from our analysis, we find that life-history traits have higher evolvability than morphology traits. This result is particularly interesting when considering the distribution of trait dimensions across the different trait types. Morphological traits are predominantly linear (65%), with smaller proportions being quadratic (5.3%) or cubic (21%). In contrast, life-history traits are mostly durations (46%) and counts (51%). These distributions suggest that apparent differences between life-history and morphological traits may reflect underlying differences in trait dimensionality rather than intrinsic differences in evolvability.

Hansen and colleagues [4] demonstrated that the residual variance (the sum of the environmental and non-additive genetic variance) is positively correlated with the additive genetic variance across individual estimates of these quantities, a pattern we replicate. They argue that the relationship exists because increases in allele frequencies are expected to increase the additive and epistatic variances [72]. Furthermore, the additive, epistatic, and dominance variances all depend on the number of loci underlying a trait. Houle [6] has also argued that the environmental variance should increase with the additive genetic variance across trait types because complex traits are probably more sensitive to both genetic and environmental perturbation.

Extending these findings, we find a strong and significant relationship between I_R and I_A across species, demonstrating that the positive relationship observed at the level of individual estimates also holds at the species level. This is consistent with our observation that interspecific variation in I_R follows the same trend as I_A—phylogenetic signal with higher estimates in plants. This suggests that factors linking the additive and residual variances operate not only within species but also across species.

Heritability

We find significant variation in heritability across species, with the upper and lower quartiles varying by nearly 2-fold. This variation could reflect differences in additive, non-additive, or environmental variance. Across species, we find no evidence of a relationship between I_A and h², suggesting that the interspecific variation in h² is unlikely to be driven primarily by differences in the additive genetic variance. In contrast, the mean-scaled residual variance (I_R) varies significantly between species and is strongly positively correlated with I_A and negatively correlated with h² across both estimates and species. This pattern suggests, as proposed by Young and Postma [9], that the observed between-species variation in h² is likely influenced more by differences in residual variance than by differences in the additive genetic variance.

As with evolvability, we find only one variable that explains some of the variation in h², this is ploidy level. Polyploids exhibit significantly lower h² than diploids, which may reflect increased opportunities for epistasis and other non-additive genetic interactions in polyploids. Some of this epistatic variance may contribute to the additive genetic variance depending on the experimental design, but overall, it is consistent with the idea that higher levels of non-additive variation obscure the relationship between I_A and h² at the species level.

We find evidence that estimates of h² are generally higher in the lab than in natural populations contrary to a previous meta-analysis [73]. However, we have many more estimates and probably increased power because we have modeled many fixed effects that influence h². The difference between lab and wild estimates is not likely to be due to higher V_A in lab populations, as I_A is not significantly elevated in lab populations. A commonly invoked explanation is that laboratory conditions reduce environmental variance, thereby inflating heritability estimates. However, we find no significant difference between lab and natural populations in I_R. This suggests that the higher heritabilities observed in laboratory studies cannot be attributed simply to reduced residual variance under controlled environments and may reflect more subtle differences between laboratory and natural settings that affect how genetic and environmental variation contribute to phenotypic variance.

Consistent with previous literature, we find that morphological traits have higher heritability than other trait types [9,4]. In contrast, count traits exhibit lower h² despite having high I_A. This pattern mirrors our findings that count traits have significantly higher I_R (Table E in S2 Appendix). Because heritability is sensitive to the residual variance whereas I_A is not, elevated I_R in count traits can lead to reduced h² even when evolvability is high.

Summary

We find substantial interspecific variation in evolvability (I_A), with clear phylogenetic signal indicating that closely related species tend to have more similar levels of evolvability, with plants standing out with particularly high levels. While we also find significant interspecific variation in heritability (h²), and the absence of a correlation with I_A suggests that this variation might be due to variation in the residual variance. Although we find significant variation in the additive genetic variance between species, we are unable to explain this. Nevertheless, our results suggest that some species have a substantially greater capacity to evolve.