Insights into the Deep Phylogeny and Novel Divergence Time Estimation of Patellogastropoda from Complete Mitogenomes

To further understand the origin and evolution of Patellogastropoda, we determined the mitochondrial genome sequence of Cellana toreuma, and compared its mitogenome characteristics with the other four limpets of Nacellidae. The ratio of Ka and Ks indicated that these Nacellidae species were suffering a purifying selection, with exception of the atp6 gene. The gene sequence is basically consistent among families, while there are great differences among Lottidae species. According to the mitogenome sequences of selected gastropod species, we reconstructed a new phylogenetic tree with two methods. The data complement the mitogenome database of limpets and is a favorable research tool for the phylogenetic analysis of Gastropoda. It is found that there is a long-branch attraction (LBA) artefact in the family Lottiidae of Patellogastropoda. Therefore, the Patellogastropoda was separated by Heterobranchia, and Lottiidae is located at the root of the whole phylogenetic tree. Furthermore, we constructed the divergence time tree according to the Bayesian method and discussed the internal historical dynamics, and divergence differences among the main lineages of 12 Patellogastropoda under an uncorrelated relaxed molecular clock. In turn, we made a more comprehensive discussion on the divergence time of limpets at the molecular level.


Introduction
Mitochondria are circular double-membrane semiautonomous organelle, which exist in the cells of most eukaryotic species. It has an independent and complete mitochondrial genome. They originate from an endosymbiotic α-proteobacterium and usually provide chemical energy sources through oxidative phosphorylation [1][2][3][4]. Because mutations affecting mitochondrial function are related to aging and disease, it also has certain biomedical significance [5,6]. The advantages of mitochondria are that their evolution rate is faster than that of nuclear genes in most species. Each cell has multiple copies of the mitochondrial genome and higher A-T content [7]. Moreover, it also has the characteristics of conservation gene function [8]. Mitogenome recombination is a common process in protists and plants [9,10]. Single genes may affect the progress of species relationships because of their different evolutionary rates [11]. Thus, the complete mitogenome is considered significant in population genetics and phylogeny, as well as an important tool for an in-depth understanding of gastropod phylogenetics.
Patellogastropoda, as an archaic mollusk, has caused concern in the scientific community. Due to their important ecological status and biodiversity, it is often researched in

Sample Collection, Identification, and DNA Extraction
Cellana toreuma wild specimens were collected from Qingdao of Shandong Peninsula in the Yellow Sea (October 2019; E 120 • 45, N 36 • 07). The specimens were deposited in absolute ethyl alcohol. The specimens were preliminarily identified through the published taxonomic book [46], and we consulted morphology experts from the marine biology museum of Zhejiang Ocean University. Fresh samples were immediately placed in absolute ethyl alcohol to ensure their quality. Using the rapid salting-out method, we extracted the genomic DNA from the adductor muscle [47]. The quality was determined by 1% agarose gel electrophoresis and stored in −20 • C for sequencing. We selected the best quality DNA from the six samples for the next-generation sequencing. All animal experiments were conducted under the guidance approved by the Animal Research and Ethics Committee of Zhejiang Ocean University.

Mitogenomes Sequencing, Assembly, and Annotation
Mitogenome sequencing of C. toreuma by the Illumina HiSeq X Ten platform was used to conduct high-throughput sequencing; this work was carried out by Origingene Bio-pharm Technology Co., Ltd. (Shanghai, China). The preliminary results showed that a sequencing library set with an average insert size of 400 bp was generated, and each library had about 10 Gb of the raw data. After that, it was necessary to delete contaminated reads and low-quality sequence fragments. The de novo assembled separate clean readings of the sequence via the NOVOPlasty software (https://github.com/ndierckx/NOVOPlasty (accessed on 26 May 2021)) [48].
The mitochondrial genome of C. toreuma was annotated and analyzed based on invertebrate genetic code by the MITOS web server (http://mitos2.bioinf.uni-leipzig.de/index.py (accessed on 30 May 2021)) [49]. We also referred to the uploaded mitogenome sequence of other Nacellidae species to ensure the accuracy of start and stop codons and gene sequences of the species in our study. The circular mitogenome visualization of C. toreuma was completed through the common CGView server (http://stothard.afns.ualberta.ca/ cgview_server/index.html (accessed on 30 May 2021)) [50].

Sequence Analyses of Mitogenomes
The nucleotide composition of the whole mitogenome, PCGs, rRNA, tRNA genes, and A-T content were analyzed by MEGA 7.0 [51]. Meanwhile, it was also determined necessary to study the codon usage and the relative synonymous codon usage (RSCU) of PCGs. Then, the base skew values were calculated using the formulas at A-T skew = (A − T)/(A + T) and G-C skew = (G − C)/(G + C) [52]. In addition, we selected DnaSP6.0 [53] to analyze the nonsynonymous (Ka) and synonymous (Ks) substitutions rates of mitogenomes in Nacellidae species to study their evolutionary adaptation. Of which other species in Nacellidae were downloaded from the GenBank database of NCBI (National Center for Biotechnology Information, https://www.ncbi.nlm.nih.gov/ (accessed on 17 March 2022)).

Phylogenetic Inference
To determine the phylogenetic position of the Patellogastropoda species in gastropods, phylogenetic analyses were performed based on the 13 protein-coding genes (PCGs) of the mitogenomes. A total of 87 sequences were downloaded from GenBank (https://www.ncbi.nlm.nih.gov/ (accessed on 17 March 2022)). Furthermore, two bivalves Donax variegatus and Donax trunculus were classified as outgroup [54] (Table 1). The software DAMBE 5.3.19 [55] was used to adjust the nucleotide sequence of each PCGs, and the substitution saturation was calculated via the GTR substitution model. The sequences were aligned using ClustalW of MEGA 7.0 [51] with the default parameters. Phylogenetic relationships were reconstructed for the maximum likelihood (ML) and Bayesian inference (BI) analyses with IQ-TREE [56] and MrBayes v3.2 [57]. The GTR + F + I + G4 model was chosen according to BIC and was the best fit selected in the ML methods analyses; we used ModelFinder to determine the best model data [58]. We reconstructed the consensus tree and used 1000 bootstrap replicates in ultrafast likelihood bootstrap. The BI analysis selected the best substitution GTR + I + G model under the AIC by MrModeltest 2.3 [59]. The first burn-in 25% of the trees were discarded, and two Markov chain Monte Carlo (MCMC) of simultaneous were operated for 2,000,000 generations. Sampling occurred every 1000 ultrafast bootstrap replicates to determine the branch support of the dataset. The phylogenetic tree was displayed using the online tool iTOL (Interactive Tree of Life, https://itol.embl.de/ (accessed on 20 March 2022)) and annotated with various datasets [60].

Divergence Time Estimation
The divergence time of the subclass Patellogastropoda species was estimated only at the nucleotide level, 13 PCGs used the Bayesian framework, and we used the uncorrelated and lognormal relaxed molecular clock model in BEAST v1.8.4 [61]. For the tree prior, we used a Yule process of speciation. Furthermore, we used two calibration points as the priors of the divergence times for calibration. The uniform distribution of the estimated divergence times was drawn by Priors [62] for fossil ages, and the 53 Mya point calibration was set as the root rate of Cellana based on the fossil of Cellana tramoserica (14-93 Mya). The 5.6 Mya point calibration was set as the root rate of Nacella based on the fossil of Nacella clypeater [62]. The final Markov chain samples every 1000 generations, discards 10% of the burn-in samples, and runs twice for 100 million generations through the TreeAnnotator v1.8.4 software in the BEAST software package. After that, we use Tracer v. 1.6 [63] to check the convergence of the chain, to ensure that the parameters of the effective sample sizes (ESSs) were greater than 200. We used the software FigTree v1.4.3 to visualize the divergence time tree [64].

General Features of Entire Mitogenome
The whole mitogenome sequence of C. toreuma was sequenced with a length of 16,260 bp (GenBank accessions: MZ329338) which is consistent with previously reported species of the four Nacellidae families, approximately 16,153 to 16,767 ( Table 1). The circular molecules are similar to other gastropods, which contain a highly variable control region and typically 37 genes including 2 ribosomal RNA genes, 13 protein-coding genes (PCGs), and 22 transfer RNA genes. Among them, a total of seventeen genes on the forward strand, including seven PCGs (cox1-3, atp8, atp6, nad3, and nad2), and ten tRNA genes (trnD, trnT, trnG, trnE, trnR, trnN, trnA, trnK, trnI, and trnS1). The other genes are encoded on the reverse strand ( Table 2). The control region was located between the trnC and trnG gene, similar to other previously reported Nacellidae species ( Figure 1) [65,66]. The genome structure of C. toreuma was identical to other Nacellidae mitogenomes, without gene rearrangement in this family. However, there was a big difference between their gene order to other families of Patellogastropoda species, whose rearrangement always brought up concern within the scientific community.  The nucleotide compositions of complete C. toreuma mitogenomes were A 28.9%, T 39.5%, G 19.9%, and C 11.7% (Table 3). Moreover, the nucleotide compositions of the mitogenomes from the other four species in the family Nacellidae, Patellogastropoda were downloaded and organized, and we compared the base compositions of these other Nacellidae species. In general, the A content of the five mitogenomes were from 26.5 to 28.9%, T 38.0 to 39.5%, G 19.9 to 22.7%, and C 11.7 to 13.9%, these base contacts only had a slight distinction. Results show that the A and T content of C. toreuma exhibited higher values than other species in the same family, with a range from 64.6% (Cellana nigrolineata) to 68.4%. The nucleotide compositions were all skewed away from C in favor of G, the G-C skews were from 0.177 (N. clypeater) to 0.283 (C. nigrolineata), and from A in favor of T, the A-T skews were negative from −0.180 (Nacella concinna) to −0.155 (C. toreuma), indicating the occurrence of more Ts than As. Table 3. Nucleotide compositions of the mitogenomes from five species in family Nacellidae of the Patellogastropoda. For the tRNA genes of C. toreuma, the length ranged from 1497-1560 bp, which had an A and T content of 69.6%, which is the highest compared to other Nacellidae species (Table 3). Due to the tRNAs of this limpet having similar values of A and T base, the A-T skew was 0, while the A-T skew of the other four species were negative. However, all the G-C skews were slightly positive from 0.103 (C. toreuma) to 0.125 (N. clypeater). For the rRNA genes, with lengths ranging from 2143 to 2222 bp, the A-T skews were positive from 0.169 (C. toreuma) to 0.237 (N. concinna), indicating a strong skew away from A. Furthermore, almost all G-C skews of these species were negative, except the new limpet C. toreuma.

Lengh (bp) A (%) T (%) G (%) C (%) A+T (%) AT-Skew GC-Skew
For the PCGs, the length ranging from 11,046 to 11,286 bp, each species in this family exhibited a negative A-T skew and a positive G-C skew, with values from −0.221 to −0.201. The first codon position of PCGs was observed to be a negative A-T skew, with the A-T content reaching about 65%, and the most value being 67% in C. toreuma. The second and third codon positions of PCGs were similar, while their G-C skews had the opposite results. The mitogenomes of Nacellidae are rich in A-T, which is similar to other invertebrates. The C. toreuma mitogenomes conventional started with the initiation codon ATG or ATT and stopped with TAA or TAG.
Generally, the mitogenome of metazoans is quite compact, whereas a total of 1641 bp in 30 intergenic spacers were found in C. toreuma mitogenome, ranging from 1 to 643 bp in length. Additionally, the control region was found between trnC and trnG (Table 2). Simultaneously, three overlapping sites (totally 87 bp) are observed ranging from 2 to 47 bp, which is commonly identified in other Nacellidae.

Mitochondrial Gene Codon Usage
The amino acids of five Nacellidae species, Leu1, Phe, and Val, were most frequently utilized ranging from 7.60 to 15.11% in Phe of N. clypeater (Figure 2). The rare amino acid was concentrated in His, Gln, and Arg, most of them less than 2 with the least being Gln of C. toreuma (only 1.38%), which is similar to other Patellogastropoda as previously reported [67,68].

tRNA, rRNA, PCGs Genes, and Control Region
For the tRNA genes of C. toreuma, the length ranged from 1497-1560 bp, wh an A and T content of 69.6%, which is the highest compared to other Nacellidae (Table 3). Due to the tRNAs of this limpet having similar values of A and T base, skew was 0, while the A-T skew of the other four species were negative. However G-C skews were slightly positive from 0.103 (C. toreuma) to 0.125 (N. clypeater). rRNA genes, with lengths ranging from 2143 to 2222 bp, the A-T skews were positi 0.169 (C. toreuma) to 0.237 (N. concinna), indicating a strong skew away from A. F more, almost all G-C skews of these species were negative, except the new limp reuma.
For the PCGs, the length ranging from 11,046 to 11,286 bp, each species in thi exhibited a negative A-T skew and a positive G-C skew, with values from −0.221 to The first codon position of PCGs was observed to be a negative A-T skew, with content reaching about 65%, and the most value being 67% in C. toreuma. The seco third codon positions of PCGs were similar, while their G-C skews had the oppo sults. The mitogenomes of Nacellidae are rich in A-T, which is similar to other brates. The C. toreuma mitogenomes conventional started with the initiation codo or ATT and stopped with TAA or TAG.
Generally, the mitogenome of metazoans is quite compact, whereas a total of in 30 intergenic spacers were found in C. toreuma mitogenome, ranging from 1 to in length. Additionally, the control region was found between trnC and trnG (T Simultaneously, three overlapping sites (totally 87 bp) are observed ranging from bp, which is commonly identified in other Nacellidae.

Mitochondrial Gene Codon Usage
The amino acids of five Nacellidae species, Leu1, Phe, and Val, were most fre utilized ranging from 7.60 to 15.11% in Phe of N. clypeater ( Figure 2). The rare ami was concentrated in His, Gln, and Arg, most of them less than 2 with the least be of C. toreuma (only 1.38%), which is similar to other Patellogastropoda as previo ported [67,68].  The relative synonymous codon usage (RSCU) for 13 PCGs of these species was measured to understand the genetic codon bias of their sequenced mitogenomes. The results showed that the synonymous codon preferences are conserved among the five species, as may ascribe to their close relationships belonging to the same family. These preferences were also recognized in some other Patellogastropoda. The four most used codons for the five species sequenced are consistently UUA (Leu2), UCU (Ser2), GCU (Ala), and CCU (Pro) for five Nacellidae species (Figure 3), and the most frequent codons of them were UUA at 2.4% (Leu2) in N. magellanica. The least used codon was CUC (Leu1) in C. toreuma, which was 0.1%. cies, as may ascribe to their close relationships belonging to the same family. These pr erences were also recognized in some other Patellogastropoda. The four most used codo for the five species sequenced are consistently UUA (Leu2), UCU (Ser2), GCU (Ala), a CCU (Pro) for five Nacellidae species (Figure 3), and the most frequent codons of th were UUA at 2.4% (Leu2) in N. magellanica. The least used codon was CUC (Leu1) in toreuma, which was 0.1%.

Selective Pressure Analysis
The selection pressure analysis also used these five species (Figure 4) to measure the ratio of non-synonymous and synonymous substitutions (Ka/Ks). We aimed to investigate the evolutionary and selective pressure relation. The results showed the average Ka/Ks ranging from cox3 (0.124) to atp6 (1.106). The ratio for most PCGs was below one, indicating that the mutations were swapped by synonymous substitutions; 13 PCGs of these mitogenomes were evolving under purifying selection. The remaining atp6 gene reached 1.106, which may be due to the influence of positive selection during evolution. Among these species, the cox3 gene had little change in amino acids and the lowest Ka/Ks ratio. It is widely used as a potential molecular marker for species identification, genetic diversity, and phylogenetic analysis [69,70].

Selective Pressure Analysis
The selection pressure analysis also used these five species (Figure 4) to measure the ratio of non-synonymous and synonymous substitutions (Ka/Ks). We aimed to investigate the evolutionary and selective pressure relation. The results showed the average Ka/Ks ranging from cox3 (0.124) to atp6 (1.106). The ratio for most PCGs was below one, indicating that the mutations were swapped by synonymous substitutions; 13 PCGs of these mitogenomes were evolving under purifying selection. The remaining atp6 gene reached 1.106, which may be due to the influence of positive selection during evolution. Among these species, the cox3 gene had little change in amino acids and the lowest Ka/Ks ratio. It is widely used as a potential molecular marker for species identification, genetic diversity, and phylogenetic analysis [69,70]. The substitution saturation index of the combined dataset for 88 Patellogastropoda mitogenomes for 13 PCGs (Iss = 0.823) was significantly lower than the critical values (Iss.cSym = 0.860 or Iss.cAsym = 0.847, p = 0.000). Therefore, substitution of combined sequencing was unsaturated and suitable for phylogenetic analysis.

Gene Arrangement
According to the hypothetical gene order of ancestral gastropods, we compared the PCG gene arrangement of 12 species in four families of the subclass Patellogastropoda ( Figure 5). The results showed that the gene order of the family Nacellidae, where C. toreuma was located, was the same as that of the family Acmaeidae, and was consistent with the ancestral gene order. In Patellidae, we found that only the nad3-nad2 fragment moved, and the position was transferred from one fragment of cox1 to the other. The nad6-nad1-rrnL-rrnS-cox3 fragment is still completely preserved. In addition, the atp8-atp6-nad5-nad4-nad4l-cytb gene fragment was completely reversed. It is worth noting that the gene rearrangement rate of the family Lottiidae is very high. Among them, L. goshimai retains the nad3-nad2 gene fragment, while the nad4 and nad4l genes in L. goshima are reversed; this situation also occurs in Lottia digitalis. Moreover, the short gene fragment of rrnL-rrnS was retained in L. digitalis, and two rrnL fragments were reversed in L. goshimai. In particular, The substitution saturation index of the combined dataset for 88 Patellogastropoda mitogenomes for 13 PCGs (Iss = 0.823) was significantly lower than the critical values (Iss.cSym = 0.860 or Iss.cAsym = 0.847, p = 0.000). Therefore, substitution of combined sequencing was unsaturated and suitable for phylogenetic analysis.

Gene Arrangement
According to the hypothetical gene order of ancestral gastropods, we compared the PCG gene arrangement of 12 species in four families of the subclass Patellogastropoda ( Figure 5). The results showed that the gene order of the family Nacellidae, where C. toreuma was located, was the same as that of the family Acmaeidae, and was consistent with the ancestral gene order. In Patellidae, we found that only the nad3-nad2 fragment moved, and the position was transferred from one fragment of cox1 to the other. The nad6-nad1-rrnL-rrnS-cox3 fragment is still completely preserved. In addition, the atp8-atp6-nad5-nad4-nad4l-cytb gene fragment was completely reversed. It is worth noting that the gene rearrangement rate of the family Lottiidae is very high. Among them, L. goshimai retains the nad3-nad2 gene fragment, while the nad4 and nad4l genes in L. goshima are reversed; this situation also occurs in Lottia digitalis. Moreover, the short gene fragment of rrnL-rrnS was retained in L. digitalis, and two rrnL fragments were reversed in L. goshimai. In particular, N. fuscoviridis hardly retains the gene fragments of its ancestors; however, we can find that it has the same fragment cox1-cox3 as L. digitalis in the same family. It also has the inversion of atp6 and cox2 genes with L. goshima, which may be the reservation of unique gene fragments produced in the evolutionary process of this family. In general, the rearrangement of Lottiidae remains the focus of our research. This irregular rearrangement may be responsible for the separation of this family from the other three families in the subclass Patellogastropoda.
Genes 2022, 13, x FOR PEER REVIEW 12 of N. fuscoviridis hardly retains the gene fragments of its ancestors; however, we can find th it has the same fragment cox1-cox3 as L. digitalis in the same family. It also has the inversio of atp6 and cox2 genes with L. goshima, which may be the reservation of unique gene fra ments produced in the evolutionary process of this family. In general, the rearrangeme of Lottiidae remains the focus of our research. This irregular rearrangement may be r sponsible for the separation of this family from the other three families in the subcla Patellogastropoda.

Phylogenetic Relationship
After a long period of evolution, the species of Patellogastropoda are mainly divide into two superfamilies, Lottioidea (Gray, 1840) and Patelloidea (Rafinesque, 1815). W concatenated the alignment of 13 common PCGs from C. toreuma, as combined in 86 sp cies that represent 26 families of these 2 superfamilies in Gastropoda, (i.e., Patellogastro oda, Heterobranchia, Caenogastropoda, Neritimorpha, Neomphaliones, and Vetigastro oda). The bivalves D. variegatus and D. trunculus were set as outgroups ( Figure 6). For M and BI trees, the 88 species could be divided into seven major clades. Of which mo branches exhibited high confidence of being coincidently classified with different clade (i.e., BI: 1 posterior probability and 100% bootstraps value). Strikingly, we found that t Patellogastropoda species were divided into two branches, located on both sides of He erobranchia. The family Lottiidae of Patellogastropoda was located at the foundation position of the integral phylogenomic tree, while other families such as Patellidae, Acma idae, and Nacellidae were located between Heterobranchia and Caenogastropoda. W speculated that this is probably a result of the long-branch attraction (LBA) artefact. Th phenomenon was also discovered in the study of two limpets N. fuscoviridis and L. g

Phylogenetic Relationship
After a long period of evolution, the species of Patellogastropoda are mainly divided into two superfamilies, Lottioidea (Gray, 1840) and Patelloidea (Rafinesque, 1815). We concatenated the alignment of 13 common PCGs from C. toreuma, as combined in 86 species that represent 26 families of these 2 superfamilies in Gastropoda, (i.e., Patellogastropoda, Heterobranchia, Caenogastropoda, Neritimorpha, Neomphaliones, and Vetigastropoda). The bivalves D. variegatus and D. trunculus were set as outgroups ( Figure 6). For ML and BI trees, the 88 species could be divided into seven major clades. Of which most branches exhibited high confidence of being coincidently classified with different clades, (i.e., BI: 1 posterior probability and 100% bootstraps value). Strikingly, we found that the Patellogastropoda species were divided into two branches, located on both sides of Heterobranchia. The family Lottiidae of Patellogastropoda was located at the foundational position of the integral phylogenomic tree, while other families such as Patellidae, Acmaeidae, and Nacellidae were located between Heterobranchia and Caenogastropoda. We speculated that this is probably a result of the long-branch attraction (LBA) artefact. This phenomenon was also discovered in the study of two limpets N. fuscoviridis and L. goshimai [68]. The LBA artifact was defined as taxa with long branches that evolve rapidly in phylogenetic analysis, regardless of their real evolutionary location, which is highly misleading because the inferred false relationship is often highly statistically supported [71]. Both BI and ML analyses are generally unable to avoid this artifact due to the limitations of their model, and there is no good solution at present, which is also a place that needs to be improved for the future. Some scholars have studied and discussed the LBA artifact in previous reports, but no studies can counteract phylogenetics. We still obtain such a branch under our repeated verification in both methods. After investigating, we speculated that the LBA artifact often produces misleading results usually resulting from uneven or adequate sampling, and the mitogenome data of Patellogastropoda is insufficient and cannot be solved temporarily. This makes us more interested in the evolutionary research of the subclass Patellogastropoda, and it is imperative to continue to find new data on Patellogastropoda. After that, these clades combined at the subclass Caenogastropoda. Phylogenetic relationships of the subclasses are recovering as ((((((Vetigastropoda + Neomphaliones) + Neritimorpha) + Caenogastroopoda) + Patellogastropoda) + Heterobranchia) + Patellogastropoda). Their overall branching is basically consistent with our previous evolutionary analysis. Arquez et al. [72] used mitogenome to analyze the evolution of five subclasses of Gastropoda. The results showed that Heterobranchia and Patellogastropoda were sister groups to each other, and this branch was located at the root of the entire evolutionary tree. Caenogastropoda, Vetigastropoda, and Vetigastropoda had the same branch as our study [72]. Subsequently, Osca et al. [73] also used the mitogenome to study the phylogeny of Gastropoda species. However, their result was slightly different, Caenogastropoda and Neritimorpha were sister groups, followed by the subclass Vetigastropoda. In addition, Uribe et al. [74] obtained a similar result to Osca et al. [73]; nevertheless, they added the species of the subclass Neomphalina, which are situated between Vetigastropoda and the sister group Heterobranchia and Patellogastropoda. Sun et al. [75] reconstructed the tree of several species of gastropods using mitogenome, species from Cocculiniformia and Neomphalina have been updated and divided into the Neomphalines subclass. These results were completely consistent with the evolutionary branches of several subclasses in our study [75]. In addition, branches from the present study were similar to the results of Feng et al. [76], except for the minimal change of Caenogastropoda and Neritimorpha. The simple phylogenetic analysis of the deep-sea limpet by Li et al. [77] was the same as the results of the present study, which confirms the reliability of our predicted branches of the phylogenetic tree.
Patellogastropoda has four families that make up the whole mitogenome, Nacellidae and Acmaeidae are sister groups, followed by Patellidae. The Lottiidae, as an independent branch, is located at the outermost of this subclass. Furthermore, our research species C. toreuma and C. nigrolineata are undoubtedly sister species, both belonging to the genus Cellana. They were divided into one branch with three species of the genus Nacella in the same Nacellidae family. Phylogenetic trees of the Nacellidae families stand for (C. toreuma + C. nigrolineata) + (N. concinna + (Nacella magellanica + N. clypeater). These evolutionary branches conform to the conventional evolutionary characteristics, complying with the results of existing studies.
from Cocculiniformia and Neomphalina have been updated and divided into the Neomphalines subclass. These results were completely consistent with the evolutionary branches of several subclasses in our study [75]. In addition, branches from the present study were similar to the results of Feng et al. [76], except for the minimal change of Caenogastropoda and Neritimorpha. The simple phylogenetic analysis of the deep-sea limpet by Li et al. [77] was the same as the results of the present study, which confirms the reliability of our predicted branches of the phylogenetic tree. Figure 6. Phylogenetic tree inferred using maximum likelihood (ML) and Bayesian inference (BI) methods based on concatenated sequences of 13 PCGs from 88 gastropod mitogenomes. The sequences of two bivalves were chosen as the outgroups. The purple dots indicate C. toreuma sequenced in this study. The number at each node is the bootstrap probability.

Divergence Times
We found that only seven available species of complete Patellogastropoda mitogenomes were recorded in the Timetree (http://www.timetree.org/ (accessed on 21 March 2022)) with only three known divergent time points. These were concentrated in Nacella and Cellana during the Cenozoic Era (65 million years ago, Mya) [78]. However, the divergence time information of Bathyacmaea nipponica, Tectura paleacea, and Patella aspera were unknown, which is a very unusual situation in Gastropoda species. Our study was designed to estimate more information on the evolutionary time and understand the historical evolution and dynamics of the limpets.
Our study demonstrated that N. magellanica and N. clypeater were differentiated around 9.74 Mya, and N. concinna differentiated at 15.21 Mya (Figure 7). They are slightly earlier than the results in Timetree, but all the current differentiation occurred in the Neogene era. Nacella and Cellana were differentiated around 50.87 Mya, which is consistent with previous estimates of the divergence time. The analyses supported the supposition that the major Patellogastropoda lineages originated in the early Cretaceous and diversified in the middle and later Cretaceous; these have never been reported previously. In addition, the N. fuscoviridis family diversified to the genus Lottia in the later Cretaceous. The genera Lottia, Bathyacmaea, Patella, and Cellana were differentiated into Cenozoic Paleogene. All Nacellidae species differentiation is concentrated in the Cenozoic period. The geographical isolation in this period provided the environmental conditions for the Nacellidae differentiation, and marine sediments provided food sources for the growth of Patellogastropoda.

Conclusions
To carry out a good systematic study of limpets, a new potential species was found, and the mitogenomes of two species in the family Lottiidae were published. This provided more insights into mitogenomes in Patellogastropoda and supplemented the molecular database in such a classification group. We obtained the mitogenome sequences of C. toreuma by high-throughput sequencing with the length of 16,260 bp, similar to other limpets. Each mitogenome has the same composition and similar results of nucleotide composition for five Nacellidae. The gene order was generally uniform within families, except for the family Lottidae. Most PCGs were initiated with the ATG codon and terminated with TAA codon. For the analysis of selective pressure, we found that most of the 13 PCGs of Nacellidae were below 1, especially the cox3 gene which exhibited the lowest value that demonstrated high conservation. This indicates that PCGs were subject to purifying selection in the family, while the atp6 gene shows a high value, indicating that this gene may have been mutated in the process of evolution. The phylogenetic tree provided a further complement to the scientific classification of Patellogastropoda species. It is found that there is an LBA artifact in the family Lottiidae in Patellogastropoda, which deceived phylogenetic methods caused by the outgroup as distant or that the taxon sampling is poor. To solve this problem, we use two methods to construct evolutionary trees, but the Patellogastropoda is also divided into two branches on both sides of Heterobranchia. This study could provide the basic information for genetic characteristics, phylogenetic position, and evolution for these limpets, and provide a basis for resource management and selective breeding in aquaculture. This Cellana species was differentiated in the late Paleogene, and their evolution may be related to the geological events that changed their living environment.

Conclusions
To carry out a good systematic study of limpets, a new potential species was found, and the mitogenomes of two species in the family Lottiidae were published. This provided more insights into mitogenomes in Patellogastropoda and supplemented the molecular database in such a classification group. We obtained the mitogenome sequences of C. toreuma by high-throughput sequencing with the length of 16,260 bp, similar to other limpets. Each mitogenome has the same composition and similar results of nucleotide composition for five Nacellidae. The gene order was generally uniform within families, except for the family Lottidae. Most PCGs were initiated with the ATG codon and terminated with TAA codon. For the analysis of selective pressure, we found that most of the 13 PCGs of Nacellidae were below 1, especially the cox3 gene which exhibited the lowest value that demonstrated high conservation. This indicates that PCGs were subject to purifying selection in the family, while the atp6 gene shows a high value, indicating that this gene may have been mutated in the process of evolution. The phylogenetic tree provided a further complement to the scientific classification of Patellogastropoda species. It is found that there is an LBA artifact in the family Lottiidae in Patellogastropoda, which deceived phylogenetic methods caused by the outgroup as distant or that the taxon sampling is poor. To solve this problem, we use two methods to construct evolutionary trees, but the Patellogastropoda is also divided into two branches on both sides of Heterobranchia. This study could provide the basic information for genetic characteristics, phylogenetic position, and evolution for these limpets, and provide a basis for resource management and selective breeding in aquaculture. This Cellana species was differentiated in the late Paleogene, and their evolution may be related to the geological events that changed their living environment.
Author Contributions: J.F. and J.M. analyzed the data, wrote the paper, and prepared the figures and tables; Y.Y. conceived and designed the experiments, reviewed drafts of the paper; J.L. contributed analysis tools, reviewed drafts of the paper; K.X. collected field material and processed the samples; B.G. and X.Y. supervised and directed the work. All authors reviewed the manuscript. All authors have read and agreed to the published version of the manuscript.