Patterns and determinants of mitogenomic evolution in Bilateria

Domain: biology
Published: 2026-03-12
Editorial outcome: minor_revision

Abstract
Mitochondrial genomes of some animals remained almost unchanged architecturally for over half a billion of years, while other lineages exhibit conspecific variability. The drivers underlying this disparity remain unknown. We analyse over 10,000 bilaterian mitogenomes to test whether locomotory capacity and parasitic lifestyle can explain the variability in architectural evolutionary rates, and attempt to answer multiple other evolutionary questions. Double-stranded (genes on both strands) architecture is the most likely ancestral state for most major radiations (phylum and above). We identify over twenty transitions to single-stranded architectures, and the evidence of state-reversals (in Annelida and Mollusca). Gene order rearrangement and sequence evolution rates are positively correlated (r=0.69), and both are increased in single-stranded mitogenomes, parasites, and species with low locomotory capacity. The latter two categories exhibit increased prevalence of single-strandedness and GC skew inversions. Mitogenome size is negatively correlated with evolutionary rates, whereas the effective population size is decoupled from all tested evolutionary variables. Endotherms exhibit slower evolutionary rates than ectotherms in Bilateria, but faster in Chordata. Here, we provide evidence that purifying selection mediated by locomotory capacity and species ecology is a major driver of mitogenomic architectural evolution, but other variables are required to fully explain the observed patterns.

Introduction
Mitogenomic architectural rearrangement rates, i.e. changes in the order and strand distribution of genes encoded by the animal mitochondrial genomes, vary dramatically among bilaterian lineages, but factors underlying this disparity in architectural evolutionary rates remain poorly understood 1–3. On the conserved end of the spectrum, most species of Chordata exhibit a mitogenomic architecture that remained almost unchanged over half a billion years 1,4,5. In contrast to this, some lineages of Mollusca, Isopoda, and Nematoda exhibit within-family/genus architectural variability 2,6–9, and nematodes from the family Mermithidae even exhibit conspecific architectural variability 7. Such stasis in the architectural evolution over long time spans is an indicator that architectural rearrangements are commonly deleterious, and likely disrupting mitogenomic replication, transcription, and gene expression co-regulation 1–3. However, no known models explain this apparent disparity in the stringency of selective constraints among lineages.

Recently, we have shown that parasitism and locomotory capacity have the power to explain a major proportion of variability in mitogenomic sequence evolutionary rates across the Bilateria 10. We proposed that reduction in energy demands allows reduced (ATP) production, which results in relaxation of purifying selection pressures – i.e. increased accumulation of deleterious mutations and accelerated evolutionary rates. Reduction in locomotory capacity and high metabolic dependence on the host (“auxotrophy”) characterise many parasitic animals, resulting in disproportionately elevated evolutionary rates, particularly in endoparasites 10. Based on that finding, and the fact that isolated parasitic lineages exhibit highly unusual, “disrupted”, rapidly evolving mitogenomic architectures 7–9,11–13, we hypothesised that parasitism and locomotory capacity variability may also have the power to explain the variability in architectural evolution.

As this hypothesis posits purifying selection pressure as the central mediator, it implies that architectural rearrangement rates (commonly quantified as gene order rearrangement rates - GORR) and sequence evolution (quantified as branch length) should exhibit a positive correlation. While some studies found support for this hypothesis 9,14–18, others found no correlation in isolated lineages 1,19,20. Some studies proposed that GORR and branch length are correlated only in lineages with high rearrangement rates, and not in those with low rearrangement rates 1,14, but this hypothesis remains untested statistically.

As regards parasitism, studies of its putative association with mitogenomic architecture evolution are scarce and results inconsistent: no correlation in Diptera (Insecta) 21, inconsistent results in Hymenoptera (Insecta) 22,23, no correlation in Annelida 18, and a recent study of Platyhelminthes found decreased GORR rates in parasites 20. Aside from a few sporadic mentions of this putative association 1,24–26, no other study had attempted to test it systematically. Our own studies of mitogenomes of parasitic animals (Supplementary Table 1) indicate that isolated lineages of nematodes and crustaceans do exhibit rapidly evolving architectures 8,9,11,27, but parasitic flatworms (Neodermata) exhibit comparatively slow-evolving architectures 28. As regards the hypothesis that species with low locomotory capacity should exhibit increased GORR rates, to our knowledge, it remains untested.

Several other questions regarding the correlates and patterns of mitogenomic evolution remain open. The observation that some species exhibit disproportionately rapid architectural evolutionary rates has produced a hypothesis that mitogenomic architecture evolution is “discontinuous”; i.e. that once a rearrangement event has destabilized the mitogenomic architecture, this is often followed by a cascade of rearrangements 11. Conversely, several comparative analyses suggested that architectural evolution is often re-stabilised after a burst of rearrangements 9,18,20, indicating an “episodic” evolutionary pattern, where purifying selection pressure is restored. Overall, the dynamics of architectural evolutionary rates remains poorly understood.

Further, mitogenomes can have genes encoded on both strands ("double-stranded") or on only one strand ("single-stranded")2. Boore4 proposed that evolutionary transitions from double-stranded to single-stranded architectures may be irreversible if they are accompanied by the loss of regulatory elements on the noncoding strand (“ratchet mechanism”). The evolutionary history of mitogenomic strandedness has never been systematically studied in animals, so the total number of state switches remains unknown, and the “ratchet” hypothesis untested. Even the ancestral state for bilaterians is debated, with some studies proposing a single-stranded ancestral architecture 2,29, whereas both studies that attempted to systematically infer the ancestral bilaterian (ur-bilaterian) architecture settled on a vertebrate-like double-stranded architecture 5,30.

A correlation between single-strandedness and highly rearranged mitogenomic architectures was proposed 2, but not tested statistically. We also noticed that some of the fastest mitogenomic DNA sequence evolutionary rates are exhibited by parasitic lineages that have exclusively or predominantly single-stranded architectures, such as Neodermata (Platyhelminthes), Chromadorea (Nematoda), and Acanthocephala 10. This is not a rule, as parasitic Orthonectida and Enoplea (Nematoda) exhibit both highly elevated evolutionary rates and double-stranded architectures 2,10,31. We hypothesise that single-strandedness may be an accidental consequence of relaxed purifying selection pressures allowing a wider range of rearrangement types, thus increasing the likelihood of a rearrangement producing a single-stranded architecture, which may then become fixed via the ratchet mechanism. This implies that single-stranded lineages should exhibit elevated sequence and architectural evolutionary rates, and be more common in parasitic and low-locomotory species.

Compositional asymmetry between strands (GC/AT skew) is linked to mitogenomic architecture via asymmetrical mutational pressures 32. Metabolic rate should be positively associated with GC skew magnitude via mutagenic reactive oxygen species 33, but crustaceans exhibited a negative correlation between these two variables, mediated by purifying selection driven by locomotory capacity variability 34. This suggests the existence of a positive correlation between GC skew and both GORR and branch length, as well as increased average GC skew values in single-stranded mitogenomes, parasites, and low-locomotory lineages. Aside from the association between GC skew and GORR in Annelida 18, and elevated GC skews in parasitic crustaceans 34, the universality of these hypotheses remains untested. Certain types of architectural rearrangements may change the direction of replication and transcription, thus changing the direction of mutational pressures, and eventually cause skew inversions 35. These have been identified in several phyla, and associated with increased branch lengths 1,35–37 and increased dN/dS ratios in crustaceans 36. Along with other arguments outlined herein, this suggests that skew inversions should be more likely in lineages evolving under relaxed selection pressures, i.e. parasitic and low-locomotory species.

Mitogenomic size is another variable that exhibits pronounced variability among bilaterian species for unknown reasons. Most animal mitogenomes evolve under a strong selection for small size, putatively driven by the “race for replication”, wherein smaller mitogenomes are replicated at higher rates due to their shorter replication time 38,39. This selection pressure should be more pronounced in species with high metabolic rates due to their faster mitogenomic metabolism and higher mitogenomic replication rates 39, implying that mitogenomes should be larger in species evolving under relaxed purifying selection pressures, and mitogenome size thus positively correlated with GORR, branch length and GC skew. Indirect support for the negative correlation between size and purifying selection was observed in some geckos and frogs 40. However, increased metabolic rates were not associated with smaller mitogenomes in selected tetrapod lineages 41. Aside from the supporting evidence for positive correlation between GORR and mitogenomic size in flatworms 20, correlates of mitogenomic size have not been tested statistically on a comprehensive dataset.

The above hypothesis also implies that mitogenomic size should be affected by the thermoregulation type: endotherms should have smaller mitogenomes and reduced size variability compared to ectotherms via their higher energy expenditure driving higher metabolic rates and producing a stronger race for replication 39,42. Support for this hypothesis was found on a dataset comprising mostly vertebrates and a small number of invertebrates 39. Further, two different models proposed that ectotherms should exhibit lower sequence mutation rates: the “functional constraints” hypothesis proposes that variations in the thermic environment may restrict physiologically acceptable amino acid substitutions 42,43, and the ”metabolic rate” hypothesis links the reduction in energy expenditure in ectotherms to a reduction in metabolic rates, production of mutagenic reactive oxygen species, and finally DNA mutation rates 33,42. Indeed, there is strong evidence that ectotherms exhibit significantly lower evolutionary rates than endotherms among vertebrates 33,42–44, and there is some evidence that the effect of thermoregulatory type remains consistent when invertebrates are included 45. However, this contradicts our central premise that high energy expenditure (as in endotherms) strengthens purifying selection. Indeed, vertebrates (which comprise all endothermic animals) exhibit some of the shortest branches among all bilaterians 10. Therefore, our working hypotheses contradict the proposals of reduced evolutionary rates in ectotherms, and instead propose reduced rates in endotherms, expected to be reflected in lower branch lengths, GC skew, and GORR.

Finally, a reduction in the effective population size (Ne) weakens the purifying selection 46–48, so we hypothesised that lower Ne should be reflected in higher average branch lengths, GORR values, GC skews, and mitogenome sizes 3,10,34,47. However, previous studies produced contradictory results 34,40,49–54, so the effects of Ne on mitogenomic evolution remain debated.

Herein, we rely on a mitogenomic dataset comprising almost 11,000 bilaterian species, categorised by locomotory capacity and lifestyle and assembled for our recent study 10, to test the hypotheses concerning the correlates of mitogenomic evolution outlined in this section, and provide evidence that purifying selection mediated by locomotory capacity and species ecology is a major driver of mitogenomic architectural evolution. However, analyses indicate that other variables must be invoked to fully explain observed patterns.

Results
To test hypotheses aimed at answering open questions about the evolution of mitogenomic architectures and sequences (Table 1 and Supplementary Note 1), we analysed 10,860 bilaterian mitogenomes, among which 5.2% were single-stranded, 3.6% belonged to macroparasitic species (ectoparasitic+endoparasitic), and 11.1% to species assigned to the low locomotory capacity category (Figure 1; Supplementary Table 2; Supplementary Data 1).

Ancestral mitogenomic architectures and strand distribution of genes

Hypothesis 1: On the basis of our own reconstruction of ancestral mitogenomic architectures and previous findings 5,30, the architecture identical to the one exhibited by most chordates, including humans (only PCGs and rRNAs were considered), was selected as the most likely ancestral architecture of all bilaterian animals (hypothesis 1a; Figure 1). This also implies that double-stranded architecture and negative GC skew on the plus strand are the ancestral bilaterian states. There was some instability in the exact position of nad6 across different analyses, but its location on the minus strand was inferred with high confidence (more details in Supplementary Note 2 and Supplementary Figs. 1-2). Ancestral double-strandedness was further corroborated by the Maximum Likelihood ancestral state reconstruction for the two architecture types. This analysis included tRNA genes (eight on the minus strand in the human mitogenome), and resolved the ancestral state for Bilateria as double-stranded with a probability of 0.9999 (hypothesis 1b).

The above ur-bilaterian architecture was used as the reference for subsequent gene order rearrangement analyses, and to infer the evolutionary history of mitogenomic strandedness shifts (double-stranded / single-stranded).

Hypothesis 2: a double-stranded architecture was inferred as the most likely ancestral state for Deuterostomia/Xenambulacraria and all belonging phyla with high confidence. We inferred one architectural evolutionary transition from genes encoded on both strands to all genes encoded on a single strand (henceforth “double-to-single”) in Xenacoelomorpha and at least one in Tunicata (Chordata), likely in the common ancestor of classes Thaliacea+Ascidiacea, but not Appendicularia. Details of these analyses are provided in Table 2 and Supplementary Note 3 and Supplementary Figs. 2-38. BayesTraits reconstruction produced double-stranded as the ancestral state for the superclade Protostomia and both of its superphyla, Ecdysozoa and Lophotrochozoa, with high probability (>0.999). Within the Ecdysozoa, double-stranded was inferred as the ancestral state for six out of seven phyla included in the dataset (Kinorhyncha was the outlier), and we inferred five putative double-to-single transitions. The superphylum Lophotrochozoa exhibited by far the most rearranged mitogenomic architectures and the most complex scenario of strandedness transitions. Likely due to the highly rearranged mitogenomic architectures of most lineages, we failed to reconstruct of the ancestral architecture for this superphylum with precision, but multiple sources of evidence indicate that a double-stranded architecture was the ancestral state. These comprise most results of the ancestral architecture reconstruction (Figure S8), BayesTraits reconstruction (probability of double-stranded > 0.999), and the fact that double-stranded was the most likely ancestral state for 11 out of 13 individual phyla, with Brachiopoda and Syndermata as the only exceptions. Syndermata comprise sensu-stricto phyla Acanthocephala (all single-stranded) and Rotifera (mostly single-stranded) (Figures S9-S12). We found indications of single-stranded to double-stranded transitions in several species of Rotifera, but due to their fragmented mitogenomes, within-species strandedness variability, and strand switches comprising only tRNA genes, these are not confirmed. Another double-to-single transition was inferred in Platyhelminthes after Catenulida branched-off. Three lophotrochozoan phyla exhibited particularly complex scenarios: Bryozoa, Annelida, and Mollusca. In Annelida we identified a double-to-single evolutionary shift after Palaeoannelida branched-off, followed by a state reversal (i.e. single-to-double shift) in a single derived Sedentaria species, Sabella spalanzani (Sabellidae). Finally, Mollusca exhibited the most complex scenario, with multiple putative shifts in strandedness, comprising several putative state reversals observed in the class Bivalvia (Table 2).

GORR values distribution

Hypothesis 3: GORR values in comparison to the ancestral bilaterian architecture ranged between 0 and 12 in the complete dataset (Figures 1 and 2). The distribution of GORR values across individual phyla could be roughly classified into three basic, albeit not necessarily mutually exclusive, patterns of architectural evolution. (1) Discontinuous: in most phyla, the ancestral architecture (disregarding volatile tRNA genes) for the phylum is highly conserved in a majority of species, with exponentially elevated rearrangement rates in a small proportion of species (e.g. Arthropoda, Chordata, Nematoda, and Echinodermata). (2) Episodic: mitogenomic architecture is re-stabilised after an evolutionary period of increased levels of architectural rearrangements, commonly occurring deep in the evolutionary history (e.g. Nematoda and Arthropoda). (3) Continuous: GORR values exhibit an approximately linear distribution, characterised by the absence of statistical outliers. Most phyla exhibited a combination of episodic and discontinuous patterns, whereas Mollusca and Bryozoa were the only phyla exhibiting a continuous pattern. In addition, Platyhelminthes exhibited a unique pattern, with parasitic Neodermata exhibiting a strongly episodic pattern (exceptionally highly rearranged, but mostly stable architecture: GORR=9), and predominantly free-living “turbellaria” exhibiting much less rearranged, but fast-evolving in an apparently continuous manner, pattern of GORR distribution (Figure 2 and Supplementary Fig. 39).

GORR values, branch length, and strand distribution of genes

Hypothesis 4: Mitogenomic architectural evolutionary rates (GORR) strongly varied among the superphyla: highly conserved in Deuterostomia (0.24), intermediate in Ecdysozoa (3.29), and high in Lophotrochozoa (6.86). In support of the working hypothesis. branch lengths followed this pattern, but with less pronounced differences (Supplementary Fig. 40). Congruent patterns across the two variables were also observed across the individual phyla, with the slowest evolution in Chordata, Hemichordata, and Echinodermata, and fastest in Platyhelminthes, Orthonectida, and Mollusca/Acanthocephala (GORR/branch length respectively) (Figure 1 and 2, Supplementary Table 2). As a result, the correlation between GORR and branch length was strong and positive in the overall dataset (hypothesis 4a; r=0.69; p<0.001; Figure 3A), and branch length variability explained 37 - 40% of the GORR variability in Bilateria (Table 3). As Chordata exhibited a unique evolutionary pattern, and represented a majority of the dataset, major hypotheses were additionally tested with Chordata removed from the dataset: the overall conclusions remained the same, although the correlation value between GORR and branch length was slightly reduced (0.57). Correlation between GORR and branch length was also positive (albeit mostly weak) in all individual phyla aside from Mollusca (-0.08) (Figure 3), as well as in mitogenomes with low rearrangement rates (GORR<7; r=0.69, R2=42-45%; rejecting the hypothesis 4b), and in the dataset split along the lines of architectural type: single-stranded (r=0.86, R2=66-70%) and double-stranded (r=0.70, R2=47-49%). Contrary to this, highly rearranged mitogenomes (GORR≥7) exhibited a very weak negative correlation (-0.18) between GORR and branch length, reflected in low explanatory power (R2=3%), but the pattern was restored when Mollusca were removed from this dataset (r=0.71, R2=51%; further details in Table 3; Supplementary Note 4, Figs. 41 - 44, and Data 2).

Hypothesis 5: Next, we tested the hypothesis that mitogenomes with all genes encoded on the same strand are exhibit signatures of relaxed purifying selection pressures compared to double-stranded mitogenomes. In agreement with the working hypothesis 5a, GORR values were significantly higher in single-stranded mitogenomes across the Bilateria (Figure 4A; p<0.001, details in Supplementary Data 3). At the level of individual phyla, results were inconsistent, with Chordata and Arthropoda exhibiting a congruent pattern, but Annelida and Nematoda exhibiting an inverted GORR pattern (Figure 4A). In support of the hypothesis 5b, single-stranded mitogenomes exhibited significantly higher sequence evolutionary rates across the Bilateria (Figure 4B; p<0.001). At the level of individual phyla, results were only partially consistent, with single-stranded lineages exhibiting significantly longer branches in Chordata and Mollusca, nonsignificantly longer in Annelida and Arthropoda, and Nematoda exhibiting an inverted pattern (Figure 4B). Both hypotheses were also confirmed when Chordata were removed from the dataset (Supplementary Fig. 45).

As regards GC skews, mitogenomes with inverted (positive) values were observed in almost all phyla, aside from three phyla represented by only a few specimens: Entoprocta (n=2), Hemichordata (n=6), and Nematomorpha (n=1) (Figure 4C and Supplementary Fig. 46). GC skews were predominantly negative (the ancestral state) in double-stranded mitogenomes, and predominantly positive in single-stranded mitogenomes, but there were multiple exceptions from this rule, with Annelida as a major outlier (+ in ds, - in ss). In agreement with the hypothesis 5c, the absolute skew magnitude was significantly higher in single-stranded than in double-stranded mitogenomes (confirmed in both Bilateria - Figure 5A, and Chordata - Supplementary Fig. 47).

GC skews and Ne

Hypothesis 6: contrary to the working hypothesis 6a, GC skew (│value│>0.1) was negatively weakly correlated with GORR (r=-0.25), with a negligible explanatory power (R2≈3%). However, the working hypothesis was weakly-moderately supported when we separately analysed Chordata and the remaining Bilateria (both r=0.38; all p=<0.001), with higher explanatory power values (R2≈26-28% and 16-17% respectively). Hypothesis 5b was also rejected in the Bilateria dataset, where the correlation between GC skew and branch length was negligible (r=0.13, R2<1%), but supported when we separately analysed Chordata (r=0.37, R2≈20%) and the remaining Bilateria (r=0.50; R2≈28%; p=<0.001 in all cases; Supplementary Fig. 48).

Hypothesis 7: the Ne dataset comprised 202 species, belonging to four phyla: Chordata (182 species), Mollusca (9), Arthropoda (8), and Echinodermata (3). All were free living, but all three locomotory capacity categories were represented: 180 species in high, 14 in intermediate, and 9 in low. There were four single-stranded mitogenomes in the dataset. Contrary to the working hypotheses 7a-c, Ne was not correlated with GORR (r=0.07, p=0.30; R2=0.010), branch length (r=0.03, p=0.66; R2<0.006), and GC skew (r=-0.14, p=0.05; R2≈0.030). It was also not associated with strandedness (R2<0.030) and with skew inversions (R2<0.005) (Supplementary Fig. 49 and Data 2).

Parasitism, locomotory capacity, and thermoregulatory type

Hypothesis 8: In agreement with the working hypothesis 8a, endoparasitic species exhibited by far the most highly rearranged mitogenomic architectures, followed by ectoparasites, parasitoids, micropredators, and free-living species (Figure 5B). In pairwise comparisons, all differences were statistically significant (Tukey HSD tests; details in Supplementary Data 3), and Bayes factors were huge (Figure 6A, Figure S50). Statistical power was good and effect sizes were very large or huge for key pairwise comparisons: endoparasites vs. free-living (Cohen’s d=2.6), and ectoparasites vs. free-living (1.5) (Figure 6B). The results were replicated with Chordata removed from the dataset and across the Ecdysozoa and Lophotrochozoa (Supplementary Fig. 45). The patterns across the individual phyla were congruent in Platyhelminthes, and partially congruent in Arthropoda (elevated GORR in ectoparasites but not in endoparasites) and Nematoda (GORR was nonsignificantly increased in endoparasites) (Supplementary Fig. 51), and mixed across lower taxonomic levels (class and order): consistent in Insecta, Chromadorea, Malacostraca, Hexanauplia, and Rhabditida, but not in Arachnida, Enoplea, and Diptera (Supplementary Figs. 52 and 53). While lmekin and brms analyses indicated that the division along the lifestyle lines (parasitism) explained only around 6% of the variability in GORR in the bilaterian dataset, these values were higher once Chordata were excluded from the dataset: ≈33% in Ecdysozoa, ≈16% in Lophotrochozoa, and ≈20% in Protostomia (Table 3; Supplementary Data 2). To reduce the effect of variability in locomotory capacity, we focused only on the low locomotory capacity category divided along the lifestyle lines (hypothesis 8b). The pattern produced by the overall bilaterian dataset was only partially upheld, with the highest GORR values exhibited by endoparasites, but ectoparasites exhibited a lower value than free-living (Figure 5C). The explanatory power remained similar to the one in the full dataset (≈8%, Table 3).

Hypothesis 8c: in support of the working hypothesis, the far end of the positive GC skew distribution was marked by the strong prevalence of endoparasites: 24 species with highest skews were all endoparasites, and among the top 100 species with highest skews (GC range = 0.477 – 0.684), 93 were endoparasitic. A similar pattern was observed at the negative end of the GC skew distribution, where four highest values were exhibited by endoparasites (18 endoparasites in total in the GC-negative dataset). In partial agreement with these observations and the working hypothesis, absolute GC skews were the highest in endoparasites, but ectoparasites and parasitoids exhibited lower average values than the free-living species (Figure 5D, Supplementary Fig. S47). In support of hypothesis 8d, we found an increased likelihood (>37×) of skew inversions (positive GC skew) in macroparasitic (ectoparasitic+endoparasitic) species compared to free living (p<0.001, Fisher’s test, one-sided; Supplementary Fig. 54), with life history explaining 22% of the variability in skew direction (plus or minus).

Finally, in support of the hypothesis 8e, we found an increased likelihood of single-stranded mitogenomes in endoparasitic (>17×) and ectoparasitic (>3×) species , and decreased in other categories (<0.6× in free-living; Figure 7). Also, there was a highly significant association between life history and strandedness (p<0.001, Chi-square test of independence; Supplementary Data 3), with life history explaining 38-41% of the variability in the distribution of strand types across lineages (Table 3).

Hypothesis 9: in agreement with the working hypothesis 9a, GORR values were significantly increased in the low locomotory capacity category (Figure 5E). For the high vs. low category comparison, Bayes Factor was huge (Figure 6A, Supplementary Fig. 50), statistical power was good, and the effect size was huge (Cohen’s d=2.8; Figure 6C). The results were replicated with Chordata removed from the dataset and across all but one major taxonomic clades when both categories were present in sufficient numbers (n≥3): superphyla (Ecdysozoa+Lophotrochozoa), phyla (Chordata, Arthropoda, and Mollusca), and classes (Insecta, Malacostraca, and Arachnida) (insufficient data at the order level). The class Hexanauplia was the only major lineage where the hypothesis was rejected (Supplementary Figs. 45, 55, and 56). The explanatory power of locomotory capacity categorisation on GORR variability was ≈14% in Bilateria, ≈11% in Lophotrochozoa, and 33% in Ecdysozoa (Table 3). To remove the effect of parasitism on the locomotory capacity classification, we conducted analyses using only free-living species (hypothesis 9b). The pattern observed in the full dataset was mostly fully upheld (Figure 5F), and the explanatory power was ≈11% (Table 3). Among the ectoparasites (hypothesis 9c), those with low locomotory capacity had strongly (1.6×) but non-significantly increased GORR values (Figure 5G). In support of hypothesis 9d, GC skews were significantly increased in low-locomotory species (Figure 5H). In support of hypothesis 9e, there was an increased likelihood (>99×) of skew inversions (positive GC skew) in low-locomotory category compared to the high-locomotory category species (p<0.001, Fisher’s test, one-sided; Figure S54), and locomotory capacity explained 47% of the variability in skew direction. In support of hypothesis 9f, we found an increased likelihood (>3×) of single-stranded mitogenomes in species classified into the low locomotory capacity category, and decreased among the high-locomotory species (<0.4×, Figure 7, Supplementary Data 3). There was a highly significant association between the locomotory capacity and strandedness (p ≈ 0.0, Chi-square test of independence), with locomotory capacity explaining ≈33% of the variability in strand types across lineages (Table 3).

Hypothesis 10: in the bilaterian dataset 8,498 species were ectotherms, and 2,362 were endotherms (all Chordata and almost all free-living). Among the endotherms, >99% were classified in the high locomotory capacity category (none were low-locomotory). Among the ectotherms, ≈75% were high-locomotory, and ≈14% were low-locomotory. In agreement with the working hypotheses 10a/b, ectotherms exhibited significantly higher GORR (hypothesis 10a, Figure 5I), and branch lengths (hypothesis 10b, Figure 5J) than ectotherms, but Chordata exhibited an inverted pattern (Supplementary Fig. 45). The explanatory power of thermoregulatory type on these two variables was small: ≈2% and ≈5%, respectively (Table 3). Contrary to the working hypothesis 10c, endotherms exhibited significantly higher average GC skew value than ectotherms (both across the Bilateria and only Chordata; Figure 5K and Supplementary Fig. 47).

Mitogenome size patterns and correlates

Hypothesis 11: Aside from hypothesis 11i, all other hypotheses associated with correlates of mitogenomic size were rejected. Mitogenomes were the smallest (≈12.6 Kb) in Longidoridae (Nematoda: Enoplea), and the largest in two Bivalvia (Mollusca) species: Anadara sativa (48.2 Kb) and Dreissena rostriformis (46.7 Kb). Only eight species exhibited mitogenomes >30 Kb (5 Mollusca and 3 Arthropoda). Multiple Chordata species (Squamata and Anura) were among the largest mitogenomes (>25 Kb) (Supplementary Data 1). The correlation between GORR and mitogenome size was weakly/moderately negative in Bilateria (hypothesis 11a; r=-0.39, Figure 3B), and size explained ≈10% of the GORR variability (Table 3). In individual phyla, correlation between GORR and size was highly variable: positive in Annelida, Echinodermata, and Nematoda (weak – moderate); negative in Platyhelminthes and Mollusca (weak); and negligible in Arthropoda, Chordata, and Bilateria (Figure 3B, Supplementary Fig. 57). The correlation between mitogenome size and branch length was weakly negative across the bilaterian dataset (hypothesis 11b, r=-0.33, Figure 3C). It also varied across the individual phyla: positive in Mollusca, Nematoda, Echinodermata, and Annelida (weak – moderate); negative in Arthropoda and Platyhelminthes (negligible – weak); and negligible in Chordata and Bilateria (Figure 3C, Supplementary Fig. 57). Ne values were not correlated with size (hypothesis 11c, r=-0.09, p=0.21, R2<0.013, Supplementary Fig. 49). GC skew was also negligibly correlated with size across all tested datasets (hypothesis 11d, r≈-0.07 - 0.17; R2≈0.005 - 0.030; Table 3, Supplementary Fig. 48). Mitogenomes were the largest in free-living, and the smallest in ectoparasitic and endoparasitic lineages (hypothesis 11e, Figure 5L). Mitogenomes were larger in high-locomotory than in low-locomotory categories (hypothesis 11f, Figure 5M). Double-stranded mitogenomes were larger than single-stranded ones (hypothesis 11g, Figure 5N). Endotherms had larger mitogenomes than ectotherms (hypothesis 11h, Figure 5O). Finally, in agreement with hypothesis 11i, variability in size was lesser in endotherms than in ectotherms (SD = 537 vs. 1340, respectively, Figure 5O).

Methods
Dataset, categorisation, phylogeny, branch lengths, and Ne

Features of all mitogenomes of bilaterian animals available in the NCBI’s manually curated RefSeq database were extracted, calculated, and comparative tables generated using PhyloSuite 75,76. The dataset was categorised according to the lifestyle (endoparasites, ectoparasites, parasitoids, micropredators, and free-living species), locomotory capacity (high, intermediate, and low), and thermoregulatory type (endotherms and ectotherms). Duplicated genes were resolved using PhyloSuite following the methodology described previously 10,75. To infer branch lengths, we used the phylogeny from our previous study 10, reconstructed using concatenated amino acid sequences of 12 mitochondrial protein-coding genes (atp8 removed) of all mitogenomes in the dataset, and mtZOA+G4+F+C50 model in IQ-TREE v2.2.0.7.mix 77. We constrained the phylum-level topology according to the bilaterian phylogeny in Laumer et al.78, with Xenacoelomorpha as the sister lineage to all other bilaterians 79,80 (“Nephrozoa” topology henceforth). The molecular sequence evolution rate was defined as the root-to-tip branch lengths, which were extracted using the TreeSuite function in PhyloSuite. More details about the dataset, including the limitations of these classifications, and methodology are available in Supplementary Methods and in Jakovlić et al.10. The effective population size (Ne) analyses were extracted by overlapping the species in our dataset with those in the dataset recently compiled by Clarke et al.81. When there were multiple populations per species, we inferred average Ne values.

GC skews

To exclude the confounding factor of variability in the noncoding region size, GC skews were calculated 82 by PhyloSuite using only coding sections of the plus strand (GC skew is more pronounced than AT skew and thus more suitable for testing these hypotheses). To attempt to reduce the amount of noise in the GC skew dataset stemming from inversions in mutational pressures caused by inversions of the origin of replication 35,36 and facilitate the analysis of magnitude irrespective of skew direction, all GC skews were converted to their absolute values. To further reduce the amount of noise, the lower end of distribution (GC<0.1) was removed, as such low GC values commonly suggest ORI inversions or otherwise disrupted replication mechanisms, so they poorly reflect the strength of purifying selection pressure 36. In addition, we attempted to remove all annotation artefacts caused by the authors submitting the wrong mitogenomic strand to the GenBank via a manual inspection of the dataset (removed 65 species).

Mitogenomic architecture analyses

Mitogenomes of bilaterian animals commonly encode 37 genes (13 protein-coding genes - PCGs, 2 rRNAs, and 22 tRNAs) distributed in varying arrangements across both DNA strands, by convention named “plus” and “minus”. This arrangement is referred to as mitogenomic architecture, and it can be quantified as the gene order rearrangement rate (GORR). For the ancestral architecture inference, we relied on previous proposals of ancestral architectures for different lineages (Supplementary Table 3), and conducted the ancestral mitogenomic architecture inference using TreeREx 83 and MLGO 84 using a linearised gene order dataset from which tRNAs were removed due to their highly variable content and difficulties in annotation. To minimise the noise and reduce the size of dataset (362 species after the reduction), we used only the ancestral architecture for a phylum when it was inferred with confidence by previous studies (e.g. only the mitogenome of Homo sapiens for Chordata 5,30 and Limulus polyphemus for Arthropoda 85), and only the “basal” lineage of large phyla exhibiting high rearrangement rates (e.g. only Catenulida for Platyhelminthes 57). As the phylogeny of Bilateria remain debated, we tested two topologies, both constrained at the phylum level. Aside from the “Nephrozoa” topology, described in the first section, we also tested the topology where Xenacoelomorpha is the sister group to Ambulacraria (Hemichordata+Echinodermata) 86,87, with paraphyletic Deuterostomia 88 - named “Xenambulacraria” henceforth. As this topology implies that the common ancestor of Chordata and Xenambulacraria (Ambulacraria+Xenacoelomorpha) was the last common ancestor of all Bilateria, and characters common to these two clades would be bilaterian plesiomorphies 88, we removed the highly rearranged Protostomia (Ecdysozoa and Lophotrochozoa) 5,30, and tested a dataset comprising only the outgroups, Xenambulacraria, and Chordata, constrained according to the relationships proposed previously 86,88. The outgroups were Cnidaria, Porifera, and Placozoa. Further methodological details are provided in Supplementary Methods. The inferred ur-bilaterian mitogenomic architecture (tRNAs excluded) was then used to calculate mitogenomic GORR values using the Breakpoint Distance measure in CREX 89. As TreeREx and CREX cannot handle missing genes, for these analyses, atp8 was removed from all mitogenomes, and species with further missing PCGs or rRNAs were removed from these analyses, which included the entire phylum Chaetognatha - lacking atp6 2,90. In addition, all mitogenomes were classified into two categories according to the strand distribution of genes: double-stranded (genes encoded on both strands) and single-stranded (all genes encoded on a single strand), and then the Maximum Likelihood method in BayesTraits v4.0.1 91 was used to infer the ancestral states of strand distribution of genes. To get a better resolution for strandedness analyses, we further analysed complete GenBank mitogenomic datasets (as opposed to only RefSeq) for all poorly represented major lineages, and for all lineages with complex evolutionary histories.

Statistical analyses

As GORR values equalled 0 for more than half of the dataset, mostly due to a majority of Chordata having a perfectly conserved ancestral bilaterian architecture and being overrepresented in the dataset (>50% of all species), we recoded GORR=0 values to 0.0000001 for all statistical analyses. To further confirm that Chordata did not skew our statistical analyses, we additionally tested all major hypotheses using a bilaterian dataset from which all Chordata were removed. Outliers in GORR distributions were identified using the IQR method (InterQuartile Range), selected due to the non-normality of GORR distribution. Pairwise comparisons of between different groups were conducted using Tukey HSD tests in the R package agricolae v1.3.5 92. As the limitations of the traditional p-value-based statistics are well-documented, we further assessed the reliability of statistical analyses by employing multiple recommended alternative tests 93. Effect sizes were inferred using the R package efsize. Statistical power analyses were conducted for each effect size measure using the pwr.t2n.test method in the R package pwr. Herein, power refers to the probability of detecting a statistically significant effect if one exists, and the commonly accepted cut-off value for a good statistical power is ≥ 0.8. Bayesian hypothesis testing was conducted using BayesFactor 94. As the dataset violated the assumption of independence of data, two algorithms designed to account for phylogenetic relatedness were used to conduct multilevel linear regression analyses: linear fixed-effect models accounting for kinship implemented in the lmekin function in coxme v2.2.16 95, and phylogenetic multilevel Bayesian models implemented in brms v2.16.150 96. For these analyses, we used a matrix of phylogenetic distances extracted from the phylogenetic tree described above using the ape v5.5 package 97. The R2 value of lmekin models was calculated using the r.squaredLR method available in the MuMIn v1.46.0 package in R. The Bayesian R2 value was inferred using the bayes_r2 function of brms package in R. For brms and lmekin analyses, we log-transformed the obtained GORR data to reduce the non-normality of distribution. Pearson’s correlation coefficients (r) were inferred using R. Correlation values were classified as follows: negligible (<0.2), weak (0.2-0.4), moderate (0.4-0.6), strong (0.6-0.8), very strong (>0.8). Due to the very large number of samples in our study (>10,000), we did not rely on p-values to assess the support for the working hypotheses in correlation tests, as statistical tests are practically always highly significant on such high numbers of samples 98. Unless if specified otherwise, statistical tests were inherently two-sided.

Peer review rounds
[{"round": 1, "reviewer_comments": {"reviewer_1": "• What are the noteworthy results?\nThe manuscript presents a large-scale comparative analysis of evolutionary patterns across bilaterian mitochondrial\ngenomes. It examines correlations among diverse genomic and organismal traits, including mitochondrial genome size,\nstrand distribution of protein-coding genes, gene rearrangements, phylogenetic branch lengths, life history characteristics,\nand locomotory behavior. The study is remarkable for its scale, incorporating over 10,000 mitochondrial genomes with broad\nrepresentation across bilaterian phyla.\n• Will the work be of significance to the field and related fields?\nThis work is significant as a comprehensive compilation and statistical exploration of correlations between mitochondrial\ngenome architecture and organismal traits. However, I am concerned that an important confounding factor is not explicitly\naddressed: the inherently high mutation rates of animal mitochondrial genomes. This elevated mutation rate not only drives\nstrand-specific base composition skew but also results in codon and amino acid usage biases linked to the strand on which\nprotein-coding genes reside. For example, a switch in coding strand due to a change in replication origin orientation could\nrapidly alter the trajectory of protein evolution, thereby increasing phylogenetic branch lengths inferred from protein\nsequences. This could potentially explain some of the apparent contradictions noted in the analyses.\nInterestingly, the authors have computed GC and AT skews for all taxa but do not seem to incorporate this data into their\nmain analyses. Given the authors’ own previous work identifying origin-of-replication inversions and their consequences on\nprotein evolution (e.g., Zhang et al. 2019), integrating strand asymmetry and skew into the analyses could substantially\nstrengthen the manuscript’s conclusions.\nAdditionally, I question the assumption that mitochondrial genome size is directly linked to metabolic rate, suggesting\ninstead that it may be more closely related to mitochondrial or organismal population sizes. Another life history dimension\nworth exploring is the well-documented difference in mitochondrial protein evolutionary rates between endotherms and\nectotherms, with competing hypotheses attributing this to either metabolic demands or to fluctuating constraints requiring\noptimization across diverse temperatures in ectotherms versus more stable ranges in endotherms.\n• How does it compare to the established literature?\nThe authors cite a recent study of their own (Reference #10) that appears to analyze the same dataset of ~10,000 genomes,\nfocusing on two life history traits that are also examined in this manuscript. It would be helpful for the authors to clarify how\nthis new analysis builds upon or differs from their prior work to better establish its novelty.\nGiven the authors’ established expertise on the effects of strand skew, genome architecture, and evolutionary rates, this\nmanuscript would be considerably enhanced by explicitly analyzing correlations involving strand asymmetry alongside the\nother genomic and organismal traits.\n• Does the work support the conclusions and claims, or is additional evidence needed?\nThe statistical analyses seem to support the existence of correlations among the traits examined, although these\nrelationships often include numerous exceptions. Incorporating strand skew and replication architecture could help interpret\nsome of these exceptions.\n• Are there any flaws in the data analysis, interpretation, or conclusions? Do these prohibit publication or require revision?\nNo critical flaws are noted. However, integrating strand asymmetry data into the analyses and more explicitly linking\nstatistical patterns to biological hypotheses would improve the interpretative strength of the study.\n• Is the methodology sound? Does the work meet the expected standards in your field?\nThe methodology appears technically sound and meets the standards of the field. That said, many of the analyses are\nframed as correlations among categorical variables rather than tests of explicit, biologically motivated hypotheses.\nStrengthening this aspect would enhance the impact of the findings.\n• Is there enough detail provided in the methods for the work to be reproduced?\nThe effort to categorize life history traits across such a large number of species is commendable. However, I found some of\nthe definitions and classifications difficult to interpret. For instance, all nematodes are categorized as low locomotion (LOC),\nyet the link between this and low constitutive metabolic demand, which seems to be the underlying hypothesis, is not fully\njustified. Similarly, it is unclear how parasitic taxa with free-living stages were treated, or why some necromenic\nCaenorhabditis species are classified as free-living. Given the inherent complexity of categorizing diverse life histories, it\nwould be helpful for the authors to emphasize potential ambiguities and how these might contribute to exceptions in their\nanalyses.", "reviewer_2": "This manuscript presents a comprehensive analysis of mitogenomic architecture evolution across Bilateria, leveraging a\ndataset of over 10,000 mitogenomes. The study addresses several unresolved questions regarding gene order\nrearrangements, strandedness transitions, and correlations with lifestyle and locomotory capacity. The work is well-\nstructured, methodologically robust, and provides significant insights into mitogenomic evolution. However, some areas\nrequire clarification or additional analysis to strengthen the conclusions.\nMajor Comments\n1. Ancestral State Reconstruction\n- The inference of the ancestral bilaterian mitogenomic architecture as double-stranded is well-supported, but the instability\nin \"nad6\" positioning across analyses warrants further discussion. Could this instability reflect methodological limitations\n(e.g., alignment ambiguity or model sensitivity)? A deeper exploration of alternative topologies (e.g., Xenambulacraria vs.\nNephrozoa) and their impact on ancestral state inference would strengthen the conclusions. The authors may want include a\nsensitivity analysis to assess the robustness of the ancestral state.\n2. Single- to Double-Stranded Transitions\n- The identification of rare single-to-double-stranded transitions in Annelida and Mollusca challenges the \"ratchet\"\nhypothesis. However, the evidence for these transitions (e.g., Sabella spalanzani) is based on limited taxa. Could annotation\nartifacts or incomplete lineage sampling explain these observations?\n3. Mitogenome Size and Evolutionary Rates\n- The negative correlation between mitogenome size and evolutionary rates contradicts the \"race for replication\"\nhypothesis. The proposed explanations are plausible but speculative.\n4. Parasitism and Locomotory Capacity\n- While the model linking low locomotory capacity to relaxed purifying selection is compelling, the explanatory power of\nthese variables for GORR variability is modest (6–14%). Other factors (e.g., effective population size, recombination rates)\nmay play significant roles. A discussion of how neutral processes might interact with selection might help.\n5. Population size of mitochondrion may substantially different from the Ne of species.\n---\nMinor Comments\n1. Statistical Analyses\n- The use of phylogeny-corrected models (lmekin, brms) is appropriate, but the R² values for some comparisons (e.g., 3% for\nGORR≥7/ds) are very low. Emphasize that these results suggest additional undetected factors drive architectural evolution.\n2. Outliers\n- The exceptional conservation in Cirripedia (Hexanauplia) is intriguing. The hypothesis of large effective population size is\nuntested. Could life-history traits (e.g., sessility, larval dispersal) also contribute?\n3. Supplementary Material\n- Not all supplementary figures/tables are directly cited in the main text.", "reviewer_3": "I co-reviewed this manuscript with one of the reviewers who provided the listed reports. This is part of the Nature\nCommunications initiative to facilitate training in peer review and to provide appropriate recognition for Early Career\nResearchers who co-review manuscripts.\nI did not run the code, I only reviewed it.\nThe R code to create the figures and tables is provided in Zenodo. The R code is mostly intuitive and allows the figures to be\nreproduced, however a README is not provided with a description of each script. This is not necessary but listing out the\nrequired dependencies would help others run the code with new data.\nThe input files should be explained further in the manuscript. For example, how were organisms categorized into lifestyles\nwithin the input spreadsheet?"}, "author_response": ">>Response: Thank you for investing your time and expertise in reviewing our study. We appreciated\nyour comments and critiques - they helped us improve our study. Four major issues were raised in\nthis comment: (1) the confounding effect of high mutation rates (branch length), (2) the absence of\nanalyses of base composition skews and the impact of skew inversions on our findings, (3) the\nimpact of Ne on mitogenomic size, and (4) the impact of thermoregulatory type on the protein\nevolution. We discuss them individually below.\n(1) As regards long branches, we agree with the reviewer that they can interfere with multiple\nevolutionary analyses, but we were careful not to include any analyses that can be rendered\nmeaningless by rapidly evolving lineages in this study, e.g. time-tree reconstruction, dN/dS analyses,\nand inference of phylogenetic relationships. While we did infer a phylogeny using our dataset, we\nconstrained it at the phylum level according to a previously inferred nuclear genomic phylogeny, and\nused only branch lengths for our analyses. This also applies to our inference of the evolutionary\nhistory of strand distribution of genes, for which we also relied on multiple published phylogenies\n(as opposed to our own). Therefore, the only relevant variable is: do our branch length estimates\nreflect the actual cumulative mitogenomic sequence evolutionary rates of individual lineages? A\nsource of noise in this aspect may be substitutional saturation (multiple substitutions at the same\nsite) and across-site compositional heterogeneity, both of which can impact the inference of branch\nlengths, as well as tree topologies (Foster 2004; Lartillot and Philippe 2004; Lartillot et al. 2007;\nMoody et al. 2022). We considered this problem when designing our analyses. As ‘heterogeneous’ or\n‘mixture’ models (i.e. models accounting for compositional heterogeneity) tend to outperform\n‘homogeneous’ models in the correct estimation of branch lengths on highly saturated,\ncompositionally heterogeneous datasets (Lartillot et al. 2007; Feuda et al. 2017; Kapli et al. 2020;\nWilliams et al. 2021), in order to achieve maximum accuracy of branch length estimates, we\nemployed a ‘mixture’ evolutionary model (C50) implemented in IQ-TREE. Furthermore, to confirm\nthat results are robust to the model selection, we tested several other models, and found that\nbranch lengths were highly correlated across them (>96%) (Jakovlić et al. 2023). On this basis, we\ncan conclude that our findings are robust to the evolutionary model employed. Nevertheless, we\ncannot exclude a possibility that the model that we chose slightly underestimated branch lengths of\nfast-evolving lineages due to imperfectly modelled reverse mutations and varying numbers of silent\nmutations. We attenuated the issue of reverse mutations by relying on amino acid alignments to\ninfer branch lengths, as reverse mutations are much less likely at the amino acid level than at\nnucleotide level. A minor downside is that they may not fully capture the extent of silent nucleotide\nsubstitutions. Although this issue likely affects most lineages included in the dataset, it is reasonable\nto assume that the fastest-evolving lineages, such as Neodermata, may be disproportionately\naffected. Importantly, as there is a strong correlation between branch lengths (brl) and GORR, there\nis a pronounced overlap in fastest-evolving lineages between these two variables. For example,\nNeodermata have both the highest average sequence and architectural evolutionary rates among\nbilaterians. As “reverse mutations” are extremely rare in terms of gene rearrangements (Boore\n2006; Babbucci et al. 2014), this issue was more likely to cause an underestimation than an\noverestimation of correlation and regression values between these two variables. We observed this\nwhile conducting preliminary analyses for our previous Nature Communications study (Jakovlić et al.\n2023): the explanatory power of key variables on branch length variability was slightly increased (a\nfew % points) when we used branch lengths inferred using ‘mixture’ models, as opposed to using\nbranch lengths inferred using ‘homogeneous’ models. In addition, as the effect sizes were huge for\nthe key findings in this study, there is no reason to assume that a small level of imprecision in branch\nlength estimates can produce a type 1 error (false positive findings). Therefore, our branch length\nestimates should be sufficiently accurate for the purposes of our analyses, and minor inaccuracies\nwere highly unlikely to affect our findings in a significant manner. Moreover, if the noise in branch\nlength inference were to be completely eliminated, our hypotheses would likely be even more\nstrongly supported.\nAs a side note, another issue that the reviewer may have had in mind when pointing out the long\nbranch issue is the discrepancy in evolutionary rates between the comparatively fast-evolving\nmitogenomes (mtDNA) and the comparatively slow-evolving nuclear genomes (nDNA). The\nunderlying reasons for this discrepancy are believed to be the reduced effective population size (Ne)\nof mtDNA caused by the uniparental inheritance (along with a few other factors) and the high\nproduction of mutagenic ROS (reactive oxygen species) within the mitochondria. This discrepancy in\nevolutionary rates between the two genomic compartments would affect - and possibly invalidate -\nany comparative analysis that comprised a mixture of both mtDNA and nDNA data. However, we\nfocused strictly on mtDNA data in our study, and all bilaterian lineages share this reduced mtDNA Ne\neffect.\n(2) As regards the question of base composition skews, as the reviewer correctly pointed out, and as\nour own analyses show (Jakovlić et al. 2021b), any analysis involving GC skew magnitude will be\ninvariably hampered by multiple sources of noise. This was the reason why we originally decided to\nexclude skew analysis from our study. Following the reviewer’s comment, and as GC skews are\nindeed closely associated with mitogenomic architecture, we included them in this revised study,\nbut multiple sources of noise in this dataset must be understood and acknowledged when\ninterpreting our findings. We briefly mentioned the limitations in the revised study, but did not go\ninto details in efforts to save space. Briefly, strand switches of individual genes can strongly affect\noverall GC skew values, whereas variations in the replication mechanism and strand switches of the\norigin of replication can produce full skew inversions (e.g. from negative to positive GC skew values)\n(Hassanin et al. 2005; Jakovlić et al. 2021b) – i.e. render direct comparisons of skew values across\nlineages meaningless. We attempted to ensure comparability of GC skew values by converting them\nto their absolute values. This allowed us to analyse the skew magnitude, irrespective of direction.\nFurther, the outliers at the lower end of distribution (GC<0.1) were removed as our experience\nshows that they suggest ORI inversions, recent inversions of large mitogenomic segments, or\notherwise disrupted replication mechanisms (Jakovlić et al. 2021b). To further attempt to make sure\nthat comparison of GC skew magnitudes is meaningful across different lineages, they were\ncalculated using only coding sections of the plus strand – thus minimising the impact of the\nvariability in the noncoding regions size from the dataset. Finaly, we attempted to remove all\nannotation artefacts caused by the authors submitting the wrong mitogenomic strand to the\nGenBank via a manual inspection of the dataset (removed 65 identified species). As expected, the\nresults were rather noisy, with most hypotheses accepted, but often only partially, or with notable\noutliers. The support for working hypotheses was generally much stronger when the two datasets\n(Chordata and Bilateria-minus-Chordata) were analysed separately due to high GC skews in\nChordata. The inclusion of GC skews also produced some interesting findings: results indicate that\nthe evolutionary transition to single-strandedness largely (but not exclusively) results in a GC skew\ninversion (from minus to plus), and that skew inversions are strongly associated with parasitic\nlifestyle and low locomotory capacity. We have included these findings in the revised study and its\nsupplementary materials, and discussed their implications in the revised manuscript.\n(3) As regards the impact of the effective population size (N ) on mitogenomic size, we have tested it\ne\na few times before, and each time we found that it had a low predictive power on mitogenomic\nevolution (Jakovlić et al. 2021a, 2023, 2024). As these were not isolated outlier findings (Bazin et al.\n2006; Li et al. 2024), we initially decided to exclude N from the analyses. Following the reviewer’s\ne\ncomment, we included it in the revised study. It should be noted that molecular N estimates are\ne\nsensitive to a number of methodological and life history parameters (Lynch and Conery 2003; Palstra\nand Ruzzante 2008), so N inference is considered difficult and error-prone. In our previous studies,\ne\nwe relied on population data for mitogenomic cox1 gene to infer N , which is a methodology prone\ne\nto artefacts stemming from circular reasoning. For this study, we relied on N data from a recently\ne\npublished dataset, compiled from a range of sources (Clarke et al. 2024). Despite the different\nmethodology, once again we found that N is practically completely decoupled from all tested\ne\nvariables. We discuss the most likely underlying reasons for this observation in the revised study.\n(4) As regards the impact of thermoregulatory type on mitogenomic evolution, we also included this\nvariable in the revised study. Predictions regarding the impact of this variable are intriguingly\ncontradictory. We found partial support for our hypothesis that locomotory capacity-mediated\npurifying selection pressure outweighs the predictions of “functional constraints” and” metabolic\nrate” hypotheses – with endotherms exhibiting lower sequence and architectural evolutionary rates.\nThe explanatory power of thermoregulatory type was lower than that of parasitism and locomotory\ncapacity for these two variables. However, as in a multiple other cases, GC skew and size did not\nfollow this pattern. We discuss these findings and contradictions in the revised study.\n• How does it compare to the established literature?\nThe authors cite a recent study of their own (Reference #10) that appears to analyse the same\ndataset of ~10,000 genomes, focusing on two life history traits that are also examined in this\nmanuscript. It would be helpful for the authors to clarify how this new analysis builds upon or differs\nfrom their prior work to better establish its novelty.\nGiven the authors’ established expertise on the effects of strand skew, genome architecture, and\nevolutionary rates, this manuscript would be considerably enhanced by explicitly analysing\ncorrelations involving strand asymmetry alongside the other genomic and organismal traits.\n>>Response: As outlined in detail elsewhere in this response letter, we included a thorough analysis\nof strand asymmetry into the revised study. As regards the former study in Nature Communications\n(Jakovlić et al. 2023), it was focused strictly on the mitogenomic sequence evolution – i.e. it did not\ntouch the topic of architectural evolution (GORR and mitogenome size) at all. Initially, we designed\nthe current study to complement the first one by focusing primarily on the architectural evolution.\nAs mentioned earlier, following reviewers’ comments, we included additional variables for the\nrevised version, so the revised study mostly focuses on architectural evolution, but also touches\nupon aspects of sequence evolution not studied in our previous research. Following the reviewer’s\nobjection, we clarified this in the opening paragraph of the revised Discussion: “We offer an in-depth\nanalysis of the evolutionary forces influencing mitogenomic architecture across Bilateria. It builds upon our previous study,\nwhere we used a similar dataset to show that parasitism and locomotory capacity have the power to explain the variability\nin sequence evolution rates (Jakovlić et al. 2023), but did not include architectural rearrangement rates, size,\nthermoregulatory type, N, and GC skew among the studied variables.”\ne\n• Does the work support the conclusions and claims, or is additional evidence needed?\nThe statistical analyses seem to support the existence of correlations among the traits examined,\nalthough these relationships often include numerous exceptions. Incorporating strand skew and\nreplication architecture could help interpret some of these exceptions.\n>>Response: As mentioned in a previous answer, we have thoroughly incorporated skew analyses\ninto the revised manuscript, tested multiple hypotheses associated with skews and replication\narchitecture, and discussed the implications of our findings.\n• Are there any flaws in the data analysis, interpretation, or conclusions? Do these prohibit\npublication or require revision?\nNo critical flaws are noted. However, integrating strand asymmetry data into the analyses and more\nexplicitly linking statistical patterns to biological hypotheses would improve the interpretative\nstrength of the study.\n>>Response: We have incorporated strand asymmetry analyses and rewritten parts of the\nmanuscript to more explicitly link our analyses to biological hypotheses. Notably, we made sure that\nbiological underpinnings of hypotheses tested are clear in the Introduction. As adding them to Table\n1 would make it user-unfriendly, we added them to the overview of hypotheses in Supplementary\nfile 1: Text S1.\n• Is the methodology sound? Does the work meet the expected standards in your field?\nThe methodology appears technically sound and meets the standards of the field. That said, many of\nthe analyses are framed as correlations among categorical variables rather than tests of explicit,\nbiologically motivated hypotheses. Strengthening this aspect would enhance the impact of the\nfindings.\n>>Response: As mentioned in the previous response, the biological background of expected\nrelationships was outlined in the Introduction of the original paper. Following the reviewer’s\ncriticism, we reframed the text to make sure that biological underpinnings of hypotheses tested are\nclear in the Introduction. We also made some changes in the Discussion of the main text. Major\nchanges in this aspect were made in Supplementary file 1: Text S1, where we provided an overview\nof hypotheses, along with the biological underpinning of each one. In Table 1, we did not make these\nchanges, as the table was designed to provide readers with a quick overview of tested hypotheses,\nand thus we feel that it would lose its purpose if we added a lot of text to it. If the reviewer still feels\nthat the issue is insufficiently addressed, we would be happy to move the overview of hypotheses\nfrom Text S1 to the main file.\n• Is there enough detail provided in the methods for the work to be reproduced?\nThe effort to categorize life history traits across such a large number of species is commendable.\nHowever, I found some of the definitions and classifications difficult to interpret. For instance, all\nnematodes are categorized as low locomotion (LOC), yet the link between this and low constitutive\nmetabolic demand, which seems to be the underlying hypothesis, is not fully justified. Similarly, it is\nunclear how parasitic taxa with free-living stages were treated, or why some necromenic\nCaenorhabditis species are classified as free-living. Given the inherent complexity of categorizing\ndiverse life histories, it would be helpful for the authors to emphasize potential ambiguities and how\nthese might contribute to exceptions in their analyses.\n>>Response: These are all valid concerns. We did not discuss these issues in the original manuscript\nin detail simply because the problems with both lifestyle and locomotory capacity categorisations\nwere discussed in detail in our previous study published in Nature Communications (Jakovlić et al.\n2023), which relied on a similar dataset to test hypotheses related to sequence evolution (as\nopposed to this study, which primarily focuses on the architecture evolution). To avoid repetition,\nand as this is an open-access journal (so the previous study is freely available to all readers), we\ninitially decided to refer readers to the previous study for these explanations at the end of the first\nsection of Methods: “More details about the dataset and phylogenetic analyses in Supplementary\nfile: Text S2 and in Jakovlić et al. (2023).” Following the reviewers’ objections, we slightly rewrote\nthe discussion of these issues from previous study and added it to Supplementary file: Text S9 of the\nrevised submission. We inserted the text below (at the end of response) to help reviewers and\neditors assess our discussion of these problems. As can be seen in the text, there were multiple\nlineages whose categorisation was difficult. This is expected in biology, where most categorisations\nare doomed to be deeply flawed. This is why we opted for the strategy of maximising the number of\nspecies, as a small number of outliers is much less likely to affect statistical analyses of large\nsamples.\nOverall, Nematoda were probably the most difficult phylum in this aspect, as there are many\nlineages in this phylum that resisted categorisation. We hypothesize that this is likely to be a\nconsequence of their extremely complex evolutionary history in terms of lifestyle, with multiple\nindependent origins of parasitism (Blaxter and Koutsovoulos 2015). This appears to produce very\nsubtle differences between nominally free-living and parasitic lineages. In our previous study, we\nalso found that nematodes were a major outlier among bilaterian animals - in parasitism having a\nnonsignificant impact on branch length - and attributed it to the evolutionarily complex mix of\nparasitic and free-living life histories that exhibit many similarities between nominally parasitic and\nfree-living nematodes (Jakovlić et al. 2023). The complex evolutionary history of nematodes is also\nevident from the count of independent origins of parasitism in our dataset introduced in the\nprevious study: “Accounting for some uncertainty in evolutionary scenarios (e.g. in Nematoda and\nPlatyhelminthes), independent origins of endoparasitism were observed in Acanthocephala (1 origin),\nArthropoda (2: Pentastomida and Rhinonyssidae), Orthonectida (1), Platyhelminthes (2 or 3), and between 5\nand 7 origins in Nematoda according to the evolutionary scenario proposed previously 34. We identified 12\nputative origins of ectoparasitism in Arthropoda and 1 in Platyhelminthes (see Supplementary Note 2 and\nSupplementary Table 3).”\nWe also mentioned this in the discussion of the current study. As regards the specific example raised\nby the reviewer - necromenic Caenorhabditis species - this is another example of difficulties\nassociated with the categorisation of nematodes. However, we would still argue that necromeny of\nthese nematodes does not fit the definition of parasitism employed in our classification, outlined in\nthe previous study: “parasitism was defined as a consumer interaction in which the consumer feeds\non a single individual (the host) during at least one life-history stage, where both parasites and hosts\nbelong to the Animalia, which includes strategies employed by parasitoids, parasitic castrators,\nmacroparasites and pathogens, but excludes micropredators, brood ‘parasites’, kleptoparasites,\nsymbiotic egg predators, inquilines, non-feeding symbionts, and plant-parasites (Lafferty et al. 2015;\nWeinstein and Kuris 2016).” Looking at this definition, it is missing the word “alive”, as necromenic\ndauer juveniles of nematodes simply remain on the body of the associated animal until it dies, and\nthen feed, further develop, and propagate on the decomposing cadaver (Kiontke and Sudhaus 2006).\nWe added the missing word to the definition used in this study. As they don’t feed on the live host,\nthis could be considered a form of non-feeding symbiosis. Also, it does not differ substantially from\nphoresy (Kiontke and Sudhaus 2006), where dauer juveniles of nematodes use the associated animal\nonly as a means of transportation to carry them to a new food source (e.g. a carrion of a different\nanimal). Therefore, their classification as free-living is relatively solid. We added the discussion of\nthis problem to Text S9.\nAs regards other difficult-to-categorise lineages, five Strongyloididae (Nematoda) species were\nexcluded from the dataset because they had a life cycle that alternates between free-living and\nparasitic generations 44,45 (this was mentioned in the previous study, and is also available in Text S9).\nParasitic species with free-living life stages were and remain a problem in our classification, as many\nparasitic lineages have a life stage that takes place outside the host. In most cases, however, it is the\nvery early life stages that are free-living (such as the egg or some early larval form), and to our\nknowledge a majority of these have a very limited locomotory ability. Finally, all of them have the\nadult stage that matures within the host. As this is in stark contrast to parasitoids, where the\nlarva is parasitic, whereas the adults are free-living and highly locomotory, we ultimately decided to\ncategorise all these species as parasitic, despite the complexity of their life-cycles. We added a\nmention of this problem to the revised Text S9, and stressed that it may cause a certain amount of\nnoise in the data.\nFrom Text S9:\nCategorisation of life history strategies\nThis categorisation was introduced in detail in our previous study (Jakovlić et al. 2023). We adopted the definition of\nparasitism as a consumer interaction in which the consumer feeds on a single live individual (the host) during at least one\nlife-history stage, where both parasites and hosts belong to the Animalia. This includes strategies employed by parasitoids,\nparasitic castrators, macroparasites and pathogens, but excludes micropredators, brood parasites, kleptoparasites,\nsymbiotic egg predators, inquilines, non-feeding symbionts, and plant-parasites (Lafferty et al. 2015; Weinstein and Kuris\n2016; Jakovlić et al. 2023). The dataset was classified according to life history into five categories: endoparasites (EndoP),\nectoparasites (EctoP), parasitoids, micropredators (MP), and free-living (F). Parasites living inside the host’s body (including\nthe intestines, nasal passages, etc.) were classified as endoparasites (this category also included mesoparasites) (274\nspecies). Parasites attaching themselves to or permanently living on the surface of the host’s body were classified as\nectoparasites (117 species). As we expected endoparasites to have a lower locomotory capacity and higher metabolic\ndependence on the host, we expected their evolutionary rates to be higher than those of ectoparasites, for which we\nfound strong support before (Jakovlić et al. 2023). Parasitoids (66 species) comprise organisms that are parasitic during a\npart of their life cycle, commonly this comprises a parasitic larva stage, followed by a free-living adult stage. As adults\nrequire a locomotory capacity that does not differ from standard free-living species, they should evolve under a stronger\npurifying selection pressure than ecto- and endo-parasites, and thus exhibit lower mitogenomic evolutionary rates. To\nfurther confirm our predictions, we separately classified micropredators (186 species). These organisms resemble many\nectoparasites in their haematophagous feeding habits, but they only visit the host for feeding, and otherwise they are free-\nliving (e.g. mosquitoes). We hypothesised that their evolutionary rates should not differ from other free-living organisms.\nAll other bilaterian species were classified as free-living (10,217 species). We found support for all these hypotheses, but\nonly regarding mitogenomic sequence evolution rates, in our previous study (Jakovlić et al. 2023), but evolutionary\npatterns of mitogenomic architecture in this aspect remain untested.\nBelow, we also provide further discussion of some potentially confounding lineages, also adopted from our previous study\n(Jakovlić et al. 2023). As mentioned above, all plant parasites were included in the free-living category, although some,\nsuch as plant-parasitic nematodes, may be similar to those classified as parasitic in terms of reduced locomotory capacity\nand confinement to a single host. Due to the high overall similarity between parasitoid and parasitic castrator strategies\n(Kuris 1974), the latter were also classified as parasitoids in most cases. Exceptions were some cestodes, which exhibit a\nmix of endoparasitic and parasitic castrator strategies (Arnott et al. 2000). With respect to our hypothesis (locomotory\ncapacity, metabolic dependence on the host, and physical confinement to a single host) they are more similar to\nendoparasites than to parasitoids, so they were classified as endoparasites. We assigned botflies (Insecta: Oestridae) to\nparasitoids, because their larvae are endoparasitic, and adults are free-living, which is a typical parasitoid lifestyle.\nHowever, as adults generally do not feed during their short life, their diet resembles that of typical parasites. However,\nselection for high locomotory capacity and the absence of a confinement to a single host strongly outweigh this argument.\nA micropredator attacks more than one host (but one at a time), reducing each host's fitness by at least a small amount.\nMost micropredators are haematophagous (feeding on blood). They include annelids such as leeches, many insects (such\nas mosquitoes, tsetse flies and bed bugs), and even some vertebrates, such as lampreys and vampire bats. They often\nexhibit a mix of feeding strategies. For example, in many mosquitoes (Culicidae), typically both males and females feed\nnon-parasitically. Only females need a blood meal in order to obtain the nutrients needed to produce eggs. Oxpecker\nBuphagus erythrorhynchus (Aves) is also such an example; their relationship with mammals was originally thought to be an\nexample of mutualism, but more recent evidence suggests that oxpeckers may be facultative parasites (micropredators)\n(Plantan et al. 2013). Following the above definition, sheep ked Melophagus ovinus was classified as an ectoparasite and\nnot as a micropredator, because it is believed to spend its entire life in the wool of a single sheep (Small 2005).\nTinaminyssus melloi and Ptilonyssus chloris (Arachnida: Mesostigmata: Rhinonyssidae) were classified as endoparasites.\nNasal mites show characteristics typical of endoparasitic species: reduced shielding, reduced setation, and overall, a body\ntype that would have reduced mobility compared to free-living or ectoparasitic relatives (Hilario-Pérez and Dowling 2018),\nso they were classified as endoparasites (along with other mesoparasites). Trouessartia rubecula, a feather mite, lives on\nbirds’ feathers, but feeds on feather oils, so we did not classify it as parasitic. In general, mites were often difficult to\nclassify, as their lifestyles often exhibit characteristics of free-living, ectoparasitic, and micropredatory lineages. The\nglochidium (plural glochidia) is a microscopic larval stage of some freshwater mussels, aquatic bivalve molluscs in the\nfamilies Unionidae and Margaritiferidae, the river mussels and European freshwater pearl mussels. Glochidia have hooks\nthat enable them to attach themselves to fish (for example to the gills of a fish host species) for a period before they\ndetach, fall to the substrate, and take on the typical form of a juvenile mussel. Since a fish is active and free-swimming, this\nprocess helps distribute the mussel species to areas that it could not reach otherwise. Despite its temporary attachment to\nthe host, we did not classify any of these species as parasitic or parasitoid. We failed to find sufficiently precise life-history\ninformation for a few insect species: Canthesancus helluo (Reduviidae) and two Sarcophagidae species: Goniophyto\nhonshuensis and Miltogramma oestracea. As larvae of some Sarcophagidae are internal parasites of other insects such as\nOrthoptera, and some Reduviidae are blood-sucking ectoparasites, we were not sure how these three species should be\nclassified, so we removed them from the dataset used for statistical analyses. The categorisation of nematodes was very\ndifficult, as there are many lineages in this phylum that resist categorisation, likely due to their extremely complex\nevolutionary history in terms of lifestyle, with multiple independent origins of parasitism (Blaxter and Koutsovoulos 2015).\nThis appears to produce very subtle differences between nominally free-living and parasitic lineages. For example, we\nremoved five Strongyloididae species from the dataset because they alternate between free-living and parasitic stages\n(Harvey et al. 2000; Luong and Mathot 2019). We should note that our initial analyses indicated that they exhibit almost\nthe same branch lengths as strictly endoparasitic nematodes. Contrary to this, we did not remove necromenic\nCaenorhabditis species from the dataset because necromeny of these nematodes does not fit the definition of parasitism\nemployed in our classification, outlined at the beginning of this section: “…a consumer interaction in which the consumer\nfeeds on a single live individual (the host) during at least one life-history stage, where both parasites and hosts belong to\nthe Animalia…”. As necromenic dauer juveniles of nematodes simply remain on the body of the associated animal until it\ndies, and then feed, further develop, and propagate on the decomposing cadaver (Kiontke and Sudhaus 2006), they do not\nfeed on the live host. This could be considered a form of non-feeding symbiosis, and it does not differ substantially from\nphoresy (Kiontke and Sudhaus 2006), where dauer juveniles of nematodes use the associated animal only as a means of\ntransportation to carry them to a new food source (e.g. a carrion of a different animal).\nParasitic species with free-living life stages were a major problem in our classification, as many parasitic lineages have a life\nstage that takes place outside the host. In most cases, however, it is the early life stages that are free-living (such as the\negg or some early larval form), and to our knowledge a vast majority of these have a very limited locomotory ability.\nFinally, all of them have the adult stage that matures within the host. As this is in stark contrast to parasitoids, where the\nlarva is parasitic, whereas the adults are free-living and highly locomotory, we ultimately decided to\ncategorise all these species as parasitic, despite the complexity of their life-cycles, and the potential noise that it may introduce in our analyses.\nCategorisation of locomotory capacity\nThe categorisation of locomotory capacity was based on the “visual interaction hypothesis”, which is the best explanatory\nvariable for variability in scaling coefficients between the mass-specific metabolic rate in marine animals: high metabolic\ndemand follows strong selection for locomotory capacity for pursuit and evasion in visual prey/predators inhabiting well-lit\noceanic waters, whereas limited visibility allows for reduced locomotory capacity, reflected in low metabolic rates\n(Childress and Mickel 1985; Seibel and Drazen 2007). We divided the dataset into three locomotory capacity categories: (1)\nHigh comprised all species expected to rely on locomotion for pursuit and evasion of prey/predators. This category\ncomprised a vast majority of species (8,713); (2) Low comprised all species that have merely a rudimentary locomotory\ncapacity. This comprises sedentary species, or those capable only of minimal, very slow locomotion, that do not rely on\nlocomotion for pursuit and evasion of prey/predators. This category comprised 1,202 species, mostly from lineages such as\nnematodes, flatworms, sessile molluscs and crustaceans, parasitic crustaceans, ticks, some annelids, mites,\nAcanthocephala, Diplura, Brachiopoda, Bryozoa, Ascidiacea, Echinodermata, Entoprocta, Gastrotricha, Gnathostomulida,\nHemichordata, Kinorhyncha, Nemertea, Onychophora, Orthonectida, Phoronida, Priapulida, Tardigrada, and\nXenacoelomorpha; (3) Because locomotory capacity is a continuous trait, it is impossible to categorise species into only two\ncategories without creating a large overlap between the two categories. Therefore, to minimise the overlap between the\nLow and High locomotory capacity categories, we designed a third category: Intermediate locomotory capacity. This\ncategory comprised 945 species that would be expected to possess more than a rudimentary locomotory capacity, but also\nrely on strategies other than locomotion to evade/pursue predators/prey. Examples are certain flightless insects and\nCollembola, certain bathyal, abyssal and stygobitic crustacean and fish lineages, Chaetognatha, certain amphibian and\nreptilian lineages, some mammals, some Bivalvia and Cephalopoda, some Gastropoda, Nematomorpha, and Rotifera.\nCertain lineages were rather difficult to classify. Chaetognaths swim in short bursts using a dorso-ventral undulating body\nmotion, where their tail fin assists with propulsion and the body fins with stabilization and steering, so we classified them\ninto the Intermediate category. Among the Annelida, Sipuncula merely have a rudimentary LOC, so they were classified as\nLow. Crassiclitellata (earthworms), Hirudinida (leeches) and Eunicida are capable of more than rudimentary locomotion, so\nwe classified them as Intermediate. Sabellida are sedentary, so we classified them as Low. Terebellida and Urechidae\nappear to possess merely a rudimentary locomotory capacity, so we placed them into the Low category. Magelonidae are\nburrowing worms with limited mobility, so we classified them as Low. Among the Nemertea, some species may possess\nmore than a basic locomotory capacity, but due to limited evidence, we classified them all as Low. Phoronida, Priapulida,\nEntoprocta, and Hemichordata possess merely a rudimentary locomotory capacity, so all were classified as Low.\nOnychophora can move at up to 5 cm per second, so we classified them as Intermediate.\nReviewer #2 (Remarks to the Author):\nThis manuscript presents a comprehensive analysis of mitogenomic architecture evolution across\nBilateria, leveraging a dataset of over 10,000 mitogenomes. The study addresses several unresolved\nquestions regarding gene order rearrangements, strandedness transitions, and correlations with\nlifestyle and locomotory capacity. The work is well-structured, methodologically robust, and\nprovides significant insights into mitogenomic evolution. However, some areas require clarification\nor additional analysis to strengthen the conclusions.\n>>Response: Thank you for investing your time and expertise into in reviewing our study. We\nappreciate your comments and critiques; they helped us improve our study.\nMajor Comments\n1. Ancestral State Reconstruction\n- The inference of the ancestral bilaterian mitogenomic architecture as double-stranded is well-\nsupported, but the instability in \"nad6\" positioning across analyses warrants further discussion.\nCould this instability reflect methodological limitations (e.g., alignment ambiguity or model\nsensitivity)? A deeper exploration of alternative topologies (e.g., Xenambulacraria vs. Nephrozoa)\nand their impact on ancestral state inference would strengthen the conclusions. The authors may\nwant include a sensitivity analysis to assess the robustness of the ancestral state.\n>>Response: We have included the discussion of this problem in the original Supplementary file 1:\nText S3: “MLGO: In combination with the “Nephrozoa” topology, this algorithm produced the architecture\nexhibited by most chordates, including humans (“ur-chordate” architecture henceforth) as the ur-bilaterian\narchitecture, but with all genes encoded on the same strand, i.e. nad6 was inverted from the minus to the plus\nstrand (Figure S1). In more detail, the ur-chordate architecture is characterised by 12 PCGs and 2 rRNAs\nencoded on the “plus” strand, whereas nad6 is encoded on the “minus” strand (opposite transcription\ndirection), within the “nad5 -nad6 cytb rrnS” gene box (“-“ indicates that the gene is encoded on the minus\nstrand). MLGO produced a “nad5 nad6 cytb rrnS” gene box as the ancestral state for Bilateria. This gene (nad6)\nis on the minus strand in almost all Chordata, Hemichordata, Echinodermata, and Xenacoelomorpha\nmitogenomes, aside from a few exceptions in the latter two phyla. Given the putatively key position of\nXenacoelomorpha in this topology, we suspected that it may have skewed the reconstruction. To get a better\nresolution, we additionally analysed all 14 annotated Xenacoelomorpha mitogenomes available in the\nGenBank (Figure S2). They exhibited highly scrambled gene orders, some with all genes on a single strand, so it\nwas almost impossible to infer the ancestral position of nad6. When we tested the performance of MLGO in\ncombination with the “Xenambulacraria” topology, the ancestral chordate architecture characterised by the\n“nad5 -nad6 cytb” gene box was resolved as the ancestral state for Bilateria (Figure S1). This supports our\nsuspicions of an artefact caused by the “basal” position of Xenacoelomorpha. Nad6 was also resolved on the\nminus strand in the ur-Xenambulacraria (Xenacoelomorpha+Ambulacraria) architecture, but it was\ntranslocated to the rrnL -nad6 nad1 gene box. As regards the ur-Ambulacraria (Echinodermata+Hemichordata),\nit was also on the minus strand, but translocated to the cytb -nad6 rrnS gene box (Figure S1). TreeREx: In\ncombination with the “Nephrozoa” topology, this algorithm produced the ancestral chordate architecture\ncharacterised by the “nad5 -nad6 cytb rrnS” gene box as the ancestral state for Nephrozoa. However, nad6\nwas translocated to the “rrnS rrnL -nad6 nad1” in the ur-bilaterian architecture (Figure S1). This gene box (rrnL\n-nad6 nad1) was not observed in the ur-chordate architecture, nor in any of the Hemichordata and\nEchinodermata mitogenomes. Among the six available genera for Xenacoelomorpha, it was found only in\nXenoturbella (Figure S2). As this genus was overrepresented (6/14 mitogenomes in the dataset), we again\nsuspected that this skewed the TreeREx reconstruction (Xenacoelomorpha was the sister-lineage to all other\nbilaterians in this topology). The remaining Xenacoelomorpha species exhibited highly scrambled gene orders,\nfrom which we failed to infer the ancestral arrangement of nad6 for this phylum. When we tested the\nperformance of TreeREx in combination with the “Xenambulacraria” topology, nad6 was translocated between\ncytb and rrnS (“nad5 cytb -nad6 rrnS” gene box). An identical architecture was also produced as the ancestral\nstate for the Xenambulacraria (Xenacoelomorpha+Ambulacraria) and Ambulacraria\n(Echinodermata+Hemichordata) (Figure S1). Therefore, in both cases we rejected results inferred using the\nNephrozoa topology as likely to be artefactual, which leaves us with two candidates for the ancestral bilaterian\narchitecture, both produced using the Xenambulacraria topology, and differing only in the position of nad6: (1)\n“nad5 -nad6 cytb rrnS” (MLGO), or (2) “nad5 cytb -nad6 rrnS” (TreeREx) (Figure S1). As mitogenomic\nrearrangements are considered to be (nearly) a Dollo character, i.e. almost irreversible (Babbucci et al. 2014;\nStruck et al. 2023), the probability of homoplasy in mitogenomic arrangements is statistically very low.\nTherefore, the architecture that is found across more distantly related lineages is the most parsimonious\ncandidate for the ancestral architecture (Boore 2006; Zhang et al. 2019; Struck et al. 2023). The “cytb -nad6”\narrangement was found only in Hemichordata, whereas the “-nad6 cytb” arrangement is present in most\nChordata and Echinodermata. Neither arrangement was observed in Xenacoelomorpha. In addition, “nad6\ncytb” arrangement is present in the putative ancestral architectures of both Ecdysozoa and Lophotrochozoa\n(Protostomia), albeit with nad6 on the plus strand (Webster et al. 2006; Podsiadlowski et al. 2009;\nBernt et al. 2013). On this basis we accepted the “-nad6 cytb” arrangement as the most parsimonious candidate for the\nancestral bilaterian mitogenome. “\nWe moved this discussion to supplementary data as it is very long, and we felt that the ancestral\narchitecture reconstruction was not the key element of our study, as it was done twice before with\nthe same result, i.e. producing the PCGs+rRNAs arrangement identical to the one exhibited by the\nhuman mitogenome as the ancestral bilaterian arrangement (Lavrov and Lang 2005; Oxusoff et al.\n2018). As regards the main objectives of our study, the instability may affect two different analyses:\nstrand distribution of genes and GORR values. The first issue can be safely ignored, as double-\nstranded architecture was inferred as the ancestral state with almost 100% confidence (multiple\ntRNA genes are encoded on the minus strand in the ancestral architecture, but these were not\nincluded in our GORR analyses). As regards the second issue, GORR values, originally, we merely\nlogically assessed that the impact of this instability on the overall GORR dataset should be minimal,\nso we did not pursue this problem further. Following the reviewer’s comment, to confirm that this\ninstability in reconstruction did not affect our results, we inferred GORR values for all bilaterian\nspecies in the dataset using both the accepted ancestral architecture characterised by the “-nad6\ncytb” arrangement (named “ancestral” henceforth), and the alternative ancestral architecture\ncharacterised by the “cytb -nad6” arrangement (named “alternative” henceforth). The two GORR\ndatasets exhibited a high correlation value of 0.94 (p<0.0001). We further confirmed that both\ndatasets reproduce key findings. The correlation coefficient (Pearson’s) between GORR and branch\nlength in Bilateria (hypothesis 4a) differed only marginally between the two GORR datasets: r =\nancestral\n0.69 (p<0.001), and r = 0.65 (p<0.001). The correlation coefficient between GORR and\nalternative\nmitogenome size (hypothesis 10a) also differed only marginally between the two GORR datasets:\nr = -0.39 (p<0.001); r = -0.42 (p<0.001). We also confirmed that GORR produces\nancestral alternative alternative\nthe same pattern regarding in the lifestyle and locomotory capacity categorisation.\nThese results were added to the revised Supplementary file 1: Text S2 and Figure S45, panels A and E\n- GORR (alternative). Therefore, regardless of the instability in the ancestral architecture\nreconstruction, main conclusions of our study remain unchanged. We have slightly rewritten the\ndiscussion to clarify certain points for the revised manuscript: “Our results strongly support the hypothesis\nthat the ancestral bilaterian mitogenome was double-stranded and architecturally very similar to the typical vertebrate\nmitogenome, also found in humans (Lavrov and Lang 2005; Oxusoff et al. 2018), albeit with some instability noted in the\nspecific positions of nad6 and cytb genes. However, alternative ancestral reconstructions did not alter our fundamental\nconclusions (Supplementary file 1: Text S2 and Figure S45).” If the reviewer thinks that some or all of this\ndiscussion should be moved to the main paper, and the editor agrees, we shall happily oblige.\n2. Single- to Double-Stranded Transitions\n- The identification of rare single-to-double-stranded transitions in Annelida and Mollusca challenges\nthe \"ratchet\" hypothesis. However, the evidence for these transitions (e.g., Sabella spalanzani) is\nbased on limited taxa. Could annotation artifacts or incomplete lineage sampling explain these\nobservations?\n>>Response: Some of these issues were discussed in the original submission, and some in its\nSupplementary file 1. In attempts to keep the Discussion not exceedingly long, we decided to keep\nthe detailed explanation of these issues in the revised Supplementary file: Text S7: “As regards switches\nin strandedness, caution needs to be exercised in the identification of these events due to the possible annotation\nartefacts, questionable functionality of some tRNAs undergoing strand switches, and limited representation of some\nlineages. As a result, some of these events were not identified with certainty, and we cannot exclude a possibility that the\ncomplete evolutionary history of bilaterian animals comprises additional state transitions. For example, transitions\ninvolving only tRNA genes may be artefacts caused by recent transitions and loss of functionality of tRNA genes(s)\ntransposed to the minus strand, as we suspect may be the case in Rotifera. However, in the case of subfamilies\nBathymodiolinae and Modiolinae, a majority of mitogenomes exhibited trnT translocated to the minus strand. Such high\nlevel of conservation indicates that the event took place a long time ago, which indirectly indicates that this gene might be\nfunctional, or otherwise its sequence would have degraded over the long evolutionary time. Therefore, this may represent\na true transition in strandedness. It is much easier to make such claims with certainty in cases when transitions comprise\nPCGs, but two of three such cases were limited to a single species. There was only one mitogenome available for Scintilla\nsp., so we cannot exclude an annotation or species identification artefact. Contrary to this, all three available mitogenomes\nof Sabella spalanzani (sequenced by two different teams) exhibited the identical architecture – with 19 genes on the minus\nstrand – so there is strong evidence that this represents a true single-to-double transition. In the third case, all five\navailable Dreissenidae (Bivalvia: Imparidentia: Myida) mitogenomes (nominally three species from two genera) exhibited a\ntransition of multiple PCGs to the minus strand. As our analyses indicate that single-stranded is the ancestral state for\nImparidentia, this represents another solid evidence of a single-to-double transition...”. If the reviewer feels that\nthis should be discussed in the main study, we shall gladly oblige.\n3. Mitogenome Size and Evolutionary Rates\n- The negative correlation between mitogenome size and evolutionary rates contradicts the \"race for\nreplication\" hypothesis. The proposed explanations are plausible but speculative.\n>>Response: In response to this objection, we explained better why gene losses are the most likely\nexplanation for this unexpected relationship and provided some further examples , as well as\nstressed that this does not exclude the possibility that other variables may also play a role: “While\nthese explanations do have a sufficient power to explain the observed patterns, we cannot exclude that other variables,\nsuch as the replication mechanism, repair machinery, recombination, etc., may also play a role.”\n4. Parasitism and Locomotory Capacity\n- While the model linking low locomotory capacity to relaxed purifying selection is compelling, the\nexplanatory power of these variables for GORR variability is modest (6–14%). Other factors (e.g.,\neffective population size, recombination rates) may play significant roles. A discussion of how\nneutral processes might interact with selection might help.\n>>Response: In the early drafts of the manuscript, we included a long discussion of this specific\nproblem. However, as the manuscript was far above the length cut-off required by the journal’s\npolicies, and as the discussion involved several variables that we did not explicitly test for in our\nanalyses (effective population size – Ne, recombination, and replication\nmechanism/machinery/architecture), ultimately, we decided to remove this section from the paper.\nDuring the revision, we included Ne and thermoregulatory type among the tested variables, and\nfound that both had negligible explanatory power (1-2%). We discuss both variables in the revised\nmanuscript and supplementary data. Importantly, we emphasise that these values cannot be\nconsidered low. In addition, we discuss the outlier of Chordata/Deuterostomia, and show that\nexplanatory powers of both variables, but especially locomotory capacity (56%!), are strongly\nincreased when Protostomia are analysed separately: “First, it should be emphasised that these values cannot\nbe considered low. For example, thermoregulatory type and N both explained less more than 2% of the GORR variability in\ne\nour results, whereas body size explained 6% of the branch length variability (but limited to mammals) (Nabholz et al. 2008),\nand generation time explained ≈13% of this variability in both vertebrates and invertebrates (Nabholz et al. 2008;\nThomas et al. 2010). This indicates that additional factors modulate architectural evolution, with candidates comprising selective\nsweeps (Bazin et al. 2006), stabilisation of architectural evolution rates caused by shifts to single-strandedness (discussed\nfurther below). In addition, the predictive power of locomotory capacity and parasitism was much stronger within\nProtostomia alone (≈20% and ≈56% respectively), suggesting that the extreme architectural conservation in\nDeuterostomia, which dominated the full dataset, may be governed by different mechanistic constraints (e.g., replication\nor repair machinery).” We discuss additional details in Text S6.\n5. Population size of mitochondrion may substantially different from the Ne of species.\n>>Response: We assume that the reviewer is referring to the fact that the population size of mtDNA\ndata (mitochondrion) is expected to be much lower (around ¼) than the Ne of the nuclear genome.\nWhile this is absolutely true, this problem needs to be considered only when comparing mtDNA and\nnDNA. As we compared only mtDNA across different bilaterian lineages, all of them are affected by\nthe same phenomenon. This comprises the DUI lineages, where the inheritance is still uniparental,\neven if doubly so. Therefore, this issue should not affect any of our analyses and conclusions. As\nmentioned earlier, Ne analyses were added during the revision, and the impact of Ne on key\nvariables was nonsignificant and negligible.\n---\nMinor Comments\n1. Statistical Analyses\n- The use of phylogeny-corrected models (lmekin, brms) is appropriate, but the R² values for some\ncomparisons (e.g., 3% for GORR≥7/ds) are very low. Emphasize that these results suggest additional\nundetected factors drive architectural evolution.\n>>Response: As mentioned in a response to an earlier question, we discuss the magnitude or R2\nvalues in more detail in the revised manuscript. As regards the specific example, during the revision\nwe managed to identify the reason for the loss of correlation between branch length and GORR in\nsome subdatasets: Mollusca. The corresponding sections were rewritten thoroughly, both in Results\nand Discussion. In the results section (please see hypothesis 4), some analyses were removed, and\nnew ones, clearly showing the effect that Mollusca produce on the relationship between the two\nvariables, were added. In the discussion, we rewrote the related part as follows: “…we did not find\nsupport for the proposal (Bernt et al. 2013) that GORR and branch length correlate only in lineages with high\nrearrangement rates. The correlation was strong across Bilateria, with Mollusca being a notable exception, possibly due to\nits high recombination rate (Ladoukakis and Zouros 2001; Smith et al. 2023), resulting in the decoupling between these two\nvariables. Despite this, locomotory capacity still produced predicted effects on mitogenomic evolution in Mollusca – with\nhigh-locomotory lineages exhibiting shorter branches (Jakovlić et al. 2023) and lower GORR values.” We provide\nadditional details in Text S8: “The correlation between GORR and branch length disappeared in highly rearranged\nmitogenomes with double-stranded architectures (it was found in all other tested combinations). This appears to\ncontradict the results of Bernt et al. (2013), but results are not directly comparable as they studied a much smaller dataset\n(650 vs. 10,860 mitogenomes), and they only statistically tested the correlation on the overall dataset, without testing\nsubdatasets statistically: “Complete shuffling of the mitochondrial genome is clearly correlated with long branches (equals\nhigh substitution rate), while a moderate gene rearrangement (2–6 breakpoints) has virtually no effect.” (Bernt et al.\n2013). Notably, in the overall dataset, they obtained an identical value to the one found in our analyses (r = 0.69).” We\nalso highlighted in multiple places throughout the discussion that other factors seem to be\ncontributing to the architectural evolution. A potential candidate for other variables that\ndrive/determine the rate of architectural rearrangement rate that emerged during this revision is\nthe possibility that evolutionary transition to single-strandedness results in a stabilisation of\narchitecture, and there are indications that Deuterostomia may employ a different molecular\nmachinery, resulting in low architectural evolution rates.\n2. Outliers\n- The exceptional conservation in Cirripedia (Hexanauplia) is intriguing. The hypothesis of large effective population size is\nuntested. Could life-history traits (e.g., sessility, larval dispersal) also contribute?\n>>Response: Indeed, this is an intriguing outlier. When we first observed it in the study of skews in\ncrustaceans (Jakovlić et al. 2021a) - we feared that it may invalidate our entire working hypothesis.\nIn that light, it was comforting to discover in subsequent studies that it is not a major outlier only\namong the crustaceans, but also among all bilaterians. Larval dispersal is closely associated with our\nhypothesis of large Ne (effective population size), as larval dispersal should be positively associated\nwith Ne. If high larval dispersal rates result in high Ne values, along with high density of populations\n(locally), it may have the power to explain the absence of the expected relaxation of purifying\nselection in association with the transition to the sessile lifestyle. For this revision, we added the\nanalysis of Ne data recently compiled by Clarke et al. (Clarke et al. 2024), but there were no\nCirripedia (Hexanauplia) included in it, so this remains hypothetical. Also, we should emphasise that\nwe cannot reject the possibility that some other variable, such as more efficient mtDNA repair\nmachinery or a different replication mechanism, may contribute to this perplexing outlier. Another\nfactor that may hypothetically play a role is the duration of the larval stage: often short in parasitic\nlineages, and long in some (but not all) barnacles. This prompts a hypothesis that barnacles undergo\na strong selection for locomotory capacity during their long larval stage, whereas (most) parasites\ndon’t. As it remains speculative, we did not discuss it in the paper. In attempts to shorten the\nDiscussion, we did not go into many speculations in the main text: “An outlier where selection for\nlocomotory capacity failed to explain the pattern is Hexanauplia (Arthropoda): low-locomotory species exhibited lower\nGORR values and branch lengths (Jakovlić et al. 2023). This class comprises sessile crustaceans - barnacles (subclass\nCirripedia), which have been identified as a major outlier in the relationship between locomotory capacity and\nmitogenomic evolution among all bilaterian animals (Jakovlić et al. 2021a, 2023). For reasons that we still fail to fully\nexplain, despite their predominantly sessile lifestyle, mitogenomes of this lineage are evolving under some of the most\nstringent purifying selection pressures among all Arthropoda .” However, we expanded the discussion of this\nproblem in the Supplementary file 1: Text S8: “Hexanauplia (Arthropoda): low-locomotory species exhibited low\nGORR values (Figure S49 and S56) and branch lengths (Jakovlić et al. 2023). This can be attributed to the unique evolution\nof barnacles (subclass Cirripedia; class Hexanauplia in NCBI’s taxonomy database, and Thecostraca in WORMS). There are\nnow multiple lines of evidence across several studies (Jakovlić et al. 2021a, 2023) indicating that mitogenomes of this\nlineage have some of the slowest evolutionary rates among all Arthropoda (Figure S61-S62) and they are a major outlier\namong the Bilateria as well in relation to the expected relationship between the locomotory capacity and mitogenomic\nevolutionary rate. This indicates that they are undergoing high levels of purifying selection, or possibly exceptionally low\nmutational pressures. As mutational pressures are directly positively associated with the time that the lagging strand\nspends in a single-stranded state, replication mechanisms vary across different bilaterian lineages (Lewis et al. 2015;\nCiesielski et al. 2016; Yasukawa and Kang 2018), we should account for a possibility that Thecostraca/Hexanauplia may\nhave strongly reduced metabolic rates or even employ a replication/repair mechanism that produces lower mutational\nrates. Another factor that may explain the outlier is variability in N . Previously, we observed that N values inferred using\ne e\ncox1 data indicate that sessile crustaceans (Thecostraca/Hexanauplia) exhibited a non-significantly increased average N\ne\nvalue among the studied crustaceans (Jakovlić et al. 2021a). A factor that may explain their high N is high larval dispersal\ne\n(Minchinton and Scheibling 1991). Another factor that may contribute to increased purifying selection pressures are their\nhigh population densities, known to cause massive density-dependent mortality in some circumstances (Bertness 1989).\nElevated N and competition for high metabolic efficiency may produce stronger purifying selection pressure than\ne\npredicted by the locomotory capacity hypothesis, resulting in lower-than-expected evolutionary rates. Previously, we also\nfound that sessile species exhibited the lowest ω values in the entire crustacean mitogenomic dataset and a wide range of\nGC skews (Jakovlić et al. 2021a). There was also some evidence for a correlation between the GC skew and N variables\ne\n(r(20) = –0.26, p = 0.238)) and between the ω (i.e. dN/dS) and N variables (r(20) = –0.346, p = 0.114)) in that dataset\ne\n(Jakovlić et al. 2021a). It remains unclear whether this applies to all sessile species, and whether it is sufficient to explain\nthis outlier, so other explanations should be sought as well. More data are needed to clarify the unique evolutionary\ndynamics of this crustacean taxon.“\n3. Supplementary Material\n- Not all supplementary figures/tables are directly cited in the main text.\n>>Response: We thoroughly checked the manuscript and made sure that all supplementary\nelements are cited. Please note that in some cases we cite multiple figures, e.g. Figures S9-S12, so\nfigures S10 and S11 cannot be identified using a simple text search tool.\nReviewer #3 (Remarks to the Author):\nI co-reviewed this manuscript with one of the reviewers who provided the listed reports. This is part\nof the Nature Communications initiative to facilitate training in peer review and to provide\nappropriate recognition for Early Career Researchers who co-review manuscripts.\nReviewer #3 (Remarks on code availability):\nI did not run the code, I only reviewed it.\nThe R code to create the figures and tables is provided in Zenodo . The R code is mostly intuitive and\nallows the figures to be reproduced, however a README is not provided with a description of each\nscript. This is not necessary but listing out the required dependencies would help others run the\ncode with new data.\n>>Response: Thank you for the helpful suggestion. We have added a new file,\nreadme_for_codes.txt, which explains how to run each script to generate the corresponding figures\nand tables. This file also lists all required dependencies along with instructions for installing them. In\naddition, we have expanded the code annotations within the individual scripts to clarify the key\nsteps. The updated scripts and documentation have been uploaded to Zenodo:\nhttps://doi.org/10.5281/zenodo.17809659.\nThe input files should be explained further in the manuscript. For example, how were organisms\ncategorized into lifestyles within the input spreadsheet.\n>>Response: As explained in much more detail elsewhere in this response letter, to avoid repetition,\nwe opted to refer readers to our previous Nature Communications study, where the categorisation\nwas introduced in detail. Following the reviewers’ objections, we included a detailed discussion of\ncategorisation to the supplementary data of the revised paper (Text S9).\nThank you for your comments."}, {"round": 2, "reviewer_comments": {"reviewer_1": "The authors sufficiently addressed the concerns discussed in the review, adding substantial changes to the text and supplementary tables to better describe the biological hypotheses and novelty of this paper compared to their recent publication.\nWe also appreciate that the authors provided additional analyses and explanations related to GC skew, effective population size, thermoregulatory type, and life history categories, which we feel were important to include in this manuscript. We recommend this revised manuscript for publication.", "reviewer_2": "Besides the font size in the main figures are too small, I have no more further comments.", "reviewer_3": "I co-reviewed this manuscript with one of the reviewers who provided the listed reports. This is part of the Nature Communications initiative to facilitate training in peer review and to provide appropriate recognition for Early Career Researchers who co-review manuscripts.\nWe appreciate that the authors took our suggestion to add a README to the Zenodo to aid in reproducing the results and making the code usable for the community."}, "author_response": ">>R: Thank you for investing your time and expertise into reviewing our manuscript, as well as your positive assessment or our revisions.\n>>R: We have revised all figures and increased font as much as it was possible during this round of revision. Thank you for reviewing our manuscript.\n>>R: We thank both reviewers for investing their time and expertise into reviewing our manuscript."}]

Decision letter