1 INTRODUCTION

Solemyidae is an intriguing group of protobranch bivalves whose fossil records trace back to the early Ordovician (Pojeta, 1988) and includes, as modern solemyids, two reciprocally monophyletic genera, Acharax Dall, 1908 and Solemya s. l. Gray 1840. These are distinguished by the position of the ligament: external in Acharax while internal in Solemya s. l. (Taylor et al., 2008; Kamenev, 2009; Oliver et al., 2011; Sharma et al., 2013). Solemyid bivalves are distributed at various depths around the world, generally to inhabit in anaerobic muddy and sandy bottoms and organic falls (Conway et al., 1992; Fujiwara, 2003), and nutritionally depend on intracellular chemosynthetic symbionts (Fisher and Childress, 1986; Rodrigues et al., 2010; Walton, 2015; Fukasawa et al., 2017). Acharax generally has a deeper distribution (400–5 379 m) associated with reducing environments (Neulinger et al., 2006), while Solemya is found mainly in the continental shelf and upper slope (0–600 m) (Neulinger et al., 2006; Oliver et al., 2011), however, with further research, Solemya has also been recorded from bathyal depths (Oliver et al., 2011; Sato et al., 2013).

Because the shell characteristics of solemyids are conservative and difficult to assess (Kamenev, 2009; Oliver et al., 2011), and they lack features such as shell sculpture and hinge teeth and are uniformly covered by distinctive, thick, shiny periostracum (Taylor et al., 2008; Saether et al., 2016), little morphological and molecular work has focused on protobranch bivalves. Therefore, the taxonomy and systematic status of Solemyidae are still problematic (Walton, 2015; Fukasawa et al., 2017). Moreover, solemyids live in deep burrows (Seike et al., 2012) and are able to swim when disturbed, which hinders the collecting work.

Mitochondrial genomes (mitogenomes) have been proven to be a useful source of molecular markers for taxonomic identification, population genetics studies, phylogeography and phylogenetic analysis of bivalves at different taxonomic levels, including members of Protobranchia (Sharma et al., 2013; Sato et al., 2020). The typical metazoan mitogenome is a closed circular molecule of 12–22 kb, usually contains 37 genes: 13 protein-coding genes (PCGs), 2 rRNAs, and 22 tRNAs (Boore, 1999). With the maturity of next-generation sequencing technology, comparative analysis and phylogenetic studies of bivalves based on complete mitogenomes have recently become popular (Uribe et al., 2018; Liu et al., 2020; Rahuman et al., 2020). At present, there are more than 950 species of protobranchia recorded around the world, but only 4 mitogenomes have been sequenced (Plazzi et al., 2013; Russell et al., 2018), all of which belong to Solemyidae (Solemya elarraichensis, S. pervernicosa, S. velesiana, and S. velum), and no mitogenome of the genus Acharax has been sequenced. This means evaluation of the systematic status of Solemyidae using mitogenomes is unachievable.

Our goal is to show that mitogenomes can be used to reconstruct the phylogeny of bivalves and further demonstrate its potential for resolving the systematic position of Solemyida. We attempted to accomplish this by: (ⅰ) characterize the mitogenome organization of Acharax sp. and compared the composition, relative synonymous codon usage (RSCU), and genetic distances with other available Solemyidae mitogenomes to evaluate the variation and conservation of mitogenomes among solemyid bivalves; (ⅱ) assess the intraspecific divergence of protobranch bivalves using the COX1 gene; (ⅲ) investigate the systematic status of solemyid bivalves in Bivalvia based on phylogenetic analyses. Since the only mitogenomes available within the Soleymida are all in the genus Solemya, we will further explore phylogenetic relationships within the Solemyida by targeting 3 nuclear and 2 mitochondrial genes. Lastly, we postulate that bivalves that live in the deep-sea are likely under different selection pressures due to life in an extreme environment than bivalves of shallower depths. Anatomical and physiological studies have shown that in order to adapt to the deep-sea environment, solemyid bivalves are significantly different from shallow protobranch bivalves in a series of physiological and biochemical aspects, such as the labial palp, digestive tract, and gills (Xu, 1999). Therefore, we speculate that the extreme deep-sea environment may affect the mitogenomes of solemyid bivalves. We will demonstrate this by evaluating the amino acid preferences of mitogenomes among bivalves from different depths, we were expect to see differences in the relative percentage of non-polar and polar amino acids between shallow water and deep-sea species.

2 MATERIAL AND METHOD 2.1 Sampling, DNA extraction, PCR amplification, and sequencing

The Acharax sp. was captured from the Haima methane seep in the northwestern slope of the South China Sea at 1 390 m using a remotely operated vehicle (ROV) in 2018. After being photographed onboard, the specimen was frozen in liquid nitrogen until DNA extraction. Morphological observations indicate that the specimen belonged to the solemyid genus Acharax sp. (Taylor et al., 2008; Oliver et al., 2011).

Total genomic DNA was extracted from the muscle tissues using TIANamp Marine Animals DNA Kit (TIANGEN, China) following the manufacturer's instructions. Three nuclear genes (18S rRNA, 28S rRNA, and H3) were amplified with primer pairs (Supplementary Table S1). PCR amplifications were carried out in 25-μL volume containing 12.5-μL 2X Premix Ex Taq (TaKaRa, China), 1-μL template DNA, 1 μL of each primer and 9.5-μL ddH₂O. The cycling conditions were as follows: initial denaturation for 2 min at 94 ℃, followed by 35 cycles of denaturation at 94 ℃ for 1 min, annealing at 50 ℃ for 30 s, extension at 72 ℃ for 50 s, and a final extension at 72 ℃ for 10 min. PCR products were purified using the PCR clean-up Kit (Axygen). The purified PCR products were bidirectionally sequenced on ABI 3730XL DNA Analyzer.

2.2 Mitogenome sequencing, assembly, and annotation

Total genomic DNA was used for Illumina sequencing and Oxford Nanopore sequencing (Shanghai BIOZERON Co. Ltd.). The Illumina data was used to evaluate the complexity of the mitogenome and correct the Nanopore long reads. For Illumina sequencing, the high quality genomic DNA was fragmented to ~400 bp using the Covaris M220 system, used to construct short-insert libraries according to the instructions (TruSeqTM Nano DNA Sample Prep Kit, Illumina), and then sequenced on an Illumina NovaSeq 6000 platform with 150-bp paired-end reads length. Sequencing produced 4 217.1 MB of clean data. Clean data were obtained by using Trimmomatic 0.39 (http://www.usadellab.org/cms/index.php?page=trimmomatic). For Oxford Nanopore sequencing, gTube was used to break the genomic DNA, and long-read sequencing followed the protocol in the SQK-LSK109 genomic sequencing kit (ONT, Oxford, UK). Nanopore sequencing was performed with the GidION using the R9.4.1 flow cells.

For the mitogenome assembly, the ABySS (http://www.bcgsc.ca/platform/bioinfo/software/abyss) was used to do genome assembly with multiple-Kmer parameters and got the optimal results of the assembly. Then, canu (https://github.com/marbl/canu) was used to assemble the Nanopore corrected long reads. Finally, according to the paired-end and overlap relationships of reads, GapCloser (https://sourceforge.net/projects/soapdenovo2/files/GapCloser/) was subsequently applied to fill up the remaining local inner gaps and correct the single base polymorphism for the final assembly results. The circular mitogenome map of the Acharax sp. was drawn with Organellar Genome DRAW (Lohse et al., 2007).

The assembled mitogenome was annotated in the MITOS webserver (Bernt et al., 2013). The PCGs were analyzed with the ORF finder available on the NCBI (https://www.ncbi.nlm.nih.gov/orffinder/). The positions of the rRNA genes were delineated by the rRNAmmer webserver (https://services.healthtech.dtu.dk/service.php?RNAmmer-1.2). The tRNA genes were confirmed using the tRNAscan-SE webserver (Chan and Lowe, 2019).

2.3 Data analysis

The nucleotide composition, codon usage, and relative synonymous codon usage (RSCU) values of five Solemyidae mitogenomes (one newly sequenced and four previously published: ON023263, NC_ 034904, NC_034905, NC_034906, NC_017612) were analyzed by MEGA 7.0 (Kumar et al., 2016). Composition skew values were calculated using the following formulas: AT skew=(A–T)/(A+T) and GC skew=(G–C)/(G+C) (Perna and Kocher, 1995). The nucleotide diversity (Pi) and genetic distances (Kimura's 2-parameter) of 13 PCGs in Solemyidae mitogenomes were estimated by DnaSP 6.0 (Rozas et al., 2017) and MEGA 7.0, respectively. The rates of nonsynonymous to synonymous substitutions (Ka/Ks) of 13 PCGs of the five solemyids were also determined in DnaSP 6.0. The heterogeneities of the COX1 gene sequence divergence among the protobranch bivalves were evaluated by AliGROOVE (Kück et al., 2014).

The compositions of amino acids with different properties (nonpolar, polar uncharged, and polar charged) were compared between 18 deep-sea and 45 shallow water bivalves (Supplementary Table S2) with the PCGs of the complete mitogenomes. The pairwise differences and significance levels between the compared groups were evaluated by two-tailed t-tests in SPSS 16.0. Before analysis, the normality of the distributions of the data was assessed using the Kolmogorov-Smirnov test in SPSS 16.0. The results show that the data conform to normal distribution.

2.4 Gene order analysis

The putative molluscan ancestral pattern and the mitogenome order of Acharax sp. were pairwise compared to predict the mitogenome rearrangement events (e.g., transpositions, reverse transpositions, reversals, and tandem duplication random loss) using Common interval Rearrangement Explorer (CREx), heuristically determines mitochondrial rearrangement scenarios based on common intervals (Bernt et al., 2007).

2.5 Phylogenetic analysis

The phylogeny was performed based on the mitogenomes of Acharax sp. and 47 other bivalve species belonging to two subclasses, 12 orders (Supplementary Table S3). Two species of Polyplacophora served as outgroups. The nucleotide sequences of 13 mitochondrial PCGs were aligned using MAFFT version 7 webserver (Katoh et al., 2019) with default parameters. Ambiguously aligned regions were removed by Gblocks 0.91b (default setting) (Talavera and Castresana, 2007). Then, PhyloSuite 1.2.2 (Zhang et al., 2020b) was used to concatenate the trimmed alignments into a single dataset. The best partition schemes and optimal substitution models were inferred using PartitionFinder (Lanfear et al., 2016) (Supplementary Table S4). Phylogenetic relationships were inferred using both maximum likelihood (ML) and Bayesian inference (BI) methods. The ML tree was constructed by IQ-TREE (Nguyen et al., 2015), and branch support was assessed by ultrafast bootstrapping with 1 000 replicates (Hoang et al., 2018) as well as Bayesian-like transformation of the aLRT test (Anisimova et al., 2011). The Bayesian analysis was performed by MrBayes (Ronquist et al., 2012). Two independent runs were carried out with four Markov Monte Carlo chains for 1 000 000 generations started from a random tree, with sampling every 100 generations. The average standard deviation of split frequencies and the likelihood values were monitored, to determine whether the two runs converged onto the stationary distribution. The initial 25% of the trees were discarded as burn-in. The remaining trees were used to construct the 50% majority rule consensus tree and to estimate the Bayesian posterior probabilities (PP).

Then, phylogenetic relationships within Solemyida were assessed based on 2 mitochondrial (16S rRNA and COX1) and 3 nuclear (18S rRNA, 28S rRNA, and H3) gene sequences from 17 in-group species (Supplementary Table S5). Nucleotide sequences of each gene were aligned in batches using MAFFT Version 7 (Katoh et al., 2019), ambiguously aligned regions were deleted using Gblocks 0.91b (Talavera and Castresana, 2007); sequences were subsequently concatenated into a single alignment used to generate nexus files in PhyloSuite 1.2.2 (Zhang et al., 2020b). Maximum likelihood analysis was performed by the GTR+F+R3 substitution model in IQ-TREE 1.6.10 (Nguyen et al., 2015), and branches evaluated by ultrafast bootstrap (BP) with 1 000 replicates. All phylogenetic trees and node labels were graphically edited with iTOL (Letunic and Bork, 2007).

3 RESULT AND DISCUSSION 3.1 Mitogenome organization and characterization

After assembly, the complete mitogenome of Acharax sp. was obtained, with a length of 18 970 bp (GenBank accession number: ON023263). This length is the largest compared with other Solemya bivalves, which range from the 15 660 bp of S. velum (Plazzi et al., 2013) to the 16 554 bp of S. pervernicosa (Russell et al., 2018). Similar to most bivalve mitogenomes, the mitogenome of Acharax sp. contains 37 genes: 13 PCGs, 2 rRNA genes, and 22 tRNA genes. Eighteen genes (18 PCGs, 10 tRNAs) were encoded by the heavy (H) strand, while the remaining genes (5 PCGs, 2 rRNAs, and 12 tRNAs) were encoded by the light (L) strand (Fig. 1; Table 1).

The nucleotide composition of the Acharax sp. mitogenome is as follows: A=32.80%, C=11.30%, G=14.80%, and T=41.10%. The GC content (26.10%) was significantly lower than the AT content (73.90%) (Table 2), and this AT richness was typical in other Solemyidae species (Table 3). Additionally, for the Acharax sp. mitogenome, the GC and AT skews were 0.134 and -0.112, respectively, which are opposite to the trends found in Solemya species and suggest a bias toward G and T (Table 3). Twenty-nine intergenic regions, totaling 4 066 bp, were scattered throughout the Acharax sp. mitogenome (Table 1). The longest noncoding region (2 357 bp) located between trnE and trnQ, which was identified as the putative control region (CR) and has high AT content (86.40%). Additionally, four overlaps were found between adjacent genes in the Acharax sp. mitogenome ranging from 1 bp to 25 bp (NAD4L/NAD4, NAD6/trnP, trnL1/rrnL, and trnV/rrnS).

The lengths of PCGs, rRNAs, tRNAs, and the CRs for the present five Solemyidae species were compared (Fig. 2). The CRs showed maximum length diversification, and this variation was considered responsible for the differences in the entire mitogenome length. In addition, the control region length in Acharax was longer than that in the Solemya species, in accordance with its longest mitogenome. Thus, even though they belong to the same family, there may be large length variation among species.

3.2 Protein-coding gene

The combined length of 13 PCGs in the Acharax sp. mitogenome was 11 139 bp, accounting for 58.72% of the mitogenome (Fig. 2; Table 3). All PCGs of the present five Solemyidae mitogenomes started with ATD (ATA, ATG, and ATT) as initiation codons, and ATG was the most common (Table 4). Most PCGs ended with two conventional stop codons (TAA and TAG), while incomplete termination codons (TA/T) were detected (Table 4). TA was found in NAD5 of S. velum, and T was found in NAD1 of S. pervernicosa and CYTB of S. velum. Incomplete termination codons are common in metazoan mitogenomes, which can be corrected by post transcriptional polyadenylation (Dreyer and Steiner, 2006; Satoh et al., 2016).

Generally, the AT skew and GC skew are considered main indicators to measure the strand asymmetry and nucleotide composition (Lü et al., 2019). For the 13 PCGs of Acharax sp., the AT and GC skews were -0.203 and 0.134, respectively (Table 2), which demonstrates that A and C were more scarce than T and G. NAD3 has the lowest AT skew and highest GC skew. The equilibrium between selection pressure and mutation pressure in the process of replication and transcription might be the main reason for nucleotide skew (Wei et al., 2010). Here, NAD3 had a greater fluctuation in the AT skew and GC skew, indicating that its mutational and selection pressure might be different from the remaining PCGs.

The Acharax sp. mitogenome encodes a total of 3 699 amino acids. Leu was the most frequently used (15.49%), while Cys was the least frequently used amino acid (Supplementary Fig.S1). The use frequency for amino acids in S. elarraichensis, S. pervernicosa, S. velesiana, and S. velum were very similar (Supplementary Fig.S1).

Numerous studies have indicated that metazoan mitogenomes often have a preference toward a higher representation of A and T, which results in a subsequent bias of the corresponding amino acids (Salvato et al., 2008; Wang et al., 2016b). The AT content of PCGs in Acharax sp. was 71.5%, the third codon position was especially A+T rich (85.9%), which was similar to the four previously published Solemya bivalves (Table 3). RSCU further demonstrated nucleotide composition bias in those aforementioned Solemyidae mitogenomes. In Acharax sp. and four published Solemya mitogenomes, the RSCU values of NNA and NNU are usually > 1 (Fig. 3), suggesting a strong AT bias in the third codon. This result supports the hypothesis that in mitogenomes the AT bias of the third codon might be correlated positively with the codon usage bias (Hao et al., 2012; Li et al., 2019a; Kou et al., 2020).

3.3 Nucleotide diversity and genetic distance

In a comparison of the five Sloemyidae mitogenomes (Acharax sp., S. elarraichensis, S. pervernicosa, S. velesiana, and S. velum), the nucleotide diversity values were highly variable among the 13 PCGs (Fig. 4). In the 13 PCGs, NAD2 (Pi=0.337), ATP8 (Pi=0.319), and NAD6 (Pi=0.318) presented high level of polymorphisms; while COX1 (Pi=0.175) and COX2 (Pi=0.219) were the most conserved genes, which had relatively low diversity values.

To reveal the divergence of the PCGs among solemyid bivalves, the pairwise genetic distances were calculated (Fig. 5). The results showed that the genetic distance at the third nucleotide position was significantly higher than that at the first+second nucleotide position, suggesting that the evolution of the first and second positions was slower than that of the third position. The highest genetic distances were detected in ATP6 (1.739) and NAD2 (1.345) at the third nucleotide of codons, while explored in ATP8 (0.564) based on the first+second nucleotide position (Fig. 5). The COX1 and NAD1 showed low genetic distance in both the first+second and third codon analyses, which indicated that they might have low evolutionary rates among the solemyid bivalves, while those of NAD2, NAD5, and ATP6 were high (Fig. 5).

The AliGROOVE analysis of the COX1 gene sequences for 66 protobranch bivalves showed strongly divergent patterns among species of different orders and thus they can be used for further phylogenetic studies of Protobranchia. The interspecific genetic divergence of COX1 between Acharax sp. and congeneric species A. johnsoni was 19.6%, and the average interfamilial genetic divergence was 24.9% (Fig. 6).

3.4 Evolution rates and adaptive characters of PCGs

The value of nonsynonymous substitution (Ka)/synonymous substitution (Ks) is a common index to evaluate the selection pressure of species in molecular evolution (Zhu et al., 2017). Ka/Ks < 1, =1, and > 1 indicate purifying selection, neutrality, and positive selection, respectively (Yang and Bielawski, 2000). Mitochondria are the main place of aerobic respiration, being important for energy metabolism and are strongly functionally constrained; therefore, purifying selection is considered the dominant force in mitogenome evolution (Tomasco and Lessa, 2011). All 13 PCGs of five Solemyidae mitogenomes were under powerful purifying selection with Ka/Ks values far less than 1. The ATP8 and NAD6 genes had high Ka/Ks values (mean: 0.305, 0.235) compared to other genes, while those of the COX1 and COX3 genes were low (Fig. 7). Researchers have suggested that due to DNA repair mechanisms, low mutation rates tend to occur in highly expressed genes (Barrientos et al., 2002). That is, COX1 and COX3, which have low Ka/Ks, may have higher expression levels than other mitochondrial genes.

However, weak and/or episodic positive selection may occur against this background of strong purifying selection under extremely harsh environments. Studies have found evidence of unique evolutionary processes and adaptations of mitochondrial PCGs, such as the ATP synthase genes of the artiodactylian tribe Caprini and Tibetan loaches (Hassanin et al., 2009; Wang et al., 2016a), the cytochrome b gene in alpaca living at high altitudes (da Fonseca et al., 2008), the cytochrome c oxidase genes of Tibetan antelope and plateau pika (Xu et al., 2005; Luo et al., 2008), the NADH dehydrogenase genes of Tibetan horses and high altitude snub-nosed monkeys (Xu et al., 2007; Yu et al., 2011). Similarly, the deep-sea is a harsh environment characterized by high hydrostatic pressure, hypoxia, darkness, and low temperature (Rex, 1981). Positively selected sites were also identified in mitogenomes of deep-sea organisms (Zhang et al., 2017a, 2020a; Sun et al., 2018; Yang et al., 2021). Bivalvia is among the most diverse molluscan classes, inhabiting various aquatic environments, from shallow continental shelves to the deep sea (Bieler et al., 2013; Lemer et al., 2019). Inhabitants living in different depth ranges would suffer to varying degree from the effects exerted by hydrostatic and temperature pressures. In order to study how mitochondria adapted to the aforementioned pressures, the amino acid constitutions of PCGs were compared between bivalves from shallow water (45 species) and the deep sea (18 species). The results showed that the nonpolar amino acids in the PCGs of these deep-sea mitogenomes (56.74%) were markedly more plentiful than those from the shallow water (55.48%) (P < 0.01) (Fig. 8a); accordingly, the percentages of polar uncharged amino acids and charged amino acids were markedly lower in the deep-sea batch (P < 0.05) (Fig. 8b, c).

Previous studies noted that the the function of proteins is affected by the composition and properties of amino acids (Metpally and Reddy, 2009; Yang et al., 2015). Our research revealed significant differences in nonpolar amino acid compositions between deep-sea bivalves PCGs and those of shallow water. Studies of amino acid composition in deep-sea amphipods have also obtained similar results, that is, with the increase of depth, the proportion of nonpolar amino acids in the mitochondrial PCGs of deep-sea amphipods present an increasing trend (Li et al., 2019b). The protein products of mitochondrial coding genes are embedded in the hydrophobic lipid chains of the mitochondrial membrane, and an increase proportion of nonpolar amino acids would contribute to compact interactions between lipid chains and membrane proteins (Li et al., 2019b). Recently, the adaptive molecular mechanisms of organisms in the deep-sea environment have already come to the attention of researchers. In addition to the differences in amino acids composition and properties, amino acid substitutions of certain proteins have also been found in the mitochondria of deep-sea scale worms, vesicomyids, crabs, starfish, snailfish, and so on (Zhang et al., 2017b, 2020a; Mu et al., 2018; Wang et al., 2019; Yang et al., 2019), which provides strong evidence of adaptive changes in mitogenomes of the deep-sea organisms.

3.5 RNA gene

The classical set of 22 tRNAs was detected in Acharax sp. mitogenome, with size ranging of 64 (trnD)–70 bp (trnV) (Table 1). The total length, A+T content, and AT skew of the 22 tRNAs in Acharax sp. mitogenome were 1 478 bp, 69.8% and 0.026, respectively (Tables 2 & 3). The secondary structure of 22 tRNAs were schematized in Supplementary Fig. S2. Furthermore, many noncomplementary and T-G base pairs were detected in the stem regions of the tRNA genes, which seem to be universal phenomena in many organisms and could be repaired through post-transcriptional editing (Lavrov et al., 2000).

Two ribosomal RNA genes (rrnS and rrnL) of Acharax sp. are located on the L strand, with lengths of 893 bp (AT%=70.6%) and 1 418 bp (AT%=73.8%), respectively. rrnS is located between trnV and trnM, and rrnL is located between trnL1cta and trnV, and the pattern rrnL-trnV-rrnS has also been found in the four previously published solemyid mitogenomes (Plazzi et al., 2013; Russell et al., 2018). A comparison of the 22 tRNA and 2 rRNA genes among solemyids indicated an extremely high similarity in different species from the same family (Table 3).

3.6 Gene rearrangement

Gene order rearrangements in molluscs have been discussed in detail in numerous studies (Cunha et al., 2009; Guerra et al., 2018; Uribe et al., 2018; Zhang et al., 2021) and vary widely within and among classes and orders. However, within Protobranchia, only gene order patterns of five species have been published previously. Among them, S. elarraichensis, S. velesiana, S. velum, and Acharax sp. shared identical mitochondrial gene orders, while S. pervernicosa showed differences in the arrangement of trnT.

Compared to the molluscan ancestral pattern, at least ten genes were rearranged in Acharax sp. mitogenome. In the mitogenomes of more primitive taxa, trnD is located between COX2 and ATP8, but in Acharax sp. it was found between ATP6 and ATP8. In addition, the trnD-ATP8-ATP6-trnF-NAD5-trnH-NAD4-NAD4L-trnT-trnS2 cluster and trnP were involved in shifting between two strands (Fig. 9) (Plazzi et al., 2013). trnG and trnE have switched places, and trnG moved to the H strand. Furthermore, based on CREx analysis, gene rearrangement scenarios were deduced from the putative ancestral molluscan pattern to Acharax sp., being the result of consecutive events of transposition/reverse transposition and reversal (Fig. 9).

3.7 Phylogenetic analysis

The Protobranchia are intriguing bivalves in terms of early evolution, unique anatomy, and ecological diversification in deep sea (Zardus, 2002; Sharma et al., 2013). However, owing to lack of prominent shell characteristics and difficulties in collecting live material, the higher systematics and evolutionary history of Protobranchia have long been controversial. Recently, based on larger taxon sampling and sequence data from molecular loci, the monophyly of Protobranchia has been supported (Lemer et al., 2019; Sato et al., 2020). Here, the phylogenetic reconstruction of bivalve lineages based on nucleotide sequences of 13 PCGs by ML and BI resulted in almost identical topologies. As a representative of Protobranchia, Solemyida was monophyletic and at the most basal position (Fig. 10). Then, we performed phylogenetic analyses based on sequence data of nuclear 18S rRNA, 28S rRNA, H3, and mitochondrial 16S rRNA, and the COX1 genes of solemyids. In Solemyoidea, the two genera Acharax and Solemya formed a monophyletic clade each. Acharax sp. clustered with other Acharax bivalves with high bootstrap support values (Fig. 11).

Considering that the present taxon coverage of existing mitogenomic trees is undeniably limited, the higher systematics of Protobranchia remains less clear. Deeper phylogenomic focus on this group is urgently needed, especially among members of Nuculida and Nuculanida. Hence, more comprehensive taxon samplings and mitogenomes from Protobranchia are needed to move us closer to the goal of reconstructing the natural evolutionary history of protobranch bivalves.

4 CONCLUSION

In the present study, the complete mitogenome Acharax sp. was determined and compared with other available mitogenomes of solemyids. The mitogenome of Acharax sp. is 18 970 bp in length and encoded 37 typical mitochondrial genes. The Ka/Ks ratio analyses indicated that five Solemyidae mitogenomes were suffering powerful purifying selection with Ka/Ks values far less than 1. Compared to 45 existing shallow water equivalents, the 18 accessible mitogenomes from the deep sea, including the Acharax sp., show significantly more nonpolar amino acids in the 13 PCGs, which may be a clue of an adaptation to the deep-sea environment. Compared to the molluscan ancestral pattern, at least 10 genes were rearranged in Acharax sp. mitogenome, as a result of consecutive events of transposition/reverse transposition and reversal. Phylogenetic analysis supported that Solemyida monophyletic and Acharax sp. clustered with other Acharax bivalves with high support values. The data presented here provides new information about mitogenome organization, adaptive characters, and phylogenetic position for solemyids, which may shed light on the understanding of the mitogenomic adaptation and evolution in protobranch bivalves.

5 DATA AVAILABILITY STATEMENT

The sequence data that support the findings of this study are openly available in GenBank of NCBI at https://www.ncbi.nlm.nih.gov.

6 ACKNOWLEDGMENT

The authors thank the captain and crew of the R/V Tan Suo Yi Hao and the pilots of the HOV Shen Hai Yong Shi for providing technical support and thank the Center for Ocean Mega-Science, Chinese Academy of Sciences.

Electronic supplementary material

Supplementary material (Supplementary Tables S1–S5 and Figs.S1–S2) is available in the online version of this article at https://doi.org/10.1007/s00343-022-2217-8.