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

1. Introduction

The species of the genus Pomacea (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 Caribbean Islands [1]. They have significant cultural and economic importance, serving as a traditional food for local populations in the Amazon region since pre-Hispanic times [2,3]. Moreover, due to their protein value and easy cultivation, they offer promising opportunities for aquaculture [4]. However, in the last 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 [5].

Consequently, numerous genetic studies have focused on the identification [6,7], phylogeography [8] and genetic diversity [9] of these invasive species. Since mitochondrial DNA has proved to be very useful to understand animal systematics and evolution [10,11,12], most genetic studies in Pomacea are based on mitochondrial markers, and mitochondrial genomes are available for invasive species such as P. canaliculata [13,14], P. maculata [15] or P. diffusa Blume, 1957 [16,17]. However, despite comparative studies have been performed on these mitogenomes [18,19], no mitochondrial genome has been published for any native Pomacea species to date.

In Peruvian Amazon, where apple snails are commonly referred to as "churos", up to 20 species of Pomacea have been reported. However, their taxonomy and evolutionary relationships are poorly known [20] and only recently they have received more attention [21,22,23]. Although the phylogenetic relationships of some Peruvian species such as Pomacea aulanieri (Deville & Huppé, 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 using mitochondrial markers [22], mitochondrial genome information could offer novel insights into 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: Pomacea reevei and Pomacea aulanieri. We analyzed the gene content, nucleotide composition and secondary structure of the ribosomal RNAs (rRNAs), transfer RNAs (tRNAs) and control region, as well as the phylogenetic relationships based on a comparison with other Pomacea mitogenomes to improve our understanding of the molecular evolution in Pomacea.

2. Materials and Methods

2.1. Sample Collection and DNA Extraction

Samples of Pomacea reevei and Pomacea aulanieri were bought from local fishermen from Loreto, Peru (Napo and Huallaga basins, respectively). They were relaxed using ethanol, euthanized by thermal shock (4°C), and then fixed and preserved in 96% ethanol [24]. Genomic DNA (gDNA) was extracted from 1 cm³ of tissue using the E.Z.N.A. Mollusc DNA Kit (OMEGA Bio-tek, Norcross, GA, USA). Concentration and quality of gDNA was measured using Nanodrop Lite Spectrophotometer. DNA integrity was verified with 1% agarose gel electrophoresis. Voucher specimens were deposited in the Museo de Historia Natural of the Universidad Nacional Mayor de San Marcos, Lima, Peru (MUSM).

2.2. Genome Sequencing and Assembly

The whole genome of Pomacea reevei and P. aulanieri was 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 with FastQC (https://www.bioinformatics.babraham.ac.uk/projects/fastqc/) and processed using Fastp [25] removing reads with >40% bases with Phred quality < Q15, trimming adapter sequences with >6 bases and trimming polyG tails. Mitogenome assembly were obtained using GetOrganelle [26] with the Pomacea canaliculata mitochondrial genome (KU052865.1) as a starting reference. Bandage [27] was used to verify if the assembly graph was circular.

2.3. Annotation of Mitochondrial Genome and Predictions of Secondary Structures

The MITOS2 web server (http://mitos2.bioinf.uni-leipzig.de/index.py) [28] was used to predict the protein-coding genes (PCGs), transfer RNAs (tRNAs) and ribosomal RNAs (rRNAs). The annotation of the tRNAs was confirmed using ARWEN 1.2 (http://130.235.244.92/ARWEN/index.html) [29]. Protein-coding genes were manually corrected searching for open reading frameworks (ORFs) with the NCBI tool ORFfinder (https://www.ncbi.nlm.nih.gov/orffinder/). Secondary structure of the RNAs was inferred using R2DT (https://rnacentral.org/r2dt) [30] and manually optimized using the 12S rRNA Eukaryote reference model [31] and the 16S rRNA Mollusca reference model [32]. Potential secondary structure of the control region was predicted using the RNAfold web server (http://rna.tbi.univie.ac.at//cgi-bin/RNAWebSuite/RNAfold.cgi). All secondary structures were drawn using Inkscape 1.1 (http://www.inkscape.org/). The mitochondrial genome maps of both species were drawn using the CGView tools [33] available at the web server Proksee (https://proksee.ca/).

2.4. Comparative Analysis

The mitochondrial genome and annotation files of P. canaliculata (KJ739609.1), P. diffusa (MF373586.1), P. maculata (MF401379.1) and P. occulta Yang & Yu, 2019 (KR350466.1) were retrieved from GenBank. We performed a re-annotation of these mitochondrial genomes using MITOS2 Web Server [28] and ORFfinder to verify current annotations. Nucleotide composition was calculated using MEGA X [34] for whole mitogenomes, PCGs, tRNAs, rRNAs and control region. Compositional asymmetry was calculated using the formulas for AT-skew = (A - T)/(A + T) and GC-skew = (G - C)/(G + C) [35]. Relative synonymous codon usage (RSCU) values were calculated using MEGA X [34]. Nucleotide diversity (π) was determined using DnaSP6 [36] 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 [34]. 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 [36].

2.5. Phylogenetic Analysis

The phylogenetic analyses included 35 sequences of Caenogastropoda species and two outgroup species: Haliotis rubra (AY588938) and Aplysia californica (AY569552). Orthologs of the nucleotide sequences of the 13 PCGs, 12S rRNA, and 16S rRNA were individually aligned using Muscle 3.8 [37] implemented in Aliview 1.27 [38]. Ambiguous aligned regions were identified and removed using Gblocks 0.91b [39] under relaxed conditions in PhyloSuite v1.2.3 [40]. Individual gene alignments were concatenated using PhyloSuite v1.2.3 [40]. PartitionFinder2 v2.1.1 [41] in PhyloSuite v1.2.3 [40] was used to determine the best partitioning scheme and substitution model. Maximum Likelihood (ML) analysis was performed using IQ-TREE 2.1.2 [42] on the CIPRES Science Gateway web server [43], with 10000 ultrafast bootstrap replicates. Bayesian Inference (BI) analysis was performed using MrBayes 3.2 [44] on the CIPRES Science Gateway web server [43], with two independent runs of four Markov Chain Monte Carlo (MCMC) chains running for 2 million generations and sampling every 1000 generations with a burn-in of 25%. The resulting phylogenetic trees were visualized through FigTree v1.4.3 [45] and edited with Inskape 1.1. (http://www.inkscape.org/).

3. Results and Discussion

3.1. Organization and Structure of the Mitogenomes

The mitochondrial genomes of Pomacea reevei and P. aulanieri (GenBank accesion numbers OR253802 and OR253803, respectively) contained the typical 37 genes (13 PCGs, two rRNAs and 22 tRNAs), and a large non-coding region known as the control region (Figure 1, Table 1 and Table 2), although circular assemblies were not achieved due to repetitive sequences within this region. Previous studies in Pomacea also failed to recovered the complete control region [19] or have not reported it [18]. Most genes, including all PCGs and rRNAs, where encoded on the heavy strand (H-strand), whereas eigth tRNAs (tRNA-Met, tRNA-Tyr, tRNA-Cys, tRNA-Trp, tRNA-Gln, tRNA-Gly, tRNA-Glu, tRNA-Thr) were located on the light strand (L-strand). The lengths of both mitogenomes were similar: 15,660 bp in P. reevei and 16,096 bp in P. aulanieri. These values closely matched previous reports in other Pomacea species, where lengths between 15,516 (P. occulta) and 16,373 bp (P. diffusa) have been reported [13,14,15,16,17,18,19]. Both mitogenomes showed identical gene arrangements, which were consistent with other Ampullaridae species. Including the control region, 24 intergenic regions were identified in P. reevei, ranging from 2 to 292 bp, for a total of 714 bp (4.6% of mitogenome), whereas 25 intergenic regions were identified in P. aulanieri, ranging from 3 to 524, for a total of 1,081 bp (6.7% of mitogenome). The intergenic regions represented 4.8% of the mitogenome in P. canaliculata [13], 3.5% in P. maculata [19] and 10.1% in P. diffusa [19].

The longest overlapping region was found between ND5/tRNA-Phe in P. aulanieri (21 bp), although this overlapping was only 2 bp in P. reevei. This overlapping has also been reported in P. diffusa [18], with a similar length (20 bp) to P. aulanieri, but no in P. canaliculata, P. maculata nor P. occulta [18,19]. Both P. reevei and P. aulanieri also showed an overlapping region of 7 bp between NAD4L/NAD4 that has been reported in other Ampullariidae such as Marisa cornuaretis [46] and P. diffusa [18], as well as other Caenogastropoda [47,48,49]. The re-annotation of the mitogenomes of P. canaliculata, P. maculata and P. occulta in the present study showed that these species also had this overlapping, although it had not been previously reported [13,14,18,19]. P. aulanieri additionally showed two small overlapping regions between tRNA-Pro/ND6 (1 bp) and tRNA-Ser1/ND2 (1 bp). The overlapping between NAD1/tRNA-Pro reported for P. canaliculata, P. maculata and P. occulta [13,14,18,19] has not been found in P. reevei nor P. aulanieri. The protein-coding genes, rRNAs and tRNAs comprised 71.78%, 14.37% and 9.36% of the whole mitochondrial genome of P. reevei, respectively. Similar values were observed in P. aulanieri, where protein-coding genes, rRNAs and tRNAs represented 69.93%, 14.34% and 9.23%, respectively.

3.2. 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 set of elements of the mitogenome, the control region (77.4%) showed the highest A + T content in P. reevei, whereas in P. aulanieri, the tRNAs (71.58%) and the control region (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 control region of the first showed a positive value, whereas in the second this region had a negative value. This difference in the GC-skew of the control region can also be observed in the 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 values 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 control region, 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 control region 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 control region 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), were the whole mitogenomes and all the sets of elements showed positive values, except for the control region of P. aulanieri and P. diffusa. Perna & Kocher suggested that differences in the number and direction of AT-skew and GC-skew values reflect variations in selective pressures and nucleotide substitution processes [35]. Therefore, the positive AT-skew values of rRNAs and slightly negatives values for tRNAs comparing to the strongly negative values of PCGs could be a signal of the strong selective pressures that force PCGs to eliminate deleterious mutations.

3.3. Protein-Coding Genes and Codon Usage

The length of the PCGs ranged from 159 bp (ATP8) to 1,707 bp (NAD5) in P. reevei and the same was observed in P. aulanieri, but NAD5 was slightly longer (1,728 bp). A comparison of the length of all PCGs among reported Pomacea mitogenomes [13,14,15,16,17,18,19] revealed that the size of COI (1536 bp), COII (687 bp), ATP8 (159 bp), NAD6 (495 bp), CytB (1140 bp), NAD4L (297 bp), COIII (780 bp) and NAD3 (354 bp) is conserved among all species. The ATP6 of most Pomacea showed the same length (699 bp), with the exception of P. diffusa (714 bp). In NAD1, NAD4 and NAD5 we identified three groups: (1) P. reevei, with 939, 1,374 and 1,707 bp, respectively; (2) P. maculata, P. canaliculata and P. occulta, with 960, 1,368 and 1,710 bp, respectively; and (3) P. aulanieri and P. diffusa, with 945, 1,362 and 1,728 bp, respectively. The most conspicous variation in size was found in NAD2: 1,062 bp (P. canaliculata and P. maculata), 1,065 bp (P. occulta), 1,071 bp (P. diffusa) and 1,074 bp (P. reevei and P. aulanieri).

The start codon in all PCGs of both species was the same (ATG), similar to the reported in P. diffusa [16] and Marisa cornuarietis [46], 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 stop codon [18,19].

Excluding the stop codons, 3,734 and 3,739 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 5, Figure 6), although here, CGC (Arg) was the least frequently used codon. In both species, all amino acids have their preferred codons.

Sun et al. indicated that nucleotide compositional asymmetry in the coding strand could affect codon and aminoacid usage [50]. Therefore, we expected that based on the negative AT-skew and positive GC-skew of 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 most abundant aminoacids were Leu2 (TTN) and Phe (TTN). Yang et al. (2018b) 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.

3.4. 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 7) and P. aulanieri (Figure 8) had an identical organization, with four structural domains, although the most studied has been the Domain III as is the most conserved region allowing its successul PCR across a wide variety of taxa [51]. Here we found that the structure of the Domain III in P. reevei and P. aulanieri is similar to previous reports in other mollusks [31,52,53,54]. On the other side, the Domains I and II have been poorly studied and they were difficult to infer in both species, although Simon et al. 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 [51].

The 16S rRNA had a length of 1,340 bp in P. reevei (Figure 9 and Figure 10) and 1,348 bp in P. aulanieri (Figure 11 and Figure 12), showing the same organization, with six structural domains. These secondary structures agree with mollusk consensus structure proposed by Lydeard et al. [32], and they have the same number of stem and loops as the models inferred for the Caenogastropoda Cacozeliana lacertina and Paracrostoma paludiformis (available at https://crw-site.chemistry.gatech.edu/). Smith & Bond [55] 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 y 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 folded into 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 13), whereas in P. aulanieri only tRNA-Ser2 showed this loop (Figure 14). Yang et al. [19] reported a similar loop in the tRNA-Ser1 of P. maculata, but Yang et al. [18] indicated that all tRNAs in P. canaliculata, P. maculata and P. diffusa had the typical clover leaf structure. However, after inferring with ARWEN 1.2 [29] the secondary structure of the tRNA of P. canaliculata (KJ739609.1), P. maculata (MF401379.1), P. occulta (KR350466.1) and P. diffusa (MF373586.1) 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.

Figure 13[M1] . Secondary structures of the tRNA genes in the mitogenome of P. reevei.

3.5. Control Region

The incomplete control region of the mitochondrial genome of both species was located between tRNA-Phe and COIII and showed all the specific characteristics for the control region [56,57]: (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 control region 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) [46]. Brauer et al. also found interspecific differences in the control region length of Conus, with the size of this region in C. consors (698 bp) being five-fold longer than in C. textile or C. borgesi [57]. The control region in C. quercinus is even longer (943 bp) [49]. Within Pomacea, species from P. bridgesii clade (e.g. P. aulanieri) seem to have a longer control region than those of P. canaliculata clade (e.g. P. reevei).

A repetitive unit of 12 bp (ACATACATACAT) was manually identified in the control region of P. reevei, with eigth repetitions on H-strand, and three repetitions on L-strand. A repeat unit was also found in the control region of P. aulanieri, with a length of 23 bp (ATATAATCTCTATATGTGTATGG) and six copies on H-strand. Different repetitive units have been reported in Pomacea: five copies of GATACTATAATATAAA (16 bp) in P. occulta [15]; 13 copies on H-strand and 11 repeats on L-strand of AAGATACTATAATATA (16 bp) in P. maculata [19]; 11 repeats of TAAGATATAAAGAAACTAAGAGA (23 bp) in P. canaliculata [13]; and 19 copies on H-strand and 18 repeats on L-strand of ATCTATACATAC (12 bp) in P. diffusa [16]. Maynard et al. found many AT repeats in Haliotis rubra mitogenome, suggesting that if the number of AT repeats were polymorphic, it could become a marker for individual typing [58]. The variations observed in the number and motif of the repetitive units of the mitochondrial control region of Pomacea suggest it would be valuable to evaluate alternative strategies to resolve this repetitive region, such as long read sequencing [59].

The inferred secondary structure of the control region of P. reevei and P. aulanieri showed differences in the number and organization of stems and loops between both species (Figure 15). Similar differences have been found between the control region of Conus consors, C. borgesi and C. textile [57]. It has been suggested that stem-loop structures has a role to start replication in animal mitochondria [60], but they have been poorly studied in molluscan mitochondrial genome.

3.6. Genetic Variation among Pomacea

The nucleotide diversity along the whole mitogenome was evaluated for six Pomacea species (P. reevei, P. aulanieri, P. canaliculata, P. maculata, P. occulta and P. difusa) using a sliding window approach (Figure 16). 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 control region displayed the highest value of the whole mitogenome (0.435). Castellana et al. found among Vertebrates that COI, COII and COIII were the most conserved mitochondrial genes, contrasting with ATP8 and NAD genes that were the most variables, whereas CytB and ATP6 showed intermediate values [61]. Peretolchina et al. compared four mitogenomes among Caenogastropoda family Baicaliidae nad found little variability among COI, recommending instead COII and COIII for phylogenetic studies and CytB and NAD genes for population studies [62]. Here in Pomacea, we found similar patterns, with ATP8, CytB and NAD genes being the most variables, and COI, COII and COIII as the most conserved. Fonseca et al. suggested that the phylogenies inferred from NAD5 are the most topologically concordant with the one produced by the whole mitogenome [63]. Although we have no 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, 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 previous study in the three invasive Pomacea species [18].

3.7. Phylogenetic Analysis

The phylogenetic trees recovered by Maximum Likelihood (ML) and Bayesian Inference (BI) showed an almost identical topology, except for Costapex balwinae within Neogastropoda (Figure 18 and Figure S1). The inference generated with BI exhibited higher support values than ML, and it was used to represent the evolutionary relationships of Caenogastropoda (Figure 18), whereas the ML tree is shown in Figure S1. Most clades on Figure 18 showed a high support (BI: 1.00, ML: >90), with the exception of Ampullarioidea + Viviparoidea (BI: 0.83, ML:61), that have been classified with Obscurella hidalgoi (Cyclophoroidea) within Architaenioglossa, although the monophyly of this group has not been recovered using mitogenomes [64].

Although mitogenome sequence KY008698.1 is identied in GenBank as P. bridgesii, it was published as P. diffusa [17], and posteriorly analyzed with this species name [18,19]. We compared this sequence with P. diffusa MF373586.1, and found that the only difference was on the sizes of the control region, being longer in KY008698.1 due to the higher number of repetitive units. Therefore, we concluded that both represent the mitogenome of P. diffusa.

The phylogenetic relationships within Ampullariidae were recovered with the highest support values, both in the BI and ML trees, with Pomacea recovered as monophyletic. Within Pomacea, we recovered two main clades: (1) P. canaliculata, P. maculata, P. occulta and P. reevei; and (2) P. diffusa and P. aulanieri. Hayes et al. identified four groups within Pomacea based on nuclear and mitchondrial markers: (1) P. canaliculata clade, (2) P. bridgesii clade, (3) Effusa clade, and (4) Flagellata clade [65]. Posteriorly, P. maculata was differentiated within P. canaliculata clade [66]; and Ramírez et al. [22] recovered P. reevei (as Pomacea sp. 1) within P. canaliculata clade, and P. aulanieri within P. bridgesii clade. These phylogenetic relationships are confirmed here based on mitochondrial genomes, and agree with previous results [18,19]. So far, recovered Pomacea mitogenomes have been restricted to P. canaliculata and P. bridgesii clades, but no mitogenome is available for Effusa or Flagellata clades.

Along this study, we have found many differences between the mitogenome structure in Pomacea such as the presence and length of overlapping regions, PCGs lengths, COIII start codon, tRNA length range, secondary structure of tRNAs for Serine and ribosomal RNAs length that are summarized in Table 5. These differences are concordant with the topology of the molecular phylogeny (Figure 18), although P. reevei of the P. canaliculata clade showed some similitudes to P. bridgesii clade species such as P. aulanieri.

4. Conclusions

The comparative analysis of the mitogenomes of P. reevei, P. aulanieri and previously published mitogenomes in Pomacea, has revealed differences in the organization and structure of the mitogenomes, 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.