Variability of mitochondrial ORFans hints at possible differences in the system of doubly uniparental inheritance of mitochondria among families of freshwater mussels (Bivalvia: Unionida)

Background Supernumerary ORFan genes (i.e., open reading frames without obvious homology to other genes) are present in the mitochondrial genomes of gonochoric freshwater mussels (Bivalvia: Unionida) showing doubly uniparental inheritance (DUI) of mitochondria. DUI is a system in which distinct female-transmitted and male-transmitted mitotypes coexist in a single species. In families Unionidae and Margaritiferidae, the transition from dioecy to hermaphroditism and the loss of DUI appear to be linked, and this event seems to affect the integrity of the ORFan genes. These observations led to the hypothesis that the ORFans have a role in DUI and/or sex determination. Complete mitochondrial genome sequences are however scarce for most families of freshwater mussels, therefore hindering a clear localization of DUI in the various lineages and a comprehensive understanding of the influence of the ORFans on DUI and sexual systems. Therefore, we sequenced and characterized eleven new mitogenomes from poorly sampled freshwater mussel families to gather information on the evolution and variability of the ORFan genes and their protein products. Results We obtained ten complete plus one almost complete mitogenome sequence from ten representative species (gonochoric and hermaphroditic) of families Margaritiferidae, Hyriidae, Mulleriidae, and Iridinidae. ORFan genes are present only in DUI species from Margaritiferidae and Hyriidae, while non-DUI species from Hyriidae, Iridinidae, and Mulleriidae lack them completely, independently of their sexual system. Comparisons among the proteins translated from the newly characterized ORFans and already known ones provide evidence of conserved structures, as well as family-specific features. Conclusions The ORFan proteins show a comparable organization of secondary structures among different families of freshwater mussels, which supports a conserved physiological role, but also have distinctive family-specific features. Given this latter observation and the fact that the ORFans can be either highly mutated or completely absent in species that secondarily lost DUI depending on their respective family, we hypothesize that some aspects of the connection among ORFans, sexual systems, and DUI may differ in the various lineages of unionids.


Background
Many species of gonochoric bivalve molluscs from four different orders (Unionida, Mytilida, Venerida, Nuculanida) possess a peculiar mode of mitochondrial transmission called doubly uniparental inheritance, or DUI, which is particularly well-documented (~80 species) in freshwater mussels of the order Unionida [1,2]. DUI basically consists of a cyclic parent-specific inheritance of two distinct mitochondrial (mt) genomes (or mtDNA), where females segregate only the so called F ("female-transmitted") mtDNA in their eggs and males segregate only the other mitotype, named M ("male-transmitted"), in their sperm. The zygote is heteroplasmic but, depending on its sexual development in the subsequent stages, i.e. whether it will become a female or a male as an adult, an individual will transmit only one of these two types of mt genomes to the next generation [1]. DUI seems to be strictly associated with the gonochoristic sexual system of a particular species, as it was discovered that four species of unionids and one species of margaritiferid each appear to have lost independently their M mtDNA during the transition from dioecy to hermaphroditism [3]. Notice that, however, it is not known if the loss of the M type is perfectly contemporaneous with the switch to hermaphroditism. All these mentioned obligate hermaphroditic species now retain only a modified version of the F genome that is called the H genome (for "hermaphrodite") [3]. In an attempt to understand if and how DUI and sex determination are connected, genomic studies have highlighted the following features of DUI in freshwater mussel mtDNAs: (1) a high level of sequence divergence between F and M (up to~40% of difference in their nucleotide sequences, resulting into an almost 50% difference in the amino acid sequence of their encoded proteins [4,5]); (2) the presence of a 3′-elongation of the cox2 gene in the M mtDNA (compared to that found in most metazoans) [6]; (3) the presence in both mt genomes of "ORFan" genes, i.e., genes without obvious homology or function [7], named F-orf in the F mtDNA and M-orf in the M mtDNA, respectively [8]. An ORFan is found also in the H mtDNAs mentioned above: it is a highly mutated F-orf and is named Horf for distinction [3]. It is worth noting that these features of DUI-positive freshwater mussel mtDNAs also appear in practically all DUI bivalves outside Unionida, although with many variations and combinations (e.g., F vs M divergence can be lower; elongated cox2 genes can appear in the F instead of the M, or be absent from both; and mtDNAspecific ORFans can be present in only one of the two mtDNAs or duplicated in a same mt genome [1,5,[9][10][11][12]). Additional coding sequences are sometimes found in the mtDNA of bivalves (with or without DUI) and even other molluscan species [9,[13][14][15][16][17][18][19][20][21][22]. Usually they are identifiable as duplicated standard mitochondrial protein coding genes, with a variable degree of similarity to the original sequence. In fact, the two copies can evolve in tandem but, in other cases, the only recognizable parts consist of small segments encoding functional domains of the original mitochondrial genes. In absence of further functional studies, this leads us to hypothesize that the primary function of these genes, if indeed they retain functionality, might be related to oxidative phosphorylation.
Functional studies focusing on the ORFan genes F-orf and M-orf in DUI bivalves provided evidence they are transcribed and translated into proteins (from here on respectively indicated as F-ORF and M-ORF) located inside and outside the mitochondria [3,10,[23][24][25][26]. Moreover, extensive bioinformatic studies on both the gene sequences and their translated proteins produced evidence for two options for the origin of the ORFans, which may be the result of either (1) the insertion of viral sequences into the host mt genome [23,27] or (2) the duplication and subsequent modification of extant mitochondrial genes and sequences [5,28]. Analyses aimed at understanding their origin by predicting the functions of their protein products broadly converged on similar patterns in all species considered (freshwater mussels and others). F-ORFs and M-ORFs, for example, were predicted to interact with nucleic acids and/or membranes (for signaling and/or interactions with the immune system), and M-ORFs, in particular, were also predicted to interact with the cytoskeleton and to have a role in the ubiquitination processes [27,28]. The high variability of the ORFans among distantly related families and orders, however, questions their homology in all DUI bivalves. To have a better understanding of ORFans evolution and the functions of their encoded proteins, efforts should focus on a taxon in which DUI is widespread, such as the Unionida [2]. Until a few years ago, only F, M, and H mt genomes from the family Unionidae and a few F and H genomes of Margaritiferidae were available, but recently mtDNAs from families Hyriidae, Iridinidae, Mulleriidae, and the first margaritiferid M mtDNAs have been published [5,29]. Given the hypothesized, although still untested, link between ORFans and sexual systems in freshwater mussels [3], the sequencing of these new mt genomes allowed examination of the evolution of ORFan genes, DUI, and sexual systems in a phylogenetic context [5]. It was suggested that DUI may have been present as an ancestral state before the radiation of the order Unionida, and that some ORFans have been partially or totally purged from the remaining mtDNAs of some lineages that may have lost DUI in their early stages of radiation (Iridinidae and Mulleriidae). However, the ORFans have been maintained in the other families that regularly show DUI (Hyriidae, Margaritiferidae, Unionidae), and in each of these taxa, the mt genomes, especially the M, show their own family-specific peculiarities. For example, in margaritiferid M mtDNAs the M-orf is duplicated (one copy, M-orf1, appears to be homologous to the M-orf of Unionidae and Hyriidae, while the second copy, M-orf2, is specific to Margaritiferidae only), whereas the hyriid M-orf is much longer than those of Margaritiferidae and Unionidae [5].
In this study, we present eleven new mt genomes from ten species of freshwater mussels, with or without DUI and with different sexual systems (for each species, references describing DUI status and/or reproductive modes on which we relied for this study are given): Chambardia rubens (Lamarck, 1819) [30,31] [37] for Margaritiferidae. First, we characterize the overall structure of the new mt genomes and highlight their unique features, some of which are described for the first time, and we build a phylogeny of freshwater mussels using these and other genomes. Then we identify new Forf and M-orf genes and describe the new F-ORF and M-ORF proteins by comparing them to a set of already published sequences, showing how, despite having evolved different three-dimensional configurations, they share some key features. Finally, considering our findings, we discuss whether the DUI system works the same way in all DUI freshwater mussels or if there may be family-specific differences, as well as the modifications occurring in the mtDNAs after DUI is lost.

Sequencing, assembly, and general features of the new mt genomes
We obtained a complete sequence for ten of the eleven new mt genomes; the M mtDNA of W. carteri showed a sequencing gap in a non-coding region between trnV and trnH. All sequences were deposited in GenBank under the accession numbers MK761136-46 (Table 1). A summary of their length and a comparison of their gene order are given in Table 2. In summary, they present the same general features recently highlighted by [5] (for families Iridinidae, Mulleriidae, Hyriidae, Margaritiferidae) and [29] (for Margaritiferidae). As is typical for DUI freshwater mussels, the cox2 gene carried by W. carteri M mt genome is longer compared to its F counterpart (respectively 1329 bp and 693 bp) and to those of non-DUI species [6] (Table 3). Apart from the occasional species-specific difference in length of some noncoding regions, particularly in Hyriidae (i.e., between atp8 and nad4L in D. suavidicus, and between trnV and trnH in W. carteri M mtDNA) and Iridinidae (1049 bp between nad5 and trnF in C. rubens, compared to the 23-76 bp of the other species mtDNAs), the most notable features lie in the presence/absence of ORFans, on which we will focus.

Phylogeny of freshwater mussel mt genomes
A Bayesian inference analysis was performed with MrBayes [38] on 12 protein coding genes and their respective protein sequences extracted from the new 11 mt genomes, from 51 additional mtDNAs of freshwater mussels available in GenBank (25 F, 22 M, and other 15 mtDNAs from non-DUI species; Additional file 1: Table S1), and from three outgroup For each species for which at least one mt genome sequence was obtained in this study (either complete or largely complete), the respective family, and the provenance of the specimen(s) used for the sequencing (columns 'Country', 'Latitude', 'Longitude') are indicated, as well as the respective sexual system. For each mtDNA sequenced, its type is specified: non-DUI., mtDNA of dioecious or hermaphroditic species with no evidence of DUI presence or secondary loss; F and M, female-and male-transmitted mtDNAs of dioecious DUI species, respectively. The only available mtDNA of Pseudunio auricularius, a dioecious species for which no evidence of DUI has been produced yet (no M mtDNA sequence available yet), which carries an F-orf, has been considered as F for simplicity    In the nucleotide-based tree ( Fig. 1), freshwater mussel mtDNAs form a monophyletic group divided in two main branches: one containing all M mtDNAs and another containing all female-transmitted ones, from both DUI and non-DUI taxa. Relationships among DUI species in these two main branches are maintained in most cases, i.e. the phylogeny of F mt genomes mirrors that of M ones. In the few cases where this situation does not occur, either the nodes usually have posterior probabilities < 1.000 (e.g., see the different relative position of Aculamprotula tortuosa mt genomes in the Unionidae clades) or the number of F and M genomes for a taxon is different (e.g., Margaritiferidae). In both of the two main branches, family Hyriidae is sister group to Margaritiferidae and Unionidae, which always form reciprocally sister groups. For Hyriidae clades, W. carteri mt genomes are always sister to the respective Echyridella menziesii ones. In Margaritiferidae, P. auricularius female-transmitted mtDNA is sister to Pseudunio marocanus F mtDNA. Iridinidae and Mulleriidae form a single branch sister to Hyriidae + Margaritiferidae + Unionidae female-transmitted mtDNAs, where Mulleriidae form a monophyletic clade but Iridinidae do not: C. rubens is sister to all Mulleriidae (node posterior probability = 0.983), while the other iridinid Mutela dubia is sister to all. Inside Mulleriidae, the dioecious [32] A. elongata is recovered as distantly related to the congeneric hermaphrodite [39][40][41] Anodontites trapesialis.
The main difference of the protein-based tree ( Fig. 2) with the nucleotide-based one is the position of the M mtDNAs clade: here, it is sister to a clade containing N. margaritacea mtDNA as sister to all female-transmitted mt genomes of freshwater mussels. Minor differences with the nucleotide tree are as follows. In the femaletransmitted mtDNAs clade: for Hyriidae, W. carteri

Search and annotation of ORFan genes
To search for new ORFan genes, a set of 35,032 open reading frames (ORFs) ( Table 4) was extracted from the new mt genomes in Table 1 and the additional 51 mtDNAs of freshwater mussels from GenBank (Additional file 1: Table S1), for a total of 62 mt genomes analyzed. Using as a criterion of choice the level of similarity between the known ORFan proteins and the translated proteins of the nucleotide sequences in the ORFs set, we found with the HMMER [43] suite of programs 25 F-orfs, 5 H-orfs, and 26 M-orfs, three of which are M-orf2 from Margaritiferidae species [5] and one appears to be a recent duplication specific to the unionid S. woodiana (we named the two copies M-orfa and Morfb; see also [12]) (Additional file 1: Tables S2 and S3). Range, mean, and standard deviation (SD) of the cox2 gene length in egg-and sperm-transmitted mtDNAs in different freshwater mussel families. Spermtransmitted mtDNAs are present only in families showing DUI. Uniparental transmission of mtDNA through female gametes is assumed in Iridinidae and Mulleriidae, which do not show DUI. In Mulleriidae, the partially sequenced cox2 gene from Anodontites trapesialis was excluded from this analysis; however, it does not show the typical 3′ elongation of M cox2 genes [5] The nucleotide sequences of these 56 ORFans are listed in Additional file 2. In summary, single ORFan protein sequences used as seeds mainly recognized homologous ORFan sequences, with few non-ORFan hits (as in the case for E. menziesii M-ORF which recognized some full, in-frame nad4L sequence) (Additional file 1: Table S2). The F-ORF Hidden Markov Model (HMM) profile we used recognized, in addition to all sequences from which it was assembled, also the H-ORFs (Additional file 1: Table S3). The M-ORF HMM profiles recognized all the M-ORFs forming them, plus, with lower scores and Evalues, some protein translated from ORFs overlapping trnD, atp8 (in two cases the hit comprised the whole inframe sequence of its protein), nad6, and nad4L genes  Table S3). The putative proteins from the AtraUR219 and AtraUR2218 ORFans of A. trapesialis, located between atp8 and nad4L of this species mt genome (Fig. 3) (which have been proposed to have originated from duplication and divergence of atp8 and might be related to M-orfs of DUI species [5]), have no apparent homologs in other species. However, some ORFs from other species were retrieved that overlap the same region where the two A. trapesialis ORFans are located and from which are hypothesized to have originated: among these hits, for example, four include either a segment of the atp8 gene or its full in-frame sequence (Additional file 1: Table S4).
In short, as shown also in Fig. 3, we found that: (1) all DUI species of freshwater mussels analyzed carry a F-orf in their F mtDNA and at least one M-orf in their M; (2) secondarily hermaphrodite (i.e., that switched from gonochorism to hermaphroditism) species of Unionidae and Margaritiferidae that lost DUI always possess a Horf [3]; (3) species that do not show evidence of DUI (i.e., no evidence of heteroplasmy) from families Iridinidae, Mulleriidae, and Hyriidae have none of these ORFans. The only mtDNA we retrieved from P. auricularius presents a standard F-orf, and given that this species is dioecious, it is plausible it will be revealed as a DUI species, and therefore we treated this genome as F mtDNA.

Sequence-based analyses of ORFan protein products
The amino acid composition of F-ORFs appears to be rather homogeneous among families, with no clear differences and, although more variable, a general trend is observable also in the M-ORFs (Fig. 4).  (Fig. 5), but form a separate branch in the ML tree (which has, however, low resolution) (Fig. 6). We also attempted to add AtraUR219 and AtraUR2218 sequences to the M-ORF alignment for the ML reconstruction, but this disrupted the clustering of M-ORFs, especially for the more numerous Unionidae (not shown).

Tertiary structure prediction of ORFan proteins
Currently, there are no data derived from crystallographic studies on the ORFan proteins, nor established For each taxonomic group of freshwater mussels, the number of species and mtDNAs considered in this study are given, as well as the number and kind of ORFans (F-, M-, or H-orfs) retrieved from them. mtDNA types: F and M, female-and male-transmitted mtDNAs of dioecious DUI species, respectively; H, mtDNA of stably hermaphroditic species that lost DUI sensu [3]; non-DUI, mtDNA of dioecious and hermaphroditic non-DUI species. mtDNAs of dioecious species carrying an F-orf but for which no evidence of DUI has been produced yet, or of dioecious DUI species that can show occasional hermaphroditism, were considered as F for simplicity structural similarities with known proteins in databases, that may guide bioinformatic analyses aimed at predicting the ORFan proteins folding. Therefore, for completeness, we decided to show the best results we obtained with I-Tasser [45] regardless of the C-scores assigned to the models: C-scores are usually comprised between − 5 and 2, therefore higher values indicate higher confidence of the model. The predicted tertiary structure of some selected ORFan proteins may appear to be highly different at first sight, but similarities among proteins of the same kind can be recognized (Figs. 7 and 8). A common feature between the F-ORFs of Hyriidae and Unionidae is the presence of two antiparallel helices separated by a loop. This conformation is not found in the F-ORF of Margaritiferidae, in which only one small helix is predicted (preceded by a small beta strand in Cumberlandia monodonta and Pseudunio marocanus but not in Margaritifera margaritifera). Hyriidae M-ORF, Margaritiferidae M-ORF1 and M-ORF2, and S. woodiana M-ORFa all share the presence of three antiparallel helices in their N-terminus. The three helices have the same relative orientation, but the third one is in front of the first two in Hyriidae and behind them in Margaritiferidae and S. woodiana (and it is also much smaller in this species).
The portion beyond this third helix varies for each species, but, for example, Hyriidae M-ORFs are similar in this part of the protein and are clearly discernible from those of Margaritiferidae, and M-ORF1s and M-ORF2s of this family are again distinguishable between them. Cumberlandia monodonta M-ORF1 structure is less defined compared to the homologous proteins from M. margaritifera and P. marocanus, and in its M-ORF2, the third helix appears to be on the same plane as the other two. The three Unionidae M-ORFs examined have extremely divergent configurations, and no obvious similarities can be recognized among them. The two S. woodiana M-ORFs, most probably the product of a duplication event specific to this species [12], do not resemble one another. AtraUR219 protein is constituted by a short N-terminal beta strand, two helices crossing each other and connected by a simple loop, and a small C-terminal beta strand. AtraUR2218 protein is very short (22 aa) and it is predicted to be only a single helix.

Three-dimensional alignments of ORFan proteins
Summarized statistics for pairwise three-dimensional (3D) alignments of the ORFan protein models are shown in   Table 1 and Additional file 1: Table S1. Standard mitochondrial genes are in grey, while ORFan genes (see the main text for a complete description of these genes) are colored following this code: green, Anodontites trapesialis specific ORFans; blue, M-orfs; pink, F-orfs; light pink, Horfs. Genes are pointed according to their relative direction on the mtDNAs. tRNA genes are indicated with the one-letter code of their respective amino acid. Dotted lines represent the segments between the two regions, which are not indicated for simplicity. Sinanodonta woodiana annotation is based on [12] and on the current study each single freshwater mussel family, the distances (expressed as root mean square deviation of distances) between Cα and between Cß atoms (respectively the α- of sequences the alignment could successfully produce a tree but, because of the low number of sequences available, the trees did not show any informative clear-cut clustering patterns that could split the proteins into, for example, family-specific (e.g., Hyriidae vs Margaritiferidae vs Unionidae) or kind-specific (e.g., Margaritiferidae M-ORF1 vs M-ORF2) patterns as in the sequence-based analyses.

Discussion
The deep relationships among freshwater mussel families have long been debated [46] but, although mt genomes from the sixth family Etheriidae are at present still not available, the phylogenies presented here support a sister group relationship between Hyriidae and Margaritiferidae + Unionidae, and an equally strict relationship between Mulleriidae and Iridinidae, although not well resolved: comparable family-level topologies were also found by other recent studies [5,42,46] using different analytical methods and/or taxa. Such topology supports an early classification by [47] that splits freshwater mussels in two superfamilies, Unionoidea (Unionidae + Margaritiferidae + Hyriidae) and Etherioidea (Mulleriidae + Iridinidae + Etheriidae). In the light of our and the other mentioned results [5,42,46], and as properly discussed by [46], the separation in these two major taxa better reflects the monophyly of shared characters among families than others that introduce a third superfamily, Hyrioidea, for Hyriidae only (as in [5]). It is worth noting how the female-transmitted mt genomes of Unionoidea, the superfamily where DUI is common, form a single clade sister to one containing mtDNAs of the non-DUI superfamily Etherioidea. This, however, does not imply that DUI was present only in the common ancestor of the Unionoidea. Indeed, the variable position of the M mtDNAs clade suggests the presence of DUI either in the last common ancestor of all freshwater mussels (Fig. 1) or even earlier, before the split between orders Unionida and Trigoniida, represented by the species N. margaritacea (Fig. 2). This is because when speciation occurs after DUI appears, F and M genomes evolve according to a "sex-associated" phylogenetic pattern [48] in two distinct clades and, inside these two clades, the relationships among mtDNAs of the various species are the same. Our phylogenies, therefore, suggest that (1) DUI was lost by Iridinidae and Mulleriidae, as well as by the South American lineage of Hyriidae (as discussed in detail below), and (2) that at least the last common ancestor of all Unionida had DUI. Comparable phylogenies were already retrieved with different methods and taxa [5,42,46] and, although the mt genomes sequenced in this study largely follow already described architectures [5,29] (Fig. 3), we observed new interesting features that help us reconstruct the evolution of freshwater mussel mt genomes and their relationship with DUI. Starting with family Mulleriidae, the four new mtDNAs have no trace of rearrangements that are reminiscent of the ORFans AtraUR219 and AtraUR2218 found in A. trapesialis ( Table 2, Fig. 3), which were hypothesized to be early versions of M-orfs [5]. The structure of Mulleriidae mt genomes, therefore, makes them much more similar to Iridinidae mtDNAs, which also lack additional ORFans between atp8 and nad4L (Table 2, Fig. 3). This can be interpreted as evidence for negative selection against the rise of novel coding sequences in both families Iridinidae and Mulleriidae. It is however notable how the two mentioned ORFans are present only in A. trapesialis, sister species to all other Mulleriidae in our phylogenies (Figs. 1, 2 and 3). Until further studies, this might indicate an independent and relatively recent genomic rearrangement in A. trapesialis giving rise to its two ORFans: therefore, future works investigating ORFans evolution should consider the possibility that both AtraUR219 and AtraUR2218 may be relicts of a speciesspecific duplication event, unrelated to M-orfs of other DUI families and possibly non-functional. Also, the selection against new sequences may not be correlated to the sexual Fig. 6 Unrooted maximum likelihood (ML) trees for F-ORF and M-ORF proteins of freshwater mussels. Color code for each family are indicated inside the panels. Bootstrap values are indicated at each node system of these species, as both dioecious (A. elongata, F. fossiculifera, C. rubens) and hermaphroditic (M. parchappi) ones (Table 1) seem subject to it. For Iridinidae, C. rubens mtDNA also confirms that the non-coding region between nad5 and trnF is unusually large compared to other freshwater mussel families [5]: whether it is a control region containing regulative motifs not found in other freshwater mussels will be a matter for future studies. Finally, the relationships among the examined Iridinidae and Mulleriidae species obtained from our phylogenies do not help in solving the history of their mt genome architectures. The two families, although strictly related, do not form separate monophyletic clades, and the two congeneric species A. trapesialis (hermaphrodite [39][40][41]) and A. elongata (dioecious [32]) are distantly related (Figs. 1 and 2): whether this situation calls for taxonomic revisions or not should be a matter for further ad hoc studies.
The new Margaritiferidae female-transmitted mt genome of P. auricularius has an F-orf, as previously described for this taxon [3,5,29] (Table 2, Fig. 3), and is strictly related to F mt genomes: further studies are surely needed to fully confirm the presence of DUI in this species, but the available evidence points to this direction. The re-analysis of published Unionidae mt genomes allowed us to confirm a recent species-specific duplication of the M-orf in S. woodiana M mtDNA, as noted by [12] (Fig. 3). We are unable to say what effects (if any) this mutation may have caused to S. woodiana DUI system, but our results suggest that the M-orfa (immediately upstream of nad4L) is the one more similar to those of other Unionidae, while M-orfb (immediately upstream of trnD) appears different, although still recognizable as an M-orf. This feature of S. woodiana and the previously described rearrangements in A. trapesialis support the idea that the atp8-nad4L region in the mtDNAs of freshwater mussels is a hotspot for significant rearrangements and gene duplications, therefore offering additional support to the hypothesis for which the

Sinanodonta woodiana (b)
169 aa, C-score: -3.85 Fig. 8 3D models of the proteins encoded by Anodontites trapesialis ORFans and of representative M-ORF proteins of DUI freshwater mussels. The models shown are the first of the top five predicted by I-TASSER for each sequence. Number of amino acids (aa) of each protein and C-score of the models are indicated under the relative species names. C-scores are usually comprised between − 5 and 2: higher values indicate higher confidence of the model. The color shading of each protein goes from the blue of the N-terminus to the red of the C-terminus unionid M-orf may have originated from a duplication of atp8 [5,28]. The mtDNAs of the only dioecious species of Hyriidae in this study showing DUI, W. carteri, are comparable in all aspects to those of the previously sequenced DUI species E. menziesii [5], a strictly related species with the same Australasian distribution. In particular, both the M-orf and the cox2 in the M mtDNA are confirmed to be longer in this family compared to the others ( Table 3). The other three Hyriidae species, C. ambigua, D. suavidicus, P. obliquus (all from South America and always forming a single monophyletic clade; Table 1, Figs. 1 and 2), did not show evidence of DUI and, contrary to DUI gonochoric and DUI-less hermaphroditic unionids and margaritiferids, their F-like mtDNAs do not possess any F-orf or H-orf, and we did not find evidence of their translocation in the unassigned regions of these mt genomes (Table 2, Fig. 3). Even if the sex determination system is known only for C. ambigua, a gonochoristic DUI-less species, we can propose a working hypothesis stating that, in Hyriidae, losing DUI (1) may always cause the complete disappearance of the F-orf in the former F mtDNA (as in C. ambigua, D. suavidicus, P. obliquus), and (2) may not affect the gonochoristic sexual system (as in C. ambigua). In this family, similarly to Mulleriidae and Iridinidae which are hypothesized to have lost DUI in the early stages of their radiation [5], it seems therefore that the relationship among DUI, presence of ORFans, and gonochorism may be somewhat different compared to Unionidae and Margaritiferidae, which retain a H-orf in their mt genomes after losing DUI [3]. The parallelism between Hyriidae and Mulleriidae + Iridinidae may, however, lead to another hypothesis: we can see that among all Hyriidae species studied until now, only the Australasian ones (E. menziesii, W. carteri) show DUI, while the Neotropical ones (C. ambigua, D. suavidicus, P. obliquus) do not (Table 1). This may hint that the last common ancestor for these two lineages had DUI, which was lost (together with the ORFan in the remaining F-like mtDNA) only in the South American lineage during its radiation. Considering the current information from all families of freshwater mussels, we can speculate that once DUI and the M mtDNA are lost by a species, the ORFan  in the remaining F mtDNA (i.e. the F-orf) no longer plays a role in the DUI mechanism and gradually disappears. First, it may start to accumulate mutations and degenerate (like the H-orf in Unionidae and Margaritiferidae [3]) and then, given enough time, it completely disappears from the mtDNA without leaving recognizable traces (as in the South American Hyriidae, which no longer carry traces of the F-orf; Table 2, Fig. 3).
Despite very similar amino acid compositions (Fig. 4), sequence-based analyses of ORFan proteins managed in most cases to distinguish families (Figs. 5 and 6), with F-ORFs giving better resolution compared to M-ORFs (as in the ML analysis, for example). Following this lead, we explored for the first time the total putative 3D folding of freshwater mussels ORFan proteins (as secondary structures and other features have been already thoroughly characterized [27,28,49]), to search for patterns that could help unravel their evolutionary history. Indeed, even if we sampled only a few representative species, the predicted 3D foldings demonstrate how in each family the F-and M-ORFs have their own peculiar shape (Figs. 7 and 8). The tertiary structure of a protein is influenced by its amino acid sequence, and as a consequence, the proteins of closely related species and/or with a common origin exhibit comparable shapes. In Margaritiferidae, for example, all M-ORF1 and M-ORF2 proteins fold in a similar way, and this finding could also hint at a full functionality of the M-ORF2 protein in DUI margaritiferids, possibly with a physiological role comparable to that of M-ORF1 and other M-ORFs. In contrast, in the case of Unionidae M-ORFs (Fig. 8), the proteins exhibit a range of very divergent foldings; this could be due to the method used, the phylogenetic distance among species (see their positions in the trees in Figs. 1 and 2), or the exclusive evolutionary history of each protein (as we discussed above for S. woodiana M-ORFs). A broader sampling from each family and a comparison among different methods of tertiary structure prediction in the future could help to better define these ambiguous foldings, as well as to improve the distance analyses that lacked resolution in the current study.
Nonetheless, despite the technical difficulties, our explorative study presents evidence of possibly conserved features among ORFan proteins of the same kind, such as the relative arrangement of certain helices in F-ORFs and M-ORFs. These structural features, together with properties already characterized (e.g., [5,27,28,49]) and others yet to be discovered, will lead us to give a precise physiological role to the ORFan proteins and their respective genes (like those already hypothesized before, for example [28]). With further study, we might also be able to answer the longstanding questions about the relationships among the ORFan genes, the sex determination system, and the peculiar mitochondrial inheritance mode of freshwater mussels and of all other bivalves showing DUI [3]. However, given the rather different length and shape of Hyriidae M-ORFs compared to those of Unionidae and Margaritiferidae, and the fact that Hyriidae DUIless species are not always hermaphroditic (as C. ambigua) and their mtDNA does not possess any For H-orf (C. ambigua, D. suavidicus, P. obliquus), as opposed to what occurs in Margaritiferidae and Unionidae [3], we suspect that the nature of the link among ORFans and sexual system may not be exactly the same in all DUI freshwater mussels. On the other hand, it is also possible that long evolutionary times in the absence of DUI, as well as various ecological pressures [50], may shape freshwater mussel mt genomes and/or sexual systems, leading to the situations described for the first time in this study (i.e., no ORFan genes in both hermaphroditic and gonochoric species). To answer these questions, a broader sampling and in vivo studies on the ORFan proteins will be needed.

Conclusion
In this study we produced ten entire, plus one almost complete, new mtDNA sequences of freshwater mussel species with or without DUI from still poorly sampled families (Iridinidae, Mulleriidae, Hyriidae, and Margaritiferidae). Besides being a useful basis for future sequencing efforts, much needed for this group of endangered animals [51], we provided new data on the mitochondrial ORFan genes of freshwater mussels, whose origin is still unknown and whose function and conservation are likely related to their sex determination system [3]. We observed that the rearrangements occurring in mitochondrial genomes of species and lineages that secondarily lost DUI (i.e., that lost their ancestral M mtDNA and retain only a mutated F) are not always the same, and that losing DUI is not always linked to a switch to hermaphroditism. By analyzing the 3D structures of their translated proteins, we also evidenced common characteristics and similarities among them, hinting at conserved physiological roles of F-orf and M-orf genes in all DUI lineages of freshwater mussels, as well as familyspecific ones. We therefore questioned if the familyspecific structures of the ORFan proteins can influence some detail of the DUI system in different manners, so that the downstream effect of losing DUI on the sexual system of a species may vary. An alternative, but not necessarily mutually exclusive, hypothesis we propose to explain the observed differences among non-DUI lineages is that time and other factors may play an important role in reshaping both the mitochondrial genome and sexual system of a species after it loses DUI.

Sequencing of new mt genomes
Freshwater mussel species were selected across the main families within the order Unionida and to cover distinct sexual strategies, i.e. gonochorism and hermaphroditism (Table 1). W. carteri specimens were taken from the wild with permits for field and laboratory studies obtained from the Western Australian Department of Environment and Conservation under Regulation 17 of the Wildlife Conservation Act 1950 (SF007049) and Department of Fisheries under Exemption from the Fisheries Resources Management Act 1994 (1724-2010-06). Sex of the specimens was determined by observing gonad tissue smears for sexual cells and/or the demibranchs for the presence of marsupia, using a dissecting microscope. Tissue samples were excised from one specimen per species, and placed in 100% ethanol for DNA extraction: for all species, a foot clip was available for DNA extraction, while for W. carteri an additional male gonad sample from the same specimen was also used. DNA extraction followed [52]. DUI presence or absence for every species was assumed from previous studies [30][31][32][33][34][35][36][37]. The complete mitogenomes sequencing and assemblage was accomplished using the pipeline proposed by [53]. Annotations were performed using MITOS [54] with the final tRNA genes limits being rechecked with ARWEN [55]. Finally, personal scripts were developed and applied to adjust the mtDNA protein-coding limits since MITOS seems to underestimate gene length (for details, go to https://figshare.com/s/a756ef19cec8f65d506a).

Phylogenetic analyses of freshwater mussel mt genomes
The set of eleven new mt genomes (Table 1) was expanded by adding other 51 freshwater mussel mt genomes from GenBank (see Additional file 1: Table S1 for the complete list and details), for a total of 62 mtDNAs for 40 species. The mt genomes added encompass families Iridinidae, Mulleriidae, Hyriidae, Margaritiferidae, and Unionidae. For Hyriidae and Unionidae, we took the mt genomes from DUI species for which both the F and M mtDNAs were available (resulting in two mtDNAs per species), for Margaritiferidae all F and M mt genomes available, and for Margaritiferidae and Unionidae only also those from secondarily hermaphrodite species (i.e., the H mtDNAs). The final set was thus composed of 22 M, 25 F, and 15 mt genomes from non-DUI species, either secondarily hermaphrodite that lost DUI (sensu [3]) or gonochoric ones. For the purpose of the phylogenetic analyses, three outgroup species mt genomes were also used: the bivalves Neotrigonia margaritacea (Trigoniida) and Solemya velum (Solemyida), plus Chaetoderma nitidulum (Caudofoveata, Chaetodermatida) as outgroup to all bivalves (respective GenBank accession numbers: KU873118, JQ728447, EF211990). A total of 65 mt genomes was therefore considered for the phylogenies.
We extracted all protein coding gene sequences, except atp8 because of its short length and high variability, from the 65 mtDNAs and translated them with the invertebrate mitochondrial genetic code to obtain the relative protein sequences. The 12 protein sets were then aligned with M-Coffee (http://tcoffee.crg.cat/apps/ tcoffee/do:mcoffee) [56,57] using all multiple methods available, then from these protein alignments we retroaligned the codons of the respective genes using the TranslatorX server (http://translatorx.co.uk) [58]. Both protein and codon alignments were trimmed on the Gblocks server version 0.91b (http://molevol.cmima. csic.es/castresana/Gblocks_server.html) [59,60] using the option for a more stringent selection. jModelTest2 [61,62] was used to calculate, under the Bayesian Inference Criterion (BIC), the best-fit models of nucleotide substitution for the trimmed codon alignments. Finally, the two sets of trimmed alignments were concatenated respectively into a codon alignment and an amino acid alignment, with a respective length of 7914 and 2363 gapless positions. MrBayes 3.2.3 [38] was used to infer amino acid-and nucleotide-base phylogenies of the mt genomes. The analysis for both concatenated alignments consisted of two separate runs of four chains, 5, 000,000 generations, sampling every 100 trees with a burn-in of 0.1%. In the nucleotide analysis, the models retrieved with jModelTest2 were specified for each gene partition, and a '4by4' nucleotide substitution model was adopted for the whole alignment. In the amino acid analysis, a 'mixed' rate matrix was specified. Completed runs were accepted for further examination after checking that their standard deviation of split frequencies stabilized at values < 0.01 over the generations (as in [63]). jModelTest2 and MrBayes 3.2.3 were ran on the CIPRES Science Gateway (http://www.phylo.org) [64]. Trees were graphically edited with FigTree v. 1.4.3 [65].

Annotation of F-and M-orf genes
To locate the F-and M-orf genes in the DUI genomes in which they were not annotated, and at the same time validate previous annotations, we first used the EMBOSS tool getorf [66] to extract all possible ORFs ≥30 nucleotides long (i.e., coding at least 10 codons) under the invertebrate mitochondrial genetic code from the 62 freshwater mussel mtDNAs dataset described above, and then translated them into the corresponding proteins. This set of protein sequences was first searched with the HMMER tool jackhmmer [43] (10 iterations) using as seeds the F-ORF of V. ellipsiformis and the M-ORFs of E. menziesii, C. monodonta (M-ORF1) and V. ellipsiformis, in separate runs. The proteins retrieved from each run were then aligned with PSI-Coffee [56,67] (http:// tcoffee.crg.cat/apps/tcoffee/do:psicoffee; all pairwise methods selected), and the alignments used to build HMM profiles with hmmbuild [43] (options: -fast -symfrac 0 -fragthresh 0 -wnone -enone; see also [28]). These profiles were used to search again the whole set of proteins with hmmsearch [43] (−-max option active to allow maximum sensitivity) to confirm the presence of F-orfs and M-orfs previously found with jackhmmer, retrieve the known ones not recognized by jackhmmer, and search for new homologs of these genes. Finally, phmmer [43] (−-max option active) was used to search the protein set for homologs of the two proteins putatively encoded by the two ORFans AtraUR219 and AtraUR2218 in A. trapesialis (Mulleriidae) mtDNA, hypothesized to be related to M-orf genes [5].
Sequence-based analyses of proteins MEGA7 [68] was used to calculate the amino acid composition of F-ORFs, M-ORFs, and AtraUR219 and AtraUR2218 of A. trapesialis. To visualize the relationships among the already annotated and newly discovered ORFan proteins based on pairwise similarity, we ran CLANS [44]. Specifically, the CLANS program was conducted online on the MPI Bioinformatic Toolkit website [69] (https://toolkit.tuebingen.mpg.de/#/tools/clans) using the BLOSUM45 scoring matrix and BLAST HSP's Evalues up to 1E-4. The output was then run locally on the CLANS application for ≥10,000,000 rounds to obtain reliable 3D clustering data. The alignments of the F-ORFs and M-ORFs used to build their HMM models, plus an alignment of the two A. trapesialis ORFan proteins and the M-ORFs (constructed with PSI-Coffee as mentioned above), were used to build ML trees in MEGA7 [68], using 1000 bootstraps, the 'mtREV' model of substitution, uniform rates among sites, and a gap partial deletion of 95%. The trees were built unrooted because given the uncertain origin of the ORFan genes it is not possible to choose a reliable outgroup sequence.

3D structure-based analyses of proteins
We used I-TASSER [45] online (https://zhanglab.ccmb. med.umich.edu/I-TASSER/) to obtain the 3D models of the F-and M-ORF of some representative DUI species (Hyriidae: E. menziesii, W. carteri; Margaritiferidae: C. monodonta, M. margaritifera, P. marocanus; Unionidae: V. ellipsiformis, S. woodiana) plus AtraUR219 and AtraUR2218 of A. trapesialis (Mulleriidae). The most supported models (i.e., the ones with the best C-score) were then used as input for MATRAS [70] (http:// strcomp.protein.osaka-u.ac.jp/matras/) to perform pairwise and multiple 3D alignments of the proteins. The multiple alignments aimed at obtaining trees based on DRMS (root mean square deviation of Cα atoms, measured in Å) distances among them, using as a minimal set the ORFan proteins from E. menziesii, C. monodonta, and V. ellipsiformis (which have been thoroughly characterized in past studies [3,5,27,28]) and adding as much proteins as MATRAS would allow from the other species. When an I-TASSER model made MATRAS fail in producing a tree, we refined it with ModRefiner [71] (https://zhanglab.ccmb.med.umich.edu/ModRefiner/) and repeated the 3D alignment. If the refining did not succeed in improving the results, the protein was removed from the analysis.