Phylogeny of Strombidae (Gastropoda) Based on Mitochondrial Genomes

Introduction

As the largest group within Stromboidea, the family Strombidae comprises about 120 extant species (MolluscaBase, 2022a). They are widely distributed in tropical and subtropical seas all over the world, and some species possess high ecological and economic values. For example, the west Atlantic strombid Aliger gigas is considered one of the most important fishery resources in the Caribbean (Machkour-M’Rabet et al., 2021). For the reason of overfishing, A. gigas has been listed as a vulnerable commercial species in the CITES Appendix II since 1992. The Indo-Pacific strombid Conomurex luhuanus is also an exploited species in coastal countries (Ulm et al., 2019). The shells of strombids varied greatly, from small and fusiform to large and decorated with a strongly stretched outer lip (Latiolais et al., 2006). However, the morphological diversity of Strombidae made it difficult to establish a proper classification system, especially at the subgenus or genus level (Kronenberg and Vermeij, 2002).

The traditional classification divided Strombidae into several genera, within which Strombus and Lambis were the two most species-rich groups (Abbott, 1960; Abbott, 1961). Although Strombus and Lambis differ greatly in shell morphology, they possess similar characteristics of soft tissue anatomies, egg masses, and radulae (Abbott, 1961). The subsequent morphological (Stone, 2001) and molecular (Latiolais et al., 2006) analyses revealed the non-monophyly of either Strombus or Lambis and called for further revision of Strombidae taxonomy, which relied on a robust phylogenetic framework.

The complete mitochondrial genome has been proved as a useful tool in resolving phylogenetic relationships of different mollusk groups (Uribe et al., 2016; Abalde et al., 2017; Uribe et al., 2019; Irisarri et al., 2020; Wang et al., 2021). Within Stromboidea, a total of 14 mt genomes have been reported and used for phylogenetic reconstructions at different levels (Irwin et al., 2021; Jiang et al., 2019; Machkour-M’Rabet et al., 2021). For example, Machkour-M’Rabet et al. (2021) sequenced one mt genome of the Mexican Caribbean strombid A. gigas and reconstructed the mitogenomic phylogeny of Hypsogastropoda including 110 species. Irwin et al. (2021) reported seven stromboid mitogenomes and analyzed them along with other published sequences to reconstruct the caenogastropod phylogeny. Both studies recovered the inclusion of Xenophoridae (represented by Onustus exutus and Xenophora japonica) into Stromboidea, which has been supported in morphological or behavioral characteristics (Berg, 1974; Simone, 2005; Simone, 2011). By far, Stromboidea is considered to comprise six extant families: Aporrhaidae, Rostellariidae, Seraphsidae, Strombidae, Struthiolariidae, and Xenophoridae (MolluscaBase, 2022b). However, the relationship between families within Stromboidea is still controversial, as different data show different topologies (Irwin et al., 2021; Machkour-M’Rabet et al., 2021). For example, one phylogeny in Irwin et al. (2021) based on mitochondrial (13 protein-coding genes (PCGs) + 2 ribosomal RNA (rRNA) genes) and nuclear genes (18S + 28S) revealed two monophyletic clades, with the first clade including Xenophoridae + (Aporrhaidae + Struthiolariidae) and the second formed by (Seraphsidae + Rostellariidae) + Strombidae. The phylogeny by Machkour-M’Rabet et al. (2021) derived from 13 mitochondrial PCGs only recovered the first clade, whereas Seraphsidae was sister to Rostellariidae + Strombidae in the second clade. Within Strombidae, an Eastern Pacific/Atlantic clade and an Indo-West Pacific clade were biogeographically structured (Machkour-M’Rabet et al., 2021), but the phylogenetic relationships between strombid genera remained controversial. Therefore, it is still necessary to improve the phylogenetic resolution of Strombidae by increased taxon sampling.

Paleontological studies indicated that strombids originated during Cenomanian-Turonian from the aporrhaids (Stromboidea). Strombids remained at very low diversity for the rest of the Cretaceous but diversified rapidly in the early Cenozoic, whereas the aporrhaids diversity was affected in the end-Cretaceous mass extinction and declined rapidly in the early Cenozoic (Roy, 1996). The strombids achieved maximum diversity during the middle Eocene, with plenty of fossils reported from the Eocene to Pliocene (Roy, 1996; Bandel, 2007).

In the present study, the complete mitochondrial genomes of Canarium labiatum, Lentigo lentiginosus, and Euprotomus aratrum were sequenced and analyzed together with those of other stromboids published before (Table 1). Following Uribe et al. (2017), two different datasets were derived from mt genomes. The first dataset contained the deduced amino acid sequences of the 13 PCGs and the nucleotide sequences of the two rRNA genes of all available Stromboidea mitogenomes, whereas the second dataset included the nucleotide sequences of the 13 PCGs and the two rRNA genes of all available Strombidae mitogenomes. Our aims were 1) to confirm the phylogenetic positions of the three species representing three genera within Strombidae, 2) to reconstruct robust internal phylogenetic relationships of superfamily Stromboidea, and 3) to date major cladogenetic events within Strombidae.

Table 1 List of mt genomes used in the present study.

Materials and Methods

Sample Collection and DNA Extraction

The specimen of C. labiatum was collected in the intertidal zone of Wuzhizhou Island (18.18′47″N; 109.45′59″E), while the specimens of L. lentiginosus and E. aratrum were sampled in a local market of Ximaozhou Island (18°14′22″N; 109°22′42″E). Samples were deposited in 95% alcohol in the Laboratory of Economic Shellfish Genetic Breeding and Culture Technology (LESGBCT), Hainan University.

Genomic DNA was extracted from small pieces of foot tissue (about 30 mg) using TIANamp Marine Animals DNA Kit (Tiangen, Beijing, China) following the instructions. Only one specimen of each species was used for DNA extraction. The genomic DNA was visualized on 1% agarose gel for quality inspection.

DNA Sequencing and Mitogenome Assembly

Genomic DNA of three species was submitted to Novogene Company (Beijing, China) for library construction and next-generation sequencing. Two sequencing libraries with average insert sizes of approximately 300 bp were prepared and then sequenced as 150-bp paired-end runs on the Illumina HiSeq 2000 platform. Finally, about 10 Gb of raw data were generated for each library. The raw reads were cleaned using Trimmomatic v.0.39 (Bolger et al., 2014) with the following settings: ILLUMINACLIP : TruSeq3-PE.fa:2:30:10 LEADING:3 TRAILING:3 SLIDINGWINDOW:5:20 MINLEN:40. The generated clean data were imported in Geneious Prime 2021.0.1 (Kearse et al., 2012) for mitogenome assembly following Irwin et al. (2021).

Genome Annotation and Sequence Analysis

The three newly determined mitogenomes were annotated with Geneious Prime. The 13 PCGs were determined by ORF Finder (http://www.ncbi.nlm.nih.gov/orffinder) with the invertebrate mitochondrial genetic code. The secondary structure of transfer RNA (tRNA) genes was generated by MITOS Webserver and modified according to the tRNA structure provided by ARWEN (Laslett and Canbäck, 2008) in Microsoft Visio 2016. The rRNA genes were identified and annotated by comparing the MITOS results and the previously published Strombidae mitogenomes.

The nucleotide composition of the mt genomes, PCGs, rRNA, and tRNA genes were computed using MEGA X (Kumar et al., 2018). The base skew values for a given strand were calculated as AT skew = (A − T)/(A + T) and GC skew = (G − C)/(G + C), where A, T, G, and C are the occurrences of the four nucleotides. Codon usage of PCGs was estimated using MEGA X. The mitochondrial genome map was generated using CGView (Grant and Stothard, 2008).

Phylogenetic Analysis

The three newly sequenced mitogenomes were aligned along with those of other stromboids available in GenBank (Table 1). Two datasets were constructed and analyzed following Uribe et al. (2017). The first dataset was comprised of the deduced amino acid sequences of the 13 mt PCGs and the nucleotide sequences of the two rRNA genes of all the available Stromboidea mitogenomes. It was aimed to test the internal phylogenetic relationships of the superfamily Stromboidea. Charonia lampas and Bufonaria rana that belong to Tonnoidea were used as outgroups (Jiang et al., 2019). The second dataset that contained the nucleotide sequences of the 13 PCGs and the two rRNA genes of all the 11 available Strombidae mitogenomes (Table 1) was intended to reconstruct the phylogeny of the family Strombidae. Terebellum terebellum (Seraphsidae) and Varicospira cancellata (Rostellariidae) from Stromboidea were used as the outgroup.

The nucleotide sequences of the 13 PCGs were aligned separately as codons using ClustalW implemented in MEGA X. The deduced amino acid sequences of the 13 PCGs were translated from the aligned codon sequences, according to the invertebrate mitochondrial genetic code. Nucleotide sequences of the rRNA genes were aligned separately using MAFFT v7 (Katoh and Standley, 2013) with default parameters. Ambiguously aligned positions were removed using Gblocks v.0.91b (Castresana, 2000) with default parameters. Finally, the different single alignments were concatenated into a single dataset in Geneious Prime 2021.0.1. Sequences were converted into different formats for further analyses using DAMBE5 (Xia, 2013).

Phylogenetic analyses were conducted using maximum likelihood (ML) (Felsenstein, 1981) and Bayesian inference (BI) analyses (Huelsenbeck and Ronquist, 2001). ML analyses were carried out using RAxML-HPC2 on XSEDE (Stamatakis, 2006) with the rapid bootstrap algorithm and 1,000 replicates using the CIPRES Science Gateway v.3.3 (Miller et al., 2010). BI analyses were conducted with MrBayes 3.2.6 (Ronquist and Huelsenbeck, 2003), running four simultaneous Monte Carlo Markov chains (MCMC) for 10,000,000 generations, sampling every 1,000 generations, and discarding the first 25% generations as burn-in. Two independent BI runs were performed to increase the chance of adequate mixing of the Markov chains and to increase the chance of detecting a failure to converge, as determined using Tracer v1.6. The effective sample size (ESS) of all parameters was more than 200.

The best partition schemes and best-fit substitution models for the two datasets were conducted using PartitionFinder 2 (Lanfear et al., 2017), under the Bayesian information criterion (BIC). For the PCGs at the amino acid level, the partitions tested were all genes combined, all genes separated (except atp6-atp8 and nad4-nad4L), and genes grouped by enzymatic complexes (atp, cob, cox, and nad). For PCGs at the nucleotide level, the partitions tested were all genes combined, all genes separated (except atp6-atp8 and nad4-nad4L), and genes grouped by subunits. Additionally, these three partition schemes were tested considering separately the three codon positions. The rRNA genes were analyzed with two different schemes (genes grouped or separated). The best-fit substitution models of the two datasets are provided in Table S3.

Estimation of Divergence Times

The divergence times within Strombidae were estimated based on the second dataset (the nucleotide sequences of the 13 PCGs and the two rRNA genes) and modified slightly (using only T. terebellum as the outgroup), using an uncorrelated, lognormal relaxed molecular clock model in BEAST v.1.10.4 (Drummond and Rambaut, 2007). For the tree prior, a Yule process of speciation was employed. The partitions selected by PartitionFinder 2 (see above) were applied. The final Markov chain was run twice for 100 million generations, sampling every 10,000 generations, and the first 10 million were discarded as burn-in, according to the convergence of chains checked with Tracer v.1.6. The ESS of all the parameters was above 200.

The posterior distribution of the estimated divergence times was obtained by specifying two calibration points that were based on fossil records as priors for divergence times of the corresponding splits. The first calibration point was set at the origin of Strombidae. A lognormal distribution was applied, with a minimum of 89.8 Mya and a 95% upper limit of 100.5 Mya (offset, 89.8; mean, 2.5; SD, 2.8) based on the probable ancestor of strombids evolved during the Cenomanian-Turonian epoch (Roy, 1996). A second calibration point was set for the divergence of Strombus and Aliger, with a minimum of 23.0 Mya and a 95% upper limit of 28.5 Mya (exponential distribution; offset, 23.0; mean, 1.5) based on the proposed pre-Miocene common ancestry between the two groups (Machkour-M’Rabet et al., 2021).

Results and Discussion

Genome Structure, Organization, and Composition

The mitochondrial genome composition of the three strombids measured in this study is shown in Tables 2, S1, and S2. The genome lengths of L. lentiginosus, C. labiatum, and E. aratrum are 16,054, 15,843, and 16,187 bp, respectively, and the differences in genome size are mainly in the non-coding regions. They encode for 13 PCGs, 22 tRNA genes, and two rRNA genes, and their gene order is the same as the stromboid consensus gene order (Irwin et al., 2021). Among all the mt genes, eight tRNAs are encoded on the minor strand, while the others are encoded in the major strand (Figure 1, represented by the mt genome of L. lentiginosus). The AT content, AT skew, and GC skew values of the whole mitochondrial genomes are shown in Table 3. The AT content of the three species varied from 64.1% to 67.9%, which indicated a high A + T bias. The negative AT skew and positive GC skew of L. lentiginosus and C. labiatum suggested the nucleotide compositions of major strands were skewed from A toward T and insignificantly skewed from C toward G, which has also been revealed in other mollusk taxa (Williams et al., 2014; Zhan et al., 2018; Sun et al., 2020). However, the GC skew of E. aratrum showed a value of −0.13, which is opposite to the other two species but also reported in mollusk species (Cheng et al., 2013; Sun et al., 2018; Yang et al., 2020).

Table 2 Gene annotations of the complete mt genome of Lentigo lentiginosus.

Figure 1 Mitochondrial genome map of the mt genomes of strombids represented by Lentigo lentiginosus.

Table 3 List of AT content, AT skew, and GC skew of Lentigo lentiginosus (Llen), Canarium labiatum (Clab), and Euprotomus aratrum (Eara).

Protein-Coding Genes, Transfer RNA, and Ribosomal RNA Genes

The AT content and AT skew of the PCGs were similar in the mt genomes of L. lentiginosus and C. labiatum, but the GC skew of E. aratrum showed the opposite asymmetry with the other two species (Table 3). All PCGs of the three species started with the conventional initiation codon ATG and ended with the complete stop codons TAA and TAG. In the mt genome of L. lentiginosus, eight PCGs ended with TAA and five with TAG (Table 2), while in C. labiatum, 10 PCGs used TAA, and three ended with TAG (Table S1). In E. aratrum, all PCGs ended with TAA except for nad4l, which used TAG as the stop codon (Table S2). Codon usage of PCGs is shown in Table 4. All three mt genomes had 3,740 codons (excluding the 13 stop codons), among which UUA (Leu) was the most frequently used codon (Table 4 and Figure 2). The least chosen codon CGC (Arg) was observed in C. labiatum and L. lentiginosus, which has also been found in other gastropod mt genomes (Yang et al., 2020), whereas in E. aratrum, the least chosen codon was GCG (Ala). Several codons (UUA, AUU, UUU, GCU, AUA, UCU, GUU, GGA, and UAU) were detected more frequently used than others, indicating a synonymous codon usage bias in strombid mt genomes. These preferred codons, which were found ending in A or U, resulted in a strong A + T bias at the third codon position and contributed to the increase of high A + T content in the whole mitochondrial genome. The synonymous codon usage bias might be caused by mutational bias alone or by both mutation bias and natural selection (Wei et al., 2014).

Table 4 Codon and relative synonymous codon usage (RSCU) of 13 protein-coding genes (PCGs) in the mt genomes of Lentigo lentiginosus, Canarium labiatum, and Euprotomus aratrum.

Figure 2 Relative synonymous codon usage (RSCU) of mitochondrial genomes for Lentigo lentiginosus, Canarium labiatum, and Euprotomus aratrum.

The average AT contents of tRNA genes were 66.7%, 66.5%, and 64.2% in L. lentiginosus, C. labiatum, and E. aratrum, respectively (Table 3). The lengths of the tRNA genes were similar within three mt genomes, ranging from 65 to 73 bp (Tables 2, S1, and S2). All tRNA genes could be folded into typical clover-leaf secondary structures except for the trnS-AGN in three mt genomes because of the missing dihydrouracil (DHU) arms (Figure 3), which was common in metazoan mt genomes (Wolstenholme, 1992).

Figure 3 Inferred secondary structures of 22 transfer RNAs (tRNAs) from Lentigo lentiginosus.

The rrnS genes of L. lentiginosus, C. labiatum, and E. aratrum were 987, 978, and 1,003 bp in length, with AT contents of 66.5%, 67.6%, and 63.9%, respectively. In contrast, rrnL genes were 1,394, 1,382, and 1,373 bp, with AT contents of 70.0%, 70.1%, and 68.4%, respectively. All rRNAs of three mt genomes showed a positive AT skew, but those of L. lentiginosus and C. labiatum were substantially GC skewed compared with E. aratrum, which showed a week GC skew value (Table 3).

Phylogenetic Relationship

The phylogeny of Stromboidea was reconstructed based on the first dataset (5,644 positions in length) using probabilistic methods (Figure 4). The best partition scheme for the amino acid sequences of PCGs was combining genes by subunits, while for the nucleotide sequences of rRNA genes, the best scheme was combining rrnL and rrnS genes. Both ML (−lnL = 45,826.59) and BI (−lnL = 45,699.12 for run1; −lnL = 45,700.08 for run2) arrived at identical topologies (Figure 4).

Figure 4 Phylogenetic relationship of Stromboidea based on the combined amino acid sequences of 13 mitochondrial protein-coding genes (PCGs) plus nucleotide sequences of 2 ribosomal RNA (rRNA) genes: the reconstructed maximum likelihood (ML) phylogram using Charonia lampas and Bufonaria rana as outgroup is shown. The first number at each node is bootstrap proportion (BP) of maximum likelihood (ML) analyses, and the second number is Bayesian posterior probability (PP). Nodal with maximum statistical support (BP = 100; PP = 1) is marked with a solid red circle. BP values of ML under 60 are not shown in the tree.

The reconstructed phylogeny recovered two fully supported clades (namely, Clade I and Clade II) within Stromboidea as identified before (Irwin et al., 2021). Clade I was comprised of Rostellariidae + (Seraphsidae + Strombidae) (Figure 4). The affinity of Clade I has been supported by the morphological similarities in the possession of a “stromboid notch” in the outer lip of the shell and the position of the eye on the end of peduncles and a diminished cephalic tentacle that arises from the middle to the end on that peduncle, inconsistent with the members of Clade II where the eye is located at the base of the cephalic tentacle (Simone, 2005; Maxwell et al., 2019). Based on the morphological synapomorphies, a new crown clade Neostromboidea was constructed to separate families of Clade I from Clade II (Maxwell et al., 2019). Within Clade I, the sister group of Seraphsidae and Strombidae was moderately supported (Figure 4). However, the internal relationships within Clade I could be affected by gene choice and the phylogenetic method. According to Irwin et al. (2021), analyses based on nuclear genes only recovered Seraphsidae as the sister group of Strombidae, consistent with our phylogeny based on amino acid sequences of PCGs plus nucleotide sequences of rRNA genes. Nevertheless, analyses using nucleotide sequences of PCGs recovered the sister relationship between Rostellariidae and Strombidae, while combined sequences of mitochondrial and nuclear genes recovered Strombidae + (Rostellariidae + Seraphsidae) (Irwin et al., 2021). This result indicated that not only gene choice but data type could influence the tree topologies within Stromboidea.

Within Clade II, Xenophoridae was sister to Struthiolariidae + Aporrhaidae (Figure 4). The close relationship between Xenophoridae and Stromboidea has been disclosed in previous studies based on molecular phylogenies (Irwin et al., 2021; Machkour-M’Rabet et al., 2021) and behavioral (Berg, 1974) and morphological (Simone, 2005; Simone, 2011) traits. The topology within Clade II presented here was consistent with that of Irwin et al. (2021) and Machkour-M’Rabet et al. (2021) derived from nucleotide sequences of PCGs but different from that resulting from nuclear genes only (represented as (Xenophoridae + Aporrhaidae) + Struthiolariidae). However, none of these molecular phylogenies was in accordance with the topology derived from more than 100 morphological characteristics (e.g., shell form, eye location, and radular length), which revealed Xenophoridae as sister to Seraphsidae + Strombidae (Simone, 2005).

The final matrix of the second dataset was 13,249 positions in length. According to the BIC, the best partition scheme for the PCGs was the one combining genes by subunits but analyzing each codon position separately (Table S3). For the rRNA genes, the best partition scheme was the one combining together rrnL and rrnS genes. Both ML (−lnL = 80,868.59) and BI (−lnL = 78,968.49 for run 1; −lnL = 78,969.77 for run 2) analyses arrived at identical topologies (Figure 5). The reconstructed phylogeny derived from the second dataset was different from the first dataset within Strombidae. In order to test if this contradiction resulted from data type or outgroup selection, the PCG sequences of the first and second datasets were converted to nucleotide (dataset three; nucleotide sequences of 13 PCGs and two rRNA genes of all available Stromboidea mitogenomes) and amino acid (dataset four; amino acid sequences of 13 PCGs and nucleotide sequences of two rRNA genes of all available Strombidae mitogenomes) levels, respectively. The reconstructed phylogeny based on dataset three revealed the same topologies within Strombidae (Figure S1) as dataset two (Figure 5), the phylogenetic relationships of Strombidae derived from dataset one (Figure 4), and dataset four (Figure S2) showed several inconsistencies, but they still differed from those based on nucleotide datasets (Figures 5, S1). These results suggested that the incongruent topologies generated by the first and second datasets were mainly due to data type. In addition, the outgroup selection could influence the topologies of amino acid datasets. The contradicted phylogenies within Strombidae between different data types have also been observed in Irwin et al. (2021). They were likely caused by low levels of variation in the amino acids at this hierarchical taxonomic level (Uribe et al., 2017). In order to maximize phylogenetic information, the topologies derived from dataset two were considered and discussed in the following.

Figure 5 Phylogenetic relationship of Strombidae based on nucleotide sequences of 13 mitochondrial protein-coding genes (PCGs) and 2 ribosomal RNA (rRNA) genes: the reconstructed maximum likelihood (ML) phylogram using Terebellum terebellum and Varicospira cancellata as outgroup is shown. The first number at each node is bootstrap proportion (BP) of maximum likelihood (ML) analyses, and the second number is Bayesian posterior probability (PP). Nodal with maximum statistical support (BP = 100; PP = 1) is marked with a solid red circle. The branches in purple indicate Eastern Pacific/Atlantic strombids, while those in green represent Indo-Pacific lineages.

This study increases the taxon sampling to more than 30% of the genera of Strombidae (11/30). Within Strombidae, a total of three lineages were recognized (Figure 5). The lineage represented by the single species L. lentiginosus was firstly branching off. According to Abbott (1960), Lentigo was originally defined as a subgenus under Strombus, with five species from both Indo-Pacific and Eastern Pacific/Atlantic regions assigned to this group. This allocation, however, was proved as arbitrary due to the different conchological characteristics between the Indo-Pacific species (L. lentiginosus and Lentigo pipus) and the Eastern Pacific/Atlantic members (Strombus granulatus and Strombus latus) (Kronenberg and Vermeij, 2002). The molecular phylogeny by Latiolais et al. (2006) suggested that the genus Strombus sensu Aboott was non-monophyletic, even though only the Eastern Pacific/Atlantic members S. granulatus and S. latus were included. The current Lentigo was treated as one valid genus to which the Indo-Pacific species L. lentiginosus and L. pipus were assigned (MolluscaBase, 2022c). Represented by the type species L. lentiginosus, the phylogenetic position of Lentigo in the present study was fully supported (Figure 5).

The second lineage was comprised of E. aratrum + (A. gigas and Strombus pugilis) (Figure 5). Machkour-M’Rabet et al. (2021) revealed that Strombidae was grouped by two biogeographically structured clades, which was contradicted by the inclusion of Indo-Pacific E. aratrum within the Eastern Pacific/Atlantic clade formed by A. gigas + S. pugilis (Figure 5). Previous studies also supported that the family Strombidae could not be split into Eastern Pacific/Atlantic versus Indo-Pacific regions (Latiolais et al., 2006; Dekkers, 2008). The close affinity of A. gigas and S. pugilis has been supported by Irwin et al. (2021) and Machkour-M’Rabet et al. (2021), and it corresponded to the Atlantic and West Pacific endemic species with large size (Bandel, 2007). All the Eastern Pacific/Atlantic strombids, according to Latiolais et al. (2006), formed a monophyletic clade whose ancestor might originate from the European Tethys (Maxwell et al., 2020). The position of E. aratrum within this lineage needs to be further confirmed with a complete taxon sampling.

The remaining seven Indo-Pacific species comprised the third lineage (Figure 5). Their phylogeny in the present study is consistent with that of Machkour-M’Rabet et al. (2021) but differs from that of Irwin et al. (2021), according to the relative positions of C. luhuanus. The different topologies were attributed to the different datasets used for phylogenetic reconstructions since Irwin et al. (2021) not only included mitochondrial genes but also added two nuclear genes 18S and 28S. Based on the phylogeny of Latiolais et al. (2006), Dekkers (2008) concluded that Strombidae could be divided into two major clades, with one clade corresponding to medium- to large-sized shells and decorated with course sculpture with knobs usually on the shoulder and another with smooth and mostly small shells. Those shell morphological characteristics did not correspond with the present phylogeny. However, the topology within the third lineage still reflected morphological similarities to some extent. The shell of C. luhuanus, which was firstly branching off, was characterized by a distinct conic shape with a depressed spire and simple outer lip (Abbott, 1960; Simone, 2005; Bandel, 2007). Within the remaining species, both Harpago chiragra and Lambis lambis showed a spider-like shell (Bandel, 2007). The sister group formed by Tridentarius dentatus and C. labiatum was consistent with a small-sized and slender fusiform, while the rest clade grouped by Laevistrombus canarium and Dolomena variabilis possessed large body whorl (Bandel, 2007).

Divergence Times

The estimated divergence times within Strombidae in the present study (Figure 6) were much more recent compared with those of Jiang et al. (2019) and Machkour-M’Rabet et al. (2021). This contradiction could be explained by the employment of older calibration points based on the fossils, which were outside Stromboidea in previous studies. The diversifications within Strombidae were dated from about 35.0 to 4.6 Mya (Figure 6), with most speciation events falling in the periods after 26 Mya. This result supported the hypothesis that the present-day levels of strombid diversity were achieved by the Miocene (Roy, 1996). The divergence between the Eastern Pacific/Atlantic A. gigas + S. pugilis and the Indo-Pacific E. aratrum was dated to 26 Mya (Figure 6), similar to the divergence time (about 24 Mya) between the corresponding geographical groups of marine gastropod nassariids (Caenogastropoda: Nassariidae) (Yang et al., 2021). It was inferred that the ancestor of Eastern Pacific/Atlantic strombids originated from the European Tethys (Maxwell et al., 2020). The estimated origin time of the Eastern Pacific/Atlantic species A. gigas and S. pugilis fell in the late Oligocene. During this period, the onset of deep-water circulation between the Arctic and North Atlantic Oceans (from about 35 to 15 Mya; Davies et al., 2001) may have constituted a strong barrier to marine taxa dispersal between North America and Eurasia. In this regard, the Western Atlantic strombids were separated from the European Tethys group. They might cross into the Eastern Pacific in advance of the closure of the Isthmus of Panama and form the current distribution pattern. During the Oligocene–Miocene transition, there was a rapid expansion of the Antarctica ice sheet and large-scale climate fluctuations attributed to the global cooling condition. Since the ancestor of Strombidae originated in the tropics (Roy, 1996), this global cooling event is likely to have favored species that could be adapted to a subtropical environment. These global climate oscillations might have triggered diversification events leading to the extant Eastern Pacific/Atlantic strombids. The origins of most extant Indo-Pacific species within the third lineage were dated to the period from 16.6 to 4.6 Mya, falling in the epoch when there was a sustained global cooling event that started about 12 Mya and culminated about 5.4 Mya (Herbert et al., 2016). A trend toward cooler conditions in the Indo-Pacific region might act as another diversification event that has favored these extant species to adapt to a subtropical environment. This result therefore suggests another correlation between speciation events and glacial climate change as revealed in other marine organisms (Davis et al., 2016).

Figure 6 Divergences time estimations for the Strombidae using Bayesian relaxed dating methods (BEAST) based on nucleotide sequences of 13 mitochondrial protein-coding genes (PCGs) and 2 ribosomal RNA (rRNA) genes: dates (and credibility intervals) are in millions of years, and horizontal bars represent 95% credibility intervals of relevant nodes. Calibration constraints are indicated with an asterisk on the corresponding nodes.

Conclusion

The complete mitochondrial genomes of L. lentiginosus, E. aratrum, and C. labiatum were similar to those of other strombids published before in genome size, gene order, and nucleotide composition. The reconstructed mitogenomic phylogeny recovered two monophyletic clades corresponding to different morphological synapomorphies and supported the inclusion of Xenophoridae within Stromboidea. The internal phylogenetic relationships of Strombidae also reflected several morphological similarities to some extent. Furthermore, the present study calls for an increased taxon sampling to resolve the controversial phylogenetic positions within Strombidae. The divergence time estimation suggests that the strombid diversifications might be caused by paleoclimate change.

Data Availability Statement

The datasets presented in this study can be found in online repositories. The names of the repository/repositories and accession number(s) can be found below: https://www.ncbi.nlm.nih.gov/genbank/, ON310809 https://www.ncbi.nlm.nih.gov/genbank/, ON310803 https://www.ncbi.nlm.nih.gov/genbank/, ON310804.

Author Contributions

FL: data curation, formal analysis, investigation, and writing—original draft preparation. JZ: data curation and investigation. QM: data curation and investigation. ZG: software and supervision. AW: supervision and funding acquisition. YY: data curation, formal analysis, methodology, funding acquisition, supervision, and writing—review and editing. CL: funding acquisition, supervision, and writing—review and editing. All authors listed have made a substantial, direct, and intellectual contribution to the work and approved it for publication.

Funding

This work was supported by the National Key Research and Development Program of China (2019YFD0901301), the Key Research and Development Project of Hainan Province (ZDYF2021SHFZ269), and the Starting Research Fund from the Hainan University is (RZ2100003081).

Conflict of Interest

The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.

Publisher’s Note

All claims expressed in this article are solely those of the authors and do not necessarily represent those of their affiliated organizations, or those of the publisher, the editors and the reviewers. Any product that may be evaluated in this article, or claim that may be made by its manufacturer, is not guaranteed or endorsed by the publisher.

Supplementary Material

The Supplementary Material for this article can be found online at: https://www.frontiersin.org/articles/10.3389/fmars.2022.930910/full#supplementary-material

Supplementary Table 1 | Gene annotations of the complete mt genome of Canarium labiatum.

Supplementary Table 2 | Gene annotations of the complete mt genome of Euprotomus aratrum.

Supplementary Table 3 | Best fit partitions and substitution models.

References