1. Introduction

The variation in genome size among animals is remarkable, spanning four orders of magnitude, from 0.02 Gb in the nematode Pratylenchus coffeae to 120 Gb in the marbled lungfish Protopterus aethiopicus (Gregory 2023). Understanding why such a huge variation occurs is a long-standing question in evolutionary biology. It is now established that genome size is not related to organismal complexity (C-value enigma) or the number of coding genes in eukaryotes. Rather, variation in DNA content depends on the amount of non-coding DNA such as Transposable Elements (TEs), introns and pseudogenes (Elliott and Gregory 2015; Kidwell 2002; Lynch et al. 2011). However, the evolutionary mechanisms leading certain lineages to inflate their genome size or to maintain streamlined genomes are still debated (Galtier 2024).

The various hypotheses that have been proposed can be divided into adaptive and non-adaptive. Adaptive theories such as the nucleoskeletal (Cavalier-Smith 1978) and the nucleotypic one (Gregory and Hebert 1999) consider genome size to be mainly composed of “indifferent” DNA (Graur et al. 2015), whose bulk is indirectly selected as a consequence of its effect on nuclear and cellular volume. Cellular phenotypes, in turn, are thought to influence the fitness of organisms by affecting traits such as cell division rate (Bennett and Riley 1997), metabolic rate (Vinogradov 1995; Vinogradov 1997) and developmental time and complexity (Jockusch 1997; Olmo et al. 1989; Gregory 2002). On the other hand, non-adaptive theories emphasize the importance of the neutral processes of mutation and genetic drift in determining genome size (Lynch and Conery 2003; Petrov 2001). In particular, the concepts originally proposed by Lynch and Conery were later on formalized within the framework of the Mutational Hazard Hypothesis (MHH) (Lynch 2007).

The MHH posits that tolerance to the accumulation of non-coding DNA depends on its mutational liability, which is minimized when the mutation rate is low and the effective population size is high (N_e, which is inversely proportional to the intensity of genetic drift). The fundamental assumption is that most of such extra DNA is mildly deleterious and its fate in the population depends on the interplay between selection and genetic drift: a newly emerged nearly-neutral allele with a given negative selective effect (s) should be effectively removed from the population when N_e is high (|N_es| >> 1), while it should have approximately the same chance of fixation of a neutral allele when N_e is low (|N_es| ∼ 1 or lower than 1) as drift will be the predominant force (Ohta 1992).

Because of the pervasiveness of TEs and their generally neutral or slightly deleterious effect (Arkhipova 2018), their dynamics in response to changing N_e are of particular interest in the context of the MHH. TE insertions are expected to drift to fixation as neutral alleles and enlarge genome size in organisms with low N_e, while the genomes of organisms with large N_e are expected to remain streamlined as emerging TEs should be efficiently removed by purifying selection (Lynch 2007). The MHH is highly popular because it is based on general principles of population genetics and proposes a unifying explanation for the evolution of complex traits of genome architecture without recourse to specific molecular processes. Indeed, the studies supporting the MHH are based on phylogenetically very diverse datasets as the goal of the theory is to explain broad patterns of complexity emergence and variation (Lynch and Conery 2003; Yi and Streelman 2005; Yi 2006). Nevertheless, other authors pointed out that the application of the MHH to such distantly related taxa could suffer from confounding factors intervening across organisms with very different biologies (Charlesworth and Barton 2004; Daubin and Moran 2004), thus raising the question whether N_e could explain genome size variation patterns at finer scales. Additionally, potential issues of robustness of the original dataset to phylogenetic control have been raised (Lynch 2011; Whitney et al. 2011; Whitney and Garland 2010).

Recent attempts have been made to assess the impact of increased genetic drift on the genomic TE content and genome size increase across closely related species with similar biological characteristics, employing both genetic diversity data and life-history traits (LHTs) as predictors of N_e. While some studies do not find any evidence supporting the role of N_e in genome size and TE content variation (Bast et al. 2016; Kapheim et al. 2015; Mackintosh et al. 2019; Roddy et al. 2021; Yang et al. 2024), others do (Chak et al. 2021; Cui et al. 2019; Lefébure et al. 2017; Mérel et al. 2021; Merel et al. 2024). Thus, a univocal conclusion can hardly be drawn from such studies. The MHH has thus been investigated at either very wide - from prokaryotes to multicellular eukaryotes - or narrow phylogenetic scales (i.e., inter-genera or inter-population). However, no study spanning across an exhaustive set of distantly related taxa and relying on a phylogenetic framework has yet been performed to our knowledge. Such an approach would allow a systematic test of the effect of N_e variation on long-term patterns of genome and TE expansion at a wider evolutionary scale, while controlling for the effect of phylogenetic inertia. Synonymous genetic diversity is commonly used to inform patterns of N_e (Lynch and Conery 2003; Romiguier et al. 2014; Buffalo 2021; Lynch et al. 2023). However, while its insights are limited to the age of current alleles (mostly less than 10N_e generations in diploid organisms), complex genomic features likely have much deeper origins. Comparative measures of divergence, like the genome-wide ratio of non-synonymous substitution rate to synonymous substitution rate (dN/dS), can quantify the level of genetic drift acting on protein-coding sequences since the last common ancestor. This approach accounts for processes occurring over a longer time scale than those responsible for genetic diversity. In fact, polymorphism might reflect relatively recent population size fluctuations (Daubin and Moran 2004) and might even diverge from indices of long-term N_e if the selection-drift equilibrium is not reestablished (Lefébure et al. 2017; Müller et al. 2022). We therefore adopt dN/dS as an index of long-term N_e, as it is more likely to reflect the evolutionary lapse during which the deep changes in genome size and TE content that we are investigating occurred (Whitney et al. 2011; Whitney and Garland 2010).

In this study, we took advantage of 3,214 public metazoan reference genomes and C-value records (i.e., the haploid DNA content of a nucleus) to estimate genome sizes. A subset of 807 species including birds, mammals, ray-finned fishes, insects and molluscs was selected to test the predictions of the MHH, especially through the relationship among the level of drift, genome size and genomic TE content. A phylogeny was computed with metazoan-conserved genes. N_e was accounted for by the dN/dS and by life-history traits (LHTs) when available; TE contents were estimated de novo from read data. Controlling for phylogenetic non-independence, we (1) assessed the contribution of TE quantity to genome size differences and (2) evaluated the efficacy of N_e proxies in explaining genome size and TE content variation, across the whole dataset and within specific clades.

2. Results

2.1. Selection of high-quality assemblies

The reference genomes of the 3,214 metazoan species were downloaded via the NCBI genome database (Supplemental Table S1). From these, the genomes with contig N50 ≥ 50 kb and available C-value record were employed for genome size estimation (see 2.2, 4.2). Based on the assembly quality (contig N50 ≥ 50 kb), the completeness of metazoan core genes (complete BUSCO orthologs ≥ 70%) (Manni et al. 2021) and the raw sequencing data availability, a dataset of 930 genomes was then retained for downstream analyses. For reliable estimation of substitution rates, this dataset was further downsized to 807 representative genomes as species-poor, deep-branching taxa were excluded (Figure 1; Supplemental Table S2).

Figure 1.

2.2. Genome size estimation

For the selected genomes, the genome size records available for 465 species (Supplemental Table S3) show that the assembly size is strongly positively correlated with the C-value (Pearson’s r = 0.97, p-value < 0.001). This indicates the reliability of the use of assemblies to estimate genome size. Although a non-linear model is not statistically better than the Weighted Least Squares (WLS), assembly sizes tend to underestimate genome size in comparison to C-values, an effect becoming more and more evident as genome size increases (Supplemental Figure S1). According to the metadata that we could retrieve, whether a genome was assembled using long read or uniquely short read data does not affect the slope of the WLS with assembly size as independent variable (T-test: long reads - p-value = 0.88; short reads - p-value = 0.87). Because it is not affected by sequencing biases, C-value was used when available; otherwise, the predicted C-value according to WLS was employed. The value chosen as genome size estimation is reported for each species in Supplemental Tables S1, S2.

2.3. Transposable elements and genome size variation

Repeat content of a subset of 29 dipteran genomes was previously estimated using EarlGrey v1.3 (Baril et al. 2022) and a wrapper around dnaPipeTE (Goubert et al. 2015), an assembly-based and a read-based pipeline, respectively. dnaPipeTE leverages the sampling of reads at low-coverage to perform de novo assembly of TE consensus sequences: this approach has the advantage of being unbiased by repeat sequences potentially missing from genome assemblies. The results of the two methods overall agree with each other across the scanned genomes (Pearson’s r = 0.88, p-value < 0.001), with the most notable difference being the proportion of unknown elements, generally higher in dnaPipeTE estimations (Wilcoxon signed-rank test, p-value < 0.001; Supplemental Table S4; Supplemental Figure S2). We therefore mined the remaining genomes with dnaPipeTE which is much less computationally intensive. Repeat content could only be estimated for 672 species over the 807 representative genomes (Supplemental Table S2): for the remaining 135, the pipeline could not be run because of unsuitable reads (e.g. only long reads available or too low coverage).

A very strong positive correlation between TE content and genome size is found both across the whole dataset and within taxa (Figure 2A; Supplemental Figure S3; Table 1; Supplemental Table S5). A notable exception concerns the avian clade that deviates from this pattern: the range of TE content is wider than the one of genome size compared to the other clades (Figure 2A), resulting in a weaker power of TEs in explaining genome size variation in this group (Table 1; Supplemental Table S5).

Figure 2.

Table 1.

Alternatively to the impact of TEs on genome size, we investigated whether whole or partial genome duplications could be major factors in genome size variation among animals. BUSCO Duplicated score has indeed a slightly positive effect on genome size, which is however much weaker than that of TEs (Slope = 6.639 · 10⁻⁹, adjusted-R² = 0.022, p-value < 0.001). Of the 24 species with more than 30% of duplicated BUSCO genes, 13 include sturgeon, salmonids and cyprinids, known to have undergone whole genome duplication (Du et al. 2020; Li and Guo 2020; Lien et al. 2016), and five are dipteran species, where gene duplications are common (Ruzzante et al. 2019). In general, TEs appear as the main factor influencing genome size variation across species.

2.4. dN/dS and life history traits as proxies of effective population size

Intensity of effective selection acting on species can be informed by the dN/dS ratio: a dN/dS close to 1 accounts for more frequent accumulation of mildly deleterious mutations over time due to small N_e, while a dN/dS close to zero is associated with a stronger effect of purifying selection. We therefore employed this parameter as a genomic indicator of N_e. We compiled several LHTs from different sources (see 4.6) to cross-check our estimations of dN/dS.

We estimated dN/dS with a mapping method (hereafter referred to as Bio++ dN/dS; Dutheil et al. 2006; Guéguen et al. 2013), and with a bayesian approach using Coevol (hereafter referred to as Coevol dN/dS; Lartillot and Poujol 2011). The two metrics are reported in Supplemental Table S2. Note that for Coevol, we report both the results relative to dN/dS at terminal branches (Table 2) and the correlations inferred by the model (Table 3; Supplemental Table S5).

Table 2.

Table 3.

As expected, Bio++ dN/dS scales positively with body mass and longevity under Phylogenetic Independent Contrasts (PIC) (Table 2; Supplemental Figure S4). dN/dS estimation on the trimmed phylogeny deprived of short and long branches results in a stronger correlation with LHTs, suggesting that short branches might contribute to dN/dS fluctuations (Table 2; Supplemental Figure S5). Coevol dN/dS values are however highly concordant with Bio++ dN/dS (Supplemental Figure S6) and scale positively with body mass and longevity, as well (Table 2).

As for Coevol reconstruction, dN/dS covaries as expected with most of the tested LHTs: dN/dS scales positively with body length, longevity, mass, sexual maturity and depth range in fishes (Table 3); the same is found in mammals, in addition to negative correlations with population density and metabolic rate; in birds, mass, metabolic rate and sexual maturity correlate in the same way with dN/dS, although this is consistently observed only for the GC3-rich gene set (Table 3). Based on the available traits, the two kinds of N_e proxies analysed here correspond in general to each other, supporting dN/dS as a meaningful indicator of long-term effective population size, at least for vertebrate clades. Results are not reported for molluscs and insects, as none and very few records of LHTs (seven species with at least one trait) were available, respectively.

2.5. dN/dS does not predict genome size and overall TE content across metazoans

If increased genetic drift leads to TE expansions, a positive relationship between dN/dS and TE content, and more broadly with genome size, should be observed. However, we find no statistical support for this relationship across all species in the PIC analysis. Similarly, no association is found when short and long branches are removed (Figure 2B-C; Supplemental Figure S7; Table 2). Contrary to our expectations, Coevol dN/dS scales negatively with genome size across the whole dataset (Slope = −0.287, adjusted-R² = 0.004, p-value = 0.039) and within insects (Slope = −1.241, adjusted-R² = 0.026, p-value = 0.018).

Surprisingly, different patterns are observed relative to TE content, as a negative correlation with Coevol dN/dS is detected across all species (Slope = −0.903, adjusted-R² = 0.004, p-value = 0.050) and within Mammalia (Slope = −2.113, adjusted-R² = 0.063, p-value = 0.001), while a positive correlation is found within Actinopteri (Slope = 3.407, adjusted-R² = 0.046, p-value = 0.013) (Table 2). Therefore, the two only significant positive PIC correlations found for birds and fishes are contrasted by results with an opposite trend in other groups. However, such correlations are slightly significant and their explained variance is extremely low.

Overall, we find no evidence for a recursive association of dN/dS with genome size and TE content across the analysed animal taxa as an effect of long-term N_e variation. PIC analysis without the 24 species with more than 30% of duplicated BUSCO genes produced similar results with dN/dS as independent variable (genome size: slope = 0.234, adjusted-R² = 0.002, p-value = 0.102; TE content: slope = 0.819, adjusted-R² = 0.004, p-value = 0.054; Recent TE content: slope = 2.002, adjusted-R² = 0.012, p-value = 0.003), indicating that genomic duplications have a negligible effect on the missing link between dN/dS and genome size.

2.6. Population size and genome size: a complex relationship across clades

Although no strong signal is found across the full dataset using PIC, different trends within different clades are suggested by both PIC and Coevol approaches.

Coevol infers a negative effect of dN/dS on genome size in insects (GC3-poor: CC = −0.330, PP = 0.03; GC3-rich: CC = −0.110, PP = 0.24) and on TE content in mammals (GC3-poor: CC = −0.220, PP = 0.08; GC3-rich: CC = −0.249, PP = 0.06), and a positive correlation (even though below significance threshold) with TEs in fishes (GC3-poor: CC = 0.192, PP = 0.88; GC3-rich: CC = 0.167, PP = 0.84). Additionally, a negative correlation with TE content is found in birds for the GC3-rich geneset, while a positive – yet not significant – trend is found using the GC3-poor geneset (GC3-poor: CC = 0.065, PP = 0.73; GC3-rich: CC = −0.195, PP = 0.03) (Table 3). Even though available for fewer species, LHTs partially support these trends for vertebrates: Actinopteri display a positive correlation between longevity and recent TE content (see 2.7). Instead, mammalian TE content correlates positively with metabolic rate and population density, and negatively with body length, mass, sexual maturity, age at first birth and longevity (Supplemental Table S5). Within Aves, Coevol predicts opposite results for genome size and TE content: genome size associates positively with longevity and mass, as well as with dN/dS, and negatively with metabolic rate, while TEs correlate positively with metabolic rate (and negatively with dN/dS in one case) (Supplemental Table S5). In summary, N_e seems to negatively affect TE content in fishes, and positively in mammals. Importantly, genome size correlations seem to follow the same trends of TE content in these groups, although correlations are weaker and mostly non-significant. In the case of birds, genome size seems rather to be explained by N_e as expected by MHH but not TE content, which instead might have the opposite trend.

2.7. dN/dS weakly correlates with the recent TE content

The global TE content integrally reflects a long history of TE insertions and deletions. To have a glance at the dynamics of TEs on an evolutionary time comparable to that of the level of drift estimated using dN/dS, we additionally examined the quantity of the youngest elements. From the overall TE insertions, we estimated a recent TE content, defined by reads with less than 5% of divergence from consensus, and included it among the traits to model with PIC and Coevol. In PIC analysis, the variation of recent TE content weakly associates with Bio++ dN/dS across the full dataset (Slope = 1.963, adjusted-R² = 0.012, p-value = 0.003) and in mammals (Slope = 4.115, adjusted-R² = 0.024, p-value = 0.031). On the contrary, using Coevol dN/dS this correlation is found in Aves (Slope = 4.982, adjusted-R² = 0.028, p-value = 0.006) and Actinopteri (Slope = 4.365, adjusted-R² = 0.061, p-value = 0.005), but the opposite is detected in mammals (Slope = −3.151, adjusted-R² = 0.024, p-value = 0.032) (Table 2). In agreement with PIC, Coevol reconstruction retrieves a positive association of recent TE content with dN/dS and longevity only in fishes, and a relationship opposite to expectations with dN/dS (CC = −0.296, PP = 0.06) and LHTs in mammals. In contrast with PIC results, a negative correlation between recent TE content and dN/dS is found for birds using the GC3-rich genesets (CC = −0.240, PP = 0.01) (Table 3; Supplemental Table S5).

On the whole, only a very weak positive effect of dN/dS on recent TE insertions is observed across all species. However, considering again the taxa separately, clade-specific patterns emerge: a negative association between population size proxies and recent TE content is jointly found by the two methods only in fishes. Conversely, mammals show a positive correlation between recent TE content and population size proxies. Therefore, the coevolution patterns between population size and recent TE content are consistent with the scenarios depicted by genome size and overall TE content in the corresponding clades.

3. Discussion

Our results demonstrate the absence of a negative relationship between genome size and effective population size across a large dataset of animals, in contrast to the prediction of the MHH (Lynch and Conery 2003; Lynch 2007). Rather, our results highlight heterogeneous patterns within clades, with no consistent response of genome size and TE dynamics to N_e variations.

3.1. Assembly size underrate genome size as genomes grow bigger

Assembly size is commonly used as a measure of genome size. However, the difficulty in assembling repetitive regions generally have it underestimate the actual genome size, in particular when only short reads are employed (Benham et al. 2023; Blommaert 2020; Peona et al. 2018). Consequently, methods that directly measure C-value such as flow cytometry (Doležel and Greilhuber 2010) and Feulgen densitometry (Hardie et al. 2002) are normally preferred as they do not rely on sequence data. Despite applying quality criteria to the assembly, the relationship between assembly size and genome size might still be questioned. However, we show that assembly size can overall approximate genome size quite well and, probably because we removed lower quality assemblies, no effect related to read type (short Illumina vs long ONT/PacBio) was detected. This suggests that the selected assemblies can mostly provide a reliable measurement of genome size, and thus a quasi-exhaustive view of the genome architecture. On the other hand, because a gap with C-value is still present, we integrated this metric to correct assembly size estimations to their “expected C-values”. Similarly, the use of dnaPipeTE allowed us to quantify the repeat content without relying on assembly completeness. In summary, we extensively controlled for the effect of data quality on results and employed methods to minimize it.

3.2. Transposable elements are major contributors to genome size variation

Genome size and TE content have already been reported as tightly linked in eukaryotes (Elliott and Gregory 2015; Kidwell 2002), arthropods (Sproul et al. 2023; Wu and Lu 2019) and vertebrates (Chalopin et al. 2015). Our results are consistent with this perspective across all the animal species analysed, as well as at the level of ray-finned fishes (Reinar et al. 2023), insects (Heckenhauer et al. 2022; Mérel et al. 2024; Petersen et al. 2019; Sessegolo et al. 2016), mammals (Osmanski et al. 2023), and molluscs (Martelossi et al. 2023), strongly indicating TEs as major drivers of genome size variation in metazoans. It should be noted however that we mainly focused on some vertebrate groups and insects, while leaving out many animal taxa with fewer genomic resources currently available including much of the animal tree of life, such as most molluscs, annelids, sponges, cnidarians and nematodes. Even for better studied vertebrates, our datasets are far from comprehensive. For instance, the genomes of squamate reptiles are relatively stable in size but show a high variability in repeat content (Castoe et al. 2011; Pasquesi et al. 2018). A similar case is represented by bird genomes where, according to our analysis and consistently with other studies (Ji et al. 2022; Kapusta and Suh 2017), repeat content has a lower capacity to explain size compared to other clades. This could be due to satellites, whose contribution to genome size can be highly variable (Flynn et al. 2020; Pasquesi et al. 2018; Peona et al. 2023). However, this type of repeat should have been identified by dnaPipeTE. Indeed, the remarkable conservation of genome sizes in this taxon has prompted interpretations involving further mechanisms (see discussion below).

Another way for genomes to grow involves genomic duplications. Although a high proportion of duplicated BUSCO genes may indicate a low haplotype resolution of the assembly, many species with a high duplication score in our dataset correspond to documented duplication cases, suggesting that such BUSCO statistics may provide an insight into this biological process. However, the contribution of duplicated genes to genome size is minimal compared to the one of TEs, and removing species with high duplication scores does not affect our results: this implies that duplication does not impact genome size strongly enough to explain the lack of correlation with dN/dS. Across the animal species considered here, the activity of TEs is therefore a preponderant mechanism of DNA gain, and their evolutionary dynamics appear of prime importance in driving genome size variation.

3.3. Reduced selection efficacy is not associated with increased genome size and TE content

Our dN/dS calculation included several filtering steps by branch length and topology: indeed, selecting markers by such criteria appears to be an essential step to reconcile estimations with different methodologies (M Bastian, F Bénitière, A Marino, pers. comm.). In addition, our analyses resulted to be robust to species pruning by deviant branch lengths. Müller et al. (2022) showed that recent N_e fluctuations might perturb the expected correlation between long- and short-term estimates of N_e. According to the nearly neutral theory, alleles that start at low relative frequencies have a mean fixation of ∼4Ne generations, under the implicit assumption of constant N_e (Kimura and Ohta 1969). This implies that dN/dS, that accounts for the accumulation of substitutions over time, has a weaker sensitivity to short-term changes in N_e compared to estimates based on polymorphism (Müller et al. 2022). Additionally, inferences on simulated and empirical data showed that N_e changes along branches could be captured and generally recapitulated by dN/dS and LHTs in a framework similar to that of Coevol (Latrille et al. 2021). Accordingly, dN/dS assessments by Bio++ and Coevol are highly concordant between each other and with LHTs. Taken together, our results point at the dN/dS found with the two methods as reliable proxies of long-term N_e.

If TEs are ascribable to nearly-neutral mutations, a negative effect of N_e on TE expansion, and consequently on genome size – equivalent to a positive association with dN/dS – is expected. However, no such correlation is observed across the sampled species. To account for a shallower phylogenetic scale, we isolated recently active elements and at the same time explored the same relationships within each clade. Indeed, while the selective effect of elements might be slightly negative as long as they are active, TEs accumulated over long periods of time might be subject to changing dynamics: in the latter case, the pace of sequence erosion could be in the long run independent of drift and lead to different trends of TE retention and degradation in different lineages. Extracting recent elements should thus allow us to have a glimpse of the latest TE colonizations. A positive scaling between the quantity of young TEs and dN/dS found in some cases indicates that relatively recent expansions of TEs could be subject to a more effective negative selection. However, this trend is always very weak and often summarizes that of full TE content within clades. A potential limit of this analysis lies in the application of the same similarity threshold to all species to delimit recent elements. While this is not problematic when comparing species that recently split apart (e.g., Yang et al. 2024), some noise might be introduced at large scale, as the quantity of young repeats that evolved on the same time scale can vary according to the mutation rate and generation time of a species.

Interestingly, the correlation patterns between population size proxies and genomic traits emerging within single clades are distinct and sometimes opposite to the expectations of MHH. Mammals display a negative correlation of dN/dS with TE content, a pattern that is uniformly confirmed by LHTs. Not only does this result corroborate previous findings of no relationship between N_e and genome size in mammals (Roddy et al. 2021), but supports a correlation opposite to the predictions of the MHH. On the other hand, the positive effect of dN/dS observed on TE content in ray-finned fishes might lend support to a role of drift on genome size in this clade: this result is also consistent with a previous study which found a negative scaling between genome size and heterozygosity in this group (Yi and Streelman 2005). In birds, population size seems to negatively affect genome size and positively the TE content, a decoupling that is however not surprising given the higher variation of TE load compared to the restricted genome size range. Contrasting signals from the two genomic traits have already been observed by Ji et al. (2022) who also reported a positive correlation between assembly size and mass, but a negative correlation between TE content and generation time. It is possible for the very narrow variation of the avian genome sizes to impair the detection of consistent signals. An “accordion” dynamic has been proposed whereby higher TE loads are paralleled by equally strong deletional pressures, which could contribute to the maintenance of remarkably small and constant genome sizes in birds, in spite of ongoing TE activity (Kapusta et al. 2017; Kapusta and Suh 2017). Finally, the diffused evidence for a positive and a negative correlation of genome size with body mass and metabolic rate, respectively, is also compatible with the adaptationist perspective of powered flight indirectly maintaining small genome sizes in birds as a consequence of the metabolic needs (Wright et al. 2014; Zhang and Edwards 2012). In insects, dN/dS scales negatively with genome size, but never with TE content. As eusociality appears to bring about selection relaxation (Imrit et al. 2020; Kapheim et al. 2015; Weyna and Romiguier 2021), several studies explored the link between N_e and genome size in this taxon by focussing on social complexity as a proxy, but with contrasting outcomes: Mikhailova et al. (2023) find bigger genome size associated with eusociality in Hymenoptera, but the opposite trend in Blattodea; in contrast and partially in accordance with our findings, Kapheim et al. (2015) and Koshikawa et al. (2008) report less abundant TEs in eusocial hymenopters and smaller genomes in eusocial termites, respectively. While the approximation of N_e based on dN/dS should allow for a quantification of selection efficacy in wider terms than sociality traits, the investigated evolutionary scale might hold an important role in the outcome of such analyses. First, and specifically relative to insects, genome size seems to be subject to different evolutionary pressures – either selective or neutral – within different insect orders (Cong et al. 2022), implying that increased drift might not necessarily produce the same effect on genome size across all insect groups. More generally, the five defined clades cover quite different time scales: insects and molluscs have much more ancient origins than mammals and birds, and such distant groups also evolve at very different evolutionary rates, making it difficult to characterize the evolution of their traits on the same evolutionary scale. Nevertheless, the results are still valuable in highlighting the absence of drift effect on genome size in the long-term evolution of such broad groups, in contrast to previous work focusing at the population level or on recently diverged species (Cui et al. 2019; Mérel et al. 2021; Yang et al. 2024). At the same time, as noted by Mérel et al. (2024), comparing very distantly related species - as the insect and molluscan species of our dataset - might overshadow any relationship between genome size and N_e, either due to dN/dS predicting power being weakened by branch saturation, deep N_e fluctuations not being detected by our methods, or to additional factors affecting long-term genome size evolution.

3.4. Do lineage-specific TE dynamics affect genome size evolution?

Our findings do not support the quantity of non-coding DNA being driven in a nearly-neutral fashion by genetic drift. Moreover and in agreement with previous analyses (Pasquesi et al. 2018), we find that TE activity can, under comparable drift effects, give place to lineage-specific outcomes that mostly do not seem to depend on effective population size. These results contrast with those of other large-scale analyses which instead support the predictions of the drift-barrier hypothesis for a general impact of N_e on other genomic features, notably mutation rate (Bergeron et al. 2023; Lynch et al. 2023; Wang and Obbard 2023) and splicing accuracy (Bénitière et al. 2024). To put this in perspective, it should be emphasized that, in the framework of the MHH, the success of nearly-neutral alleles not only depends on N_e but also on the liability to mutation of non-coding DNA (Lynch et al. 2011). Overall, we studied N_e variation without accounting for the mutagenic burden posed by TEs in different lineages, inherently assuming the same distribution of selective effect and a constant activity of TEs in all species and among TE insertions. However, it is known that TEs are subject to waves of activity rather than a uniform pace of transposition (Arkhipova 2018). Moreover, given the broad phylogenetic scale of our dataset, it is likely for different levels of mutational hazard to be acting across genomes due to different “host-parasite” dynamics in different animal groups (Ågren and Wright 2011). For instance, while the big genomes of salamanders are not related to small N_e, the low synonymous substitution rates and low degree of deletions due to ectopic recombination suggest weak mutational hazard of TEs that possibly contributes to the maintenance of genomic gigantism in this group (Frahry et al. 2015; Mohlhenrich and Mueller 2016; Rios-Carlos et al. 2024). Additionally, lineage-specific TE dynamics themselves might underlie different genomic architectures: for example, mammalian genomes are generally characterized by one preponderant type of active element and by a long-term retention of old TEs (Osmanski et al. 2023; Sotero-Caio et al. 2017), as in human where a very small proportion of active elements (<0.05%) is unlikely to impose a mutation rate causing genome size variation (Mills et al. 2007). Conversely, squamate and teleost fish genomes are smaller and characterized by several, simultaneously active and less abundant TE types (Duvernell et al. 2004; Furano et al. 2004; Novick et al. 2009; Pasquesi et al. 2018; Volff et al. 2003). These different patterns of genomic organization seem overall associated with different rates of elements’ turnover (Blass et al. 2012; Lavoie et al. 2013; Novick et al. 2009; Volff et al. 2003). All such variables might alter the selective effect and differentiate TEs from gradually and constantly evolving alleles, eventually contributing to the lack of association between long-term N_e and genome size. Finally, Kapusta et al. (2017) showed that large-scale deletions can be as important as DNA gain in determining genome size, thus questioning the assumption of the rate of elements insertion being greater than their removal rate (Lynch 2007). This implies that the contribution of TEs constitutes just one side of the coin and that deletion bias could also drive the divergence of genome size across lineages, as suggested by several studies linking negatively deletion rates with genome size (Frahry et al. 2015; Ji et al. 2022; Kapusta et al. 2017; Wang et al. 2014).

3.5. Perspectives

There is some evidence of selection acting against TEs proliferation. In Anolis lizards, the ability of TEs to reach fixation can vary between populations of the same species according to population size (Ruggiero et al. 2017; Tollis and Boissinot 2013). Furthermore, a negative effect of N_e was found on genome size and TE content expansion at the inter-population level in Drosophila suzukii (Mérel et al. 2021) and at the interspecies level in fruit flies (Mérel et al. 2024), asellid isopods (Lefébure et al. 2017) and killifishes (Cui et al. 2019), supporting the role of genetic drift in determining recent differences in genome size among closely related animal species. Given the very different taxonomic scale of such works and ours, and with the perspective of lineage-specific interaction between genome and genomic parasites in mind, our negative results for the MHH at metazoan scale are not incompatible with an effect of N_e on genome size within specific clades. In a nutshell, although an increase in genetic drift seems to lead to the short-term accumulation of transposable elements, this effect is not visible in the long-term, suggesting that it fades over time. A general mechanism of selection preventing the proliferation of non-coding DNA and TEs in animals might exist but its effect be detectable only at a sufficiently short evolutionary time. In this sense, the lack of evidence for MHH in other clade-specific studies might be due to the methodological challenges of either estimating a suitable marker of N_e or investigating too distantly related lineages. Moreover, the contrasting outcomes of such studies might reflect an actual variability in the selective effect of TEs not compatible with a general selection mechanism. Further reducing the phylogenetic scale under study and systematically exploring the effect of N_e within independent biological systems could therefore provide an alternative way to test the impact of drift, while removing the confounding effects due to different genomic backgrounds.

4. Methods

4.1. Dataset

All the metazoan reference assemblies available as of November 14th 2021 were used, except for insect genomes which were drawn from Sproul et al. (2023), for a total of 3,214 assemblies. For each assembly, quality metrics were computed with Quast 5.0.2 (Gurevich et al. 2013) and genome completeness was assessed with BUSCO 5.2.2 using the 954 markers of the metazoa_odb10 geneset (Supplemental Table S1). Availability of raw reads was verified with SRA Explorer (https://github.com/ewels/sra-explorer). All the assemblies with either a contig N50 smaller than 50 kb, less than 70% of complete BUSCO orthologs, or without available reads were excluded from TE and dN/dS analyses. The subdivision into Actinopteri (N = 148), Aves (N = 260), Insecta (N = 189), Mammalia (N = 182), Mollusca (N = 28) was adopted to perform alignments, phylogenies, dN/dS estimation with Bio++ and Coevol runs (see below).

4.2. Genome size estimation

Assembly sizes and C-values were jointly used to estimate genome size. C-values measured by either flow cytometry (FCM), Feulgen densitometry (FD) or Feulgen image analysis densitometry (FIA) were collected from https://www.genomesize.com/ (last accessed 6 october 2022) for all available species of our initial dataset with contig N50 ≥ 50 kb, totalling 465 measurements for 365 species (Supplemental Table S3). To assign a unique C-value, when multiple values were present for one species, the most recent one was retained and, if dates were the same, the average was used. For all the other species having contig N50 ≥ 50 kb but with no available C-value record, genome size was calculated as an expected C-value predicted from a WLS where the 465 FCM, FD and FIA estimations were the dependent variables and assembly sizes were the independent variables (for details see https://github.com/albmarino/Meta-analysis_scripts). Out of all the records in this training dataset for genome size, 93 correspond to ray-finned fishes, 93 to mammals, 92 to birds, 106 to insects, and 9 to molluscs, overall mirroring the taxa represented in the final dataset. For the purpose of our analysis, C-values were used for the species for which such data were available, while the expected C-value was used as genome size estimation in all the other cases, regardless of the type of sequencing data used for the assembly (Supplemental Tables S1, S3).

4.3. Gene alignment

The 954 annotated single-copy BUSCO genes were aligned with the pipeline OMM_MACSE 11.05 using MACSE 2.06 (Ranwez et al. 2018; Scornavacca et al. 2019). Alignments were performed separately within each clade - Actinopteri, Aves, Insecta, Mammalia and Mollusca.

4.4. Phylogeny

Phylogenies were computed separately for each clade with IQ-TREE 1.6.12 (Nguyen et al. 2015). JTT+F+R10 substitution model was selected with ModelFinder (-m MFP option) (Kalyaanamoorthy et al. 2017). For reasons of computing power and time, we have reconstructed the phylogenies of each clade independently and then grouped them together to create a single complete phylogenetic tree (see below). The same set of 107 concatenated BUSCO amino-acid sequences was used to calculate all the phylogenies. However, since this produced a spurious relationship in the mammalian tree with paraphyly of primates, an alternative set of randomly selected 100 genes was used instead for the phylogeny of Mammalia (Supplemental Table S6). Each phylogeny was rooted using an outgroup species belonging to its respective sister clade: the outgroup sequences were added to gene alignments with the enrichAlignment function from MACSE; the outgroup + clade gene alignments were concatenated and used to recompute the outgroup + clade phylogeny taking into account the previously computed tree topology of the clade. The outgroups were then removed, and rooted clade phylogenies were merged together manually using the tree editor program Baobab (Dutheil and Galtier 2002). Finally, 50 top-shared genes across all species were chosen among the set of 107 genes (Supplemental Table S6) to recalculate the branch lengths of the whole dataset phylogeny: with the MACSE program alignTwoProfiles the nucleotide gene alignment of one clade was joined to the one of its respective sister clade until achievement of the whole dataset alignment. Branch lengths were then estimated based on the 50-genes concatenate and the tree topology (for the detailed workflow, see https://github.com/albmarino/Meta-analysis_scripts).

4.5. dN/dS estimation

When genetic drift is strong, slightly deleterious mutations are more likely to reach fixation than under conditions of high N_e and more effective selection (Ohta 1992). The genome-wide fixation rate of non-synonymous mutations is expected to drive the dN/dS ratio due to nearly-neutral mutations responding to different N_e: a higher dN/dS accounts for more frequent accumulation of mildly deleterious mutations over time due to small N_e, while lower dN/dS is associated with a stronger effect of selection against slightly deleterious non-synonymous mutations due to high N_e (James et al. 2016; Romiguier et al. 2014; Weyna and Romiguier 2021; Woolfit and Bromham 2005). This is also supported at the polymorphism level, with higher pN/pS and accumulation of slightly deleterious mutations in smaller populations (Leroy et al. 2021; Dussex et al. 2023).

Before dN/dS calculation, sequences with more than 10% of their length occupied by insertions were preemptively removed from BUSCO alignments. Estimation of dN/dS on either very long or short terminal branches might lead to loss of accuracy due to a higher variance of substitution rates or to branch saturation, respectively. Furthermore, shared polymorphism can be captured in the substitution rates when closely related species are compared, and further contribute to bias dN/dS (Mugal et al. 2020). To correct for such issues, genes with deviant topology were identified and removed from every clade with PhylteR using default parameters (Comte et al. 2023). Moreover, genes exhibiting branch lengths shorter than 0.001, for which dN/dS could have a large variance, were also not integrated in the dN/dS calculation of a species.

We then used bppml and mapnh from the Bio++ libraries to estimate dN/dS on terminal branches (Dutheil et al. 2006; Guéguen et al. 2013; Guéguen and Duret 2018; Romiguier et al. 2012). bppml calculates the parameters under a homogenous codon model YN98 (F3X4). Next, mapnh maps substitutions along the tree branches and estimates dN and dS. More precisely, the substitution rate is given by the number of substitutions mapped according to the model parameters normalized by the number of substitutions of the same category (i.e. synonymous, non-synonymous) that would occur under the same neutral model (Bolívar et al. 2019). Therefore, dN and dS are calculated for each species as follows:

Where n is the number of genes, K is the substitutions count as mapped by the substitution model calculated for the gene, O is the substitutions count as mapped under the same neutral model, and l is the branch length of the given species for that gene. In addition to the gene filtering described above, Bio++ dN/dS was recalculated on a reduced dataset where the longest (> 1) and shortest (< 0.01) branches were removed, in order to ensure that substitution saturation and segregating polymorphism did not influence the results. Terminal branches with more than 1 and less than 0.01 amino-acid substitutions per site were removed, and dN/dS recalculated on the trimmed phylogenies with the same method described above.

The same metric was estimated with Coevol 1.6 (Lartillot and Poujol 2011). Coevol models the co-evolution of dN/dS and continuous traits along branches following a multivariate Brownian diffusion process, thus reducing the variance in the dN/dS of the smallest branches (Lartillot and Delsuc 2012). Bio++ dN/dS was therefore compared with dN/dS estimated by Coevol on terminal branches to verify the consistency between the two methods.

4.6. Compilation of life history traits

Life-history traits – body mass, longevity, generation time, among others – are found to be related to N_e in mammals, birds, and amniotes in general (Bolívar et al. 2019; Figuet et al. 2016; Nikolaev et al. 2007; Popadin et al. 2007). Available LHTs were assigned to the species of our dataset using information from several resources. Adult body mass, body length, maximum longevity, basal metabolic rate, age at first birth, population size and population density were assigned to mammalian species using PanTHERIA (Jones et al. 2009). For birds, body mass information was extracted from Dunning (2007). Shallow to deep depth range, longevity in the wild, body length and body mass were compiled for ray-finned fishes from www.fishbase.org using the rfishbase R package (Boettiger et al. 2012). Additionally, age at sexual maturity, adult body mass, maximum longevity and metabolic rate were extracted from AnAge (Tacutu et al. 2018), as well as body mass and metabolic rate from AnimalTraits (Herberstein et al. 2022): such data were used to complement information when missing from the databases cited above. All the retrieved LHTs and their relative source are reported in Supplemental Table S2.

4.7. TE quantification

TEs were annotated with a pipeline employing dnaPipeTE (Goubert et al. 2015) in two rounds (https://github.com/sigau/pipeline_dnapipe). Raw reads are filtered with UrQt (Modolo and Lerat 2015) or fastp (Chen et al. 2018) and undergo a first “dirty” dnaPipeTE round. The obtained dnaPipeTE contigs are mapped against a database of organellar, fungal, bacterial and archaean reference sequences with Minimap2 (Li 2018) and the original reads matching contaminant sequences removed with SAMtools (Danecek et al. 2021). Finally, the quality- and contaminant-filtered reads are used to perform a second “clean” dnaPipeTE round. dnaPipeTE was configured with the Dfam 3.5 and RepBaseRepeatMaskerEdition-20181026 repeat libraries and was run with a genome coverage of 0.25.

To verify the consistency of dnaPipeTE estimations, the dnaPipeTE-based pipeline was benchmarked on a subset of 29 dipteran species against EarlGrey 1.3, an automated pipeline performing TE annotation on genome assemblies (Baril et al. 2022). EarlGrey was configured with the same libraries as dnaPipeTE and was run with ‘metazoa’ as search term. The results of the second dnaPipeTE round were used to extrapolate the total and recent TE content, the latter being defined by all the reads below 5% of divergence from the corresponding consensus. The two contents were extracted by adapting dnaPT_landscapes.sh from the dnaPT_utils repository (https://github.com/clemgoub/dnaPT_utils) to a custom R script (https://github.com/albmarino/Meta-analysis_scripts).

4.8. Gene duplication

To account for the effect of whole genome or segmental duplications, we used the BUSCO Duplicated score: if big-scale duplication events recently took place, a higher score should be observed genome-wide even for conserved genes. As many of the genomes with BUSCO Duplicated score above 30% corresponded to reported cases of genomic duplication, we used this threshold to perform PIC analysis with and without species whose genome size is potentially more affected by duplication events.

4.9. Phylogenetic independent contrasts and Coevol reconstruction

The correlations of Bio++ and Coevol dN/dS with LHTs, as well as with genome size and TE content, were tested with PIC to correct for the covariation of traits due to the phylogenetic relatedness of species (Felsenstein 1985). PIC were performed on the whole dataset, the trimmed dataset, and within every clade with the R packages ape 5.7.1 (Paradis et al. 2004), nlme 3.1.162, and caper 1.0.1. Results were plotted with the ggplot2 3.4.2 package. Additionally, Coevol 1.6 was run to test for the coevolution of dN/dS and traits: sequence substitution processes and quantitative traits such as LHTs, genome size, TE content, are here assumed to covary along the phylogeny as a multivariate Brownian motion process. Coevol infers trait values on internal nodes and terminal branches (those used for PIC), as well as correlation coefficients and their relative posterior probabilities. Due to computational limitations, Coevol analysis was carried out separately on every clade and on a limited number of genes. As dN/dS can be influenced by GC-biased gene conversion (Bolívar et al. 2019; Ratnakumar et al. 2010), two sets of 50 genes similar in GC-content at the third codon position (GC3) were defined. Markers and species were subject to a more stringent filtering in order to have as much information as possible for each species. 16 species with less than 50% of single-copy orthologs were further filtered out from the dataset. In addition to the PhylteR step, only genes represented in at least 95% of the species of a clade were retained. From those, the 50 GC3-poorest and the 50 GC3-richest genes were chosen. Coevol was then run with both the gene sets for every clade. Convergence of the MCMC chains were checked visually by plotting the evolution of statistics. Likelihood values and correlations were estimated, running the chains for a minimum of 1,000 steps and discarding the first 400 steps as burn-in.

5. Data access

The accession numbers, SRA experiments and C-values used in this study are reported in Supplemental Table S1. C-values are available at https://www.genomesize.com/. The respective sources of LHTs are reported in Supplemental Table S2. Detailed commands and custom scripts can be found at https://github.com/albmarino/Meta-analysis_scripts.

6. Competing interest statement

The authors declare that they have no conflict of interest.

Acknowledgements

We are grateful to Laurent Duret, Nicolas Lartillot, Tristan Lefébure, Mélodie Bastian and Florian Bénitière for the helpful discussions. We thank Mélodie Bastian and Florian Bénitière for methodological exchanges to ensure to have a dN/dS accurately estimated. We thank Laurent Guéguen for his help with dN/dS analyses. Analyses were run on the servers of the Montpellier Bioinformatics Biodiversity platform (MBB).

This research was funded by the French National Research Agency (ANR-20-CE02-0008-01 “NeGA”). A CC-BY public copyright licence has been applied by the authors to the present document and will be applied to all subsequent versions up to the Author Accepted Manuscript arising from this submission, in accordance with the grant’s open access conditions.