The Complete Mitochondrial Genome Sequences of the Philomycus bilineatus (Stylommatophora: Philomycidae) and Phylogenetic Analysis

The mitochondrial genome (mitogenome) can provide information for phylogenetic analyses and evolutionary biology. We first sequenced, annotated, and characterized the mitogenome of Philomycus bilineatus in this study. The complete mitogenome was 14,347 bp in length, containing 13 protein-coding genes (PCGs), 23 transfer RNA genes, two ribosomal RNA genes, and two non-coding regions (A + T-rich region). There were 15 overlap locations and 18 intergenic spacer regions found throughout the mitogenome of P. bilineatus. The A + T content in the mitogenome was 72.11%. All PCGs used a standard ATN as a start codon, with the exception of cytochrome c oxidase 1 (cox1) and ATP synthase F0 subunit 8 (atp8) with TTG and GTG. Additionally, TAA or TAG was identified as the typical stop codon. All transfer RNA (tRNA) genes had a typical clover-leaf structure, except for trnS1 (AGC), trnS2 (TCA), and trnK (TTT). A phylogenetic analysis with another 37 species of gastropods was performed using Bayesian inference, based on the amino acid sequences of 13 mitochondrial PCGs. The results indicated that P. bilineatus shares a close ancestry with Meghimatium bilineatum. It seems more appropriate to reclassify it as Arionoidea rather than Limacoidea, as previously thought. Our research may provide a new meaningful insight into the evolution of P. bilineatus.


Introduction
Usually, the mitochondrial (mt) genome is a closed circular structure with molecules ranging from approximately from 14,000 to 18,000 bp in length. The outer ring is a heavy chain, while the inner ring is a light chain, and most genes are transcribed by the heavy chain. Mitochondria are semi-autonomous organelles that can independently encode certain proteins. Since mtDNA is expressed as matrilineal inheritance, the structure and proportion of mitochondrial DNA can better reflect the genetic characteristics of populations [1], although it does not reflect the genetic information of the male population [2,3]. It generally contains two ribosomal RNA genes (12S rRNA and 16S rRNA), 22 transfer RNA genes (tRNAs), 13 protein-coding genes (PCGs), and one typical non-coding control region with regulatory elements necessary for transcription and replication [4]. The mtDNA control region contains the origin of replication of one strand, and the origin of transcription

Specimen Collection and DNA Extraction
P. bilineatus (which is neither an endangered animal, nor a protected species, according to the International Union for Conservation of Nature [IUCN] Red List of Threatened Species) was collected from farmland vegetables in Shanghai, China, in September 2017. The experimental procedures involved in our research were conducted in accordance with international guidelines for the care and treatment of laboratory animals. The samples were placed in 100% ethanol during collection and stored at −80 • C until DNA was extracted (No. ASTM-Mo-P1241). Total DNA was extracted from the partial body tissue using the TIANamp Marine Animal DNA Kit (Tiangen Biotech, Beijing, China). The quality of the separated DNA was detected by 1.5% electrophoresis and the DNA was stored at −20 • C until the complete mitogenomes had been amplified by PCR [14].

Sequence Amplification, Assembly, and Gene Annotation
For the amplification of the P. bilineatus mitogenome, primer (Table S1) sets were designed based on the mitogenomic sequences obtained from other Limacidae species [15][16][17] (GenBank: JN619347, KJ744205, NC_035429). All PCR amplifications were performed with LA Taq DNA polymerase, using 2 × Premix KA Taq (Takara, Beijing, China) in a 50 µL reaction volume, including 20 µL sterilized distilled water, 25 µL 2 × TaqMix buffer, 2 µL DNA, and 1.5 µL of each primer (10 nM). The fragments were amplified with an initial denaturation at 94 • C for 3 min, followed by 35 cycles of denaturation at 94 • C for 30 s, annealing at 50-60 • C for 30 s, then elongating at 72 • C for 1-3 min, with a final elongation at 72 • C for 10 min after the last cycle. All PCR products were sequenced at Shanghai Sangon Biotech. The obtained sequences had 100% coverage of the PCR products.

Sequence Analysis
After the sequence had been manually checked, the mitogenome sequence was performed using the BLAST program from NCBI (https://blast.ncbi.nlm.nih.gov), then assembled using the DNAman program (LynnonBiosoft, San Ramon, CA, USA) [18]. Thirteen PCGs were first identified using an open reading frame (ORF) finder (https://www.ncbi.nlm.nih.gov/orffinder/) [19] to specify the invertebrate mitochondrial genetic code and translate it into putative proteins based on the Limacidae sequences available in GenBank. The codon usage of thirteen PCGs was computed using the MEGA 6.0 software [20]. The tRNA genes were verified using the MITOS WebServer (http://mitos.bioinf.uni-leipzig.de/index.py), using the default settings [21]. The skewness of nucleotide composition was gauged according to the following formulas: AT skew [(A − T)/(A + T)] and GC skew [(G − C)/(G + C)] [22], where the positive AT skew means that there are more As than Ts; an AT negative skew means that there are less As than Ts, and the same for the GC skew.

Phylogenetic Analysis
A total of 38 species of the 13 PCGs were used to reconstruct the phylogenetic relationship among Gastropoda, including P. bilineatus. The BioEdit software (Tom Hall, windows 10) was used to conduct multiple alignments of amino acid sequences encoded by 13 genes. The complete mitogenome of 37 species was obtained from the GenBank database (Table S2). We chose two Nudibranchia species as outgroups [23], and the mitogenomes of Cyanoplax cavema (NC026848) and Nuttallina californica (NC026849) were downloaded from NCBI [24]. The amino acid sequence of the 13 genes of all selected species was prepared in the following specific order: cox1, nad6, nad5, nad1, nad4L, cob, cox2, atp8, atp6, nad3, nad4, cox3, nad2. After organizing the data, the amino acid sequence of 13 PCGs per species was used for phylogenetic analysis, which was performed with Bayesian inference (BI) methods, using MrBayes v3.2.6 and the maximum likelihood (ML) methods (bootstrap = 1000) by IQ-TREE v1.6.10 [25,26]. PartitionFinder 2.1.1 was used to select the optimal evolutionary models for phylogenetic analysis [27]. The BI analysis was executed using MrBayes with the following command program: Mcmc ngen = 100,000; checkpoint = yes; checkfreq = 5000; sump relburnin = yes; burninfrac = 0.25; and sumt relburnin = yes; burninfrac = 0.25. A good indication of convergence was considered to be reached when the average standard deviation of split frequencies was below 0.01. All analyses converged within 100,000 generations and the resulting phylogenetic trees were directly viewed in FigTree v1.4.2 [28].

Mitogenome Sequence Structure Analysis
The complete mitochondrial genome of P. bilineatus is 14,347 bp ( Figure 1) in length. This was deposited in GenBank with the accession No. MG722906. As in other Cenogastropoda, the mitogenomes are comprised of 38 genes: 13 PCGs (cob, nad1-6, nad4L, cox1-3, atp6 and atp8), 23 tRNAs (one for each amino acid and two each for Ser, Leu, and Cys), two rRNAs (12S and 16S), and two non-coding regions ( Table 1). The PCG region is 11,133 bp long, and nine of the 13 PCGs are encoded on the H-strand (nad1, nad2, nad4, nad4L, nad5, nad6, cob, cox1, and cox2), and the remaining four PCGs are encoded on the L-strand (nad3, cox3, atp6, and atp8) ( Table 1). Such patterns of expression and gene structure of P. bilineatus follow the classical mitochondrial gene hypothesis, as do most other Gastropoda organisms [29,30]. The yellow shows the NADH dehydrogenase family gene, the pink part is the cytochrome c oxidase family gene, the red shows two different ribosomal RNAs, the green acts as two ATPase family genes, the purple represents the cytochrome c biogenesis gene, and the deep blue region is the translation region of the transport RNA.  The yellow shows the NADH dehydrogenase family gene, the pink part is the cytochrome c oxidase family gene, the red shows two different ribosomal RNAs, the green acts as two ATPase family genes, the purple represents the cytochrome c biogenesis gene, and the deep blue region is the translation region of the transport RNA.

Skewness, Overlapping, and Intergenic Spacer Regions
The mitogenome of P. bilineatus has a 182 bp overlap between genes in 15 locations, ranging from 1 to 47 bp, with the longest 47 bp overlap located between trnT and cox3. The mitogenomes of P. bilineatus contain 344 bp of intergenic spacer sequences (1 to 143 bp in length), distributed over 18 regions. The longest spacer sequence is located between 16S and cox1, which are extremely A + T rich (75.52%), and has enough palindromic structures to fold into a specific secondary structure (Figure 2), so we speculate that the structure is a CR region. The nucleotide compositions of the H-strand of P. bilineatus are as follows: A = 4688 (32.68%), T = 5658 (39.44%), G = 2038 (14.21%), C = 1963 (13.69%). The A + T rich region accounts for 72.11% of the total nucleotide composition. Compared with other Gastropod species, the enrichments are higher than those of many marine species sequenced, for example, Aplysia californica (66.34%), Homoiodoris japonica (63.47%), Sakuraeolis japonica (65.77%), Haliotis discus hannai (60.39%), and Platevindex mortoni (62.99%). In contrast, the enrichment is lower that of Achatinella mustelina (80.06%), Achatinella sowerbyana (79.11%), and Succinea putris (76.69%). Species with relatively close A + T rich regions are the following: Aegista diversifamilia (71.07%), Naesiotus nux (73.26%), Meghimatium bilineatum (71.44%), Pupilla muscorum (71.19%), Microceramus pontificus (71.94%), and Vertigo pustelina (72.23%). The highest A + T contents (75.52%) were detected in the control region, which is consistent with previous research on other slugs [31,32]. The AT content is higher than the content of GC, as is generally shown in Gastropoda mitochondrial genomes [33]. Additionally, the AT skew (−0.0938) for the P. bilineatus mitogenome is appreciably negative, reflecting a higher occurrence of Ts to As, and its GC skew (0.0803) is positive, indicating a higher content of Gs than Cs. These findings are somewhat similar to the majority of Gastropoda mitogenomes to date, except for H. japonica [34], which has a negative GC skew.

Figure 2.
Mitochondrial DNA control region secondary structure of P. bilineatus, built using the Predict a Secondary Structure Web Server [35]. The structure is colored by base-pairing probabilities. For unpaired regions, the color denotes the probability of being unpaired.

Protein-Coding Genes
The total nucleotide length of the 13 PCGs in P. bilineatus is 10,839 bp in length, and most of the PCGs start with a typical ATN codon, except for the cox1 and atp8 genes, which start with TTG and GTG, an accepted canonical mitochondrial start codon for invertebrate mitogenomes [36][37][38][39]. The stop codons (TAA and TAG) are utilized in these PCGs. In P. bilineatus, cox3, nad5, nad4L, and cox2 end with TGA, and the remaining genes (cox1, cob, atp8, atp6, nad6, nad1, nad3, nad4, and nad2) terminate with TAA. The incomplete stop (T-) codon is usually found in metazoan mitogenomes, which is presumably completed via post-transcriptional polyadenylation [40], but this does not occur in P. bilineatus. Relative synonymous codon usage values for the P. bilineatus and M. bilineatum mitogenomes are summarized in Figure 3. By comparing P. bilineatus with M. bilineatum, which is more closely related to its evolution, it was found that the 13 PCGs have full codons. In the codon distribution chart (Figure 4), three amino acids (Leu1, 445 and 441; Ile, 265 and 253; and Phe, 322 and 310) are shown to be the most common amino acids in both P. bilineatus and M. bilineatum. Furthermore, the AT skew and GC skew values of the PCGs of the two species are shown in Figure  5. Most of the values in AT skew are negative, except for the atp8 gene. In contrast, the GC skew values show a large difference, and the most significant difference is seen in the two genes nad4L and nad4, which are negative for P. bilineatus and positive for M. bilineatum, respectively. This result suggests that more Ts and Gs are present in most PCGs, which is consistent with most previous observations [41]. Mitochondrial DNA control region secondary structure of P. bilineatus, built using the Predict a Secondary Structure Web Server [35]. The structure is colored by base-pairing probabilities. For unpaired regions, the color denotes the probability of being unpaired.

Protein-Coding Genes
The total nucleotide length of the 13 PCGs in P. bilineatus is 10,839 bp in length, and most of the PCGs start with a typical ATN codon, except for the cox1 and atp8 genes, which start with TTG and GTG, an accepted canonical mitochondrial start codon for invertebrate mitogenomes [36][37][38][39]. The stop codons (TAA and TAG) are utilized in these PCGs. In P. bilineatus, cox3, nad5, nad4L, and cox2 end with TGA, and the remaining genes (cox1, cob, atp8, atp6, nad6, nad1, nad3, nad4, and nad2) terminate with TAA. The incomplete stop (T-) codon is usually found in metazoan mitogenomes, which is presumably completed via post-transcriptional polyadenylation [40], but this does not occur in P. bilineatus. Relative synonymous codon usage values for the P. bilineatus and M. bilineatum mitogenomes are summarized in Figure 3. By comparing P. bilineatus with M. bilineatum, which is more closely related to its evolution, it was found that the 13 PCGs have full codons. In the codon distribution chart (Figure 4), three amino acids (Leu1, 445 and 441; Ile, 265 and 253; and Phe, 322 and 310) are shown to be the most common amino acids in both P. bilineatus and M. bilineatum. Furthermore, the AT skew and GC skew values of the PCGs of the two species are shown in Figure 5. Most of the values in AT skew are negative, except for the atp8 gene. In contrast, the GC skew values show a large difference, and the most significant difference is seen in the two genes nad4L and nad4, which are negative for P. bilineatus and positive for M. bilineatum, respectively. This result suggests that more Ts and Gs are present in most PCGs, which is consistent with most previous observations [41].

Transfer RNA Genes
The prediction results of the 23 tRNA genes (two for Ser, Cys, and Leu and one for each of the other amino acids) in the P. bilineatus mitogenomes are shown in Figure 6. A total of 15 tRNAs are encoded by the H-strand, and the remaining eight tRNAs are encoded by the L-strand ( Table 1). The tRNAs vary in size from 47 bp (trnK) to 70 bp (trnT). This tRNA genomic architecture is identical to all other Gastropoda species examined to date [42], in addition, to an additional trnC. All the 23 tRNA genes could be folded into a typical clover secondary structure [43] with the exception of trnS1 (GCT), trnS2 (TGA), and trnK (TTT). In P. bilineatus, two trnS genes cannot be folded into typical clover-leaf secondary structures due to the lack of the dihydrouracil (DHU) stem and loop, and trnK missing the acceptor stem. There are also two genes lacking a partial sequence structure (DHU stem and loop for trnS; TψC stem and loop for trnG) in M. bilineatum. A double trnC exists with the same sequences and secondary structure and has a length of 291 bp between the two trnCs, however, this intergenic sequence was formed by the gene duplication of trnC and the predicted secondary structure did not form the characteristic structure of CR. We predict that it does not belong to CR region [44].

Transfer RNA Genes
The prediction results of the 23 tRNA genes (two for Ser, Cys, and Leu and one for each of the other amino acids) in the P. bilineatus mitogenomes are shown in Figure 6. A total of 15 tRNAs are encoded by the H-strand, and the remaining eight tRNAs are encoded by the L-strand ( Table 1). The tRNAs vary in size from 47 bp (trnK) to 70 bp (trnT). This tRNA genomic architecture is identical to all other Gastropoda species examined to date [42], in addition, to an additional trnC. All the 23 tRNA genes could be folded into a typical clover secondary structure [43] with the exception of trnS1 (GCT), trnS2 (TGA), and trnK (TTT). In P. bilineatus, two trnS genes cannot be folded into typical clover-leaf secondary structures due to the lack of the dihydrouracil (DHU) stem and loop, and trnK missing the acceptor stem. There are also two genes lacking a partial sequence structure (DHU stem and loop for trnS; TψC stem and loop for trnG) in M. bilineatum. A double trnC exists with the same sequences and secondary structure and has a length of 291 bp between the two trnCs, however, this intergenic sequence was formed by the gene duplication of trnC and the predicted secondary structure did not form the characteristic structure of CR. We predict that it does not belong to CR region [44].

Phylogenetic Analysis
We established the phylogenetic relationships between 38 species based on the amino acid sequences of the 13 PCGs using BI and ML methods. The partition finder output shows that the best model is LG + I + G for the cox1 and nad6 genes and LG + G for the rest of the genes. C. cavema and N. californica were selected as the outgroups. Three species (P. bilineatus, M. bilineatum, and Arion rufus) cluster together within Arionoidea, and P. bilineatus shares a close ancestry with M. bilineatum, which is well supported (posterior probabilities = 0.9989, bootstrap proportions = 100) by the BI and ML analysis (Figure 7). In the Stylommatophora taxon, the Arionoidea, Halicoidea, and Urocoptoidea groups are evolutionarily related and share a close ancestry. In the phylogenetic tree, it can be seen that the earliest appearance of differentiation was marine life, including Prosobranchia (H. discus

Phylogenetic Analysis
We established the phylogenetic relationships between 38 species based on the amino acid sequences of the 13 PCGs using BI and ML methods. The partition finder output shows that the best model is LG + I + G for the cox1 and nad6 genes and LG + G for the rest of the genes. C. cavema and N. californica were selected as the outgroups. Three species (P. bilineatus, M. bilineatum, and Arion rufus) cluster together within Arionoidea, and P. bilineatus shares a close ancestry with M. bilineatum, which is well supported (posterior probabilities = 0.9989, bootstrap proportions = 100) by the BI and ML analysis (Figure 7). In the Stylommatophora taxon, the Arionoidea, Halicoidea, and Urocoptoidea groups are evolutionarily related and share a close ancestry. In the phylogenetic tree, it can be seen that the earliest appearance of differentiation was marine life, including Prosobranchia (H. discus hannai, Haliotis laevigata, and Haliotis rubra). These gradually evolved into Opisthobranchia (A. californica, Aplysia kurodai), Pulmonata (P. mortoni), and even terrestrial organisms. This also provides data support for the theory that life evolved from the ocean to the land [45]. Although species P. bilineatus and M. bilineatum were previously considered to be a member of the Limacoidea superfamily, it can be found in Figure 7 that species P. bilineatus and M. bilineatum are far related to species Deroceras reticulatum of the Limacoidea, so it can be speculated that it is more appropriate to classify the species P. bilineatus and M. bilineatum into the Arionoidea superfamily. The result graph of the evolutionary tree constructed by two different software and algorithms is consistent, which can be used as data support to make the result more reliable. However, there are some differences in the values of the two analysis results at the branches, among which the most obvious difference lies in the branch nodes of the Helicoidea superfamily (Cepaea nemoralis, Helix aspersa, Cylindrus obtusus, Helicella itala).

Evolution of the Breathing Pattern
In the process of biological evolution from sea to land, the transformation of the breathing pattern plays a very important role [46][47][48]. Pulmonates are euthyneuran heterobranchs that have

Evolution of the Breathing Pattern
In the process of biological evolution from sea to land, the transformation of the breathing pattern plays a very important role [46][47][48]. Pulmonates are euthyneuran heterobranchs that have developed a lung sac from the mantle cavity and have a narrow opening called a pneumostome. The mantle cavity acts as a lung cavity in air and has also developed a secondary gill for respiring in water. A comparison of respiratory organs (Figure 8) showed that species living in the sea, such as Aplysia californica, Pleurobranchaea sp., Pleurobranchaea nobaezealandiae, and Homoiodoris japonica had gills or dendritic gills, while amphibians (Onchidium reevesii, Paraincidium reevesii, Platevindex mortoni, and Platevindex sp.) and terrestrial species evolved lung sacs, cutaneous respiration, and a lung cavity. As a transition between Pulmonata and Opisthobranchia, Peronia verruculata, distinctively, has two respiratory organs: Lung sacs and epigenetic gills. The pneumostome (or breathing pore) is a feature (the respiratory opening) of the external body anatomy of an air-breathing land slug or land snail. It is a part of the respiratory system of gastropods and is an opening on the right side of the mantle of a stylommatophoran snail or slug. Inside the mantle cavity, the animal has a highly vascularized area of tissue that functions as a lung [49]. It is speculated that it might have been lost in the process of evolution, where the gills evolved lung sacs and then a secondary gill evolved again to successfully adapt to the surrounding environment from shallow sea to land.
Genes 2019, 10,198 11 of 14 of tissue that functions as a lung [49]. It is speculated that it might have been lost in the process of evolution, where the gills evolved lung sacs and then a secondary gill evolved again to successfully adapt to the surrounding environment from shallow sea to land.

Conclusions
The first complete mitochondrial genomes of P. bilineatus for the family Philomycidae were determined and compared with other Gastropoda species. The mitogenome sequence of P. bilineatus is 14,347 bp in length. The mitogenome is composed of the typical structure of 13 PCGs, two rRNAs, 23 transfer RNA genes (double trnC gene), and two non-coding regions. Both of the complete mitogenomes of the Philomycidae species are typical circular molecules with a similar genome organization and structure as those found in other Mollusc species, but there are some differences in the ordering of genes. Similar to other invertebrate mitogenomes, most of the PCGs utilize ATN, except for cox1 (TTG) and atp8 (GTG). Additionally, all of the tRNAs can be folded into typical cloverleaf secondary structures, with the exception of trnS and trnK, which lack a DHU stem and acceptor stem, respectively. In this paper, two different analysis software, MrBayes and IQ-TREE (ML, bootstraps = 1000), were selected to construct phylogenetic trees for the selected species, and the analysis results were consistent. Through the phylogenetic analysis, the mitogenome sequences were used to resolve the higher-level relationship of Gastropoda. The results showed that P. bilineatus

Conclusions
The first complete mitochondrial genomes of P. bilineatus for the family Philomycidae were determined and compared with other Gastropoda species. The mitogenome sequence of P. bilineatus is 14,347 bp in length. The mitogenome is composed of the typical structure of 13 PCGs, two rRNAs, 23 transfer RNA genes (double trnC gene), and two non-coding regions. Both of the complete mitogenomes of the Philomycidae species are typical circular molecules with a similar genome organization and structure as those found in other Mollusc species, but there are some differences in the ordering of genes. Similar to other invertebrate mitogenomes, most of the PCGs utilize ATN, except for cox1 (TTG) and atp8 (GTG). Additionally, all of the tRNAs can be folded into typical cloverleaf secondary structures, with the exception of trnS and trnK, which lack a DHU stem and acceptor stem, respectively. In this paper, two different analysis software, MrBayes and IQ-TREE (ML, bootstraps = 1000), were selected to construct phylogenetic trees for the selected species, and the analysis results were consistent. Through the phylogenetic analysis, the mitogenome sequences were used to resolve the higher-level relationship of Gastropoda. The results showed that P. bilineatus and M. bilineatum belong to the superfamily Arionoidea, and Arionoidea is most closely related to Helicoidea and Urocoptoidea. This work at the molecular level suggests that P. bilineatus and M. bilineatus should be classified as Arionoidea, not Limacoidea, as previously thought. In addition, this paper also compared and analyzed the evolutionary trend of respiratory organs in the process of biological evolution from sea to land. Our findings provide important molecular data for further studies on the population genetics and meaningful evolutionary biology of Gastropoda. Additional taxonomic work is needed in the future.
Supplementary Materials: The following are available online at http://www.mdpi.com/2073-4425/10/3/198/s1, Table S1: List of primers used in cloning the mitochondrial genome of P. bilineatus, Table S2: List of the complete mitogenomes of the Gastropoda species.