Contrasting Modes of Mitochondrial Genome Evolution in Sister Taxa of Wood-Eating Marine Bivalves (Teredinidae and Xylophagaidae)

Abstract The bivalve families Teredinidae and Xylophagaidae include voracious consumers of wood in shallow-water and deep-water marine environments, respectively. The taxa are sister clades whose members consume wood as food with the aid of intracellular cellulolytic endosymbionts housed in their gills. This combination of adaptations is found in no other group of animals and was likely present in the common ancestor of both families. Despite these commonalities, the two families have followed dramatically different evolutionary paths with respect to anatomy, life history, and distribution. Here, we present 42 new mitochondrial genome sequences from Teredinidae and Xylophagaidae and show that distinct trajectories have also occurred in the evolution and organization of their mitochondrial genomes. Teredinidae display significantly greater rates of amino acid substitution but absolute conservation of protein-coding gene order, whereas Xylophagaidae display significantly less amino acid change but have undergone numerous and diverse changes in genome organization since their divergence from a common ancestor. As with many bivalves, these mitochondrial genomes encode 2 ribosomal RNAs, 12 protein-coding genes, and 22 tRNAs; atp8 was not detected. We further show that their phylogeny, as inferred from amino acid sequences of 12 concatenated mitochondrial protein-coding genes, is largely congruent with those inferred from their nuclear genomes based on 18S and 28S ribosomal RNA sequences. Our results provide a robust phylogenetic framework to explore the tempo and mode of mitochondrial genome evolution and offer directions for future phylogenetic and taxonomic studies of wood-boring bivalves.


Introduction
Teredinidae and Xylophagaidae ( fig. 1) are the sole members of the molluscan class Bivalvia capable of consuming wood as food, or xylotrophy. Indeed, all but two species of Teredinidae and all species of Xylophagaidae burrow nearly exclusively in wood or woody plant materials and are thought to utilize vascular plant cell wall material (lignocellulose) as a substantial source of dietary carbon (Distel 2003;Turner 1967;Voight 2015;Nishimoto et al. 2021). Moreover, all members of both families are thought to harbor intracellular cellulolytic gammaproteobacteria in their gill tissues, a feature that is unique among all animals. Teredinidae are commonly called "shipworms" because of their worm-like morphology and historical legacy of destroying unprotected wooden ships. Even today, shipworms account for billions of US dollars in damages to man-made wooden structures in marine and brackish environments (Distel 2003). Shipworms were well-known to the earliest seafarers and have altered the course of human civilizations (Stearns 1886;Distel 2003), having been implicated in the defeat of the Spanish Armada and in the disastrous conclusion of the fourth voyage of Christopher Columbus (Rayes et al. 2015). Similarly, the burrowing activity of xylophagaids is reported to have caused the failure of early submarine telecommunications cables by damaging their insulating sheaths (Jeffries 1861;Parkes and Keeble 2016) and even today remains a concern for the design of subsea umbilicals (Parkes and Keeble 2016). Despite these negative impacts, wood-boring bivalves also play beneficial economic and ecological roles in many marine environments by converting recalcitrant wood into animal biomass that is more easily consumed by a wide range of organisms (Petra Pop et al. 2017;Cragg et al. 2020). Shipworms are also considered a delicacy in several traditional cuisines (Turner 1971) and have been proposed to have economic potential as a shellfishery (Willer and Aldridge 2020).
Despite their common xylotrophic lifestyle, these two families inhabit distinct and nearly nonoverlapping distributions wherein they are exposed to very different environmental conditions. Teredinidae are common in tropical and temperate waters but are largely absent from polar zones. With a few notable exceptions, for example (Velasquez and Shipway 2018), they occur at shallow depths, typically less than 200 m. They are most common in near-shore and intertidal, mangrove, estuarine, and coastal riparian environments where floating and deposited driftwood and submerged roots and branches provide abundant and consistent sources of food and shelter. Within these environments many teredinids experience and tolerate broad ranges of salinity and temperature with strong seasonal, diurnal, and weather-driven variation (Voight 2015). Because they often inhabit floating wood and wood deposited in the intertidal zone, they may frequently experience prolonged periods of air exposure due to tidal recession. They survive such exposure by sealing their burrows, which are lined with an impervious shell-like calcareous material, using a paired set of calcareous plates called pallets. Pallets and shell-lined burrows are common to Teredinidae but are absent in Xylophagaidae, with the exception of the genus Xyloredo in which the distal end of the burrow is lined (Voight et al. 2019). Although these adaptations may prevent desiccation and death during atmospheric exposure, they introduce other significant stresses such as extended periods of anaerobiosis and accumulation of metabolic waste.
Xylophagaidae, on the other hand, are found in sunken wood on the sea floor most commonly at depths greater than 150-200 m, with a few species coming into shallow water in boreal and high latitudes (Voight 2015;Romano et al. 2020). Typically, they are exposed to the consistent low temperatures and narrow salinity ranges characteristic of deep-sea environments and experience substantially less seasonal diurnal, and storm driven variation than is Significance Among metazoans, bivalves display unusual variation in mitochondrial genome evolution, organization, and inheritance, but the factors that influence this evolutionary lability are poorly understood. Here, we present 42 new mitochondrial genome sequences from bivalve sister clades that despite their close phylogenetic relationship and uniquely shared wood-feeding habits, differ dramatically in morphology, life history, and distribution. We show that these differences are correlated with equally dramatic differences in tempo and mode of mitochondrial genome evolution, laying the groundwork for improved understanding of the complex evolutionary interactions between marine organisms, environments, and mitochondrial genomes.
commonly encountered by teredinids in shallow water. With few exceptions, for example (Jayachandran et al. 2022), Xylophagaidae have not been reported to inhabit floating wood and have only rarely been observed within the intertidal zone (Turner 2002;Voight 2016). Thus, they are typically not exposed to the atmosphere, lack adaptations for sealing their burrows, and do not typically experience the types of anaerobic and metabolic stresses encountered by Teredinidae during prolonged periods of burrow closure.
The two families also differ with respect to reproductive strategies. Although little is known about the reproductive biology of Xylophagaidae, most are thought to be broadcast spawners with planktotrophic larvae. Although many species were once thought to brood their young, small specimens commonly found within the female's burrows have more recently been shown to be dwarf males (Haga and Kase 2013) rather than larvae. Male dwarfism is often associated with sparse food availability and spatial limitation for growth (Haga and Kase 2013), features consistent with deep-water environments where wood, which serves both as food and habitat for Xylophagaidae, is scarce. In contrast, Teredinidae employ a wide range of reproductive strategies, including pseudocopulation (Shipway et al. 2020), broadcast spawning with maximized production of gametes, and larval brooding with fewer offspring and extended parental care. In most cases, these are adaptations suited to rapid utilization of a comparatively abundant but patchily distributed wood supply (MacIntosh et al. 2014). Male dwarfism is found in only one described teredinid species, Zachsia zenkewitschi, which inhabits spatially restrictive seagrass rhizomes (Shipway et al. 2016).
Finally, the two families differ substantially in anatomy and morphology. For the most part, Xylophagaidae conform to a typical bivalve body plan, wherein the internal organs are located between the anterior and posterior adductor muscles and the entire body, with the exception of the siphons in some cases, can be retracted between the shells (Voight 2015). In contrast, Teredinidae are among the most morphologically divergent bivalves (Turner 1966). During development, the body elongates and the heart, gills, gonads, and a substantial portion of the digestive system assume a position posterior to the posterior adductor muscle, bringing them permanently outside the protection of the valves.
Despite these substantial differences (summarized in fig. 2), investigations based on 18S and 28S rRNA sequences support Teredinidae and Xylophagaidae as sister families (Voight 2015;Voight et al. 2019;Romano et al. 2020), agreeing with the earlier conclusions of Purchon (Purchon 1941) and supplanting the later placement of the deep-sea wood borers as a subfamily of Pholadidae (Turner 1967(Turner , 2002. Although no formal taxonomic revision has yet been published, this view has gained wide acceptance in the literature (Voight 2015;Voight et al. 2019;Romano et al. 2020). Nonetheless, phylogenetic relationships within the two families remain unresolved.
Mitochondrial genome data have proven useful in resolving phylogenetic relationships across a wide range of metazoans, for example (Miya et al. 2001;Osigus et al. 2013;Li, Kocot, et al. 2015) including mollusks, for example (Grande et al. 2008;Mikkelsen et al. 2018;Kong et al. 2020). Although, mollusks are among the most variable of bilaterian animals with respect to mitochondrial genome size and organization (Piccinini et al. 2021), most still conform to the canonical bilaterian mitochondrial genome complement of 37 genes including 22 tRNAs, 12-13 protein-coding genes, and 2 rRNA genes. ATP synthase subunit 8 is not detected and may be absent in a number of bivalve taxa (Ghiselli et al. 2021). The most obvious variations in length are due to the presence of noncoding regions (Formaggioni et al. 2021), which, in some bivalves, account for mitochondrial genomes 2-3 times of that of the typical 15-17 kb bilaterian size (Smith and Snyder   Williams et al. 2017;Kong et al. 2020). Gene order, orientation, and location with respect to the heavy and light strands also show considerable variation in mollusks (Piccinini et al. 2021). However, the degree of variation in genome length, gene order, and strandedness, as well as nucleotide composition, appears to be largely cladespecific. In the case of the order Myida, which includes economically important invasive and wood-boring species, genomic resources are still scarce in comparison to other major bivalve clades. Prior to this study, only one complete mitochondrial genome assembly had been published for Myida in the NCBI Reference Sequence Database (O'Leary et al. 2016;Wilson et al. 2016) with no representation of the superfamily Pholadomyoidea, which includes most wood-boring species.
Because Xylophagaidae and Teredinidae present an unusual combination of recent common ancestry and unique shared feeding strategy, juxtaposed with profound differences in morphology, reproductive strategies, life history, and bathymetric distribution, we ask whether similarly stark contrasts are also observed in their mitochondrial genome organization and evolution. To this end, we explore the mitochondrial genomes from a variety of taxa representing both Teredinidae and Xylophagaidae.

Results and Discussion
Mitochondrial genome sequences were determined for 42 bivalve specimens, including 26 Teredinidae from 15 locations and 16 Xylophagaidae from 8 locations (table 1). A single unique mitochondrial genome was recovered from each specimen examined. Although no evidence of distinct sex-specific mitochondrial lineages was detected, for example, doubly uniparental inheritance, this cannot be definitively ruled out as sex of individual specimens was not determined and reproductive organs were not specifically sampled.
All mitochondrial genomes described herein contain, at minimum, 2 ribosomal RNA (rRNA) genes, 22 tRNA genes, and 12 of the 13 canonical protein-coding genes commonly present in animal mitogenomes; atp8 was not detected. All genes are encoded on the same strand and share the same orientation. The genomes range in size from 14,450 to 18,624 bp (table 2), similar to most bilaterian mtDNA genomes. No significant difference was detected in mitochondrial genome size between the two families. Much of the observed length variation is associated with noncoding regions, which on average comprised a significantly larger fraction of total mitochondrial genome in Xylophagaidae than in Teredinidae (supplementary fig. S1, Supplementary Material online). However, significant differences were observed in GC-content and in evolutionary rates between teredinids and xylophagaids using paired Nototeredo knoxi for which bootstrap support was low in all three analyses. The site homogeneous model applied to the unpartitioned supermatrix produced a tree most similar to previously published phylogenetic analyses based on large (28S) and small (18S) subunit nuclear-encoded rRNA sequences and the cytochrome c oxidase I gene (Distel et al. 2011;Voight et al. 2019;Romano et al. 2020). For example, in agreement with previous reports, the tree presented in figure 4 supports the nonmonophyly of the Teredinidae subfamilies Teredininae and Bankiinae, showing Lyrodus and Teredo (Teredininae) nested within Bankiinae, which is in turn nested within a clade containing Neoteredo and Bactronophorus (Teredininae). Additionally, in agreement with previously published results, these analyses support the monophyly of Xyloredo and of Xylophaga as well as the divergence between these and Xylonora, a genus recently erected to remedy the nonmonophyly of Xylophaga (Voight et al. 2019;Romano et al. 2020). However, the monophyly of Xylonora is not well supported in our analyses. The most notable difference between the analyses presented here and those reported previously (Distel et al. 2011;Borges et al. 2022) is the basal position of the node connecting Kuphus (Kuphiinae) to other Teredinidae, which was previously reported to be nested within Teredininae. As Kuphus is arguably the most morphologically derived and physiologically distinct member of Teredinidae, as well as one of only two members of the family identified to date that harbor chemoautotrophic endosymbionts (Distel et al. 2017;Altamia et al. 2020), the phylogenetic position of Kuphiinae is important with respect to understanding the origins of symbiosis in teredinids. The tree is also consistent with the previously reported sister relationship between the families Teredinidae and Xylophagaidae (Distel et al. 2011). We caution, however, that the analyses presented here lack sufficient taxonomic representation both within the two families and among their closest relatives to confidently resolve the basal branching order of either Teredinidae or Xylophagaidae, or to provide additional support for the proposed sister relationship between these taxa. Based on the remarkable economic impact, ecological importance, historical influence, and potential biotechnological and medicinal

Gene Order and Copy Number
Surprisingly, despite their close phylogenetic relationship, the two families differ dramatically with respect to conservation of gene order. Among all Teredinidae examined, the order of protein-coding genes in mitochondrial genomes is conserved ( fig. 4), but variation was observed among tRNA positions (supplementary fig. S6, Supplementary Material online). In 7 of 14 taxa, 1-4 tRNA genes (trn L1, M, N, Q, R, and V ) differ in position from the consensus order. This is consistent with the observation that rearrangements and duplications involving tRNAs occur more frequently than those observed among protein-coding genes (Ghiselli et al. 2021). Interestingly, in one species, Neoteredo reynei, there is an apparent duplication of trnW and its consensus neighbor, cox3. E values and quality scores determined by MITOS2 suggest the degradation of one copy of each duplicated gene.
In sharp contrast, six distinct and highly divergent patterns of protein-coding gene order were observed among eight putative xylophagaid species examined ( fig. 4). In addition, variations in tRNA gene order are far more common in Xylophagaidae than in Teredinidae, as are gene duplications (supplementary fig. S6, Supplementary Material online). Indeed, differences in tRNA gene order were apparent between all putative species examined and were observed even within a single nominal species, Xylophaga oregona. Apparent tRNA gene duplications were observed in three of eight putative species examined, involving trnS1, H, D, Y, Q, and M, with putative copy numbers ranging from 2 to 5 per genome. Similarly multiple tRNA gene copies have been reported on other bivalve mitochondrial genomes, for example, in bathymodiolin mussels (Zhang et al. 2021). Moreover, in one putative species, represented by xylophagaid specimens E23 and E81, there is an apparent tandem duplication of the rrnS and trnM genes. Although tandem duplication is thought to play an important role in mitochondrial genome rearrangement, few animal mitochondrial genomes contain duplicated copies of proteincoding and ribosomal genes, likely due to selection for the maintenance of cytonuclear stoichiometry (Ghiselli et al. 2021). Thus, the evidence of gene duplications presented here warrant further investigation.
Interestingly, among all pairwise comparisons within and between the two families, the most similar arrangements of mitochondrial protein-coding genes are observed between all Teredinidae and Xylophaga washingtona, X. oregona and three of five Xyloredo mitogenomes examined ( fig. 4). These two protein-coding gene arrangements differ only by the position of the cox2 gene, and each can be transformed into the other with a single transposition. This observation, combined with the complete conservation of protein-coding gene order in Teredinidae, suggest the parsimonious conclusion that a similar protein-coding gene arrangement likely occurred in the most recent common ancestor of the two families. Based on this assumption, credible scenarios for evolution of gene order can be inferred by analysis of common intervals using Crex (Bernt et al. 2007). For simplicity, the scenario depicted in figure 4 considers only protein-coding gene order and proposes an ancestral gene order identical to that found in all Teredinidae. Although resolving the complex history of mitochondrial gene rearrangement in wood-eating taxa requires additional taxon sampling and is beyond the scope of this investigation, the scenario presented in figure 4 serves to demonstrate the large number and diverse types of rearrangements required to reconcile the numerous and divergent gene orders observed in this investigation. Notably, protein-coding gene rearrangements are found even among closely related members of the genus Xyloredo, with some members sharing the same gene order as X. washingtona and X. oregona and others displaying a highly divergent gene order.
In addition to greater conservation of gene order, species of Teredinidae also displayed greater conservation with respect to the predicted locations of replication origins (table 3). Within Teredinidae, the predicted locations of both the heavy-strand origins (OriH) and light-strand origins (OriL) were conserved with respect to protein-coding genes but not with respect to tRNA genes. In all species the predicted OriH was located between cob and rrnL. However, the predicted OriH locations differed with respect to the locations of trnM and trnC which, in all but one species also fall between these same protein-coding genes. The predicted OriH in all but three taxa (Bactronophorus thoracities, Dicyathifer mannii and Lithoredo abatanica), falls downstream of cob but precedes trnM and trnC. In B. thoracities, trnM is transposed to a position upstream of cob and the OriH falls between cob and trnC. In L. abatanica the predicted OriH falls beween trnM and trnC, and in D. mannii it falls between trnM and rrnL. Similarly, the predicted light-strand origin (OriL) was located between cox1 and nad4 in all species but differed with respect to the position of trnT, which also falls between these two protein-coding genes. The OriL is predicted to fall in the cox1-trnT IGR in all teredinid mitogenomes except D. mannii and L. abatanica in which it falls downstream of trnT. Thus, the locations of both origins of replication appear to be invariant with respect to protein-coding genes but variable with respect to tRNA encoding genes.
In contrast, among the eight species Xylophagaidae examined, six distinct locations were predicted for the OriH with respect to flanking protein-coding genes. Interestingly, among xylophagaid species that share the same protein-coding gene order (Xylophaga oregona, X. washingtona and Xyloredo nooi) the predicted OriH locations were conserved with respect to flanking protein-coding genes. The position of OriL could not be confidently predicted in the examined xylophagaid species.

Amino Acid Substitution Rates
In contrast to the rapid evolution of gene order in Xylophagaidae as compared Teredinidae, the estimated average rate of amino acid substitution (1.53 substitutions/site), determined as a function of tip-to-root branch lengths in the concatenated supermatrix tree (figure 4), is significantly lower (P = 6.7 × 10 −9 , t-test) than in Teredinidae (2.18 substitutions/site; fig. 3). This constitutes a difference of approximately 1.4-fold. The higher rates of amino acid substitution in Teredinidae may be tied to adaptation to anaerobic stress. Unlike xylophagaids which are adapted to deep-sea environments where oxygen concentrations are comparatively consistent  3.-Violin plots comparing genome size, GC-content, and evolutionary (amino acid substitution) rate between Teredinidae and Xylophagaidae. Statistical significance of differences between Teredinidae and Xylophagaidae with respect to each genome property was evaluated using paired t-tests with P values adjusted by Bonferroni methods (P values shown above). (a) Genome size; (b) GC content; and (c) evolutionary (amino acid substitution) rates as measured by the tip-to-root distance. Note that GC-content and amino acid substitution rate, but not genome size, are significantly different at a P value threshold of <0.01.
FIG. 4.-Phylogenetic relationships and variation in mitochondrial gene arrangement among species of Xylophagaidae and Teredinidae. Phylogenetic tree inferred by ML (single-site homogeneous model, unpartitioned, IQ-Tree 1.6.7) based on 4,135 unambiguously aligned amino acid positions selected using Gblocks from an alignment of 12 concatenated mitochondrial protein-coding genes. Bootstrap proportions, less than 100 are indicated at the nodes; where no numerical value is indicated, the bootstrap proportion = 100. Protein-coding gene arrangements associated with each species are depicted on the right. Putative origins of replication are indicated by thick vertical bars. A proposed scenario for the evolutionary history of gene rearrangement, determined by analysis of common intervals using CRex, is superimposed on the phylogenetic tree with branches color-coded according to extant and hypothetical ancestral gene orders. Symbols indicate rearrangements required to transform each hypothetical ancestral gene order to the order observed in the extant species. "/" indicates transposition, "\" indicates reverse transposition, "<" indicates reversal, "*" indicates tandem duplication with random loss (TDRL). and predictable, teredinids are adapted to shallow water, intertidal environments, and floating wood, where exposure to the atmosphere and subsequent desiccation is a significant threat. To combat desiccation, teredinids seal their burrows with their paired pallets, allowing some species to survive for weeks in wood removed from water, but also exposing them to extended periods of anoxia (Lane et al. 1955).
Also interesting is the observation that within Teredinidae, the average amino acid substitution rate (2.34 substitutions/site) estimated for the clade containing the genera Bactronophorus, Neoteredo, Bankia, Teredo and Lyrodus is significantly higher (P = 3.5 × 10 −8 , Kruskal-Wallis H-test) than that estimated for the remaining teredinid mitogenomes (1.84 substitutions/ site) (supplementary fig. S7, Supplementary Material online). Several unusual reproductive strategies are observed within this clade. Although the details of reproductive behavior are not known in Bactronophorus and Neoteredo, internal fertilization is observed in Bankia, Teredo, and Lyrodus. In Bankia species, pseudocopulation is known to occur, in which sperm is delivered from the exhalent siphon of one individual to the inhalant siphon of a neighbor (Shipway et al. 2020). After internal fertilization and before release to the environment, species of Teredo and Lyrodus brood their young to late larval stages in pouches located on the dorsal side of the gill. Pseudocopulation and larval brooding favor more localized reproduction and more limited larval dispersal. These taxa are also generalists that inhabit many types of floating wood, as opposed to other teredinids that tend to inhabit less mobile habitats such as mangrove roots, sulfidic sediments (Distel et al. 2017) and limestone riverbanks (Shipway, Distel, et al. 2019;. Finally, all of these taxa are hermaphroditic and at least some are known to be capable of self-fertilization (Eckelbarger and Reish 1972). These strategies are likely adaptive for utilization of wood that is often patchily distributed in marine environments. They promote efficient localized settlement on and consumption of floating wood that may be rafted great distances from other wood sources (Treneman, Borges, et al. 2018;Treneman, Carlton, et al. 2018). Additionally, these strategies allow few or even single larvae to initiate new populations on wood islands. However, these same adaptations may also result in genetic bottlenecks resulting from the frequent transport and isolation of small founding populations, which in turn can promote the fixation of mutations. Indeed, this notion is supported by the observed prevalence of cryptic species in some of these taxa (Borges and Merckelbach 2018).
Unfortunately, relatively little is known about reproductive strategies in Xylophagaidae. However, a recent study used dynamic energy budget modeling to propose a different strategy for adaptation to life on sparse and patchy wood islands in Xylonora atlantica, involving rapid sexual maturation at small adult body size, a long-lived larval dispersal stage and high larval survival rates (Gaudron et al. 2021). If common to other Xylophagaidae, this strategy might be less prone to geographic isolation and genetic bottlenecks, leading to greater effective population sizes and more efficient purifying selection.

Conclusions
Identifying causal relationships between the functional and physiological adaptations of taxa and the observed differences in their mitochondrial genome organization and evolution is notoriously difficult (Ghiselli and Milani 2020;Ghiselli et al. 2021). Here, we show that two bivalve sister clades, Teredinidae and Xylophagaidae, have followed very different evolutionary trajectories, one leading to thriving existence in comparatively warm productive surface waters and diverse coastal habitats and the second to considerable success on the comparatively cold, relatively uniform and oligotrophic deep-sea floor. Teredinidae includes very diverse species that are known to thrive in a broad and variable range of environmental conditions and that span a wide range of life histories, feeding and reproductive strategies, and physiologies. Unfortunately, the parallel properties with respect to Xylophagaidae are less well-known, making generalization more difficult. Although we can infer differences based on their distinct respiratory demands, reproductive strategies, bathymetric ranges and the environmental variables that characterize these distributions, distinguishing between these alternatives, or other explanations will require more taxon sampling and better understanding of wood-boring bivalve phylogeny and ecology. Perhaps the most interesting feature of this data set is the contrasting patterns observed in the two families with respect to rates of sequence evolution and rates of genome rearrangement, phenomena which previously have been proposed to be positively correlated (Mortz et al. 2021). While we cannot currently explain these differences, the presented data provide fertile ground for exploration of the environmental, biological, and molecular mechanisms that shape the tempo and mode of mitochondrial evolution.

Specimen Collection and DNA Extraction
Specimens were collected by a variety of methods ranging from collection by hand, to dredging, and deployment and recovery of wooden substrates (table 1). Specimens were frozen at −80°C or preserved in 80-100% nondenatured ethanol following collection unless specified otherwise. For Xylophagaidae, siphon tissue was dissected from each bivalve, and total genomic DNA was extracted using the DNeasy Blood & Tissue Kit (Qiagen) according to the manufacturer's protocols. Tissue samples collected in Alabama and Florida were preserved in 0.25 M EDTA, pH 8.0 (Sharpe et al. 2020 (Joshi and Fass 2011), merged, and converted to FASTA files. Merged FASTA files were assembled using IDBA_ud v2.0 (Peng et al. 2012) using default parameters. Teredo sp. and Bankia sp. gill metagenomes were sequenced using Illumina MiSeq. The raw reads were assembled using either the metaspades pipeline of SPAdes (version 3.11.1) (Bankevich et al. 2012) or IDBA-UD (version 2) (Peng et al. 2012). Before assembly, raw reads were merged using BBMerge (v9.02) (Bushnell et al. 2017). Nonmerged reads were filtered and trimmed using FaQCs (Version 1.34) (Wang et al. 2014). Mitochondrial contigs were identified using TBLASTN (Altschul et al. 1997) and the previously published bivalve mitochondrial genome from Mya arenaria (Wilson et al. 2016) as bait. For T. bartschi, whole genome shotgun (WGS) sequencing was conducted at the New York Genome Center on an Illumina HiSeqX (2 × 150 bp). Library preparation utilized a TruSeq PCR-free kit (450 bp). The mitochondrial genome was bioinformatically extracted from the WGS run using MitoFinder v1.4 (Allio et al. 2020). MitoFinder employed MEGAHIT (Li, Liu, et al. 2015) for mitogenome assembly and tRNAscan-SE (Chan and Lowe 2019) for tRNA annotation. Annotation of the mitochondrial genomes was conducted initially with MITOS2 web server  with default settings and the invertebrate genetic code (i.e., NCBI translation table 5) for mitochondria, followed by manual genome annotation of start and stop positions of each gene using Artemis (Rutherford et al. 2000). The nucleotide base composition across the complete mitochondrial genome, protein-coding gene sequences and 3rd codon position were calculated using Artemis.

Phylogenetic Methods
Forty-two specimens were included in the phylogenetic analysis (table 1), including 26 Teredinidae and 16 Xylophagaidae. M. arenaria (Myida), Meretrix petechialis (Venerida) and Solecurtus divaricatus (Cardiida) were selected as the outgroup for our phylogenetic analysis based on availability in the NCBI RefSeq (O'Leary et al. 2016) database and the current understanding of bivalve evolutionary relationships (Combosch et al. 2017;Lemer et al. 2019). Prior to alignment, nucleotide sequences were translated to amino acids using invertebrate mitochondrial genetic translation code as implemented in Artemis. Each proteincoding gene was individually aligned in MAFFT v7.2.3 (Katoh and Standley 2013) followed by manual correction. Each gene was trimmed using Gblocks (Talavera and Castresana 2007) to discard ambiguously aligned sites with default parameters. Protein-coding gene alignments were then concatenated into a final supermatrix dataset using FASconCAT (Kuck and Meusemann 2010).
For the supermatrix constructed above, we used three different approaches to infer the shipworm phylogeny: 1) the concatenation (i.e., supermatrix) approach with a single site-homogeneous model or partition, 2) the concatenation approach with data-partitioning by gene, and 3) the concatenation approach with a site-heterogenous C20 model to account for amino acid compositional heterogeneity. Maximum-likelihood (ML) phylogenetic analysis was conducted using IQ-Tree 1.6.7 (Nguyen et al. 2015). Prior to ML analyses, ModelFinder (Kalyaanamoorthy et al. 2017) was used to evaluate best substitution models for each gene partition. Nodal support for ML analyses was evaluated with 1,000 ultrafast bootstrapping replicates. Note that only bootstrap support ≥95 should be considered as a strong support for a given bipartition (Minh et al. 2013).

Mitochondrial Genome Properties
We focused our analyses on the comparisons of five mitochondrial genome properties between Xylophagaidae and Teredinidae. Specifically, for a given taxon, 1) amino acid substitution rate was estimated as the distance from the most recent common ancestor shared by Xylophagaidae and Teredinidae to each tip on the concatenation-based supermatrix ML tree; 2) GC-content was calculated as the percentage of G and C nucleotides in the complete genome; 3) gene order was determined based on analysis of gene annotations predicted using MITOS 2.0 (Donath et al. 2019); 4) genome size was determined as the number of base pairs per genome, and 5) percent coding sequence was estimated as the number of nucleotides in coding sequences divided by the total genome size 100×. To determine whether there was a significant difference of properties between Xylophagaidae and Teredinidae, we conducted a paired t-test with value adjusted by Bonferroni methods using R package rstatix (Team 2021). Branch length differences among subtrees were examined using the Kruskal-Wallis H-test. Scenarios for potential gene rearrangements were explored based on analysis of common intervals using Crex (Bernt et al. 2007) as implemented on the MITOS2 server. The DNA Skew Graphing tool (Thomas et al. 2007), available online via the Viral Bioinformatics Research Centre (https://4virology.net/), was used to search representative mitochondrial genomes for abrupt changes in base composition bias that are characteristic of both the heavy-strand origin of replication (OriH) and light-strand origin of replication (OriL). After locating putative origins of replication, we utilized the UNAFold web server (Zuker 2003) to locate stable stemloop configurations containing characteristic T-rich loops (see supplementary figs. S2 and S3, Supplementary Material online; also (Brugler and France 2008) for a list of features typically associated with origins of replication and the application of the "DNA Walker" (Lobry 1996) graphing option to locate putative Oris within a mitogenome).

Supplementary Material
Supplementary data are available at Genome Biology and Evolution online. Professor Diego Veras Wilke for the encouragement, exchange of ideas, sampling logistics, and laboratory support.
Collection Permits: Collections in Brazil were performed under the following: SISBIO license number 48388 to A.T.-S., and genetic resources were accessed under the authorization of the Brazilian National System for the Management of Genetic Heritage and Associated Traditional Knowledge (SisGen permit number A2F0DA0) to the National Institute of Metrology Standardization and Industrial Quality (Inmetro, Instituto Nacional de