Comparative Mitogenome Analysis of Two Native Apple Snail Species (Ampullariidae, Pomacea) from Peruvian Amazon

Apple snails of the genus Pomacea Perry, 1810 (Mollusca: Caenogastropoda: Ampullariidae) are native to the Neotropics and exhibit high species diversity, holding cultural and ecological significance as an important protein source in Peru. However, most genetic studies in Pomacea have focused mostly on invasive species, especially in Southeast Asia, where they are considered important pests. In this study, we assembled and annotated the mitochondrial genomes of two Pomacea species native to the Peruvian Amazon: Pomacea reevei Ampuero & Ramírez, 2023 and Pomacea aulanieri (Deville & Hupé, 1850). The mitogenomes of P. reevei and P. aulanieri comprise 15,660 and 16,096 bp, respectively, and contain the typical 37 genes of the animal mitochondria with a large control region of 292 bp in P. reevei and 524 bp in P. aulanieri—which fall within the range of what is currently known in Pomacea. Comparisons with previously published mitogenomes in Pomacea revealed differences in the overlapping of adjacent genes, the size of certain protein-coding genes (PCGs) and the secondary structure of some tRNAs that are consistent with the phylogenetic relationships between these species. These findings provide valuable insights into the systematics and genomics of the genus Pomacea.


Introduction
The members of the genus Pomacea Perry, 1810 (Mollusca: Canogastropoda: Ampullariidae) are commonly known as apple snails and can be found in freshwater habitats such as small ponds, lakes, swamps and canals, extending throughout the Neotropical Region, from Argentina to southeast US and the Caribbean Islands [1].However, in recent decades, certain species such as Pomacea canaliculata (Lamarck, 1822) or Pomacea maculata Perry, 1810 have gained growing importance as invasive species, becoming pests of high economic impact on natural habitats and agricultural areas, particularly in China and Southeast Asia [2].
Given the utility of mitochondrial DNA in elucidating animal systematics and evolutionary patterns [3][4][5], these studies have been primarily based on mitochondrial markers.Therefore, mitochondrial genomes have been made available for invasive species such as P. canaliculata [6,7], P. maculata [8] and Pomacea diffusa Blume, 1957 [9,10].Comparative analyses between Pomacea mitochondrial genomes have provided new insights into their genetic diversity [11,12].However, in order to gain a comprehensive understanding of the genus, we need to extend this type of study to native Pomacea species by sequencing and comparing their mitochondrial genomes.
Hayes et al. [13] identified four large clades within Pomacea based on three nuclear genes, as well as the mitochondrial markers COI and 16S rRNA: the widespread P. canaliculata group, including similar species such as P. maculata [14] and Pomacea occulta Yang & Yu, 2019 [15]; the Pomacea bridgesii group, including P. diffusa; and the more restricted Effusa and Flagellata clades.
In the Peruvian Amazon, where apple snails are commonly referred to as "churos", nearly 20 Pomacea species have been reported.However, our understanding of their taxonomy and evolutionary relationships remains limited [16], gaining more attention only in recent times [17][18][19].While the phylogenetic relationships of certain Peruvian species such as Pomacea aulanieri (Deville & Hupé, 1850) of the P. bridgesii clade and the newly described Pomacea reevei Ampuero & Ramírez 2023 of the P. canaliculata clade have been recently elucidated by employing mitochondrial markers [18], the integration of mitochondrial genome information has the potential to offer fresh perspectives on their evolution and phylogenetics.Therefore, the aim of this study was to assemble and annotate the first mitochondrial genome of two native species of Pomacea from the Peruvian Amazon: P. reevei and P. aulanieri.

Sample Collection and DNA Extraction
Samples of P. reevei and P. aulanieri were bought from local fishermen from Loreto, Peru (Napo and Huallaga basins, respectively).The identification of the specimens was based on shell morphology and anatomical characteristics and corroborated using molecular data [18].Furthermore, P. reevei has recently been identified and described as a new species [19].They were relaxed using ethanol, euthanized through thermal shock (4 • C) and then fixed and preserved in 96% ethanol [20].Genomic DNA (gDNA) was extracted from 1 cm 3 of tissue using the E.Z.N.A. Mollusc DNA Kit (OMEGA Bio-tek, Norcross, GA, USA).Concentration and quality of gDNA were measured using the Nanodrop Lite Spectrophotometer.DNA integrity was verified using 1% agarose gel electrophoresis.Voucher specimens (Figure S1) were deposited in the Museo de Historia Natural of the Universidad Nacional Mayor de San Marcos, Lima, Peru.

Genome Sequencing and Assembly
The whole genomes of P. reevei and P. aulanieri were sequenced by Macrogen Inc. (Seoul, South Korea) using an Illumina NovaSeq 6000 platform.Before sequencing, gDNA integrity was verified using TapeStation gDNA Screen Tape on the 4200 TapeStation System.Approximately 5 Gb of raw data from 150 bp paired-end reads were generated for each sample.Raw data quality was evaluated using FastQC (https://www.bioinformatics.babraham.ac.uk/projects/fastqc/, accessed on 15 October 2022) and processed using Fastp [21] by removing reads with >40% bases with Phred quality < Q15, trimming adapter sequences with >6 bases and trimming polyG tails.Mitogenome assemblies were obtained using GetOrganelle [22] with the Pomacea canaliculata mitochondrial genome (KU052865.1)as a starting reference.Bandage [23] was used to verify if the assembly graph was circular.

Comparative Analysis
Nucleotide composition was calculated using MEGA X [31] for complete mitogenomes, PCGs, tRNAs, rRNAs and CR.Compositional asymmetry was calculated using the formulas for AT-skew = (A − T)/(A + T) and GC-skew = (G − C)/(G + C) [32].Relative synonymous codon usage (RSCU) values were calculated using MEGA X [31].Nucleotide diversity (π) was determined using DnaSP6 [33] with a sliding window analysis of 200 bp and a step size of 25 bp.Overall mean p-distance was calculated for each PCG using MEGA X [31].The nonsynonymous nucleotide substitutions per nonsynonymous site (Ka), the synonymous nucleotide substitutions per synonymous site (Ks) and the ratio of nucleotide nonsynonymous to synonymous substitutions (Ka/Ks) were calculated for each PCG using DnaSP6 [33].

Nucleotide Composition
The base composition of the mitogenomes of P. reevei (30.10%A, 44.88% T, 13.14% C, 15.87% G) and P. aulanieri (29.48% A, 39.37% T, 14.83% C, 16.32% G) was similar (Table 3), showing a high content of A + T (70.98% in P. reevei and 68.85% in P. aulanieri), a negative AT-skew (−0.15 in P. reevei and −0.14 in P. aulanieri) and a positive GC-skew (0.09 in P. reevei and 0.05 in P. aulanieri).Comparing the different sets of elements of the mitogenome, the CR (77.4%) showed the highest A + T content in P. reevei, whereas in P. aulanieri, the tRNAs (71.58%) and the CR (71.56%) showed the highest values.In both mitogenomes, the lowest A+T content was found in the PCGs.The PCGs also showed the lowest AT-skew, whereas the rRNAs had the highest AT-skew in both mitogenomes.Interestingly, although the GC-skew values were positive for the whole mitogenome, PCGs, rRNAs and tRNAs of P. reevei and P. aulanieri, the CR of the former showed a positive value, whereas in the latter, this region had a negative value.This difference in the GC-skew of the CR can also be observed in Figure 1.The comparison of the A + T content (Figure 2), between P. reevei, P. aulanieri and other species of Pomacea, showed that the A + T content of P. reevei (70.98%) was similar to species such as P. canaliculata (71.69%),P. maculata (72.31%) and P. occulta (71.93%), whereas the A+T content of P. aulanieri (68.85%) was the lowest among the compared species and similar to P. diffusa (69.49%).A similar pattern was found after comparing the A+T content of the PCGs, rRNAs, tRNAs and especially the CR, where P. reevei (77.40%),P. canaliculata (79.39%),P. maculata (85.82%) and P. occulta (86.57%) of the P. canaliculata clade showed high values, in contrast to P. aulanieri (71.56%) and P. diffusa (73.70%) of the P. bridgesii clade, although these values could be affected by the recovered length of the CR in each species.Negative AT-skew values were observed for the whole mitogenome, protein-coding genes and tRNAs of all Pomacea, whereas positive values were observed in the rRNAs (Figure 3).Although the CR of P. reevei and P. aulanieri showed a negative AT-skew, in other species, this value was positive (e.g., P. canaliculata or P. diffusa).A similar pattern was observed in the GC-skew (Figure 4), where the whole mitogenomes and most of the features (PCGs, rRNAs and tRNAs) showed positive values, except for the CR of P. aulanieri and P. diffusa.Perna and Kocher [32] suggested that differences in the number and direction of AT-skew and GC-skew values reflect variations in selective pressures and nucleotide substitution processes.Additionally, the positive AT-skew values of rRNAs and slightly negative values for tRNAs contrast with the strongly negative values of PCGs.Therefore, this finding may be a signal of the strong selective pressures that force PCGs to eliminate deleterious mutations.

Protein-Coding Genes and Codon Usage
The length of the PCGs ranged from 159 bp (ATP8) to 1707 bp (NAD5) in P. reevei and the same was observed in P. aulanieri, but NAD5 was slightly longer (1728 bp).A comparison of the length of all PCGs among reported Pomacea mitogenomes [6-12] revealed that
The start codon in all PCGs of both species was the same (ATG), similar to reports from P. diffusa [9] and Marisa cornuarietis [43], although in P. canaliculata, P. maculata and P. occulta, the start codon of COIII was ATA.The most common stop codon was TAA, found in 12 proteins of P. reevei and 9 of P. aulanieri.Only ATP6 in P. reevei used TAG, whereas ATP8, ATP6, NAD3 and NAD2 ended with this stop codon in P. aulanieri.Most Pomacea species also showed a preference for TAA as a stop codon [11,12].
Excluding the stop codons, 3734 and 3739 codons were identified in the mitogenome of P. reevei and P. aulanieri, respectively.The codon usage pattern analysis of the PCGs in P. reevei (Table 4, Figure 5) showed that the most frequent amino acids were Leu (16.6%),Ser (14.18%) and Phe (9.1%), and the least used were Cys (1.18%) and Arg (1.5%).The most frequently detected codons in P. reevei were UUA (Leu2) and UUU (Phe), whereas the least common codon was CGG (Arg).Similar values were observed in P. aulanieri (Table S1, Figure S2), although here, CGC (Arg) was the least frequently used codon.In both species, all amino acids have their preferred codons.

Ribosomal and Transfer RNA Genes
Similar to other Pomacea species, the mitochondrial genome of P. reevei and P. aulan contained two ribosomal RNA.The 12S rRNA had a length of 909 bp in P. reevei and 9 bp in P. aulanieri, whereas the length of 16S rRNA was 1340 bp in P. reevei and 1348 bp P. aulanieri.Based on the 12S rRNA and 16S rRNA length in Pomacea, three groups can Sun et al. [47] indicated that nucleotide compositional asymmetry in the coding strand could affect codon and amino acid usage.Therefore, we expected that based on the negative AT-skew and positive GC-skew of the Pomacea mitogenome, amino acids coded by GT-rich codons might be more frequently used than AC-rich codons.We corroborated this in both species (P.reevei and P. aulanieri), since the most abundant amino acids were Leu2 (TTN) and Phe (TTN).Yang et al. [11] found that in Pomacea, codons ending in A or T were more frequent than those ending in C or G.This pattern was found in this study, since in P. reevei and P. aulanieri codons ending in T (47.2% and 43.9%, respectively) and A (35.4% and 33.3%, respectively) collectively represent more than half of all codons.

Ribosomal and Transfer RNA Genes
Similar to other Pomacea species, the mitochondrial genome of P. reevei and P. aulanieri contained two ribosomal RNA.The 12S rRNA had a length of 909 bp in P. reevei and 960 bp in P. aulanieri, whereas the length of 16S rRNA was 1340 bp in P. reevei and 1348 bp in P. aulanieri.Based on the 12S rRNA and 16S rRNA length in Pomacea, three groups can be identified: (1) P. reevei, with a length of 909 and 1340 bp, respectively; (2) P. maculata, P. canaliculata and P. occulta, with a range of lengths of 929-936 and 1331-1336 bp, respectively; and (3) P. aulanieri and P. diffusa, with a range of lengths of 952-958 and 1346-1348 bp, respectively.
Secondary structures of ribosomal RNAs have been poorly explored among mollusks and the structures presented here for P. reevei and P. aulanieri are the first inferred within Pomacea.The secondary structure of the 12S rRNA of P. reevei (Figure 6) and P. aulanieri (Figure S3) had an identical organization of four structural domains, with Domain III as the most conserved region, allowing its successful PCR across a wide variety of taxa [48].Here, we found that the structure of Domain III in P. reevei and P. aulanieri is similar to previous reports in other mollusks [27,[49][50][51].On the other hand, Domains I and II have been poorly studied and they were difficult to infer in both species.Simon et al. [48] indicated that these regions showed less homoplasy than Domains III and IV so they could be more phylogenetically informative and particularly useful for resolving ancient evolutionary relationships.
The 16S rRNA had a length of 1340 bp in P. reevei (Figures 7 and 8) and 1348 bp in P. aulanieri (Figures S4 and S5), showing the same organization, with six structural domains.These secondary structures agree with the mollusk consensus structure proposed by Lydeard et al. [28], and they have the same number of stems and loops as the models inferred for the Caenogastropoda Cacozeliana lacertina (now Cacozeliana granarium (Kiener, 1841)) and Paracrostoma paludiformis (currently Brotia paludiformis (Solem, 1966)) (available at https://crw-site.chemistry.gatech.edu/,accessed on 10 March 2023).Smith and Bond [52] observed that among spiders, the secondary structure of 16S rRNA could be phylogenetically useful at higher taxonomic levels, whereas the loops, as the most variable regions, could be used at family level relationships.No study has evaluated the phylogenetic signal of the secondary structure of the 16S rRNA among Caenogastropoda.
The length of the tRNAs varied from 62 to 73 bp in P. reevei (Table 1), and from 62 to 74 bp in P. aulanieri (Table 2), a wider range than in P. canaliculata, P. maculata and P. occulta (64-70 bp), but smaller than in P. diffusa (62-76 bp).Comparing the length variation along the whole set of 22 tRNAs in Pomacea, only tRNA-Asp had the same length in all species, whereas the most variable tRNAs were tRNA-Ile (70-76 bp), tRNA-Ala (68-74 bp) and tRNA-Thr (66-71 bp).
Most tRNAs showed the typical clover leaf secondary structure, with four arms and loops.However, both tRNA for Serine in P. reevei showed a large loop instead of the typical stem-loop structure in the DHU domain (Figure 9), whereas in P. aulanieri, only tRNA-Ser2 showed this loop (Figure S6).Yang et al. [12] reported a similar loop in the tRNA-Ser1 of P. maculata, but Yang et al. [11] indicated that all tRNAs in P. canaliculata, P. maculata and P. diffusa had the typical clover leaf structure.However, after inferring the secondary structure of the tRNA of P. canaliculata (KJ739609.1),P. maculata (MF401379.1),P. occulta (KR350466.1)and P. diffusa (MF373586.1)using ARWEN 1.2 [25], we found that in the first three species, the tRNA-Ser1 had a loop in the DHU domain, whereas in P. diffusa, tRNA-Ser2 had a loop in the same domain.Therefore, both P. aulanieri and P. diffusa of the P. bridgesii clade had a loop in tRNA-Ser2, whereas P. reevei, P. canaliculata, P. maculata and P. occulta of the P. canaliculata clade had a loop in tRNA-Ser1, although P. reevei also had a loop in tRNA-Ser2.

Control Region
The largest intergenic region of the mitochondrial genome of both species was located between tRNA-Phe and COIII and showed all the specific characteristics for the CR [53,54]: (i) represents the longest intergenic region; (ii) shows a high A+T content; (iii) has a potential secondary structure; and (iv) contains repetitive elements.We retrieved a CR length of 294 bp in P. reevei and 524 bp in P. aulanieri.Both are within the reported range in Pomacea, from 141 bp (P.occulta) to 806 bp (P.diffusa).In another Ampullariid, Marisa cornuarietis, this region is considerably smaller (63 bp) [43].Brauer et al. [54] also found interspecific differences in the CR length of Conus, with the size of this region in Conus consors Sowerby I, 1833 (698 bp) being five-fold longer than in Conus textile Linnaeus, 1758 or Conus borgesi Trovão, 1979.The CR in Conus quercinus Lightfoot, 1786 is even longer (943 bp) [46].Within Pomacea, species from the P. bridgesii clade (e.g., P. aulanieri) seem to have a longer CR than those of the P. canaliculata clade (e.g., P. reevei).

Control Region
The largest intergenic region of the mitochondrial genome of both species was located between tRNA-Phe and COIII and showed all the specific characteristics for the CR [53,54]: (i) represents the longest intergenic region; (ii) shows a high A+T content; (iii) has A repetitive unit of 12 bp (ACATACATACAT) was manually identified in the CR of P. reevei, with eight repetitions on the H-strand and three repetitions on the L-strand.A repeat unit was also found in the CR of P. aulanieri, with a length of 23 bp (ATATAATCTC-TATATGTGTATGG) and six copies on the H-strand.Different repetitive units have been reported in Pomacea: 5 copies of GATACTATAATATAAA (16 bp) in P. occulta [8]; 13 copies on the H-strand and 11 repeats on the L-strand of AAGATACTATAATATA (16 bp) in P. maculata [12]; 11 repeats of TAAGATATAAAGAAACTAAGAGA (23 bp) in P. canaliculata [7]; and 19 copies on the H-strand and 18 repeats on the L-strand of ATCTATACATAC (12 bp) in P. diffusa [9].Maynard et al. [55] found many AT repeats in the Haliotis rubra mitogenome, suggesting that if the number of AT repeats were polymorphic, it could become a marker for individual typing.The variations observed in the number and motif of the repetitive units of the mitochondrial CR of Pomacea suggest that it would be valuable to evaluate alternative strategies to resolve this repetitive region, such as long read sequencing [56].
The inferred secondary structure of the CR of P. reevei and P. aulanieri showed differences in the number and organization of stems and loops between both species (Figure 10).Similar differences have been found between the CR of Conus consors, C. borgesi and C. textile [54].It has been suggested that stem-loop structures may start replication in animal mitochondria [57], but they have been poorly studied in molluscan mitochondrial genomes.

Genetic Variation among Pomacea
The plot of sequence variation exhibited a variable nucleotide diversity along Pomacea mitogenomes (Figure 11).The overall mean p-distance calculated for each protein-coding gene ranged from 0.127 (COI) to 0.244 (ATP8), with other genes such as NAD2 (0.213), NAD6 (0.206), NAD5 (0.198) or NAD1 (0.190) also showing elevated values.In contrast, the ribosomal RNAs showed lower variability values: 0.139 (12S rRNA) and 0.136 (16S rRNA), whereas the CR displayed the highest value of the whole mitogenome (0.435).Castellana et al. [58] found among vertebrates that COI, COII and COIII were the most conserved mitochondrial genes, contrasting with the ATP8 and NAD genes that were the most variable, whereas CytB and ATP6 showed intermediate values.Peretolchina et al. [59] compared four mitogenomes among Caenogastropoda family Baicaliidae and found little variability among COI, recommending instead COII and COIII for phylogenetic studies and CytB and NAD genes for population studies.Here, in Pomacea, we found similar patterns, with ATP8, CytB and NAD genes being the most variable, and COI, COII and COIII as the most conserved.Fonseca et al. [60] suggested that the phylogenies inferred from NAD5 are the most topologically concordant with the one produced by the whole mitogenome.Although we have not tested this for Pomacea, the high variability of this gene suggests it could be useful for phylogenetic studies, along with other genes such as CytB or NAD1.We found in Pomacea that all PCGs had Ka/Ks values below 1 (Figure S7), indicating that they were evolved under purifying selection.Among the 13 PCGs, COI showed the slowest evolutionary rate, whereas NAD6 had the fastest rate, which is concordant with a previous study in the three invasive Pomacea species [11].
In this study, we have found many differences in the mitogenome structure in Pomacea such as the presence and length of overlapping regions, PCG lengths, COIII start codons, tRNA length range, secondary structures of tRNAs for Serine and ribosomal RNA lengths, summarized in Table 5.These differences are concordant with the topology of the molecular phylogeny (Figure 12), although P. reevei of the P. canaliculata clade showed some similarities to P. bridgesii clade species such as P. aulanieri.

Genetic Variation among Pomacea
The plot of sequence variation exhibited a variable nucleotide diversity along Pomacea mitogenomes (Figure 11).The overall mean p-distance calculated for each protein-coding gene ranged from 0.127 (COI) to 0.244 (ATP8), with other genes such as NAD2 (0.213), NAD6 (0.206), NAD5 (0.198) or NAD1 (0.190) also showing elevated values.In contrast, the ribosomal RNAs showed lower variability values: 0.139 (12S rRNA) and 0.136 (16S rRNA), whereas the CR displayed the highest value of the whole mitogenome (0.435).Castellana et al. [58] found among vertebrates that COI, COII and COIII were the most conserved mitochondrial genes, contrasting with the ATP8 and NAD genes that were the most variable, whereas CytB and ATP6 showed intermediate values.Peretolchina et al. [59] compared four mitogenomes among Caenogastropoda family Baicaliidae and found little variability among COI, recommending instead COII and COIII for phylogenetic studies and CytB and NAD genes for population studies.Here, in Pomacea, we found similar patterns, with ATP8, CytB and NAD genes being the most variable, and COI, COII and COIII as the most conserved.Fonseca et al. [60] suggested that the phylogenies inferred from NAD5 are the most topologically concordant with the one produced by the whole mitogenome.Although we have not tested this for Pomacea, the high variability of this gene suggests it could be useful for phylogenetic studies, along with other genes such as CytB or NAD1.We found in Pomacea that all PCGs had Ka/Ks values below 1 (Figure S7), indicating that they were evolved under purifying selection.Among the 13 PCGs, COI showed the slowest evolutionary rate, whereas NAD6 had the fastest rate, which is concordant with a previous study in the three invasive Pomacea species [11].In this study, we have found many differences in the mitogenome structure in Pomacea such as the presence and length of overlapping regions, PCG lengths, COIII start codons, tRNA length range, secondary structures of tRNAs for Serine and ribosomal RNA

Phylogenetic Analysis
The phylogenetic trees recovered using Maximum Likelihood (ML) and Bayesian Inference (BI) showed an almost identical topology, except for Costapex baldwinae

Phylogenetic Analysis
The phylogenetic trees recovered using Maximum Likelihood (ML) and Bayesian Inference (BI) showed an almost identical topology, except for Costapex baldwinae Harasewych, Uribe & Fedosov, 2020 within Neogastropoda (Figures 12 and S8).The inference generated with BI exhibited higher support values than ML, and it was used to represent the evolutionary relationships of Caenogastropoda (Figure 12), whereas the ML tree is shown in Figure S8.Most clades on Figure 12 showed a high support (BI: 1, ML: >90), with the exception of Ampullarioidea + Viviparoidea (BI: 0.83, ML:61), which have been classified with Obscurella hidalgoi (Crosse, 1864) (Cyclophoroidea) within Architaenioglossa, although the monophyly of this group has not been recovered using mitogenomes [61].
Although mitogenome sequence KY008698.1 is identified in GenBank as P. bridgesii (Reeve, 1856), it was published as P. diffusa [10] and posteriorly analyzed with this species name [11,12].We compared this sequence with P. diffusa MF373586.1 and found that the only difference was on the sizes of the CR, being longer in KY008698.1 due to the higher number of repetitive units.Therefore, we concluded that both sequences represent the mitogenome of P. diffusa.

Conclusions
Mitochondrial genomes in Pomacea are highly variable compared to other snail genera.We found striking variations in CR length between P. reevei and P. aulanieri.Moreover, compared with previously published mitogenomes, we found differences in the organization and structure of features such as the PCGs and tRNAs, consistent with the phylogenetic relationships between the P. canaliculata and P. bridgesii clades.A similar comparison with other Pomacea clades, as well as other Ampullariidae genera, would provide new insights into the evolution of these freshwater snails.
We show that comparative analysis of mitochondrial sequences can help us to better understand the evolutionary relationships within the Pomacea genus.Thus, we want to highlight the importance of generating more whole-genome sequence data in order to represent the remaining Pomacea species.

Supplementary Materials:
The following supporting information can be downloaded at https: //www.mdpi.com/article/10.3390/genes14091769/s1, Figure S1 [29].Figure S6: Secondary structures of the tRNA genes in the mitogenome of P. aulanieri.Figure S7: Evolutionary rate of the Pomacea mitogenomes published in GenBank and those from this work.Ka is the nonsynonymous nucleotide substitutions per nonsynonymous site, Ks is the synonymous nucleotide substitutions per synonymous site and Ka/Ks is ratio of nucleotide nonsynonymous to synonymous substitutions.Figure S8: Phylogenetic relationships of Caenogastropoda species based on protein-coding genes and ribosomal RNAs inferred using Maximum Likelihood (ML).Numbers on branches are bootstrap values.Pomacea species sequenced in this study are marked with an asterisk.Table S1: Codon count and relative synonymous codon usage in the mitochondrial genome of P. aulanieri.The asterisk (*) in the table indicates the stop codon.

Figure 6 .
Figure 6.Secondary structure of the 12S rRNA of P. reevei.Domains are indicated with Roman numbers and helices are numbered following Hickson et al. [27].

Figure 6 .
Figure 6.Secondary structure of the 12S rRNA of P. reevei.Domains are indicated with Roman numbers and helices are numbered following Hickson et al. [27].

Figure 7 .
Figure 7. Secondary structure of the I-III domains of the 16S rRNA of P. reevei.Domains are indicated with Roman numbers and helices are numbered following Wuys et al. [29].

Figure 7 .
Figure 7. Secondary structure of the I-III domains of the 16S rRNA of P. reevei.Domains are indicated with Roman numbers and helices are numbered following Wuys et al. [29].

Figure 8 .
Figure 8. Secondary structure of the IV-VI domains of the 16S rRNA of P. reevei.Domains are indicated with Roman numbers and helices are numbered following Wuys et al. [29].

Figure 8 .
Figure 8. Secondary structure of the IV-VI domains of the 16S rRNA of P. reevei.Domains are indicated with Roman numbers and helices are numbered following Wuys et al. [29].

Figure 9 .
Figure 9. Secondary structures of the tRNA genes in the mitogenome of P. reevei.

Figure 9 .
Figure 9. Secondary structures of the tRNA genes in the mitogenome of P. reevei.

Figure 10 .
Figure 10.Secondary structure of the control region of P. reevei and P. aulanieri.

Figure 10 .
Figure 10.Secondary structure of the control region of P. reevei and P. aulanieri.

Figure 11 .
Figure 11.Nucleotide diversity (π) among the Pomacea mitogenomes available in GenBank and the mitogenomes included in this study.

Figure 11 .
Figure 11.Nucleotide diversity (π) among the Pomacea mitogenomes available in GenBank and the mitogenomes included in this study.

Figure 12 .
Figure 12.Phylogenetic relationships of Caenogastropoda species based on protein-coding genes and ribosomal RNAs, inferred with Bayesian Inference.Numbers on branches are Maximum Likelihood bootstrap values (left) and posterior probabilities from BI (right).Pomacea species sequenced in this study are marked with an asterisk.

Figure 12 .
Figure 12.Phylogenetic relationships of Caenogastropoda species based on protein-coding genes and ribosomal RNAs, inferred with Bayesian Inference.Numbers on branches are Maximum Likelihood bootstrap values (left) and posterior probabilities from BI (right).Pomacea species sequenced in this study are marked with an asterisk.
Figure S3: Secondary structure of the 12S rRNA of P. aulanieri.Domains are indicated with Roman numbers and helices are numbered following Hickson et al. [27].Figure S4: Secondary structure of the I-III domains of the 16S rRNA of P. aulanieri.Domains are indicated with Roman numbers and helices are numbered following Wuys et al. [29].Figure S5: Secondary structure of the IV-VI domains of the 16S rRNA of P. aulanieri.Domains are indicated with Roman numbers and helices are numbered following Wuys et al.

Table 1 .
Organization of the mitochondrial genome of P. reevei.

Table 2 .
Organization of the mitochondrial genome of P. aulanieri.

Table 3 .
Nucleotide composition of the mitogenome of P. reevei and P. aulanieri.

Table 4 .
Codon count and relative synonymous codon usage in the mitochondrial genome of P. reevei.The asterisk (*) in the table indicates the stop codon.

Table 5 .
Comparison of mitogenome features in Pomacea, where (+) indicates presence and (−) indicates absence of the feature.