# Relaxed Purifying Selection is Associated with an Accumulation of Transposable Elements in Flies

Vincent Mérel, Théo Tricou, Nelly Burlet, Annabelle Haudry

► To cite this version:

Vincent Mérel, Théo Tricou, Nelly Burlet, Annabelle Haudry. Relaxed Purifying Selection is Associated with an Accumulation of Transposable Elements in Flies. Molecular Biology and Evolution, 2025, 42 (6), 10.1093/molbev/msaf111. hal-05186701

HAL Id: hal-05186701

https://hal.science/hal-05186701v1

Submitted on 25 Jul 2025

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
RELAXED PURIFYING SELECTION IS ASSOCIATED WITH AN
ACCUMULATION OF TRANSPOSABLE ELEMENTS IN FLIES

Vincent Mérel¹,², Théo Tricou¹, Nelly Burlet¹ & Annabelle Haudry¹,*

¹ Universite Claude Bernard Lyon1, LBBE UMR 5558, CNRS, Villeurbanne, France

² Department of Ecology and Evolution, University of Lausanne, Lausanne, Switzerland

* corresponding author: annabelle.haudry@univ-lyon1.fr

**running title:** Genetic drift leads to TE accumulation in flies
ABSTRACT

Although the mechanisms driving genome size evolution are not yet fully understood, one potentially important factor is the dynamics of the accumulation of transposable elements (TEs). Since most TEs are neutral or slightly deleterious, a negative correlation between genome size and the efficacy of selection is expected. However, previous empirical studies on closely related species with distinct life history traits (thought to undergo different selective regimes) have yielded inconsistent results. Here, we perform the first large-scale analysis of the effect of genetic drift on genome size evolution, without any prior assumption on the amount of genetic drift. We reconstructed a phylogeny based on whole-genome data (2,242 genes) for 77 Drosophilid species to examine correlations between genome size, TE content, and the efficacy of selection (using *dN/dS* ratios of non-synonymous to synonymous divergence). Using an integrative approach that controls for shared evolutionary history, we reveal that genome-wide *dN/dS* are strongly positively correlated with genome size and TE content, particularly in GC-poor genes. This study suggests the critical importance of controlling for heterogeneity in base composition when estimating *dN/dS*. Furthermore, we emphasize that the lack of evidence for TE accumulation due to increased genetic drift in several previous studies may be due to a secondary effect of changes in life history traits (i.e. asexuality) on TE dynamics. In conclusion, this work provides evidence for TE proliferation in fly genomes when purifying selection is reduced, shedding new light on the role of TEs and genetic drift in the evolution of genome architecture.

**1 Introduction**

Eukaryotic species exhibit a surprisingly wide range in genome size: the genome of the parasitic nematode *Pratylenchus coffeae* is approximately 0.019 Gbp (Gregory, 2005), whereas that of the fork fern *Tmesipteris oblanceolata* has recently been estimated to reach 160 Gbp (Fernández et al., 2024). The accumulation of repeated elements is thought to be the main driver of differences in genome size between Eukaryotes. In particular, the amount of transposable elements (TEs) is strongly positively correlated with genome size in various taxa (Biémont and Vieira, 2006; Tenaillon et al., 2010; Elliott and Gregory, 2015). TEs are selfish genetic elements (Orgel and Crick, 1980), mostly neutral or deleterious to their host, that proliferate, copying and pasting themselves to other genomic locations.

TEs may disrupt genes, interfere with mRNA expression levels and splicing, induce chromosomal rearrangements via ectopic recombination, and provide ‘excess’ DNA that represents a mutational liability because it increases the potential and target size for harmful mutations (Lynch and Conery, 2003; Mérel et al., 2020). Due to the deleterious effects of TEs, purifying selection is expected to prevent TE fixation and thus maintain small genomes. The effective population size (*Ne*), which affects the efficacy of selection (*Ne* \* *s*, where *s* denotes the selection coefficient), has been proposed as a key parameter controlling TE accumulation in the non-adaptive hypothesis ((Ohta, 1992), later popularised by (Lynch and Conery, 2003)). Since the efficacy of selection is positively correlated with *Ne*, TEs are expected to accumulate in lineages where *Ne* is reduced. Accordingly, Lynch and Conery (2003) found a negative correlation between synonymous diversity (an estimate of 4*Ne*μ, μ representing the mutation rate) and genome size across species. However, Lynch and Conery (2003) findings appear to be largely explained by shared common evolutionary history
36

37

(Whitney et al., 2010), and rely on the comparison of very distant species with drastically different biology, which could
act as a confounding factor (Charlesworth and Barton, 2004).

38

39

40

41

42

43

44

45

46

47

48

49

50

51

52

The distribution of quantitative traits among species is not independent of the length of the shared history: closely related species tend to resemble each other more than expected by chance due to phylogenetic inertia (Felsenstein, 1985). As a result, it is necessary to account for relatedness and the resulting dependencies in comparative analyses. Approaches that explicitly account for the phylogenetic signal when assessing correlations between traits rely on stochastic models of trait evolution — typically assuming that the traits jointly evolve according to a multivariate Brownian process. The covariance matrix of this multivariate process is then estimated either by maximum likelihood (e.g. phylogenetic generalised least squares PGLS, (Martins and Hansen, 1997)) or by Bayesian inference (Pagel and Harvey, 1988), accounting for phylogenetic dependencies in the estimation of the correlations. Methods at the intersection between the classical comparative method and molecular phylogenetics (Lartillot and Poujol, 2011; Lartillot, 2014) have suggested that a substantial amount of information about ancestral traits could be obtained from the evolution of genetic sequences. The fundamental idea is as follows: (1) detailed patterns of genetic sequence evolution can be fairly accurately inferred along a phylogeny; (2) the variation in evolutionary rate along the branches of the tree can be inferred; (3) more fundamentally, these evolutionary rates and parameters can be seen as continuous traits and modelled as Brownian processes. All this can be achieved in a single Bayesian probabilistic model, combining the molecular phylogenetic tree and the comparative method dimensions of the question.

53

54

55

56

57

58

59

60

61

62

63

64

65

66

67

68

Beyond the genetic diversity used in Lynch and Conery (2003), a variety of traits could be related to *Ne*. According to the nearly neutral theory of evolution, smaller *Ne* should lead to an accumulation of slightly deleterious mutations, resulting in a higher ratio of non-synonymous to synonymous divergence, or *dN/dS*, which serves as a measure of the efficacy of selection on protein sequence (Castillo et al., 2011; Bolívar et al., 2019). Several life history traits have also been previously reported to affect *Ne* (e.g. maturity or longevity in mammals (Lartillot and Delsuc, 2012) or mating system in plants (Pollak, 1987)). Therefore, it is possible to test the effect of a reduction in *Ne* on genome architecture using comparisons among closely related species with distinct life history traits. Furthermore, population genetic theory predicts that selection on codon usage should be much weaker than selection acting on the protein sequence (typically on the order of $4Ne \* s \approx 1$, where *s* denotes the selection coefficient in favour of optimal codons). Since weakly deleterious mutations are most susceptible to fixation when *Ne* is reduced (McVean and Charlesworth, 1999), codon usage bias metrics may provide alternative estimates of the efficacy of selection. While some fruit flies with large *Ne* show higher codon usage bias than their relatives (Haddrill et al., 2010), Kessler and Dean (2014) found no direct relation between codon usage bias and *Ne* variation in mammals. Indeed, codon usage bias results from a complex balance between mutation patterns, selection, recombination (via GC-biased gene conversion) and genetic drift (Parvathy et al., 2022). It thus remains an open question whether codon usage metrics can reflect *Ne* variation in large populations.

69

70

71

72

73

To date, comparative studies of closely related taxa testing for the expected negative relationship between genome size/TE content and *Ne* have reported conflicting results. Some have shown an accumulation of TEs resulting from a reduction of *Ne*: in ectomycorrhizal fungi in symbiotic fungi (Hess et al., 2014); in subterranean *Aselloidea* (Lefébure et al., 2017); in polyploid plants (Baduel et al., 2019); in eusocial snapping shrimps and termites (Korb et al., 2015; Chak et al., 2021); in invasive populations of *D. suzukii* (Mérel et al., 2021)—, others found no significant excess in TE
74

75

76

77

abundance in species with reduced *Ne* compared to their close relatives —in asexual daphnia (Schaack et al., 2010; Jiang et al., 2017), primroses (Ågren et al., 2015), arthropods (Bast et al., 2016) or yeast populations (Bast et al., 2019); in selfing *Caenorhabditis* nematodes (Fierst et al., 2015); in eusocial bees (Kapheim et al., 2015), or hymenoptera (Ardila-Garcia et al., 2010).

78

79

80

81

82

83

84

85

The main drivers of repeatome expansion in fly genomes remain to be identified, especially with regard to the respective roles of genetic drift and selection. *Drosophila* flies present a three-fold variation in genome size, reflecting the variation in TE content (Sessegolo et al., 2016). Contrary to theoretical expectations, a study of 12 distant species suggested that greater levels of purifying selection were associated with greater abundance of euchromatic TE, using *dN/dS* as a proxy for *Ne* (Castillo et al., 2011). This finding left the mechanism underlying genome size variation in this clade unexplained. Analyses by Sessegolo et al. (2016) revealed that strong phylogenetic inertia partially explains genome size differences, with a sampling bias towards species with small genome closely related to *D. melanogaster*, which could easily be corrected with a more balanced sample across the Drosophilid phylogeny.

86

87

88

89

Here, we investigate genome size, TE content and the efficacy of selection in 82 fly lineages from 77 species. Our study includes historical *dN/dS* reconstruction along with ancestral genome size and TE content, and tests for intrinsic correlations between these variables using a Bayesian method that accounts for phylogenetic relationships. In addition, we tested for a potential correlation between *dN/dS* and codon usage bias measures.

90

# 2 Results

91

## 2.1 Forty novel fly assemblies

92

93

94

95

96

97

98

99

100

101

102

103

104

105

106

107

108

109

To study the evolution of genome size, TE content, and effective population size in flies, we collated Next Generation Sequencing (NGS) data for 84 Drosophilid lineages. This dataset includes 79 species from the genera *Drosophila* and *Zaprionus*, with a particular focus on the subgenera Drosophila and Sophophora, and spans 25–40 million years of divergence (Obbard et al., 2012) (supplementary table S1). The dataset consists of 44 publicly available samples and 40 new samples that were carefully chosen to obtain a balanced set of species along the phylogeny (at least three species were sampled per Drosophila group or Sophophora subgroup based on taxonomy, supplementary tables S2 and S3). For each of the 40 new samples, we constructed genomic assemblies using short-insert paired-end reads. The mean N50 of the assemblies obtained is 20Kb (from 3Kb to 92Kb, supplementary table S2), with on average 89% of the 3,285 dipteran BUSCO genes complete (supplementary figure 1). Two assemblies with a high percentage of duplicated BUSCO genes (≥ 20) were not considered further (*D. algonquin* and *D. tristis*). The 40 new assemblies include 14 species for which, to our knowledge, no genome assembly was available so far: *D. affinis*, *D. helvetica*, *D. imaii*, *D. lucipennis*, *D. lutescens*, *D. madeirensis*, *D. mercatorum*, *D. microlabis*, *D. mimetica*, *D. pallidosa*, *D. paralutea*, *D. prostipennis*, and *D. tsukubaensis*, according to NCBI on July 2, 2024. For publicly available and new assemblies, genes were predicted (using AUGUSTUS). Although some assemblies were fragmented, the number of predicted genes was generally close to the 13,986 annotated coding genes of *D. melanogaster* (according to FlyBase release 112.10, supplementary tables S2 and S3). This suggests that most of the coding regions were assembled. Furthermore, it is important to note that the estimates of genome size or TE content are independent of the quality of the assemblies (see the Methods section for details).
## 2.2 *Drosophila* phylogeny reconstruction

We used a *de novo* clustering approach to identify single-copy orthologous genes present in at least 50% of the genomes (i.e. 41 genomes). The resulting 2,242 orthologous genes were aligned and then concatenated to reconstruct a phylogenetic tree, using *Scaptodrosophila lebanonensis* as an outgroup (figure 1). This Maximum Likelihood (ML) phylogeny reports highly confident bootstrap support of the nodes (only four nodes with UFBoot <100) and identifies nine main clades, consistent with the taxonomic Drosophila and Sophophora groups.

This tree was compared to an ASTRAL supertree, i.e. an estimation of the true species tree from the 2,242 unrooted gene trees. Our ASTRAL supertree topology is discordant with the ML topology in only four relationships (corresponding to a normalized Robinson Foulds of 0.0127), three of them being associated with the least supported nodes in the ML tree. First, *D. micromelanica* is found at the root of the *melanica* group in the ASTRAL tree, while it is *D. melanica* in the ML tree (UFBoot = 65, figure 1). Second, *D. hydei* is sister to the species group of *D. navojoa*, *D. arizonae*, *D. mojavensis* in the ASTRAL tree, while it is also sister to *D. mercatorum* lines in the ML tree (UFBoot = 100). Third, *D. microlabis* is sister to the obscura subgroup only in the ASTRAL tree, but to the *obscura* and the *subobscura* subgroups in the ML tree (UFBoot = 97). Last, *D. takahashii* is at the root of the whole *takahashii* subgroup in the ASTRAL tree, but sister to *D. lutescens* and *D. paralutea* in the ML tree (UFBoot = 91).

We compared the ML topology presented in figure 1 with the Suvorov' ML phylogenomic tree (Suvorov et al., 2022) using the 55 species shared between the two datasets (supplementary figure 2). Only two incongruencies were found, resulting in a normalized Robinson Foulds of 0.07692 distance metric between trees (Robinson and Foulds, 1981): (1) the sister species of *D. simulans* is *D. sechellia* in our topology while *D. simulans* groups with *D. mauritiana* in the Suvorov tree. (2) *D. ercepeae* groups in the *montium* subgroup (*melanogaster* group) while it branches at the basis of the *ananassae* subgroup (*melanogaster* group) in the Suvorov tree (Suvorov et al., 2022), as well as in Van der Linde's tree (Van Der Linde et al., 2010).

## 2.3 Strong phylogenetic inertia on genome size and TE content

We retrieved flow cytometry estimates of genome size for 77 of the 82 samples. 54 of these estimates come from the literature (Hjelmen et al., 2019; Gregory, 2005), and the remaining 23 from this study. The genomic size ranges from 136.92 Mb for *D. busckii* to 332.52 Mb for *D. virilis*, with a mean of 208.7 Mb (supplementary table S4, figure 1). We found a Pagel's lambda very close to one for this trait, indicating a strong effect of shared ancestry (Pagel's λ = 0.95, p = 4.2·10⁻⁶).

We estimated repeat content in all 77 samples using dnaPipeTE (figure 1, supplementary table S5 and S6)(Goubert et al., 2015). Briefly, dnaPipeTE reconstructs TE sequences from the assembly of a subsample of reads (a depth less than 1X ensures that only repeated sequences can be assembled), annotates these sequences by homology to the RepBase database (Kapitonov and Jurka, 2008), and then estimates the proportion of reads mapping on each category of repeat. This proportion can be converted into a number of base pairs based on the estimate of genome size. We found that the total amount of TE ranges from 3.48 Mb for *D. busckii* to 120.42 Mb for *D. suzukii*. The category with the highest average abundance is LTR elements (supplementary table S6).
<figure><img src="image_2.png" /></figure>

**Figure 1: Phylogenetic tree figuring genome size and TE abundances.** The tree was reconstructed using 2,242 orthologous genes and a maximum likelihood approach. Bootstraps values are indicated at the node when different from 100. Taxonomic groups according to Schoch et al. (2020) are indicated on the right. Genome size (Mb) was estimated using flow cytometry, and genomic repeat content (Mb) was estimated from short reads using dnaPipeTE (Goubert et al., 2015).
We compared our estimates of the abundances of TEs with those obtained in Kim et al. (2021) where the authors estimated the repeat content in long-read assemblies of 101 Drosophilid species using homology to sequences in online databases. For the 36 species common to the two studies, we found significant correlations for the abundance of all TE types ($R^2$ ranging from 0.40 to 0.46, $p_{adj} < 0.05$, supplementary figure 3). Note that for this comparison 1) we considered retrotransposons (estimated in Kim et al. (2021)) to correspond to the sum of LINEs and LTR elements, 2) RC to correspond to Helitrons (Kapitonov and Jurka, 2008), and 3) for SINE, we were unable to test any correlation as this category was not present in Kim et al. (2021). Our results suggest a very strong phylogenetic inertia on TE content. Pagel's $\lambda$ indicates phylogenetic signal for TEs as a whole (Pagel's $\lambda = 1.00$, $p_{adj} = 1.64 \cdot 10^{-5}$), but also for DNA, Helitron, LINEs and LTR elements separately. Only for SINEs the results were not significant.

Despite phylogenetic inertia acting on genome size and TE abundances, both variables were found to vary not only within subgenus (Sophophora and Drosophila), but also between groups that make up these subgenera. In particular, in the melanogaster group (Sophophora), *D. suzukii* experienced an increase in genome size associated with an accumulation of diverse repeated elements, especially since its divergence from *D. biarmipes*, as previously reported by Sessegolo et al. (2016). Furthermore, in the obscura group (Sophophora), *D. guanche* (n=2) and *D. subobscura* (n=2) lineages exhibit less TEs than close relatives (e.g. *D. obscura* (n=2) or *D. imaii* (n=1)).

## 2.4 Genome size correlates with abundances of various repeat types

We tested pairwise correlations between genome size and abundances of the different TEs using two different methods, both controlling for phylogenetic inertia (table 1). The first, which relies on ancov (Lartillot, 2014), performs a Bayesian reconstruction of traits to estimate the correlation between all of them in a single test. The second corresponds to multiple tests of Phylogenetic Generalised Least Squares (PGLS, Martins and Hansen (1997)), adjusted for multiple testing. We compared the ancov and PGLS results: there is an almost perfect correlation between the correlation coefficient returned by ancov and PGLS (supplementary figure 4). Eleven traits were found to be significantly correlated with both methods, three private to ancov, and two to PGLS. Overall, we found genome size to be positively correlated with the total amount of TEs (ancov correlation coefficients of 0.62). However, the contribution to genome size variation was different between TE categories, with LINE, DNA, and LTR as main contributors (table 1). We also found significant correlations between the abundances of various categories of transposable elements. The strongest, between DNA and LINEs, had a correlation coefficient of 0.45.

## 2.5 Genome size and TE content correlate with $dN/dS$

We further assessed whether variations in genome size and TE content could be explained by changes in effective population size $Ne$ using $dN/dS$ as a proxy. We used coevol to jointly assess correlations between genome size, TE content and estimates of the efficacy of selection (Lartillot and Poujol, 2011). Given a tree and an alignment, coevol models the evolution of mentioned variables assuming brownian motion parametrized by a covariance matrix. All parameters are estimated in a Bayesian framework using Monte Carlo Markov Chains. The posterior probability associated with each correlation can be interpreted as statistical support. A posterior probability close to 1 corresponds to strong statistical support for a positive correlation. Conversely, a posterior probability close to 0 corresponds to strong statistical support for a negative correlation.
<table><caption>Table 1: Pairwise correlations between genome size and TE abundances. Correlation coefficients were estimated using ancov. Cells in light gray correspond to significant correlations (posterior probability (pp) < 0.025 or pp > 0.975). Significant PGLS are presented in bold (Benjamini-Hochberg *p*- adjusted < 0.05).</caption><thead><tr><th></th><th>GS</th><th>DNA</th><th>Helitron</th><th>LINE</th><th>LTR</th><th>SINE</th><th>TEs</th></tr></thead><tbody><tr><th>GS</th><td></td><td><strong>0.55</strong></td><td>0.19</td><td><strong>0.56</strong></td><td><strong>0.38</strong></td><td>0.03</td><td><strong>0.61</strong></td></tr><tr><th>DNA</th><td>-</td><td></td><td>0.21</td><td><strong>0.45</strong></td><td><strong>0.31</strong></td><td>0.28</td><td><strong>0.58</strong></td></tr><tr><th>Helitron</th><td>-</td><td>-</td><td></td><td><strong>0.31</strong></td><td>-0.10</td><td>-0.14</td><td><strong>0.21</strong></td></tr><tr><th>LINE</th><td>-</td><td>-</td><td>-</td><td></td><td><strong>0.22</strong></td><td>0.10</td><td><strong>0.61</strong></td></tr><tr><th>LTR</th><td>-</td><td>-</td><td>-</td><td>-</td><td></td><td>0.14</td><td><strong>0.83</strong></td></tr><tr><th>SINE</th><td>-</td><td>-</td><td>-</td><td>-</td><td>-</td><td></td><td>0.17</td></tr></tbody></table>

<table><caption>Table 2: Correlation coefficients between genome size, TE abundances and *dN/dS*. Correlation coefficients were estimated using coevol. Asterisks indicate strength of support of the posterior probability to be different than 0 (** - pp>0.975 or pp<0.025 - very strong support, * - pp>0.95 or pp<0.05 - strong support).</caption><thead><tr><th></th><th></th><th>GS</th><th>DNA</th><th>Helitron</th><th>LINE</th><th>LTR</th><th>SINE</th><th>TEs</th></tr></thead><tbody><tr><td rowspan="3">dN/dS</td><td>Random genes</td><td>0.31*</td><td>0.04</td><td>0.29*</td><td>0.40**</td><td>0.34*</td><td>-0.02</td><td>0.41**</td></tr><tr><td>GC3-poor genes</td><td>0.48**</td><td>0.37**</td><td>0.37**</td><td>0.41**</td><td>0.40**</td><td>-0.18</td><td>0.50**</td></tr><tr><td>GC3-rich genes</td><td>0.25</td><td>0.20</td><td>0.25</td><td>0.18</td><td>0.35*</td><td>0.02</td><td>0.37**</td></tr></tbody></table>

In flies, G- and C-ending codons are favoured by selection on codon usage as they are less prone to translation errors (Wu et al., 2024). Consequently, G- and C-ending codons are preferentially used, with a stronger bias in highly expressed genes compared to those that are lowly expressed (Duret and Mouchiroud, 1999; Vicario et al., 2007). It is therefore expected that genes would differ in base composition at the third-codon position according to their level of expression, with highly expressed genes being enriched in GC3, which is defined as the GC content at the third codon position. In addition, a neutral force favouring G and C alleles could interfere with the intensity of selection measured by *dN/dS*, with varying intensities between genes of different recombination rates, and between species (Jackson and Charlesworth, 2021). To control for possible differences in evolutionary rates between genes with different levels of expression, we sampled genes with varying levels of GC3. *dN/dS* were estimated on alignments of three concatenates of 50 genes each: one composed of randomly sampled genes, one composed of genes with the lowest GC3, and one composed of genes with the highest GC3 (see supplementary figure 5 for GC3 distribution).

The abundance of TEs as a whole was consistently found to be positively correlated with *dN/dS* with correlation coefficients ranging from 0.37 to 0.50. Only the abundance of SINE was never found to be significantly correlated with *dN/dS*. We used the same approach to test for a possible relationship between two codon usage bias metrics (SCUO and ENC') and TE content. Our analysis reveals no correlation between SCUO nor ENC' and genome size/TE abundances (supplementary table S7).

We found a significant positive correlation between *dN/dS* and the TE content for each concatenate (table 2). The use of *dN/dS* as an estimate of *Ne* suggests a predominant role of the effective population size on TE abundance, however with different results according to the GC3. For concatenates composed of random and GC3-poor genes, we found five and six significant correlations, respectively, while only two correlations were significant for GC3-rich genes.
<figure><img src="image_3.png" /><figcaption>Figure 2: **A.** Per branch dN/dS estimated using coevol for genes in the GC3-poor and the GC3-rich concatenates. **B.** Positive correlation between gene median GC3 and gene variance in GC3 across species. Each dot represents values for a gene and the red line the linear relationship ($\rho_{spearman} = 0.255, p < 2.2 * 10^{-16}$).</figcaption></figure>

**2.6 Effect of intragenomic variation in GC3**

The stronger correlation between *dN/dS* and the abundance of TEs for GC3-poor genes compared to GC3-rich genes (table 2) raises the need to further investigate the possible interference between GC3 and *dN/dS* that may affect the strength of the correlation. Purifying selection is expected to be stronger in highly expressed, GC3-rich genes, leading to a reduced *dN/dS*. Additionally, selection on codon usage acting on these genes violates the assumed neutrality of *dS*.

While the estimated *dN/dS* is significantly higher for GC3-poor genes (mean *dN/dS* = 0.096) than for GC3-rich genes (mean *dN/dS* = 0.048, Wilcoxonranktest = 23345, *p* < 2.2*10<sup>-16</sup>, figure 2A), GC3-rich genes show greater variation in *dN/dS*. The variation in *dN/dS* is significantly higher between branches for the GC3-rich concatenate than for the GC3-poor one (0.0034 vs 0.0017, *F* = 0.51, *df* = 162, *p* = 9.34 * 10<sup>-6</sup>, figure 2A). When considering each gene individually, there is a significant positive correlation between the variation of GC3 between species and the median GC3 (*ρ*<sub>spearman</sub> = 0.255, *p* < 2.2 * 10<sup>-16</sup>, figure 2B) with a median GC3 ranging from 0.31 to 0.88). This suggests that selection on codon bias, which is stronger in highly expressed genes with higher GC3, drives differences in GC3 across species.

Beyond the evolutionary differences between genes with extreme GC3, we tested for a linear effect of gene GC3 on the strength of the correlation between *dN/dS* and the size/TE content of the genome. We repeated the coevol analysis using 26 concatenates of 50 genes each. This included 16 concatenates of homogeneous GC3 that span the GC3 distribution, alongside ten concatenates composed of genes with random GC3. The latter were added to assess the potential effects of merging genes that have varying GC3 content within the concatenates.
A significant correlation between genome size and *dN/dS* was found in five (out of ten) GC3-random concatenates and in seven (out of 16) GC3-homogeneous concatenates (supplementary figure 5). The lowest GC3 concatenate exhibited a particularly elevated correlation coefficient (*R* = 0.48), compared to the other GC3-homogeneous concatenates (*R* = 0.27 on average, with a significant *R* ranging from 0.30 to 0.35). Nevertheless, we did not find any negative correlation between the concatenate GC3 and the correlation coefficient (*t* = -1.40, *p*<sub>adj</sub> = 0.32).

We found more significant correlations with *dN/dS* when considering TE abundance as a whole than when considering genome size: five for GC3-random concatenates (out of ten) and eleven (out of 16) for GC3-homogeneous ones (figure 3). Once again, the lowest GC3-concatenate had a particularly elevated correlation coefficient (*R* = 0.50) compared to the other GC3-homogeneous (on average *R* = 0.33), and the correlation between the concatenate GC3 and the correlation coefficient was not significant (*t* = -1.51, *p*<sub>adj</sub> = 0.32). Investigating TE abundances per type, LINEs, LTR elements, and Helitrons were found to give similar results to TEs and genome size with at least eleven concatenates (out of 26) supporting a positive impact of relaxed purifying selection on abundances (figure 3). However, no concatenate supported such an impact for SINEs, and only two did for LINEs. Overall, no significant linear correlation was found between the correlation coefficients and the concatenate GC3 (*ρ*<sub>spearman</sub>, *p*<sub>adj</sub> > 0.05).

For validation, we compared the branch *dN/dS* of coevol, which evaluates *dN/dS* together with correlations throughout the tree, with estimates from MapNH, a substitution mapping algorithm (Minin and Suchard, 2008; Dutheil et al., 2012; Guéguen and Duret, 2018). Estimates from the two methods were positively and significantly correlated for each concatenate with *R*² ranging from 0.19 to 0.85 (average *R*² = 0.52, *p*<sub>adj</sub> < 0.05, supplementary figure 8).

# 3 Discussion

We performed comparative genomics on 82 Drosophilid lineages. To do so, we generated *de novo* assemblies for 40 *Drosophila* species, 14 of which did not have their genomes assembled before. Based on a clustering approach of predicted genes, we identified more than 2,200 single-copy orthologous genes. This allowed the reconstruction of a robust phylogeny that is highly congruent with previously established ones (Van Der Linde et al., 2010; Suvorov et al., 2022). In addition, we estimated genome size using flow cytometry for 23 lineages. These new genomes and genomic data provide valuable resources for species related to *D. melanogaster*, and open up possibilities for further investigations.

## 3.1 *Drosophila* phylogeny

The ML topology we report here is highly congruent with previous *Drosophila* tree reconstructions of Suvorov et al. (2022) and Van Der Linde et al. (2010), with only two incongruencies. First, the sister species of *D. simulans* is *D. sechellia* in our topology while *D. simulans* groups with *D. mauritiana* in the Suvorov tree. However, Suvorov and co-authors also identified *D. sechellia* as a sister species of *D. simulans* on their ASTRAL topology and suggested that this is the most likely due to a focus on low-recombining regions, which are less prone to incomplete lineage sorting. Second, *D. ercepeae* branches in the *montium* subgroup (*melanogaster* group) while it branches at the basis of the *ananassae* subgroup (*melanogaster* group) in the Suvorov tree, as well as in Van der Linde's tree. Because both Suvorov and the present trees are based on thousands of genes (some most likely in common), it does not seem reasonable to suggest horizontal transfer events to explain the discrepancy. This change in the place of *D. ercepeae* could be due
<figure><img src="image_4.png" /><figcaption>Figure 3: Correlation coefficients between TE type abundances and *dN/dS* as a function of GC3. The dots indicate the median GC3 for each concatenate. Correlation coefficients were estimated using coevol. Strongly and very strongly supported positive correlations are indicated in light and dark red respectively (pp>0.975 and pp>0.95).</figcaption></figure>
to contamination (mislabeled fly species for example) or a subsample effect. Indeed, our topology was compared to that of Suvorov et al. (2022) based on shared species, representing subsamples of original datasets; however each alignment trimming parameters and model of substitution are chosen for whole datasets. If subsampling plays a role, the discrepancy in the placement of *D. ercepeae* would indicate that this species has a complex evolutionary history.

## 3.2 TE content evolution mostly explains genome size variation within flies

Overall, genome size and TE content are strongly positively correlated within flies, mostly due to LINEs, DNA, and LTR elements. As LINEs and DNA elements are generally shorter than LTR elements in *Drosophila* (see Mérel et al. (2020) for a review), their predominant role in genome size evolution likely reflects greater variation in copy number. However, for LINEs and LTR elements, the correlations with genome size may also be blurred by a slight underestimation of elevated content. For species where Kim et al. (2021) found retrotransposons to be particularly abundant in long-read assemblies, our estimates were lower (for example, *D. paulistorum* with 71 Mb versus 36 Mb). This discrepancy likely reflects the difficulty of accurately assembling repeats from a subsample of short reads when they are too numerous and divergent. Importantly, it may also negatively affect the moderate relationship between genome size and TE abundance ($R \sim 0.6$).

TE content appears to be the main driver of genome composition at a rather short evolutionary scale (here ~ 60 My), as it was previously reported for Eukaryotes (Kidwell, 2002). As an exception, we estimated a rather low TE content for *D. virilis* (40Mb, i.e. an average value for a fly), while *D. virilis* has the largest genome reported so far in *Drosophila* species. As the genome of *D. virilis* consists of more than 40% of a few related 7-bp satellites (Flynn et al., 2020), the additional DNA accumulated in its genome does not correspond to TEs and is therefore not represented here (as we specifically chose to focus on TEs).

Our results show a strong effect of shared ancestry on both genome size and TE content, indicating that closely related species tend to exhibit more similar repeat content (and consequently, genome size) compared to distant species, in agreement with Sessigolo et al. (2016). Yet, in addition to strong phylogenetic inertia, we found that genome size and TE content vary within the Sophophora and Drosophila subgenera, as well as within the groups that constitute these subgenera (e.g. melanogaster or obscura groups). Such interspecific variation in genomic traits suggests that our dataset provides a valuable opportunity to test whether the variation in TE content (and hence in genome size) could be associated with differences in the efficacy of selection.

## 3.3 Which genomic trait best reflects variation in genetic drift?

We used a phylogenetic approach integrating ancestral trait reconstruction to estimate the relationship between genomic traits expected to be related to *Ne* and genome size and TE content. Using randomly sampled genes, we found a significant positive correlation between TE content and *dN/dS* ratio, while no significant association was detected between genomic composition (TE content or genome size) and codon usage bias statistics (ENC' or SCUO). The lack of correlation questions the use of selection on codon usage metrics as proxies of *Ne*. The nearly neutral theory of molecular evolution predicts that the efficacy of selection on codon usage bias depends on the effective population size of a species (Kimura, 1983; Charlesworth, 2009), suggesting that the variation in the efficacy of selection between close relatives may result from the variation in *Ne* since their divergence. In agreement with this prediction, Galtier et al.
(2018) found evidence for selection on codon usage only in large-*Ne* species of animals, but not in small-*Ne* ones. Evidence for selection on codon usage has previously been reported in *D. melanogaster* (Shields et al., 1988), but it remains unclear whether effective selection on codon usage acts with the same intensity in other fly species. Finally, interspecific differences in *Ne* and selection on codon usage may be weakly and inconsistently correlated in flies, as found in mammals (Kessler and Dean, 2014).

## 3.4 About the correlation between $dN/dS$ and TE content

Estimates of $dN/dS$ represent branch-specific estimates that reflect the efficacy of selection acting along the respective branch in the phylogeny. Length of the branches depend on sampled species relatedness, ranging from thousands to millions of years (Obbard et al., 2012). TE content represents a longer timescale, as it is a cumulative measure resulting from TE accumulation and deletion over millions of years, alongside with more recent activity observed in closely related species (Fonseca et al., 2019; Mérel et al., 2021). Nevertheless, using coevol, which models ancestral TE content and variation in both $dN/dS$ and TE content throughout the tree, we found substantial evidence for a positive correlation between the $dN/dS$ ratio and the TE content.

The signal for a positive correlation between $dN/dS$ and TE content varied according to the GC3 of the gene set. It tends to be more frequent when using concatenates of genes with similar GC3 compared to concatenates of random genes. This may be related to the fact that heterogeneity in base composition within a concatenate is expected to alter the codon frequency estimate (which will not be correct at the gene scale) and decrease the accuracy of the $dN/dS$ estimation. More importantly, correlation coefficients were higher for GC3-poor genes than for GC3-rich genes. To interpret this result, it is important to note that the GC3 of a gene is not independent of the intensity of purifying selection acting on it. Highly expressed genes are expected to be enriched in optimal codons due to selection on codon usage, and optimal codons are mainly G- or C-ending in *Drosophila* (Duret and Mouchiroud, 1999). Therefore, GC3-rich genes are expected to be more expressed and evolve under stronger selection than GC3-poor genes. While small-*Ne* species may prevent the fixation of non synonymous deleterious mutations in highly expressed genes (enriched in GC3-rich genes), they may show an accumulation of slightly deleterious mutations in less expressed genes (enriched in GC3-poor genes). In support of this hypothesis, we found lower $dN/dS$ for GC3-rich genes compared to GC3-poor genes.

Furthermore, the positive correlation found between the variance in GC3 across species and the median GC3 of the gene suggests that selection on codon bias, which is stronger in highly expressed genes (i.e. GC3-rich genes), partly drives differences in GC3 across species. Highly expressed genes evolve under selection for codon usage bias, at least in *D. melanogaster* (Duret and Mouchiroud, 1999), indicating that the hypothesis of neutrality of *dS* is likely violated in those genes (for the species experiencing selection on codon bias). Therefore, we may question the use of $dN/dS$ estimated on GC3-rich genes to provide an accurate proxy of *Ne*.

However, we did not detect a direct relationship between the concatenate GC3 and the correlation coefficient between $dN/dS$ and TE content: additional analyses based on more single-copy orthologs should confirm the link between GC3 and $dN/dS$. It would also be very powerful to compare our estimates of levels of drift with alternative proxies relying on diversity at 4-fold degenerated sites or short intron diversity. However, the acquisition of such data represents a large effort of sampling and population-level sequencing of natural populations, as they require polymorphism data for all the species of the dataset.
<table><caption>Table 3: Literature review of genome size and TE content changes resulting from transitions in life history traits associated with a reduction in *Ne*. Studies presented on lightgray rows did not find any evidence for accumulation of TEs, or increase in genome size, in lineages with reduced *Ne*. On the contrary, darker gray rows report studies where accumulation of TEs or increase in genome size has been identified in lineages with reduced *Ne*.</caption><thead><tr><th>trait <sup>1</sup></th><th>organism (sample size)</th><th>references</th><th>divergence time</th></tr></thead><tbody><tr><td>asexuality</td><td>Daphnia (n=40)</td><td>Schaack et al. (2010)<br>Jiang et al. (2017)</td><td>a few thousand years</td></tr><tr><td colspan="4">Sexual isolates exhibit more TEs than asexual ones.</td></tr><tr><td>asexuality</td><td>evening primroses (n=3/30)</td><td>Ågren et al. (2015)</td><td>5.5-500 Ky</td></tr><tr><td colspan="4">Genome size variation not explained by reproductive mode (nore by TE content, but simple repeats).</td></tr><tr><td>asexuality</td><td>arthropods (n=10)</td><td>Bast et al. (2016)</td><td>from 22 years to 10 My</td></tr><tr><td colspan="4">No evidence for increased TE load in genomes of asexual as compared to sexual lineages.</td></tr><tr><td>asexuality</td><td>yeast <sup>2</sup> (n=12)</td><td>Bast et al. (2019)</td><td>990 generations</td></tr><tr><td colspan="4">TE loads decrease rapidly under asexual reproduction.</td></tr><tr><td>asexuality</td><td>stick insects (n=10)</td><td>Jaron et al. (2022)</td><td>less than 5 My</td></tr><tr><td colspan="4">Asexual species did not show increased TE accumulation, but very low TE activity detected.</td></tr><tr><td>eusociality & parasitism</td><td>Hymenoptera (n=131)</td><td>Ardila-Garcia et al. (2010)</td><td>more than 100 My <sup>3</sup></td></tr><tr><td colspan="4">Eusocial and parasitoid species have smaller genomes than solitary and non parasitoid species.</td></tr><tr><td>eusociality</td><td>bees (n=10)</td><td>Kapheim et al. (2015)</td><td>20-80 My</td></tr><tr><td colspan="4">Reduced diversity and abundance of TEs with increasing social complexity.</td></tr><tr><td>eusociality</td><td>termites (n=22)</td><td>Koshikawa et al. (2008)</td><td>~200 My <sup>4</sup></td></tr><tr><td colspan="4">Eusocial species have smaller genomes than their noneusocial relatives.</td></tr><tr><td></td><td>animals (n=807)</td><td>Marino et al. (2024)*</td><td>~ 700 My <sup>5</sup></td></tr><tr><td colspan="4">No evidence for increased drift resulting in an accumulation of TEs across animals.</td></tr><tr><td>eusociality</td><td>termites (n=2)</td><td>Korb et al. (2015)</td><td>~150 My <sup>6</sup></td></tr><tr><td colspan="4">Increased TE abundance and genome size with social complexity.</td></tr><tr><td>eusociality</td><td>snapping shrimps (n=33)</td><td>Chak et al. (2021)</td><td>~5-7 My</td></tr><tr><td colspan="4">Eusocial species have larger genomes with more TEs than noneusocial species.</td></tr><tr><td>parasitism</td><td>Amanita fungi (n=7)</td><td>Hess et al. (2014)</td><td>see Wolfe et al. 2012</td></tr><tr><td colspan="4">Ectomycorrhizal fungi's genomes bear elevated TE contents compared to asymbiotic ones.</td></tr><tr><td>polyploidy</td><td>polyploid plants (n~300)</td><td>Baduel et al. (2019)</td><td>~60 ky</td></tr><tr><td colspan="4">TE over-accumulation after whole genome duplication as a consequence of relaxed purifying selection</td></tr><tr><td>subterranean lifestyle</td><td>Asellidae (n=22)</td><td>Lefébure et al. (2017)*</td><td>50-250 My</td></tr><tr><td colspan="4">dN/dS increases in parallel to an important increase in genome size attributed to repeated elements</td></tr><tr><td>invasiveness</td><td>D. suzukii (n=22) <sup>2</sup></td><td>Mérel et al. (2021) **</td><td>&lt;12 years</td></tr><tr><td colspan="4">Increase of the TE load in invasive populations correlated with a reduced level of genetic diversity</td></tr><tr><td></td><td>drosophila flies (n=82)</td><td>this study *</td><td>&lt;60 My</td></tr><tr><td colspan="4">TE accumulation associated with increase in dN/dS</td></tr></tbody></table>

<footnote n="1">¹ Trait associated with a reduction in *Ne*.</footnote>

<footnote n="2">² Populations of one species analyzed, results rely on polymorphism data.</footnote>

<footnote n="3">³⁻⁶ Divergence times based on ³Johnson et al. (2013); Elsik et al. (2016); Danforth et al. (2013),⁴ Danforth et al. (2013),⁵ Kumar et al. (2017)⁶ Bourguignon et al. (2014).</footnote>

<footnote>* Studies which tested correlation between genome size variation and TE content with dN/dS as estimate of *Ne*.</footnote>

<footnote>** Studies which tested correlation between genome size variation and TE content with genetic diversity *theta* as estimate of *Ne*.</footnote>
Comparative studies have previously been conducted among closely related species with different *Ne* associated with changes in various life history traits, yielding contrasting results (reported in table 3). Most studies that did not find any evidence of TE accumulation associated with a reduction in *Ne* compared asexual lines or species with sexual ones. However, if a transition to asexuality is expected to induce a drastic reduction in *Ne*, it has other evolutionary implications. In particular, Charlesworth and Langley (1986) proposed that TE transposition rates should decrease due to within lineage transmission in non-recombining genomes. In a nutshell, the reduction in transposition in asexuals would override the trend for TE accumulation due to increased genetic drift. The contrasting results among studies comparing eusocial species with non-eusocial relatives are more complex to interpret. Eusocial species are expected to harbour reduced *Ne*, and very high levels of recombination have been found in social insects (Wilfert et al., 2007). Therefore, TE accumulation is expected in eusocial species, as found in Korb et al. (2015); Chak et al. (2021), but not in Koshikawa et al. (2008); Ardila-Garcia et al. (2010); Kapheim et al. (2015). One possible explanation might lie in the differences of evolutionary scale. Indeed, compared to bees, snapping shrimps have diverged much more recently, which might have blurred the signal due to confounding factors. That could also explain why no significant correlation was found between *dN/dS* and TE content at the much divergent scales, within insects, mammals, or animals (Marino et al., 2024). An example of a factor that could prevent the observation of TE accumulation in species with small *Ne* is the evolution of different epigenetic TE-silencing mechanisms. In addition, if the synonymous divergence is too large, saturation could occur, leading to an underestimation of the *dN/dS* ratio and a less accurate proxy of *Ne*. Alternatively, according to the rate of horizontal transfers of TEs that occur in some species, it is possible that the TE content could be better explained by different invasion histories or variation in defense mechanisms than by the intensity of genetic drift. We acknowledge that the total TE content is a simple measure resulting from complex interactions of various mechanisms.

In flies, we found a significant positive correlation between the TE content and the *dN/dS* ratio, that is, an accumulation of TEs in species with lower *Ne* estimated through their efficacy of selection. A similar relationship between TE load and genomic estimates of *Ne* was found in subterranean Asellidae (Lefébure et al., 2017), and at an intraspecific level within *D. suzukii*, where invasive populations exhibited reduced genetic diversity and increased TE content (Mérel et al., 2021).

## 3.5 Conclusion

Genome size and TE abundance are positively correlated within flies. We found strong support for a relationship between these two variables and a long-term estimate of *Ne* (*dN/dS*), suggesting a substantial effect of relaxed purifying selection on TE proliferation. Support for this correlation was dependent on the set of genes used, which likely reflects different evolutionary rates between genes. In addition, our results suggest that GC heterogeneity should be controlled for to improve estimates of *dN/dS*. As our knowledge of flies progressively extends beyond *D. melanogaster*, we should soon be able to integrate new variables, such as local recombination rate or gene expression level into interspecific comparative studies, allowing for a better understanding of the complex relationship between relaxed purifying selection and molecular traits.
364 **4 Methods**

365 **4.1 Fly sampling**

366 Genome size, TE content, and the efficacy of selection were studied for 77 drosophilid species (*73 Drosophila*, three *Zaprionus*, and one *Scaptodrosophila* - supplementary table S1). DNA sequencing data used to infer TE content and the efficacy of selection were retrieved from public databases for 43 species (supplementary tables S2, S3 and S5).
367
368
369 Flow cytometry estimates were obtained from public databases for 59 species (supplementary table S4). For cases where sequencing data and/or flow cytometry estimates were not available, flies were collected from stock centers and research labs (supplementary table S8).
370
371

372 The strains were ordered to obtain lines from a single female that was sampled in the wild as recently as possible.
373 Since the study lines were derived from relatively recently sampled females, the genome fixed in these lines should be representative of a genome from a natural population. In addition, isofemale lines maintained in the lab with small population sizes should be homozygous, facilitating assembly.
374
375

376 Flies were bred in the laboratory under optimal living conditions in order to amplify them. Adults were frozen after sex determination. To sequence only homologous chromosomes, we chose to sequence female individuals (males are heterogametic). For each species, 60 females were frozen to constitute a sufficient stock of biological material. In parallel, for the estimation of genome size by flow cytometry, the heads of 50 females per species were frozen.
377
378
379

380 **4.2 Genome size estimations**

381 Genome size estimates were collected from the Animal Genome Size database (Gregory, 2005), retrieved from Hjelmen et al. (2019), or evaluated using flow cytometry (supplementary table S4). One to sixteen heads were crushed in buffer with heads from a control (either *D. melanogaster*, *D. simulans*, or *D. virilis*). After the addition of propidiumiodid, results were analyzed using the BD Accuri™ C6 Plus Flow Cytometer. Five to ten replicates were performed per fly line, unless biological material was lacking (only one estimate was possible for *D. tristis*, two for *D. algonquin*). 382
383
384
385

386 **4.3 DNA extraction and sequencing**

387 The genome of 40 lines was sequenced for assembly (referred to as newly built assemblies, supplementary table S6). To obtain a sufficient amount of genetic material for sequencing, DNA was extracted from 20 frozen females for each line using a Phenol-Chloroform protocol. DNA treated with RNAses, was checked on gel and then assayed. Sequencing was performed with TruSeq gDNA libraries (inserts of approximately 500 bp) and HiSeq Illumina paired-ends (2x125 pb) (supplementary table S2). 388
389
390
391

392 **4.4 Genome assemblies**

393 For *de novo* genomes, raw reads were filtered and trimmed for quality using UrQt (Modolo and Lerat, 2015). Two 394 alternative assembly strategies were tested. The first consists in progressive depth assemblies using IDBA-UD (Peng 395 et al., 2012), to reduce errors in high-depth regions. The second strategy consists in several steps: i) contig assembly 396 using Ray-2.3.1 (Boisvert et al., 2010), with k-mer size varying from 25 to 45, ii) quality estimation using N50, longest 397 contig size, number of contigs >100kb of the 21 obtained contig assemblies, iii) scaffolding of the three best contig
assemblies using SOAP de novo (Li et al., 2010), with k-mer size ranging from 59 to 73 (with an increment of 2), iv) selection of the best assembly based on N50, longest scaffold size, number of scaffolds >100kb and coverage portion of the genome size. For each of the 40 newly built genomes, we therefore obtained two different assemblies (one for each strategy). Genomes' completeness was estimated using Diptera single-copy orthologous gene sets from BUSCO v4 (Seppey et al., 2019): we kept the genome assembly maximizing the number of single-copy core genes found (30 were IDBA-UD and eight RAY+SOAP strategy - supplementary table S2). Publicly available genome assemblies analyzed in this study are listed in the supplementary table S3, along with the source, taxon and assembly ID.

## 4.5 Coding sequences identification & Phylogeny

In newly built and public genome assemblies, repeat sequences were identified and masked using RepeatMasker (Rep-Base database) (Smit et al., 2013). Genes were predicted using AUGUSTUS v3.1, pre-trained on the *D. melanogaster* genome (Stanke and Waack, 2003). The number of protein-coding genes predicted on each assembly (both newly built and public) is provided in the supplementary tables S2 and S3.

Prior to gene clustering, we removed isoform variants from the *D. melanogaster* annotation by keeping only the longest one. We clustered amino acid sequences using MMseqs2 (Steinegger and Söding, 2017), with default parameters. For each cluster, we built an HMM profile using hmmbuild from the HMMER package (Eddy, 2011). Then, we compared all HMM profiles for homology using the program hmmsearch from the same package. We merged homologous families with a minimum of 30% overlap using SiliX (Miele et al., 2011). For each family we removed sequences that were smaller than 50% of the median length.

For the phylogenetic reconstruction, we retained only single-copy ortholog groups composed of at least 50% of the 82 lineages, resulting in a dataset of 2,242 orthologs. Each group of orthologs was aligned at the codon level using MACSE v2.03 (Ranwez et al., 2011)(options -local_realign_init 1 -local_realign_dec 1 -max_refine_iter 25). Sites shared by fewer than 50% of the sequences were removed with Gblocks 0.91b (Castresana, 2000) (options -b1=[Number of sequences] -b2=[Number of sequences] -b3=10 -b4=2 -b5=h). We inferred a maximum likelihood phylogenetic tree from the concatenation of all amino acid alignments using IQ-Tree version 1.6.9 (Nguyen et al., 2015), for each partition the most likely model was estimated. We computed 1,000 ultrafast bootstraps (UFBoot) and 1,000 SH-like approximate likelihood ratio tests (options -bb 1000 -alrt 1000 -bnni). Additionally, we reconstructed the species tree using a supertree approach. For this purpose, all the 2,242 gene trees were individually reconstructed using IQ-Tree, with the most likely model estimated for each gene. We then used ASTRAL-III (Zhang et al., 2018) to infer the species supertree from the individual gene trees, while accounting for gene tree discordance. The ML tree was used for further comparative genomics analyses.

## 4.6 Estimations of TE content

TE content was evaluated from private and public datasets of paired-end reads using dnaPipeTE (supplementary table S5, Goubert et al. 2015). dnaPipeTE reconstructs consensus of repeated sequences by assembling a subsample of reads with coverage less than 1X. Only a subsample of reads (representing a sequencing depth < 1X) is used to ensure that only repetitive sequences are assembled. Prior to this analysis, low quality positions were removed from reads and read sizes were standardized to 100bp using fastx-toolkit (fastq质量_filter -q 20 -p 80; fastx_trimmer -l 99; v0.0.13).
Note that for *Drosophila mauritiana* read size was only 75bp. Putative bacterial reads were filtered using kraken2 with NCBI bacterial genomes (v2.1.2: Wood et al. (2019)). dnaPipeTE was applied to the processed reads as follows: -genomeCoverage 0.15, -sample_number 2 (v1.3.1, Goubert et al. (2015)).

434

435

436

437

## 4.7 Phylogenetic inertia

438

439

440

441

442

Phylogenetic inertia was assessed using Pagel's λ (phytools v1.5, Revell 2012). Because the underlying method assumes that branch lengths correspond to shared evolutionary time, the phylogenetic tree was converted to ultrametric using LSD (To et al., 2016)(-c -a 0 -z 1). LSD uses least-squares criteria to perform a relative dating of nodes, assuming that the "root date" is 0, and the "leaf date" is 1. Note that LSD constrains ancestral nodes to be older than their daughter nodes, and that species replicates were removed and values log-transformed prior to calculations.

443

## 4.8 Correlations between genome size and TE abundances

444

445

446

447

448

449

To evaluate correlations between log-transformed genome size and TE abundances while controlling for the phylogenetic signal, we used the software ancov (Lartillot 2014, see https://github.com/bayesiancook/coevol for code and documentation). Ancov is a Bayesian reimplementation of the classical method of independent contrasts, based on a Kalman filtering algorithm, which estimates correlations between quantitative traits. Additionally, we performed Phylogenetic Generalised Least Squares (PGLS) analysis using the caper package (comparative.data(vcv=TRUE), pglss(lambda = 1.0, kappa = 1.0, delta = 1.0), Orme et al. 2013).

450

## 4.9 Strength of selection

451

452

453

454

455

456

457

458

Molecular traits were estimated from concatenates of 50 genes constructed from the 2,242 single copy orthologs. Before concatenation, nucleotide alignments were filtered using HMMCleaner (Di Franco et al., 2019) and Gblocks, applying the same parameters as before but retaining only sites shared by 90% of the sequences. To account for intragenomic differences in base composition, we defined the GC3 of each gene as its median GC3 across species. We first used three concatenates: one including random genes, one including the 50 genes of lowest GC3, and one with the 50 genes of highest GC3. The analysis was then repeated using ten concatenates of random genes, and 16 concatenates of similar GC3. The latter were built to span the overall GC3 distribution without overlap between concatenates (i.e. a gene can not be within two concatenates).

459

460

461

462

463

464

465

Codon Usage Bias (CUB) was evaluated using the coRdon package (Elek et al., 2021). Metrics were first calculated per sequence, and then the median was computed for each concatenate. We chose two CUB metrics, the Effective Number of Codons that explicitly accounts for GC bias (ENC') (Novembre, 2002) and the Synonymous Codon Usage Score (SCUO) (Wan et al., 2004). ENC' is a statistic based on Pearson's $\chi^2$ statistics. It describes the departure of the observed codon usage from an expected distribution derived from background nucleotide composition. SCUO is an index of deviation from a uniform distribution based on Shannon entropy. Both metrics are not strongly correlated (Liu et al., 2018).

466

467

468

For each gene concatenate, two different estimations of $dN/dS$ were performed using coding sequence alignments. First, coevol was used to reconstruct $dN/dS$ evolution across the tree using a codon substitution model (Lartillot and Poujol, 2011) (see also Correlations between traits). Second, bio++ libraries were used to estimate a per-branch $dN/dS$
469

470

471

472

(Dutheil and Boussau, 2008; Guéguen et al., 2013). Using the previously inferred tree topology, the YN98 (F3X4) model was employed to retrieve the most likely branch lengths, codon frequencies at the root, and substitution model parameters. Then, MapNH was used to estimate *dN* and *dS* (Minin and Suchard, 2008; Dutheil et al., 2012; Guéguen and Duret, 2018).

473

## 4.10 Correlations between traits

474

475

476

477

478

479

480

Coevol was used to assess the correlations between genome size/TE abundances and estimates of the strength of selection. Coevol models *dN/dS* evolution given a tree and alignments, but also the evolution of other traits while assessing correlations between traits (Lartillot and Poujol, 2011). Variables are evolving assuming Brownian motion parametrized by a covariance matrix. All parameters are estimated in Bayesian framework using Monte Carlo Markov chains. The advantage of this approach is its integrative nature: in a single test, the correlation between all traits is evaluated across the entire phylogenetic history. We used coevol to evaluate the correlation between genome size, TE content, codon usage bias metrics, and *dN/dS*.

481

## 4.11 Statistical analysis

482

483

All descriptive and inferential statistics were performed using R (R Core Team, 2021). P-values were adjusted for multiple testing using the Benjamini-Hochberg procedure (Benjamini and Hochberg, 1995).

484

# 5 Data access

485

486

487

488

Details about the 82 drosophilid lineages are provided in supplementary table S1. The accession numbers of the public DNA sequencing data are specified supplementary table S2,3,5. Raw sequencing data has been uploaded to the European Bioinformatics Institute ENA under the study accession number PRJEB83007. All scripts are available at https://github.com/vpymerel/GND.

489

# 6 Competing interest statement

490

The authors declare no competing interests.

491

# 7 Acknowledgments

492

493

494

495

496

This work was performed using the computing facilities of the CC LBBE/PRABI. A.V.H is supported by the Agence Nationale de la Recherche, grant ANR-20-CE02-0008 / NEGA and the CNRS APEGE DynET. We are thankful to M. Bastian, T. Latrille and N. Lartillot for their helpful advice using ancov/coevol and discussions about phylogenetic models jointly estimating *Ne* and traits. We also thank Florian Bénitière, Laurent Duret, Julien Joseph, Carina Mugal, three anonymous reviewers and the editor for their valuable comments on the manuscript.
497 References

498 Ågren, J. A., Greiner, S., Johnson, M. T., and Wright, S. I. (2015). No evidence that sex and transposable elements drive genome size variation in evening primroses. *Evolution*, 69(4):1053–1062.

500 Ardila-Garcia, A., Umphrey, G., and Gregory, T. (2010). An expansion of the genome size dataset for the insect order hymenoptera, with a first test of parasitism and eusociality as possible constraints. *Insect Molecular Biology*, 19(3):337–346.

503 Baduel, P., Quadrana, L., Hunter, B., Bomblies, K., and Colot, V. (2019). Relaxed purifying selection in autopolyploids drives transposable element over-accumulation which provides variants for local adaptation. *Nature Communications*, 10(1).

506 Bast, J., Jaron, K. S., Schuseil, D., Roze, D., and Schwander, T. (2019). Asexual reproduction reduces transposable element load in experimental yeast populations. *Elife*, 8:e48548.

508 Bast, J., Schaefer, I., Schwander, T., Maraun, M., Scheu, S., and Kraaijeveld, K. (2016). No accumulation of transposable elements in asexual arthropods. *Molecular biology and evolution*, 33(3):697–706.

510 Benjamini, Y. and Hochberg, Y. (1995). Controlling the false discovery rate: A practical and powerful approach to multiple testing. *Journal of the Royal Statistical Society: Series B (Methodological)*, 57(1):289–300.

512 Biémont, C. and Vieira, C. (2006). Junk DNA as an evolutionary force. *Nature*, 443(7111):521–524.

513 Boisvert, S., Laviolette, F., and Corbeil, J. (2010). Ray: simultaneous assembly of reads from a mix of high-throughput sequencing technologies. *Journal of Computational Biology: A Journal of Computational Molecular Cell Biology*, 17(11):1519–1533.

516 Bolívar, P., Guéguen, L., Duret, L., Ellegren, H., and Mugal, C. F. (2019). GC-biased gene conversion conceals the prediction of the nearly neutral theory in avian genomes. *Genome Biology*, 20(1):5.

518 Bourguignon, T., Lo, N., Cameron, S. L., Šobotník, J., Hayashi, Y., Shigenobu, S., Watanabe, D., Roisin, Y., Miura, T., et al. (2014). The evolutionary history of termites as inferred from 66 mitochondrial genomes. *Molecular biology and evolution*, 32(2):406–421.

521 Castillo, D. M., Mell, J. C., Box, K. S., and Blumenstiel, J. P. (2011). Molecular evolution under increasing transposable element burden in Drosophila: A speed limit on the evolutionary arms race. *BMC Evolutionary Biology*, 11(1):258.

523 Castresana, J. (2000). Selection of Conserved Blocks from Multiple Alignments for Their Use in Phylogenetic Analysis. *Molecular Biology and Evolution*, 17(4):540–552.

525 Chak, S. T., Harris, S. E., Hultgren, K. M., Jeffery, N. W., and Rubenstein, D. R. (2021). Eusociality in snapping shrimps is associated with larger genomes and an accumulation of transposable elements. *Proceedings of the National Academy of Sciences*, 118(24):e2025051118.

528 Charlesworth, B. (2009). Effective population size and patterns of molecular evolution and variation. *Nature Reviews Genetics*, 10(3):195–205.
530 Charlesworth, B. and Barton, N. (2004). Genome size: does bigger mean worse? *Current Biology*, 14(6):R233–R235.

531 Charlesworth, B. and Langley, C. (1986). The evolution of self-regulated transposition of transposable elements. *Genetics*, 112(2):359–383.

533 Danforth, B. N., Cardinal, S., Praz, C., Almeida, E. A., and Michez, D. (2013). The impact of molecular data on our understanding of bee phylogeny and evolution. *Annual review of Entomology*, 58:57–78.

535 Di Franco, A., Poujol, R., Baurain, D., and Philippe, H. (2019). Evaluating the usefulness of alignment filtering methods to reduce the impact of errors on evolutionary inferences. *BMC Evolutionary Biology*, 19(1):1–17.

537 Duret, L. and Mouchiroud, D. (1999). Expression pattern and, surprisingly, gene length shape codon usage in caenorhabditis, drosophila, and arabidopsis. *Proceedings of the National Academy of Sciences*, 96(8):4482–4487.

539 Dutheil, J. and Boussau, B. (2008). Non-homogeneous models of sequence evolution in the Bio++ suite of libraries and programs. *BMC Evolutionary Biology*, 8(1):255.

541 Dutheil, J. Y., Galtier, N., Romiguier, J., Douzery, E. J. P., Ranwez, V., and Boussau, B. (2012). Efficient selection of branch-specific models of sequence evolution. *Molecular Biology and Evolution*, 29(7):1861–1874.

543 Eddy, S. R. (2011). Accelerated Profile HMM Searches. *PLOS Computational Biology*, 7(10):e1002195. Publisher: Public Library of Science.

545 Elek, A., Kuzman, M., and Vlahovicek, K. (2021). *coRdon: Codon Usage Analysis and Prediction of Gene Expressivity*. R package version 1.12.0.

547 Elliott, T. A. and Gregory, T. R. (2015). What's in a genome? the c-value enigma and the evolution of eukaryotic genome content. *Philosophical Transactions of the Royal Society B: Biological Sciences*, 370(1678):20140331.

549 Elsik, C. G., Tayal, A., Diesh, C. M., Unni, D. R., Emery, M. L., Nguyen, H. N., and Hagen, D. E. (2016). Hymenoptera genome database: integrating genome annotations in hymenopteramine. *Nucleic acids research*, 44(D1):D793–D800.

551 Felsenstein, J. (1985). Phylogenies and the comparative method. *The American Naturalist*, 125(1):1–15.

552 Fernández, P., Amice, R., Bruy, D., Christenhusz, M. J., Leitch, I. J., Leitch, A. L., Pokorny, L., Hidalgo, O., and Pellicer, J. (2024). A 160 gbp fork fern genome shatters size record for eukaryotes. *iScience*.

554 Fierst, J. L., Willis, J. H., Thomas, C. G., Wang, W., Reynolds, R. M., Ahearne, T. E., Cutter, A. D., and Phillips, P. C. (2015). Reproductive mode and the evolution of genome size and structure in caenorhabditis nematodes. *PLoS genetics*, 11(6):e1005323.

557 Flynn, J. M., Long, M., Wing, R. A., and Clark, A. G. (2020). Evolutionary dynamics of abundant 7-bp satellites in the genome of drosophila virilis. *Molecular biology and evolution*, 37(5):1362–1375.

559 Fonseca, P. M., Moura, R. D., Wallau, G. L., and Loreto, E. L. (2019). The mobilome of drosophila incompta, a flower-breeding species: comparison of transposable element landscapes among generalist and specialist flies. *Chromosome Research*, 27:203–219.
562

563

564

Galtier, N., Roux, C., Rousselle, M., Romiguier, J., Figuet, E., Glémin, S., Bierne, N., and Duret, L. (2018). Codon usage bias in animals: disentangling the effects of natural selection, effective population size, and gc-biased gene conversion. *Molecular biology and evolution*, 35(5):1092-1103.

565

566

567

568

Goubert, C., Modolo, L., Vieira, C., ValienteMoro, C., Mavingui, P., and Boulesteix, M. (2015). De Novo Assembly and Annotation of the Asian Tiger Mosquito (Aedes albopictus) Repeatome with dnaPipeTE from Raw Genomic Reads and Comparative Analysis with the Yellow Fever Mosquito (Aedes aegypti). *Genome Biology and Evolution*, 7(4):1192-1205.

569

570

Gregory, T. R. (2005). Animal genome size database. http://www.genomesize.com/search.php. Accessed: 2021-06-07.

571

572

Guéguen, L. and Duret, L. (2018). Unbiased Estimate of Synonymous and Nonsynonymous Substitution Rates with Nonstationary Base Composition. *Molecular Biology and Evolution*, 35(3):734-742.

573

574

575

Guéguen, L., Gaillard, S., Boussau, B., Gouy, M., Groussin, M., Rochette, N. C., Bigot, T., Fournier, D., Pouyet, F., Cahais, V., et al. (2013). Bio++: efficient extensible libraries and tools for computational molecular evolution. *Molecular Biology and Evolution*, 30(8):1745-1750.

576

577

Haddrill, P. R., Loewe, L., and Charlesworth, B. (2010). Estimating the parameters of selection on nonsynonymous mutations in drosophila pseudoobscura and d. miranda. *Genetics*, 185(4):1381-1396.

578

579

Hess, J., Skrede, I., Wolfe, B. E., LaButti, K., Ohm, R. A., Grigoriev, I. V., and Pringle, A. (2014). Transposable element dynamics among asymbiotic and ectomycorrhizal amanita fungi. *Genome Biology and Evolution*, 6(7):1564-1578.

580

581

582

Hjelmen, C. E., Blackmon, H., Holmes, V. R., Burrus, C. G., and Johnston, J. S. (2019). Genome Size Evolution Differs Between Drosophila Subgenera with Striking Differences in Male and Female Genome Size in Sophophora. *G3: Genes, Genomes, Genetics*, 9(10):3167-3179. Publisher: G3: Genes, Genomes, Genetics Section: Investigations.

583

584

Jackson, B. and Charlesworth, B. (2021). Evidence for a force favoring gc over at at short intronic sites in drosophila simulans and drosophila melanogaster. *G3*, 11(9):jkab240.

585

586

587

Jaron, K. S., Parker, D. J., Anselmetti, Y., Tran Van, P., Bast, J., Dumas, Z., Figuet, E., François, C. M., Hayward, K., Rossier, V., et al. (2022). Convergent consequences of parthenogenesis on stick insect genomes. *Science advances*, 8(8):eabg3842.

588

589

Jiang, X., Tang, H., Ye, Z., and Lynch, M. (2017). Insertion polymorphisms of mobile genetic elements in sexual and asexual populations of daphnia pulex. *Genome biology and evolution*, 9(2):362-374.

590

591

Johnson, B. R., Borowiec, M. L., Chiu, J. C., Lee, E. K., Atallah, J., and Ward, P. S. (2013). Phylogenomics resolves evolutionary relationships among ants, bees, and wasps. *Current Biology*, 23(20):2058-2062.

592

593

594

Kapheim, K. M., Pan, H., Li, C., Salzberg, S. L., Puiu, D., Magoc, T., Robertson, H. M., Hudson, M. E., Venkat, A., Fischman, B. J., et al. (2015). Genomic signatures of evolutionary transitions from solitary to group living. *Science*, 348(6239):1139-1143.
595

596

Kapitonov, V. V. and Jurka, J. (2008). A universal classification of eukaryotic transposable elements implemented in Repbase. *Nature Reviews. Genetics*, 9(5):411–412; author reply 414.

597

598

Kessler, M. D. and Dean, M. D. (2014). Effective population size does not predict codon usage bias in mammals. *Ecology and evolution*, 4(20):3887–3900.

599

Kidwell, M. G. (2002). Transposable elements and the evolution of genome size in eukaryotes. *Genetica*, 115:49–63.

600

601

602

Kim, B. Y., Wang, J. R., Miller, D. E., Barmina, O., Delaney, E., Thompson, A., Comeault, A. A., Peede, D., D'Agostino, E. R., Pelaez, J., et al. (2021). Highly contiguous assemblies of 101 drosophilid genomes. *eLife*, 10:e66405. Publisher: eLife Sciences Publications, Ltd.

603

Kimura, M. (1983). *The neutral theory of molecular evolution*. Cambridge University Press.

604

605

Korb, J., Poulsen, M., Hu, H., Li, C., Boomsma, J. J., Zhang, G., and Liebig, J. (2015). A genomic comparison of two termites with different social complexity. *Frontiers in Genetics*, 6:9.

606

607

Koshikawa, S., Miyazaki, S., Cornette, R., Matsumoto, T., and Miura, T. (2008). Genome size of termites (insecta, dictyoptera, isoptera) and wood roaches (insecta, dictyoptera, cryptocercidae). *Naturwissenschaften*, 95:859–867.

608

609

Kumar, S., Stecher, G., Suleski, M., and Hedges, S. B. (2017). Timetree: a resource for timelines, timetrees, and divergence times. *Molecular biology and evolution*, 34(7):1812–1819.

610

611

Lartillot, N. (2014). A phylogenetic Kalman filter for ancestral trait reconstruction using molecular data. *Bioinformatics*, 30(4):488–496.

612

613

Lartillot, N. and Delsuc, F. (2012). Joint reconstruction of divergence times and life-history evolution in placental mammals using a phylogenetic covariance model. *Evolution*, 66(6):1773–1787.

614

615

Lartillot, N. and Poujol, R. (2011). A phylogenetic model for investigating correlated evolution of substitution rates and continuous phenotypic characters. *Molecular Biology and Evolution*, 28(1):729–744.

616

617

618

Lefébure, T., Morvan, C., Malard, F., François, C., Konecny-Dupré, L., Guéguen, L., Weiss-Gayet, M., Seguin-Orlando, A., Ermini, L., Sarkissian, C. D., et al. (2017). Less effective selection leads to larger genomes. *Genome Research*, 27(6):1016–1028.

619

620

Li, R., Zhu, H., Ruan, J., Qian, W., Fang, X., Shi, Z., Li, Y., Li, S., Shan, G., Kristiansen, K., et al. (2010). De novo assembly of human genomes with massively parallel short read sequencing. *Genome Research*, 20(2):265–272.

621

622

Liu, S. S., Hockenberry, A. J., Jewett, M. C., and Amaral, L. A. N. (2018). A novel framework for evaluating the performance of codon usage bias metrics. *Journal of the Royal Society Interface*, 15(138):20170667.

623

Lynch, M. and Conery, J. S. (2003). The Origins of Genome Complexity. *Science*, 302(5649):1401–1404.

624

625

Marino, A., Debaecker, G., Fiston-Lavier, A.-S., Haudry, A., and Nabholz, B. (2024). Effective population size does not explain long-term variation in genome size and transposable element content in animals. *eLife*.

626

627

Martins, E. P. and Hansen, T. F. (1997). Phylogenies and the comparative method: a general approach to incorporating phylogenetic information into the analysis of interspecific data. *The American Naturalist*, 149(4):646–667.
628

629

McVean, G. A. and Charlesworth, B. (1999). A population genetic model for the evolution of synonymous codon usage: patterns and predictions. *Genetics Research*, 74(2):145-158.

630

631

Mérel, V., Boulesteix, M., Fablet, M., and Vieira, C. (2020). Transposable elements in drosophila. *Mobile DNA*, 11:1-20.

632

633

634

Mérel, V., Gibert, P., Buch, I., Rodriguez Rada, V., Estoup, A., Gautier, M., Fablet, M., Boulesteix, M., and Vieira, C. (2021). The worldwide invasion of drosophila suzukii is accompanied by a large increase of transposable element load and a small number of putatively adaptive insertions. *Molecular biology and evolution*, 38(10):4252-4267.

635

636

Miele, V., Penel, S., and Duret, L. (2011). Ultra-fast sequence clustering from similarity networks with SiLiX. *BMC Bioinformatics*, 12(1):116.

637

638

Minin, V. N. and Suchard, M. A. (2008). Fast, accurate and simulation-free stochastic mapping. *Philosophical Transactions of the Royal Society B: Biological Sciences*, 363(1512):3985-3995.

639

640

Modolo, L. and Lerat, E. (2015). UrQt: an efficient software for the Unsupervised Quality trimming of NGS data. *BMC Bioinformatics*, 16(1):137.

641

642

Nguyen, L.-T., Schmidt, H. A., von Haeseler, A., and Minh, B. Q. (2015). IQ-TREE: a fast and effective stochastic algorithm for estimating maximum-likelihood phylogenies. *Molecular Biology and Evolution*, 32(1):268-274.

643

644

Novembre, J. A. (2002). Accounting for background nucleotide composition when measuring codon usage bias. *Molecular Biology and Evolution*, 19(8):1390-1394.

645

646

647

Obbard, D. J., Maclennan, J., Kim, K.-W., Rambaut, A., O'Grady, P. M., and Jiggins, F. M. (2012). Estimating Divergence Dates and Substitution Rates in the Drosophila Phylogeny. *Molecular Biology and Evolution*, 29(11):3459-3473.

648

649

Ohta, T. (1992). The nearly neutral theory of molecular evolution. *Annual review of ecology and systematics*, pages 263-286.

650

Orgel, L. E. and Crick, F. H. C. (1980). Selfish DNA: the ultimate parasite. *Nature*, 284(5757):604-607.

651

652

Orme, D., Freckleton, R., Thomas, G., Petzoldt, T., Fritz, S., Isaac, N., and Pearse, W. (2013). The caper package: comparative analysis of phylogenetics and evolution in r. *R package version*, 5(2):1-36.

653

654

Pagel, M. D. and Harvey, P. H. (1988). Recent developments in the analysis of comparative data. *The Quarterly Review of Biology*, 63(4):413-440.

655

Parvathy, S. T., Udayasuriyan, V., and Bhadana, V. (2022). Codon usage bias. *Molecular biology reports*, 49(1):539-565.

656

657

Peng, Y., Leung, H. C. M., Yiu, S. M., and Chin, F. Y. L. (2012). IDBA-UD: a de novo assembler for single-cell and metagenomic sequencing data with highly uneven depth. *Bioinformatics*, 28(11):1420-1428.

658

659

660

Pollak, E. (1987). On the theory of partially inbreeding finite populations. i. partial selfing. *Genetics*, 117(2):353-360.

R Core Team (2021). *R: A Language and Environment for Statistical Computing*. R Foundation for Statistical Computing, Vienna, Austria.
661 Ranwez, V., Harispe, S., Delsuc, F., and Douzery, E. J. P. (2011). MACSE: Multiple Alignment of Coding SEquences
662 Accounting for Frameshifts and Stop Codons. *PLOS ONE*, 6(9):e22594. Publisher: Public Library of Science.

663 Revell, L. J. (2012). phytools: An R package for phylogenetic comparative biology (and other things). *Methods in Ecology and Evolution*, 3:217-223.

665 Robinson, D. F. and Foulds, L. R. (1981). Comparison of phylogenetic trees. *Mathematical Biosciences*, 53(1):131-147.

666 Schaack, S., Pritham, E. J., Wolf, A., and Lynch, M. (2010). Dna transposon dynamics in populations of daphnia pulex
667 with and without sex. *Proceedings of the Royal Society B: Biological Sciences*, 277(1692):2381-2387.

668 Schoch, C. L., Ciufo, S., Domrachev, M., Hotton, C. L., Kannan, S., Khovanskaya, R., Leipe, D., Mcveigh, R.,
669 O'Neill, K., Robbertse, B., et al. (2020). NCBI Taxonomy: a comprehensive update on curation, resources and tools.
670 *Database: The Journal of Biological Databases and Curation*, 2020:baaa062.

671 Seppey, M., Manni, M., and Zdobnov, E. M. (2019). BUSCO: Assessing Genome Assembly and Annotation Complete-
672 ness. *Methods in Molecular Biology* (Clifton, N.J.), 1962:227-245.

673 Sessegolo, C., Burlet, N., and Haudry, A. (2016). Strong phylogenetic inertia on genome size and transposable element
674 content among 26 species of flies. *Biology Letters*, 12(8):20160407.

675 Shields, D. C., Sharp, P. M., Higgins, D. G., and Wright, F. (1988). " silent" sites in drosophila genes are not neutral:
676 evidence of selection among synonymous codons. *Molecular biology and evolution*, 5(6):704-716.

677 Smit, A., Hubley, R., and Green, P. (2013). RepeatMasker. http://www.repeatmasker.org/. Accessed: 2021-06-
678 07.

679 Stanke, M. and Waack, S. (2003). Gene prediction with a hidden Markov model and a new intron submodel. *Bioinfor-
680 matics*, 19(suppl_2):ii215-ii225.

681 Steinegger, M. and Söding, J. (2017). MMseqs2 enables sensitive protein sequence searching for the analysis of massive
682 data sets. *Nature Biotechnology*, 35(11):1026-1028.

683 Suvorov, A., Kim, B. Y., Wang, J., Armstrong, E. E., Peede, D., D'Agostino, E. R. R., Price, D. K., Waddell, P., Lang,
684 M., Courtier-Orgogozo, V., et al. (2022). Widespread introgression across a phylogeny of 155 drosophila genomes.
685 32(1):111-123.e5.

686 Tenaillon, M. I., Hollister, J. D., and Gaut, B. S. (2010). A triptych of the evolution of plant transposable elements.
687 *Trends in Plant Science*, 15(8):471-478.

688 To, T.-H., Jung, M., Lycett, S., and Gascuel, O. (2016). Fast dating using least-squares criteria and algorithms.
689 *Systematic Biology*, 65(1):82-97.

690 Van Der Linde, K., Houle, D., Spicer, G. S., and Steppan, S. J. (2010). A supermatrix-based molecular phylogeny of
691 the family drosophilidae. *Genetics research*, 92(1):25-38.

692 Vicario, S., Moriyama, E. N., and Powell, J. R. (2007). Codon usage in twelve species of drosophila. 7(1):226.
693

694

Wan, X.-F., Xu, D., Kleinhofs, A., and Zhou, J. (2004). Quantitative relationship between synonymous codon usage bias and GC composition across unicellular genomes. *BMC Evolutionary Biology*, 4:19.

695

696

Whitney, K. D., Garland, T., and Jr (2010). Did Genetic Drift Drive Increases in Genome Complexity? *PLoS Genetics*, 6(8).

697

698

Wilfert, L., Gadau, J., and Schmid-Hempel, P. (2007). Variation in genomic recombination rates among animal taxa and the case of social insects. *Heredity*, 98(4):189-197.

699

700

Wolfe, B. E., Tulloss, R. E., and Pringle, A. (2012). The irreversible loss of a decomposition pathway marks the single origin of an ectomycorrhizal symbiosis. *PLOS ONE*, 7(7):e39597.

701

702

Wood, D. E., Lu, J., and Langmead, B. (2019). Improved metagenomic analysis with Kraken 2. *Genome Biology*, 30(1):257.

703

704

Wu, X., Xu, M., Yang, J.-R., and Lu, J. (2024). Genome-wide impact of codon usage bias on translation optimization in drosophila melanogaster. *Nature Communications*, 15(1):8329.

705

706

Zhang, C., Rabiee, M., Sayyari, E., and Mirarab, S. (2018). ASTRAL-III: polynomial time species tree reconstruction from partially resolved gene trees. *BMC Bioinformatics*, 19(6):153.