The first complete mitogenome of Acharax sp.(Protobranchia, Solemyida, Solemyidae): comparisons with other Solemyidae bivalves and deep-sea adaptive characteristics*

Mei YANG, Jixing SUI, Xinzheng LI,3,4,**

1 Department of Marine Organism Taxonomy & Phylogeny, Institute of Oceanology, Chinese Academy of Sciences, Qingdao 266071, China

2 Center for Ocean Mega-Science, Chinese Academy of Sciences, Qingdao 266071, China

3 University of Chinese Academy of Sciences, Beijing 100049, China

4 Laboratory for Marine Biology and Biotechnology, Qingdao National Laboratory for Marine Science and Technology,Qingdao 266237, China

Abstract Solemyidae is an ancient group of protobranch bivalves, and most solemyids are symbiotic with chemoautotrophic and gill-hosted bacteria, enabling them to survive in unusual habitats such as deepsea chemosynthetic environments.However, evolution of the mitogenomes in this family and their phylogenetic relationships remain poorly understood.The complete mitogenome of Acharax sp.was determined and compared with other available mitogenomes of solemyids.The mitogenome of Acharax sp.is 18 970 bp in length and consists of 13 protein-coding genes, 2 ribosomal RNA genes, and 22 transfer RNA genes.The gene arrangement was identical to those of other sequenced solemyids.For the present five mitogenomes of Solemyidae species, all protein-coding genes were initiated with the typical ATD (ATA, ATG, and ATT) codon and terminated with the TAA/TAG codon.Furthermore, the incomplete termination codon was detected.The Ka/Ks ratio analyses indicated that 13 protein-coding genes of five Solemyidae mitogenomes suffered strong purifying selection.Compared to 45 existing shallow water equivalents, the 18 available mitogenomes from the deep-sea, including the Acharax sp.in this study,show significantly more nonpolar amino acids in the 13 protein-coding genes, which indicates the adaptation to the deep-sea environment.The phylogenetic tree based on 48 Bivalvia complete mitogenomes provided further information to support the scientific classification of protobranchs.The relationships among Solemyidae were assessed based on 2 mitochondrial (16S rRNA and COX1) and 3 nuclear (18S rRNA, 28S rRNA, and histone H3) gene sequences from 17 in-group species.The two genera Acharax and Solemya formed a monophyletic clade each, and Acharax sp.clustered with previously reported Acharax bivalves with high support values.

Keyword: Solemyidae; Acharax; mitogenome; deep sea; adaptation; phylogeny

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,AcharaxDall, 1908 andSolemyas.l.Gray 1840.These are distinguished by the position of the ligament: external inAcharaxwhile internal inSolemyas.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).Acharaxgenerally has a deeper distribution (400–5 379 m) associated with reducing environments (Neulinger et al., 2006),whileSolemyais found mainly in the continental shelf and upper slope (0–600 m) (Neulinger et al.,2006; Oliver et al., 2011), however, with further research,Solemyahas 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(Solemyaelarraichensis,S.pervernicosa,S.velesiana,andS.velum), and no mitogenome of the genusAcharaxhas 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 ofAcharaxsp.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 theCOX1gene; (ⅲ)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 genusSolemya, 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

TheAcharaxsp.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 genusAcharaxsp.(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(18SrRNA,28SrRNA, andH3) 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 ddH2O.The cycling conditions were as follows: initial denaturation for 2 min at 94 °C,followed by 35 cycles of denaturation at 94 °C for 1 min, annealing at 50 °C for 30 s, extension at 72 °C for 50 s, and a final extension at 72 °C 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 (TruSeqTMNano 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 theAcharaxsp.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 theCOX1gene 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-tailedttests 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 ofAcharaxsp.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 ofAcharaxsp.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 IQTREE (Nguyen et al., 2015), and branch support was assessed by ultrafast bootstrapping with 1 000 replicates (Hoang et al., 2018) as well as Bayesianlike 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 (16SrRNAandCOX1) and 3 nuclear (18SrRNA,28SrRNA,andH3) 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 ofAcharaxsp.was obtained, with a length of 18 970 bp (GenBank accession number: ON023263).This length is the largest compared with otherSolemyabivalves, which range from the 15 660 bp ofS.velum(Plazzi et al., 2013) to the 16 554 bp ofS.pervernicosa(Russell et al., 2018).Similar to most bivalve mitogenomes, the mitogenome ofAcharaxsp.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 theAcharaxsp.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 theAcharaxsp.mitogenome, the GC and AT skews were 0.134 and -0.112, respectively, which are opposite to the trends found inSolemyaspecies and suggest a bias toward G and T (Table 3).Twentynine intergenic regions, totaling 4 066 bp, were scattered throughout theAcharaxsp.mitogenome(Table 1).The longest noncoding region (2 357 bp)located betweentrnEandtrnQ, 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 theAcharaxsp.mitogenome ranging from 1 bp to 25 bp(NAD4L/NAD4,NAD6/trnP,trnL1/rrnL, andtrnV/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 inAcharaxwas longer than that in theSolemyaspecies, 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

Fig.1 The organization of the mitogenome of Acharax sp.

The combined length of 13 PCGs in theAcharaxsp.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 inNAD5ofS.velum, and T was found inNAD1ofS.pervernicosaandCYTBofS.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 ofAcharaxsp., 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.NAD3has 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,NAD3had a greater fluctuation in the AT skew and GC skew, indicating that its mutational and selection pressure might be different from the remaining PCGs.

TheAcharaxsp.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 inS.elarraichensis,S.pervernicosa,S.velesiana, andS.velumwere 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 inAcharaxsp.was 71.5%, the third codon position was especially A+T rich (85.9%),which was similar to the four previously publishedSolemyabivalves (Table 3).RSCU further demonstrated nucleotide composition bias in those aforementioned Solemyidae mitogenomes.InAcharaxsp.and four publishedSolemyamitogenomes, 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).

Table 1 Annotation of the mitogenome of Acharax sp.

Table 2 Base composition of the Acharax sp.mitogenome

Table 3 Genomic features of solemyids mitogenomes

Fig.2 The length of protein-coding genes (PCGs), tRNAs, rRNAs,and control regions among 5 Solemyidae mitogenomes

3.3 Nucleotide diversity and genetic distance

In a comparison of the five Sloemyidae mitogenomes (Acharaxsp.,S.elarraichensis,S.pervernicosa,S.velesiana, andS.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), andNAD6(Pi=0.318)presented high level of polymorphisms; whileCOX1(Pi=0.175) andCOX2(Pi=0.219) were the most conserved genes, which had relatively low diversity values.

Table 4 Start and stop codons of the Solemyidae mitogenomes

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 inATP6(1.739) andNAD2(1.345) at the third nucleotide of codons, while explored inATP8(0.564) based on the first+second nucleotide position (Fig.5).TheCOX1andNAD1showed 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 ofNAD2,NAD5, andATP6were high (Fig.5).

The AliGROOVE analysis of theCOX1gene 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 ofCOX1betweenAcharaxsp.and congeneric speciesA.johnsoniwas 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.TheATP8andNAD6genes had high Ka/Ks values (mean: 0.305, 0.235)compared to other genes, while those of theCOX1andCOX3genes 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,COX1andCOX3, 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 deepsea 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.

Fig.4 Nucleotide diversity (Pi) of 13 PCGs of the present five Solemyidae bivalves

Fig.5 The pair genetic distances of 13 PCGs among five Solemyidae bivalves

3.5 RNA gene

Fig.6 AliGROOVE analysis of COX1 gene sequences for 66 protobranch bivalves considering their nucleotide composition

The classical set of 22 tRNAs was detected inAcharaxsp.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 inAcharaxsp.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 (rrnSandrrnL) ofAcharaxsp.are located on the L strand, with lengths of 893 bp (AT%=70.6%) and 1 418 bp (AT%=73.8%),respectively.rrnSis located betweentrnVandtrnM,andrrnLis located betweentrnL1ctaandtrnV, and the patternrrnL-trnV-rrnShas 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

Fig.7 The Ka/Ks ratios for each PCG in pairwise mitogenome of Acharax sp.(Asp), S.elarraichensis (Se), S.pervernicosa(Sp), S.velesiana (Sva), and S.velum (Svm)

Fig.8 Statistics of amino acids composition in mitochondrial PCGs of bivalves from different habitats

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, andAcharaxsp.shared identical mitochondrial gene orders, whileS.pervernicosashowed differences in the arrangement oftrnT.

Compared to the molluscan ancestral pattern, at least ten genes were rearranged inAcharaxsp.mitogenome.In the mitogenomes of more primitive taxa,trnDis located betweenCOX2andATP8, but inAcharaxsp.it was found betweenATP6andATP8.In addition, thetrnD-ATP8-ATP6-trnFNAD5-trnH-NAD4-NAD4L-trnT-trnS2cluster andtrnPwere involved in shifting between two strands(Fig.9) (Plazzi et al., 2013).trnGandtrnEhave switched places, andtrnGmoved to the H strand.Furthermore, based on CREx analysis, gene rearrangement scenarios were deduced from the putative ancestral molluscan pattern toAcharaxsp.,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 nuclear18SrRNA,28SrRNA,H3,and mitochondrial16SrRNA, and theCOX1genes of solemyids.In Solemyoidea, the two generaAcharaxandSolemyaformed a monophyletic clade each.Acharaxsp.clustered with otherAcharaxbivalves with high bootstrap support values (Fig.11).

Fig.9 The rearrangement scenarios deduced from the molluscan ancestral pattern to Acharax sp.by CREx analysis

Fig.10 Phylogenetic tree deduced from the nucleotide sequences of 13 PCGs by Bayesian and ML methods

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 mitogenomeAcharaxsp.was determined and compared with other available mitogenomes of solemyids.The mitogenome ofAcharaxsp.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 theAcharaxsp., 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 inAcharaxsp.mitogenome, as a result of consecutive events of transposition/reverse transposition and reversal.Phylogenetic analysis supported that Solemyida monophyletic andAcharaxsp.clustered with otherAcharaxbivalves 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.

Fig.11 Phylogenetic relationships of Solemyida by the ML analysis of mitochondrial (16S rRNA+COX1) and nuclear (18S rRNA+28S rRNA+H3) sequences

5 DATA AVAILABILITY STATEMENT

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

6 ACKNOWLEDGMENT

The authors thank the captain and crew of the R/VTanSuoYiHaoand the pilots of the HOVShen HaiYongShifor providing technical support and thank the Center for Ocean Mega-Science, Chinese Academy of Sciences.