The samples for whole-genome sequencing of S. rivularis and N. chinensis were collected from the coastal area of Quanzhou in Fujian Province and the Qiantang River in Zhejiang Province, respectively, following the relevant guidelines and regulations. A total of ten individuals each of S. rivularis and N. chinensis were sampled, with average shell length of 55.98 ± 3.47 mm and 45.41 ± 2.74 mm, respectively. All specimens were preserved in 85% ethanol as voucher specimens. These specimens were deposited at Zhejiang Key Laboratory of Aquatic Germplasm Resource, Zhejiang Wanli University, Ningbo, China.

Raw genome reads were acquired through both Illumina HiSeq sequencing and PacBio Sequel IIe third-generation sequencing (unpublished), and assembled for the mitogenomes of these two species. Initially, a de novo mitogenome assembly was carried out with SPAdes v3.9.0 after filtering the unqualified reads by Trimmomatic v. 0.39 (Bankevich et al. 2012; Bolger et al. 2014). The scaffold sequences were then obtained by extending the contigs using SSPACE. The assembly quality was evaluated by GetOrganelle software (Jin et al. 2020). Finally, the MitoZ program was used to annotate the protein-coding genes (PCGs), two ribosomal RNAs (rRNAs) and transfer RNAs (tRNAs) (Meng et al. 2019).

The content and proportion of nucleotide bases were analyzed by MEGA 11. The base skew values were calculated according to the formulae: AT-skew = (A − T) / (A + T) and GC-skew = (G − C) / (G + C). The RSCU of the two mitogenomes was counted using PhyloSuite v1.2.3.

To explore the evolutionary relationship of S. rivularis and N. chinensis, the published mitogenome sequences of Solenoidea and Hiatelloidea were retrieved from GenBank, and Solecurtus divaricatus was selected as the outgroup (Table 1). The phylogenetic analysis was performed using PhyloSuite software (Zhang et al. 2020). First, using an invertebrate mitochondrial code table, MAFFT was used to independently align 12 protein-coding genes. The ATP8 gene was excluded due to its deletion in the majority of mollusks. Poorly aligned regions of the sequences were pruned by Gblocks under default parameters. The resulting alignments were then concatenated and transferred to ModelFinder for the best model prediction. Phylogenetic trees were estimated through maximum likelihood (ML) and Bayesian inference (BI) methods. The ML phylogenetic tree was generated using IQ-Tree with 1000 bootstrap replicates. The BI analyses were performed by MrBayes 3.2.6 with Markov Chain Monte Carlo (MCMC) for 5000,000 generations. The first 25% of trees were discarded as burn-in and the sampling was terminated when the convergence value was less than 0.01. The iTOL tool was exploited to visualize the phylogenetic tree (https://itol.embl.de/).

In addition, the most plausible gene order rearrangement events that might have occurred between Pharidae and Solenidae were reconstructed by pairwise comparisons of mitogenomes through the Common Interval Rearrangement Explorer (CREx) (Bernt et al. 2007).

The branch-site model was used to analyze the selection pressure on 12 PCGs of razor clams in the PAML package. In this model, N. chinensis was marked as the foreground branch to investigate the evolutionary adaptation between freshwater and marine species. The null model (model = 2, Nssites = 2, fix_omega = 1, omega = 1) and alternative model (model = 2, Nssites = 2, fix_omega = 0, omega = 2) were compared by likelihood ratio test (LRT). Subsequently, P-values were calculated through the chi-square distribution. Then, the posterior probability of the amino acid sites under positive selection was calculated according to the Bayesian empirical Bayes (BEB) method. The inference of positively selected sites was based on a posterior probability of greater than 95%.

The lengths of S. rivularis and N. chinensis mitogenomes were 17,159 bp and 15,957 bp, respectively (Fig. 1A). Both mitogenomes contain 12 PCGs, 22 tRNAs, and 2 rRNAs, all of which were located on the heavy chain. The ATP8 gene was missing in this two mitogenomes. Their composition was similar to that of other species in Pharidae, indicating a certain degree of conservation in this family. The detailed genes information was shown in Table 2. The base composition of S. rivularis and N. chinensis mitogenomes was displayed in Table 3 with AT contents of 66.80% and 71.50%, respectively, both of which exhibited an obvious AT bias. The AT content of N. chinensis was the highest among the published Adapedonta mitogenomes. In addition, the two mitogenomes all exhibited negative AT-skew and positive GC-skew, reflecting that the base composition ratios were A biased to T, and G biased to C. There were some differences in the types of start and termination codons of 12 PCGs between the two species (Table 2). Specifically, the start codons of 12 genes in S. rivularis were found to be ATN, TTG and GTG types, whereas in N. chinensis, all genes began with the codon ATN, with the exception of the ND4 gene, which used TTG as the start codon. Concerning termination codons, six genes in S. rivularis (cytb, atp6, cox3, nad4, nad3, nad1) and seven genes in N. chinensis (cytb, nad6, atp6, cox3, nad4l, nad3, nad1) were detected TAA or TAG at the sequence end. The remaining genes featured an incomplete termination codon consisting of a T that might be complemented into a complete stop codon by polyadenylation following transcription to the resultant mRNA (Ojala et al. 1981). Furthermore, the non-coding regions of the mitochondrial genomes of N. chinensis and S. rivularis account for 11.92% and 19.33%, respectively. The longest non-coding region (NCR) of N. chinensis and S. rivularis was both located between nad2 and trnK, with lengths of 443 bp and 1,639 bp respectively, which was identified as a putative control region (CR).

Figure 1. 

Maps of A the mitogenomes of S. rivularis and N. chinensis and their B RSCU.

As illustrated in Fig. 1B, the preferred codons for 22 amino acids of two species ended in A or U, consisting with the result of AT bias of the mitogenome sequence. As a consequence of the duplication of tRNA-Leu and tRNA-Ser, Leu and Ser were each encoded by six and eight codons, respectively. The most frequently used codons were UUA (Leu2), UCU (Ser2), GCU (Ala) and CCU (Pro). Compared to S. rivularis, CUG (Leu1), AUC (Ile), AAC(Asn) were utilized to a lesser extent in N. chinensis.

The 12 protein-coding genes from 15 taxa were concatenated to generate a sequence matrix of 10,806 bp. The tree topologies derived from the ML and BI analyses were largely congruent exhibiting high posterior probabilities (PP) and bootstrap support values (BS) in most nodes (Fig. 2). Phylogenetic analyses revealed that the genus Hiatella from Hiatelloidea was closely related to the superfamily Solenoidea, indicating a close evolutionary relationship between them. Additionally, both analyses strongly confirmed the monophyly of Solenoidea, which was divided into two major branches, Solenidae and Pharidae. In the family Pharidae, the genus Sinonovacula (including S. rivularis and S. constricta) was clustered alongside N. chinensis, with Cultellus attenuatus emerging as a sister group. Siliqua minima and Ensis leei were clustered in a separate cluster.

Figure 2. 

The phylogenetic trees based on concatenated 12 mitochondrial PCGs, and the gene orders of Adapedonta species. Values shown next to nodes are posterior probabilities (left) and ML bootstrap support values (right). Newly assembly mitogenomes are marked with triangles. Except for Panopea abrupta (https://inverts.wallawalla.edu) and Panopea globosa (Góngora-Gómez et al. 2016), the images of the other species are all from https://www.inaturalist.org.

The mitogenomes of Solenoidea all exhibited the identical composition of 12 PCGs, 22 tRNAs, and 2 rRNAs, except for Ensis leei, which contained an additional ATP8 gene (Fig. 3). The gene arrangement was consistent within each family, and there was a certain level of conservation in gene arrangement between Solenidae and Pharidae. A large block, rrnL-ATP6-M-rrnS-cox3, and five small blocks, L2-nad1-L1, S-nad2, nad5-cytb, I-D, Q-C were shared by both families, providing further evidence of the close lineage relationship observed in the phylogenetic analysis of this study. The CREx analysis suggested that three transposition and four tandem duplication random losses (TDRLs) might have occurred between Pharidae and Solenidae.

Figure 3. 

Putative gene rearrangement events between Pharidae and Solenidae. Green and red lines represent transposition and TDRL events, respectively, which were step by step identified by CREx.

The species of Solenoidea were selected for molecular evolution analysis, with N. chinensis designated as the foreground branch (Fig. 4). The branch-site model (BSM) in the PAML package was employed to detect positively selected genes (PSGs). As illustrated in Table 4, the substitution model A was significantly better than the neutral selection model null in nad5, indicating that this gene underwent positive selection in the foreground branch (P < 0.05). According to the BEB analysis, there were five positive selection sites in the nad5 amino acid sequences (140 A 0.509, 143 F 0.547, 144 L 0.865, 442 A 0.700, 446 F 0.620). Moreover, discrepancies were observed in the 144^th site between freshwater N. chinensis (Ala) and marine razor clams (Leu) (Fig. 5). However, the evidence for each site was somewhat inconclusive. These findings suggest that the nad5 gene may have played a pivotal role in the adaptive evolution of freshwater environments.

Figure 4. 

Phylogenetic tree of Solenoidea for selective stress analysis. The branch marked in red is the foreground branch.

Figure 5. 

The difference of the 144^th positive selected amino acid site in NAD5 of eight Solenoidea species. The 144^th site is indicated by a red frame.

The mitogenomes of S. rivularis and N. chinensis were newly assembled, with lengths of 17,159 and 15,957 bp, respectively. In compared with the previously sequenced Adapedonta mtDNA size (ranged from 15,381 bp to 19,507 bp), their mitogenome sizes were within the normal range (Zheng et al. 2010; Yuan et al. 2012b; Feng et al. 2021; Li et al. 2022). Notably, the genome size of N. chinensis was the smallest in the family Pharidae, which was associated with the variation in length of the control region. The CR is the region with the largest sequence and length variation in the mitogenome, and has the fastest evolution, which is crucial for the regulation of mitochondrial DNA replication and transcription (Wolstenholme 1992; Boore 1999). The substantial differences in the content and structure of the control region within the mollusk lineage provide valuable insights for population genetic analysis (Sasuga et al. 1999; Tomita et al. 2002; Kawashima et al. 2013). Among the published mitogenomes of Pharidae, there is a large control region between nad2 and trnK, such as S. constricta (1,492 bp), S. minima (1,371 bp), C. attenuatus (1,173 bp) and E. leei (1,101 bp) (Zheng et al. 2010; Feng et al. 2021; Li et al. 2022). In this study, S. rivularis displayed a moderately larger control region size of 1,639 bp, whereas it was only 441 bp in N. chinensis, making it a different mitogenome size in the family Pharidae. Intriguingly, a similar control region was not observed in the species of Solenidae (Yuan et al. 2012a, b). This distinction provides evidence for the taxonomic division of the subfamily of Solenoidea.

The topological tree constructed from the 12 mitochondrial PCGs sequence based on the BI and ML methods yielded consistent results, demonstrating that Solenoidea is clearly divided into Solenidae and Pharidae, which is consistent with the prior research results (Yuan et al. 2012d; Feng et al. 2021). Previously, S. rivularis was identified as a new species of Sinonovacula distinct from S. constricta based on morphological studies and a comparative analysis of COI and 16SrRNA fragments (Huang and Zhang 2007; Weng et al. 2013). In this research, this classification view was supported at the level of mitogenomes, and Sinonovacula belonged to the family Pharidae (Adapedonta: Solenoidea). In addition, N. chinensis was previously classified into Solecurtidae, whereas the results of this study demonstrated that N. chinensis and Sinonovacula are clustered together, forming a novel branch in the family Pharidae, which was consistent with the taxa in WoRMS (Liu 1979; Appeltans et al. 2012). Recently, the phylogenetic tree and molecular clock of tandem mitochondrial gene and nuclear gene (COI, 16S, 28S) revealed that Siliqua, Sinonovacula, Cultellus, and Novaculina belonged to Pharellinae (Bolotov et al. 2018b). However, Cultellus and Siliqua were categorized into the subfamily Cultellinae and Siliquinae, respectively, by Ahyong (Appeltans et al. 2012). Previously, Pharidae were divided into four subfamilies: Pharinae (Nasopharus, Pharus, Sinupharus), Cultellinae (Afrophaxas, Cultellus, Ensis, Ensiculus, Phaxas, Sinucultellus), Siliquinae (Siliqua), and Pharellinae (Novaculina, Orbicularia, Pharella, Sinonovacula). However, the present results indicate that Pharidae are divided into two clades, in which Cultellus is clustered alongside the Sinonovacula and Novaculina, while Siliqua and Ensis clustered together. These observations reflect that the current categorization of the subfamily Pharidae requires further research and refinement, particularly in combination with more species information.

Unlike stable gene arrangements of Vertebrata and Arthropoda, the gene orders of all genes within mtDNA exhibit considerable variability in every major molluscan lineage, including Cephalopoda, Bivalvia, Scaphopoda, and Monoplacophora (Rawlings et al. 2001; Dreyer and Steiner 2004; Yuan et al. 2012c; Stöger et al. 2016; Ma et al. 2023). Gene rearrangements may be caused by reverse transpositions, transpositions, inversions, and TDRL, which can provide important clues about the evolutionary history of species (Boore and Brown 1998; Serb and Lydeard 2003; Wang et al. 2021). In this paper, CREx analysis predicted that three transpositions and four TDRLs might have occurred between Pharidae and Solenidae, implying that dramatic mitogenome changes occurred during species differentiation. Moreover, the gene order illustration of Adapedonta revealed that species with a closer genetic relationship tended to share a similar gene arrangement, indicating that there is a potential relationship between evolution and gene rearrangement (Fig. 2). However, in this study, three distinct gene arrangement types were observed in the family Hiatellidae, especially in the genus Hiatella with nad3 and nad1 transpositions in terms of 12 PCGs arrangement (Fig. 2). The similar case that different gene arrangements in the same genus has also been reported in the genera Dendropoma and Crassostrea (Rawlings et al. 2010; Ren et al. 2010). Therefore, the taxonomic evolution of species cannot be substantiated exclusively through the examination of gene sequences; it also necessitates the integration of phylogenetic reconstruction.

Pharidae is a major marine family, with the exception of Novaculina, that is a relict marine-derived freshwater lineage (Annandale 1922; Bolotov et al. 2018a). The branch-site model study was used to determine whether positive selection occurs at a few places in freshwater razor clam. The results suggested that the nad5 gene underwent positive selection. NADH dehydrogenase is the initial and most substantial enzyme complex in the respiratory chain, functioning as a proton pump (da Fonseca et al. 2008). Nad2, nad4, and nad5 are considered to be the actual proton pumping devices because of their sequence homology with a class of Na⁺/H⁺ antiporters (Brandt 2006). The efficiency of the proton transfer process may be interfered by the mutation of the complex, which could be a crucial factor in adaptive evolution (Hassanin et al. 2009; Yu et al. 2011). For instance, the outcomes of positive selection sites in mussels from disparate habitats reflected that the p-value of nad4 was significant in freshwater branches and six sites were identified as positive sites with BEB analysis (> 95%), which implies that nad4 may contribute to the adaptation of Limnoperna fortunei in freshwater (Zhao et al. 2022). Significant non-synonymous changes were detected in the cytb and nad5 genes by comparing mitogenomes of panpulmonate gastropods that are distributed from marine to intertidal and terrestrial habitats (Romero et al. 2016). Therefore, the positive selection of nad5 gene in N. chinensis may be the result of the adaptive evolution of freshwater environment. Moreover, divergent selection occurred at site 144 of nad5, where the amino acids Ala and Leu were identified in the freshwater Novaculina chinensis and seven marine lineages, respectively, indicating divergent evolution exists the family Pharidae. Divergent evolution is the process by which separate species with common ancestors evolve distinct features to adapt to their unique living environment, which is one of the important mechanisms for the formation of biodiversity (Gautam 2020). However, the evidence supporting the positive selection of individual nad5 sites is insufficient. To provide more robust statistical support for the differences in evolutionary adaptation between freshwater and seawater species, it is necessary to include more freshwater razor clam sequences.

The authors have declared that no competing interests exist.

No ethical statement was reported.

This work was supported by Zhejiang Major Program of Science and Technology (2021C02069-7), Ningbo Major Project of Science and Technology (2021Z114) and National Marine Genetic Resource Center Program.

Yinghui Dong and Demin Zhang conceived and designed the experiments. Yiping Meng performed the experiments and drafted the manuscript. All authors reviewed the paper.

All of the data that support the findings of this study are available in the main text.