Skip to main content

Positive selection on panpulmonate mitogenomes provide new clues on adaptations to terrestrial life



Transitions from marine to intertidal and terrestrial habitats resulted in a significant adaptive radiation within the Panpulmonata (Gastropoda: Heterobranchia). This clade comprises several groups that invaded the land realm independently and in different time periods, e.g., Ellobioidea, Systellomatophora, and Stylommatophora. Thus, mitochondrial genomes of panpulmonate gastropods are promising to screen for adaptive molecular signatures related to land invasions.


We obtained three complete mitochondrial genomes of terrestrial panpulmonates, i.e., the ellobiid Carychium tridentatum, and the stylommatophorans Arion rufus and Helicella itala. Our dataset consisted of 50 mitogenomes comprising almost all major panpulmonate lineages. The phylogenetic tree based on mitochondrial genes supports the monophyly of the clade Panpulmonata. Terrestrial lineages were sampled from Ellobioidea (1 sp.) and Stylommatophora (9 spp.). The branch-site test of positive selection detected significant non-synonymous changes in the terrestrial branches leading to Carychium (Ellobiodea) and Stylommatophora. These convergent changes occurred in the cob and nad5 genes (OXPHOS complex III and I, respectively).


The convergence of the non-synonymous changes in cob and nad5 suggest possible ancient episodes of positive selection related to adaptations to non-marine habitats. The positively selected sites in our data are in agreement with previous results in vertebrates suggesting a general pattern of adaptation to the new metabolic requirements. The demand for energy due to the colonization of land (for example, to move and sustain the body mass in the new habitat) and the necessity to tolerate new conditions of abiotic stress may have changed the physiological constraints in the early terrestrial panpulmonates and triggered adaptations at the mitochondrial level.


The transition from water to land is a fascinating evolutionary issue. The realm change occurred several times and across different taxa, from microorganisms to lichens and green plants, and later, arthropods, mollusks, annelids and vertebrates [1]. The multiple transitions required modifications on several systems and organs previously adapted to aquatic habitats. For example, the presence of internal gas exchangers (lungs) to uptake oxygen, or different skin modifications as cuticle and keratin layers to decrease evaporation rates. Other examples include the production of novel compounds (e.g., uric acid and urea) to excrete nitrogen, and the presence of a skeleton and thick body muscles to support the body mass [2].

The vast majority of mollusks are marine, but several lineages within the clade Gastropoda have conquered freshwater and terrestrial habitats [3], e.g., Neritimorpha, Cyclophoroidea, Littorinoidea, Rissooidea, and Panpulmonata [4]. Most terrestrial mollusks belong to the clade Panpulmonata (Gastropoda: Euthyneura) [5]. The invasion of the non-marine realms in Panpulmonata occurred multiple times independently in the clades Ellobioidea, Hygrophila, Stylommatophora and Systellommatophora [68]. For this reason, panpulmonates are a suitable model to study parallel adaptations to new environments, in particular to land.

Examples from invertebrates like the Collembola (Arthropoda: Hexapoda) showed that adaptive changes during terrestrialization have probably occurred in genes related to ion transport, homeostasis, immune response and development [9]. On the other hand, several examples from vertebrates showed various molecular mechanisms of adaptation: duplication and functional diversification of keratin genes [10], expansion of genes encoding olfactory receptors to detect airborne ligands [11], and positive selection on either nuclear genes involved in the urea cycle [12], or mitochondrial genes responding to the increase of oxygen during the Devonian [13]. However, the genomic basis of the transition to the land in panpulmonates is still unknown.

The animal mitochondrial genome encodes 13 proteins that are involved in the production of almost all energy in the eukaryotic cells. These proteins belong to four of the five complexes of the oxidative phosphorylation (OXPHOS) pathway and are under high functional constraints [14]. For example, inefficiencies caused by amino acid changes in these proteins alter the OXPHOS performance, and produce reactive oxygen species (ROS), i.e., molecules that lead to cellular damage and metabolism disruption [15]. In addition, amino acid substitutions have been shown to improve aerobic capacity and adaptation to new environments [16].

Several studies have reported non-neutral changes in each of the 13 mitochondrial genes [17]. Cytochrome c oxidase genes cox1 and cox3 from the freshwater fish Poecilia spp. present substitutions involved in the adaptation to toxic (H2S-rich) environments. These substitutions trigger conformational changes that block the uptake of H2S [18]. In addition, repeated selection at the same structural amino acid location has been found in the mitochondrial complex I from other fish, rodents and snakes [19]. The changes likely impacted the biomechanical apparatus that affect the electrochemical gradient in the mitochondria [17]. Moreover, cytochrome b (cob) in whales demonstrates several signatures of positive selection, in comparison to other artiodactyls. The adaptive changes in cob have been related to changes in the demand of metabolic processes during cetacean cladogenesis and the transition from land to water habitats [20].

Mitochondrial genomes of euthyneuran gastropods represent a promising dataset to screen for adaptive signatures related to water-to-land transitions. As mentioned above, this clade contains terrestrial and intertidal panpulmonates (e.g., Stylommatophora, Systellommatophora, and Ellobiodea), and also freshwater taxa (Hygrophila) as well as other marine clades (e.g., Euopisthobranchia, Nudipleura). Here, we sequenced and annotated three new panpulmonate mitogenomes from the terrestrial ellobid Carychium tridentatum (Risso, 1826) and the stylommatophorans Arion rufus (Linnaeus, 1758) and Helicella itala (Linnaeus, 1758). We used these new mitogenomes in addition to 47 already published euthyneuran mitogenomes, to reconstruct the phylogenetic relationships of Euthyneura. Finally, we evaluate the magnitude of selective pressures that occurred on the branches leading to terrestrial taxa.

Results and discussion

Characteristics of the new panpulmonate mitogenomes

The length of the three new mitochondrial genomes from the terrestrial panpulmonates Carychium tridentatum (Ellobiidae, Ellobioidea), Arion rufus and Helicella itala (Arionidae and Hygromiidae, Stylommatophora) are 13908 bp, 14321 bp, and 13966 bp, respectively. The three mitogenomes all encode for 13 protein-coding genes (PCG), 22 tRNAs, and 2 rRNAs, as reported for most other animal mitogenomes [21]. A detailed overview of the gene annotations can be found in Table 1. Nine genes are encoded on the major strand: cox1, cox2, cob, nad1, nad2, nad4, nad4L, nad5 and nad6; while four are encoded in the minor strand: atp6, atp8, cox3 and nad3. The gene arrangement is similar to other panpulmonates [21, 22]. Basically, the coding-genes are organized as follows: cox1nad6nad5nad1nad4Lcobcox2atp8atp6nad3nad4cox3nad2. The cluster cox2-atp8-atp6 is conserved among other gastropods and cephalopods [20]. However, in A. rufus, we found the small rRNA subunit between cox2 and atp8 (cox2-rrnSatp8atp6). Furthermore, clusters trnD-trnC-trnF and trnY-trnW-trnG-trnH-trnQ-trnL2 are typical in Ellobioidea and Systellomatophora [22]. We found these clusters in the ellobiid C. tridentatum, but also in the stylommatophorans H. itala, and A. rufus, the latter with a slight modification (trnW-trnY).

Table 1 Gene features in the three new panpulmonate mitogenomes

The total length of the PCG is 10923 bp in C. tridentatum, 10935 bp in A. rufus, and 11071 bp H. itala. The GC-content of the PCG is approximately 30 %, being slightly higher in H. itala (34 %). PCG start with five different initiation codons: ATA, ATG, ATT, GTG, TTG. Finally, an AT-rich intergenic spacer between cox3 and trnI has been proposed as the potential origin of replication (POR) in other euthyneurans [22, 23]. We found the same intergenic region in each of the three new mitogenomes. The AT-mean value in the potential POR region was 83 %.

Phylogenetic analyses

Our reconstructed tree is congruent with previous comprehensive phylogenetic analyses in Euthyneura, using a combination of mitochondrial and nuclear genes [5], and phylogenomics [3, 24] (Fig. 1). Tectipleura (Euopisthobranchia + Panpulmonata) [25] are highly supported in both ML and BI analyses. The clade Nudipleura is monophyletic and within this clade, Pleurobranchidae (Berthellina) is the sister group of Nudibranchia, as proposed by Göbbeler et al. [26]. In addition, there is high support for the clade Euopisthobranchia. This clade was defined by Jörger et al. [5] reuniting the clades Umbraculoidea, Anaspidea, Runcinacea, Pteropoda and Cephalaspidea. In our topology, Anaspidea (Aplysia spp.) and the cephalaspideans Bulla, Odontoglaja, Sagaminopteron and Smaragdinella conform a monophyletic group.

Fig. 1

Phylogenetic tree of the Euthyneura based on mitochondrial genomes. Support for all branches is higher than 95/0.99 (bootstrap/posterior probability) unless otherwise indicated. Habitats are represented by colors following the scheme proposed by Schrödl [30]: marine (blue), intertidal (green), freshwater (brown), and terrestrial (red). Mixed colors indicate clades that occur in multiple habitats. Terrestrial taxa (Arion, Carychium and Helicella) with new mitogenomic data are highlighted in bold letters. Photo credits: Natural Museum Rotterdam (A. dactylomela: NMR40266, B. glabrata NMR81004, Bulla sp. NMR52020, C. tridentatum NMR80709, H. itala: NMR59501, O. celtica NMR69303, P. muscorum: NMR48018, R. balthica: NMR76454, S. pectinata NMR78115), Stephen Childs (C. magnifica)/CC BY 2.0, Samuel Chow (H. physis)/CC BY 2.0

Previous mitochondrial phylogenetic reconstructions recovered the monophyly of the former accepted clade “Opisthobranchia” and the paraphyly of “Pulmonata” [22, 23, 27]. However, the topologies derived from mitogenomics have received criticism, for long-branch attraction (LBA) artifacts affecting the topologies in Heterobranchia. In these cases, long-branched stylommatophorans were recovered closer to the root of the clade while they appeared as derived in the nuclear topologies [28]. On the other hand, recent genomic evidence rejected “Opisthobranchia” in favor of Euopisthobranchia as the sister group of Panpulmonata. Phylogenetic reconstructions based on concatenated nuclear and mitochondrial genes [5, 7, 29] as well new phylogenomic studies [3, 24] recovered the paraphyly of “Opisthobranchia”, and support for Panpulmonata.

We were aware of these rooting issues; thus, we choose members of the “Lower Heterobranchia” as outgroup taxa. Our topology recovered monophyletic Panpulmonata and Euopisthobranchia as its sister group. The clade Panpulmonata, defined by Jörger et al. [5], comprises the clades Amphiboloidea, Ellobioidea, Glacidorboidea, Hygrophila, Siphonarioidea, Stylommatophora, and Systellomatophora plus Acochlidia and Sacoglossa, previously regarded as opisthobranchs [30]. Therefore, Panpulmonata possesses an extraordinary diversity in morphology (snails, slugs and intermediate forms), and habitats (marine, intertidal, freshwater and terrestrial).

Recently, the monophyly of “Pulmonata” has been challenged [30], i.e., evidence from phylogenomics did not recover “Pulmonata” as a monophyletic group [3]. In our tree, members from Amphibolidae (Salinator) and Pyramidelloidea (Pyramidella) appear between traditional “Pulmonata” clades, favoring Panpulmonata over “Pulmonata”. Our topology supports the Amphipulmonata clade (Ellobioidea + Systellomatophora) [29], and rejects the Geophila hypothesis (Stylommatophora + Systellomatophora) [30]. Finally, the association between Stylommatophora and Hygrophila has been also found using phylogenomic analyses [24], although with a small subset of euthyneuran taxa.

Patterns of evolutionary rates

The relative evolutionary rates (RER) for amino acids were not equally distributed over the alignment (Additional file 1). The mean RER value in the nad genes was 2–3 times higher than in the cox genes. The RER for the cox1 gene were below the mean rates in our dataset, indicating a higher number of conserved sites.

Stylommatophora (y = 0.6514x; R 2 = 0.9641) along with Hygrophila (y = 0.6708x; R 2 = 0.9770) presented the highest divergence slopes in our data (Additional file 2). This means that fewer nucleotide changes produced more amino acid changes in comparison to the other clades, i. e. non-synonymous changes are more frequent in both Stylommatophora and Hygrophila. Furthermore, Stylommatophora presented the highest absolute values for both nucleotide and amino acid divergence. This result explains the presence of long branches in this clade (Fig. 1).

Several hypotheses have been proposed to explain the extreme divergence found in the mitochondrial DNA of land snails (Stylommatophora) [31]: (1) exceptionally accelerated rate of evolution, (2) haplotype groups previously differentiated in isolated refuges getting into secondary contact, (3) natural (positive) selection preserving the variation, and (4) the particular population structure in pulmonates that allowed them to preserve ancient haplotypes. Accelerated evolution of the mtDNA has been found in several species of land snails and slugs (hypothesis 1). The evolutionary rate in the mtDNA (rRNA) of the snails Euhadra and Mandarina was 10 % per Ma [32, 33], and 5.2 % in the slug Arion [34]. However, this is not a general pattern for all pulmonates, for example the evolutionary rate of Albinaria and Partula was estimated to be 1–1.2 and 2.8 % per Ma, respectively [35]. The secondary contact after allopatric divergence of haplotypes (hypothesis 2) has been found in Candidula [36] and Cepaea [37]. Moreover, introgression of mitochondrial lineages as a result of hybridization has been observed in two species of the land snail genus Trochulus [38]. The effect of natural selection (hypothesis 3) shaping the genetic diversity has been proposed before for land snails [39] although it was considered to be uncommon [40]. However, a recent study has shown that mitochondrial DNA undergoes substantial amounts of adaptive evolution, especially in mollusks [41]. The particular demographic pattern of land snails that produces highly structured populations (hypothesis 4), i.e., “islands” of isolated demes, affects the probability of reciprocal monophyly of two samples and the chance that a gene tree matches the species tree [42], and explained the persistence of ancestral polymorphisms and the extreme divergence in Achatinella [43], Systrophia [44], and Xerocrassa [45]. In addition, in the case of Hygrophila, some studies found high divergence rates in Physella [46], and Radix [47], although no clear hypothesis has been proposed to explain this pattern.

Analyses of selective pressures

Codon substitution models have been widely used to detect adaptive signatures affecting protein evolution [48]. First, we tested for the presence of positive selected codons across the alignments. All comparisons (M2a-M1a, M8-M7, M8-M8a) consistently favored positive selection models M2a and M8 in cox1, cox2, cox3, and cob (p < 0.05) (Table 2). However, the proportion of sites with ω > 1 was extremely low (Additional files 3 and 4). High ω-values in the positively selected genes can be explained by the presence of few synonymous sites affecting dN estimations. On the contrary, a context of strong negative selection could explain ω < 1 for most of the genes [49, 50].

Table 2 Site test of positive selection

Positive selection tests based on either sites or branches only, are conservative for many genes [51]. This is because the test is only significant if the average ω > 1 holds true for all sites or all branches. However, one might expect that positive selection affects only specific sites in specific branches or lineages [52]. For these reasons, we used the branch-site test of positive selection [51] focusing on the branches leading to terrestrial taxa (foreground) within Panpulmonata. The alternative model (model A) fitted significantly better than the null model (model A1) in the genes cox1, cox2, cob and nad5 (Additional file 5). From these four genes, only cob and nad5 presented an ω-ratio higher than one and positively selected codons in the foreground (two sites in cob and six in nad5) (Table 3).

Table 3 Branch-site test of positive selection on the mitochondrial genes

The branch-site model has been shown to detect ancient episodes of positive selection [53]. A potential problem of the test can be the saturation over long evolutionary times; however, simulations have shown that extreme sequence divergence does not generate false positives although it can lead to a high rate of false negatives, especially in older nodes [54]. In the same study, significant levels of positive selection (ω = 6 or 12) were detected at the radiation of bony vertebrates (Euteleostomi), approximately 400–500 Ma [55]. In our data, divergence times are lower than in the vertebrate study: Euthyneura and Panpulmonata probably diverged from their sister groups 250–350 Ma and 150–250 Ma ago, respectively; while panpulmonate clades with terrestrial taxa diverged more recently (Ellobioidea: 140–160 Ma; Stylommatophora: 100–150 Ma) [3, 5].

Convergent adaptations related to realm shifts

The evolution of the lung in early panpulmonates, probably originating from the pallial cavity of an intertidal gilled-ancestor, was the key evolutionary innovation that allowed the diversification to non-marine realms [24]. Both Ellobioidea and Stylommatophora possess lungs, although they colonized the land in different times. Land invasions in Ellobioidea occurred at least twice, one within the genus Pythia (15–25 Ma) and the other in the subfamily Carychiinae (50–100 Ma) [8]; while terrestrialization in Stylommatophora appeared to be older (100–150 Ma) [5].

Different genes appeared to be under positive selection (ω > 1) in terrestrial panpulmonate branches. While cob and nad5 are both part of the OXPHOS pathway, they belong to different complexes (complex III and I, respectively), suggesting that adaptations occurred in several molecular targets. The observed non-synonymous mutations produced similar changes in the amino acid properties, albeit in different regions of each gene. Results from TreeSAAP for both cob and nad5 showed that these changes alter the equilibrium constant (ionization of COOH) property (Fig. 2). This property may influence the protein efficiency reducing ROS production while increasing individual longevity [56]. Alterations in the equilibrium constant may have allowed organisms to better cope with abiotic stress conditions in the new hot, cold or dry habitats. For example, desiccation tolerance has been shown as a limiting factor for the invasion of dry habitats in thiarid freshwater snails [57]. Moreover, the activation of the antioxidant metabolism reducing ROS excess has been linked to desiccation tolerance in the algae Mastocarpus stellatus and Porphyra columbina occurring in the upper intertidal zone [58]. Since desiccation stress (or abiotic stress in general) is linked to metabolic activity and ROS production [59], this directly affects the invasion success of an evolutionary lineage.

Fig. 2

Detection of significant physicochemical amino acids changes using TreeSAAP. This analysis was performed on the genes that were shown to be under positive selection by PAML. Regions above the z-score of 3.09 were significantly different than assumed under neutrality. Only the highest significant category (8) is shown. This category represents a significant change in the equilibrium constant (ionization of COOH) property. a cob, b nad5. Asterisks indicate regions under positive selection according to PAML and TreeSAAP (sites 308 and 512)

The increase in metabolic efficiency has been also related with the terrestrial invasion in other animals, e. g., tetrapods, during the Devonian. Amphibians, lungfishes, and coelacanths presented significant changes in the same equilibrium constant (ionization of COOH) property suggesting an adaptation to increased oxygen levels and changing metabolic requirements [13]. These mutations affected both cob and nad5 tetrapod genes. Also, it is noteworthy to mention that in the terrestrial panpulmonate nad5 gene, the different approaches used by PAML and TreeSAAP found signatures of positive selection along with amino acid property changes in similar regions (Sites 308 and 512; Fig. 2, Table 3).

Changes at the molecular level could have also occurred in different taxa, so we compared the sites under positive selection in our data against previous studies to find sites with similar adaptive patterns (Tables 4 and 5). For cob, site 151 is located in an intermembrane domain. This site is homologous to site 158 in a previous cetacean-artiodactyl alignment (Table 4) and was revealed to be under selective pressure in cetaceans, also influencing the equilibrium constant of the cd2 intermembrane helix [20]. The authors proposed that non-synonymous changes in cob are related to increasing metabolic demands during cladogenesis. Similarly, adaptive evolution in cob has been related to the increase of energy metabolism in response to the evolution of flight in bats [60].

Table 4 Alignment section of the cob gene from Garvin et al. [17]
Table 5 Alignment section of the nad5 gene from Garvin et al. [17]

In case of nad5, site 474 is homologous to site 540 from the alignment of Garvin et al. [17] (Table 5). This site is positively selected in rorquals (Balaenopteridae) and salmons (Salmonidae). Site 474 is part of the biomechanical apparatus that generates the electrochemical gradient and it has been shown to be under positive selection in the Pacific salmon species [19] and in eutherians [14]. Also, site 474 corresponds to site 519 in subterranean rodents [49]. This codon position is positively selected only in lineages that independently colonized the subterranean niche, a habitat suggested to be energetically demanding [61]. Mutations in the NADH complex, especially in the transmembrane domains, may affect the proton pump activity of this complex [14]. These changes could facilitate the proton flow and improve the efficiency of ATP production - characteristics associated to increased energetic requirements in non-marine habitats [62].


We represent evidence of positive selection on several amino acid positions in the mitochondrial complexes I (nad5) and III (cob). These episodes of positive selection occurred in independent branches of panpulmonates with terrestrial taxa (Ellobioidea and Stylommatophora), indicating their possible role during the invasion of the land realm. Most of these sites have been shown to be also under positive selection in several other taxa. Moreover, the general pattern suggests that non-synonymous mutations in both genes are probably linked to increased or altered metabolic requirements. An increased demand for energy due to the colonization of land and the necessity to cope with different abiotic stress conditions may have changed the physiological constraints in the early terrestrial panpulmonates and triggered functional adaptations at the mitochondrial level. Future studies can take into account the predicted codons and the information on the physicochemical changes to test whether these mutations also affect protein structure and function. New genomic information from panpulmonates will most likely reveal even more genes involved in metabolic and structural processes that were key to the colonization of the terrestrial realm.


Mitogenome sequencing

Our dataset comprises 50 complete mitochondrial genomes, and is representative of all described lineages within Euthyneura, except Acochlidia where no mitogenomes have been sequenced so far. We did not consider identical mitogenomes and mitogenomes that were not verified for biological accuracy by GenBank. Accession numbers for each sample are provided in the Additional file 6.

In addition, we used DNA previously isolated from specimens of Arion rufus (Linnaeus, 1758) [63], Carychium tridentatum (Risso, 1826) [64], and Helicella itala (Linnaeus, 1758) [65] for DNA shotgun sequencing. We followed the protocol described by Feldmeyer et al. [66] with some variations: 500 ng DNA per sample was used for library preparation, following the Roche GS FLX Titanium General Library Preparation specifications. Each sample was sequenced on 1/8 of a titanium plate on a 454 sequencer. It should be noted, that approximately 100 specimens of the microgastropod C. tridentatum originating from a single locality had to be pooled to obtain 500 ng of DNA.

Mitogenome assembly and annotation

Newbler v2.0.1 (Roche) was used for contig assembly, with standard settings. Then, we subjected contigs to BlastN and BlastX searches against the mitochondrial genomes of closely related taxa. In addition, to close the remaining gaps after the assembly, we designed flanking primers using Geneious R7 [67]. Primer sequences can be found in Additional file 7. Sequence amplification was conducted using the following PCR conditions: Each 25 μL PCR mix included 1 μL (10 pmol) of each primer, 2.5 μL 10x PCR buffer, 2 μL (100 mM) MgCl2, 0.2 μL (20 mM) dNTPs, 0.3 μL Taq-polymerase (Fermentas), 1.5 μL (10 mg/mL) bovine serum albumin, 12.50 μL ddH2O and 4 μL template DNA in variable concentrations. Temperature conditions: 1 min at 95 °C, followed by 30 cycles of 30 s at 95 °C, 30 s at 52 °C and 30 s at 72 °C, and finally, 3 min at 72 °C. Visualization of PCR products was performed on a 1.4 % agarose gel. Amplicons were cleaned using the QIAquick PCR Purification Kit or the QIAquick Gel Extraction kit (Qiagen) whenever multiple bands were detected. Sanger sequencing was performed using the PCR primer pair (5 pmol) and the BigDye® Terminator v.3.1 Cycle Sequencing Kit (LifeTechnologies, Inc.) on an ABI 3730 capillary sequencer, using the facilities of the Senckenberg BiK-F Laboratory Centre, Frankfurt am Main.

Both, shotgun contigs and Sanger sequences were aligned in Geneious R7 to obtain the complete mitochondrial genomes. The complete mitogenome assemblies were annotated using the MITOS webserver [68]. This program also identified rRNA and tRNA genes. Additionally, we compared the results from MITOS to other annotation strategies like NCBI ORF Finder or Geneious R7, with similar results. Finally, we compared the new gene annotations against other panpulmonate mitochondrial genes to evaluate the length of the reading frames. Newly determined genomes were deposited into GenBank (accession numbers: KT626607, KT696545, KT696546).

Sequence alignments

The 13 protein-coding genes (PCG) were translated into amino acids in Geneious R7 using the invertebrate mitochondrial genetic code, and then aligned using the MAFFT [69]. Ambiguous aligned regions were removed using Gblocks [70]. For downstream analyses, we did not use alignments that had a length below 150 aa (~450 nt) after Gblocks trimming. Thus, nine genes were selected: atp6, cox1, cox2, cox3, cob, nad1, nad2, nad4, nad5. Then, nucleotide sequences of these genes were aligned using TranslatorX [71]. This software aligns the nucleotides by codons taking into account the information from the amino acid alignment. Gblocks was used again in the codon-based alignment with less stringent parameters to trim flanking regions and long gaps. The concatenated alignment length is 9711 nt.

Phylogenetic reconstructions

It has been shown that outgroup selection is important to conceal current hypothesis of euthyneuran phylogeny [21, 72]. Thus, we decided to follow Wägele et al. [28] choosing the “Lower Heterobranchia” clade (Hydatina, Micromelo and Pupa) as the outgroup. For tree reconstruction, we used only the first and second positions of the alignment in order to reduce saturation levels. The alignment length after removing the third position was 6474 nt. Data were partitioned by gene using the partition scheme suggested in PartitionFinder [73]. Maximum likelihood analyses were conducted in RAxML-HPC2 (8.0.9) [74, 75] implemented on XSEDE [76] (CIPRES Science Gateway). We followed the “hard and slow way” suggestions indicated in the manual and selected the best-likelihood tree after 1000 independent runs. Then, branch support was evaluated using bootstrapping with 1000 replicates, and confidence values were drawn in the best-scoring tree. Bayesian inference was conducted in MrBayes v3.2.2 [77] on XSEDE (CIPRES). Two simultaneous Monte Carlo Markov Chains (MCMC) were run, with the following parameters: eight chains of 50 million generations each, sampling every 1000 generations and a burn-in of 25 %. Tracer 1.6 [78] was used to evaluate effective sample sizes (ESS > 200). We assume that a bootstrap value of >70 % [79] and a posterior probability of > 0.95 [80] are evidence of significant nodal support.

Patterns of evolutionary rates

Relative evolutionary rates were calculated in the software MEGA6 [81], using the nucleotide and amino acid information following the procedure described by Merker et al. [82]. The rates are scaled such that the average evolutionary rate across all sites is 1. Sites with a rate lower than 1 evolve slower than the average while sites with a rate higher than 1 evolve faster than the average. The relative rates were estimated under the General Time Reversible (GTR) model (+Γ) for nucleotide sequences (complete concatenated alignment: 9711 nt.) and under the mtREV model (+Γ) for amino acid sequences. The relative rates were scaled in windows with a size of 30 for nucleotides and 10 for amino acids, using the R package zoo [83]. Finally, we compared the amino acid divergence relative to the nucleotide divergence among main clades. Pairwise nucleotide and amino acid distances were calculated under the previously described substitution models.

Analysis of selective pressures

The CODEML program from PAML v4.8a [84] was used to analyze positive selection in each mitochondrial gene. PAML estimates the omega ratio (ω = dN/dS) using maximum likelihood. The omega ratio compares non-synonymous (dN) against synonymous (dS) substitutions per site. Assuming neutrality, ω -values are equal to one; however, ω > 1 is expected if the gene undergoes adapting molecular evolution. In the latter scenario, non-synonymous mutations offer fitness advantages to the individual and have higher fixation probabilities than synonymous mutations [85].

The maximum likelihood topology was set as the guide tree. We re-estimated branch lengths on the tree using codon model M0 (one-ratio) and used them as fixed when fitting the site and branch-site models. In the site models, the ω-ratio is allowed to vary among codons in the alignment [85], while in the branch-site models, the test focuses on the so-called foreground branches [50]. Specifically, we tested branches leading to the air-breathing land snails and slugs (Stylommatophora) and to the terrestrial Carychium (Ellobioidea).

We evaluated site models using a likelihood-ratio test (LRT). First, we compared the selection model (M2a; model = 0, NSites = 2, fix_omega = 0, omega = 5) against the nearly neutral model (M1a; model = 0, NSites = 1, fix_omega = 0, omega = 1) to detect signatures of ω > 1. Then, we used models that calculate ω from a beta distribution. Model M8 (model = 0, NSites = 8, fix_omega = 0, omega = 5, 10 equal class proportions plus one class with ω > 1) was compared against either model M7 (model = 0, NSites = 7, fix_omega = 0, omega = 1, 10 equal class proportions) or model M8a (model = 0, NSites = 8, fix_omega = 0, omega = 1, 10 equal class proportions plus one class with ω = 1).

For the branch-site test we compared the model A (model = 2, NSsites = 2, fix_omega = 0, omega = 5) against the null model A (model = 2, NSsites = 2, fix_omega = 1, omega = 1). The LRT was calculated as follows, 2*(lnL H1 – lnL H0). The Bayes Empirical Bayes (BEB) algorithm implemented in CODEML was used to calculate posterior probabilities of positive selected sites.

Genes detected to be positively selected in the branch-site test were then analyzed in TreeSAAP [86]. This software identifies significant physicochemical changes comparing the distribution of observed changes inferred from a phylogenetic tree against the random distribution of changes under neutrality. The magnitude of the change is rated from 1 (most conservative) to 8 (most radical). A highly significant z-score calculated in TreeSAAP (z > 3.09, p < 0.01) indicates more non-synonymous substitutions than assumed under the neutral model [87]. We followed the suggestions from George and Blieck [13] to increase the accuracy of the test and lower the rate of false positives. Thus, we considered only the most radical changes (category 8, p < 0.01) as significant. Moreover, we focused only on 20 of the 31 amino acid properties available in the software.

Alignment comparisons with previous studies

Almost all previous work on mitochondrial molecular adaptation has been done on vertebrates. We used the amino acid alignments from a recent review by Garvin et al. [17] to evaluate whether the positive selected sites found in our study are also present in a broader taxonomic context and could present a biological function. Mitochondrial sequences of cob and nad5 reported by Garvin et al. [17] and references therein [14, 19, 20, 49] were downloaded and aligned using the global homology algorithm g-insi in MAFFT [69].


  1. 1.

    Laurin M. How Vertebrates Left Water. Berkeley: University of California Press; 2010.

    Book  Google Scholar 

  2. 2.

    Little C. The Terrestrial Invasion: An Ecophysiological Approach to the Origins of Land Animals. Cambridge: Cambridge University Press; 1990.

    Google Scholar 

  3. 3.

    Zapata F, Wilson NG, Howison M, Andrade SC, Jorger KM, Schrödl M, et al. Phylogenomic analyses of deep gastropod relationships reject Orthogastropoda. Proc Biol Sci. 2014;281(1794):20141739.

    Article  PubMed  PubMed Central  Google Scholar 

  4. 4.

    Kameda Y, Kato M. Terrestrial invasion of pomatiopsid gastropods in the heavy-snow region of the Japanese Archipelago. BMC Evol Biol. 2011;11:118.

    Article  PubMed  PubMed Central  Google Scholar 

  5. 5.

    Jörger KM, Stoger I, Kano Y, Fukuda H, Knebelsberger T, Schrödl M. On the origin of Acochlidia and other enigmatic euthyneuran gastropods, with implications for the systematics of Heterobranchia. BMC Evol Biol. 2010;10:323.

    Article  PubMed  PubMed Central  Google Scholar 

  6. 6.

    Kano Y, Neusser TP, Fukumori H, Jörger KM, Schrödl M. Sea-slug invasion of the land. Biol J Linn Soc. 2015;116(2):253–9.

    Article  Google Scholar 

  7. 7.

    Klussmann-Kolb A, Dinapoli A, Kuhn K, Streit B, Albrecht C. From sea to land and beyond--new insights into the evolution of euthyneuran Gastropoda (Mollusca). BMC Evol Biol. 2008;8:57.

    Article  PubMed  PubMed Central  Google Scholar 

  8. 8.

    Romero PE, Pfenninger M, Kano Y, Klussmann-Kolb A. Molecular phylogeny of the Ellobiidae (Gastropoda: Panpulmonata) supports independent terrestrial invasions. Mol Phylogenet Evol. 2016;97:43–54.

    Article  PubMed  Google Scholar 

  9. 9.

    Faddeeva A, Studer RA, Kraaijeveld K, Sie D, Ylstra B, Marien J, et al. Collembolan transcriptomes highlight molecular evolution of hexapods and provide clues on the adaptation to terrestrial life. PLoS One. 2015;10(6), e0130600.

    CAS  Article  PubMed  PubMed Central  Google Scholar 

  10. 10.

    Vandebergh W, Bossuyt F. Radiation and functional diversification of alpha keratins during early vertebrate evolution. Mol Biol Evol. 2012;29(3):995–1004.

    CAS  Article  PubMed  Google Scholar 

  11. 11.

    Nikaido M, Noguchi H, Nishihara H, Toyoda A, Suzuki Y, Kajitani R, et al. Coelacanth genomes reveal signatures for evolutionary transition from water to land. Genome Res. 2013;23(10):1740–8.

    CAS  Article  PubMed  PubMed Central  Google Scholar 

  12. 12.

    Amemiya CT, Alfoldi J, Lee AP, Fan S, Philippe H, Maccallum I, et al. The African coelacanth genome provides insights into tetrapod evolution. Nature. 2013;496(7445):311–6.

    CAS  Article  PubMed  PubMed Central  Google Scholar 

  13. 13.

    George D, Blieck A. Rise of the earliest tetrapods: an early Devonian origin from marine environment. PLoS One. 2011;6(7), e22136.

    CAS  Article  PubMed  PubMed Central  Google Scholar 

  14. 14.

    da Fonseca RR, Johnson WE, O’Brien SJ, Ramos MJ, Antunes A. The adaptive evolution of the mammalian mitochondrial genome. BMC Genomics. 2008;9(1):119.

    Article  PubMed  PubMed Central  Google Scholar 

  15. 15.

    Moreno-Loshuertos R, Acin-Perez R, Fernandez-Silva P, Movilla N, Perez-Martos A, de Rodriguez C, et al. Differences in reactive oxygen species production explain the phenotypes associated with common mouse mitochondrial DNA variants. Nat Genet. 2006;38(11):1261–8.

    CAS  Article  PubMed  Google Scholar 

  16. 16.

    Dalziel AC, Moyes CD, Fredriksson E, Lougheed SC. Molecular evolution of cytochrome c oxidase in high-performance fish (teleostei: Scombroidei). J Mol Evol. 2006;62(3):319–31.

    CAS  Article  PubMed  Google Scholar 

  17. 17.

    Garvin MR, Bielawski JP, Sazanov LA, Gharrett AJ. Review and meta-analysis of natural selection in mitochondrial complex I in metazoans. J Zool Syst Evol Res. 2015;53(1):1–17.

    Article  Google Scholar 

  18. 18.

    Pfenninger M, Lerp H, Tobler M, Passow C, Kelley JL, Funke E, et al. Parallel evolution of cox genes in H2S-tolerant fish as key adaptation to a toxic environment. Nat Commun. 2014;5:3873.

    CAS  Article  PubMed  Google Scholar 

  19. 19.

    Garvin MR, Bielawski JP, Gharrett AJ. Positive Darwinian selection in the piston that powers proton pumps in complex I of the mitochondria of Pacific salmon. PLoS One. 2011;6(9), e24127.

    CAS  Article  PubMed  PubMed Central  Google Scholar 

  20. 20.

    McClellan DA, Palfreyman EJ, Smith MJ, Moss JL, Christensen RG, Sailsbery JK. Physicochemical evolution and molecular adaptation of the cetacean and artiodactyl cytochrome b proteins. Mol Biol Evol. 2005;22(3):437–55.

    CAS  Article  PubMed  Google Scholar 

  21. 21.

    Stoger I, Schrodl M. Mitogenomics does not resolve deep molluscan relationships (yet?). Mol Phylogenet Evol. 2013;69(2):376–92.

    CAS  Article  PubMed  Google Scholar 

  22. 22.

    White TR, Conrad MM, Tseng R, Balayan S, Golding R, de Frias Martins AM, et al. Ten new complete mitochondrial genomes of pulmonates (Mollusca: Gastropoda) and their impact on phylogenetic relationships. BMC Evol Biol. 2011;11:295.

    CAS  Article  PubMed  PubMed Central  Google Scholar 

  23. 23.

    Grande C, Templado J, Zardoya R. Evolution of gastropod mitochondrial genome arrangements. BMC Evol Biol. 2008;8:61.

    Article  PubMed  PubMed Central  Google Scholar 

  24. 24.

    Kocot KM, Halanych KM, Krug PJ. Phylogenomics supports Panpulmonata: opisthobranch paraphyly and key evolutionary steps in a major radiation of gastropod molluscs. Mol Phylogenet Evol. 2013;69(3):764–71.

    Article  PubMed  Google Scholar 

  25. 25.

    Schrödl M, Jörger KM, Klussmann-Kolb A, Wilson NG. Bye bye “Opisthobranchia”! A review on the contribution of mesopsammin sea slugs to euthyneuran systematics. Thalassas. 2011;27:101–12.

    Google Scholar 

  26. 26.

    Göbbeler K, Klussmann-Kolb A. Out of Antarctica?--new insights into the phylogeny and biogeography of the Pleurobranchomorpha (Mollusca, Gastropoda). Mol Phylogenet Evol. 2010;55(3):996–1007.

    Article  PubMed  Google Scholar 

  27. 27.

    Medina M, Lal S, Valles Y, Takaoka TL, Dayrat BA, Boore JL, et al. Crawling through time: transition of snails to slugs dating back to the Paleozoic, based on mitochondrial phylogenomics. Mar Genomics. 2011;4(1):51–9.

    Article  PubMed  Google Scholar 

  28. 28.

    Wägele H, Klussmann-Kolb A, Verbeek E, Schrödl M. Flashback and foreshadowing-a review of the taxon Opisthobranchia. Organisms Div Evol. 2014;14(1):133–49.

    Article  Google Scholar 

  29. 29.

    Dayrat B, Conrad M, Balayan S, White TR, Albrecht C, Golding R, et al. Phylogenetic relationships and evolution of pulmonate gastropods (Mollusca): new insights from increased taxon sampling. Mol Phylogenet Evol. 2011;59(2):425–37.

    Article  PubMed  Google Scholar 

  30. 30.

    Schrödl M. Time to say “Bye-bye Pulmonata”? Spixiana. 2014;37(2):161–4.

    Google Scholar 

  31. 31.

    Thomaz D, Guiller A, Clarke B. Extreme divergence of mitochondrial DNA within species of pulmonate land snails. Proceedings of the Royal Society B-Biological Sciences. 1996;263(1368):363–8.

    CAS  Article  Google Scholar 

  32. 32.

    Hayashi M, Chiba S. Intraspecific diversity of mitochondrial DNA in the land snail Euhadra peliomphala (Bradybaenidae). Biol J Linn Soc. 2000;70(3):391–401.

    Article  Google Scholar 

  33. 33.

    Chiba S. Accelerated evolution of land snails Mandarina in the oceanic Bonin Islands: Evidence from mitochondrial DNA sequences. Evolution. 1999;53(2):460–71.

    CAS  Article  Google Scholar 

  34. 34.

    Pinceel J, Jordaens K, Backeljau T. Extreme mtDNA divergences in a terrestrial slug (Gastropoda, Pulmonata, Arionidae): accelerated evolution, allopatric divergence and secondary contact. J Evol Biol. 2005;18(5):1264–80.

    CAS  Article  PubMed  Google Scholar 

  35. 35.

    Guiller A, Coutellec-Vreto MA, Madec L, Deunff J. Evolutionary history of the land snail Helix aspersa in the Western Mediterranean: preliminary results inferred from mitochondrial DNA sequences. Mol Ecol. 2001;10(1):81–7.

    CAS  Article  PubMed  Google Scholar 

  36. 36.

    Pfenninger M, Posada D. Phylogeographic history of the land snail Candidula unifasciata (Helicellinae, Stylommatophora): fragmentation, corridor migration, and secondary contact. Evolution. 2002;56(9):1776–88.

    Article  PubMed  Google Scholar 

  37. 37.

    Davison A, Clarke B. History or current selection? A molecular analysis of ‘area effects’ in the land snail Cepaea nemoralis. Proceedings of the Royal Society B-Biological Sciences. 2000;267(1451):1399–405.

    CAS  Article  PubMed Central  Google Scholar 

  38. 38.

    Dépraz A, Hausser J, Pfenninger M. A species delimitation approach in the Trochulus sericeus/hispidus complex reveals two cryptic species within a sharp contact zone. BMC Evol Biol. 2009;9:171.

    Article  PubMed  PubMed Central  Google Scholar 

  39. 39.

    Goodacre SL. Population structure, history and gene flow in a group of closely related land snails: genetic variation in Partula from the Society Islands of the Pacific. Mol Ecol. 2002;11(1):55–68.

    Article  PubMed  Google Scholar 

  40. 40.

    Parmakelis A, Kotsakiozi P, Rand D. Animal mitochondria, positive selection and cyto-nuclear coevolution: insights from pulmonates. PLoS One. 2013;8(4), e61970.

    CAS  Article  PubMed  PubMed Central  Google Scholar 

  41. 41.

    James JE, Piganeau G, Eyre-Walker A. The rate of adaptive evolution in animal mitochondria. Mol Ecol. 2016;25(1):67–78.

    CAS  Article  PubMed  Google Scholar 

  42. 42.

    Wakeley J. The effects of subdivision on the genetic divergence of populations and species. Evolution. 2000;54(4):1092–101.

    CAS  Article  PubMed  Google Scholar 

  43. 43.

    Thacker RW, Hadfield MG. Mitochondrial Phylogeny of Extant Hawaiian Tree Snails (Achatinellinae). Mol Phylogen Evol. 2000;16(2):263–70.

    CAS  Article  Google Scholar 

  44. 44.

    Romero P, Ramirez R. Intraspecific divergence and DNA barcodes in Systrophia helicycloides (Gastropoda, Scolodontidae). Rev Peru Biol. 2011;18(2):201–8.

    Google Scholar 

  45. 45.

    Sauer J, Hausdorf B. Reconstructing the evolutionary history of the radiation of the land snail genus Xerocrassa on Crete based on mitochondrial sequences and AFLP markers. BMC Evol Biol. 2010;10:299.

    Article  PubMed  PubMed Central  Google Scholar 

  46. 46.

    Nolan JR, Bergthorsson U, Adema CM. Physella acuta: atypical mitochondrial gene order among panpulmonates (Gastropoda). J Molluscan Stud. 2014;80(4):388–99.

    Article  PubMed  PubMed Central  Google Scholar 

  47. 47.

    Pfenninger M, Cordellier M, Streit B. Comparing the efficacy of morphologic and DNA-based taxonomy in the freshwater gastropod genus Radix (Basommatophora, Pulmonata). BMC Evol Biol. 2006;6:100.

    Article  PubMed  PubMed Central  Google Scholar 

  48. 48.

    Zhai W, Nielsen R, Goldman N, Yang Z. Looking for Darwin in genomic sequences--validity and success of statistical methods. Mol Biol Evol. 2012;29(10):2889–93.

    CAS  Article  PubMed  Google Scholar 

  49. 49.

    Tomasco IH, Lessa EP. The evolution of mitochondrial genomes in subterranean caviomorph rodents: adaptation against a background of purifying selection. Mol Phylogenet Evol. 2011;61(1):64–70.

    Article  PubMed  Google Scholar 

  50. 50.

    Zhang J, Nielsen R, Yang Z. Evaluation of an improved branch-site likelihood method for detecting positive selection at the molecular level. Mol Biol Evol. 2005;22(12):2472–9.

    CAS  Article  PubMed  Google Scholar 

  51. 51.

    Yang Z, dos Reis M. Statistical properties of the branch-site test of positive selection. Mol Biol Evol. 2011;28(3):1217–28.

    CAS  Article  PubMed  Google Scholar 

  52. 52.

    Yang Z. Molecular Evolution: A Statistical Approach. Oxford: Oxford University Press; 2014.

    Book  Google Scholar 

  53. 53.

    Young JN, Rickaby RE, Kapralov MV, Filatov DA. Adaptive signals in algal Rubisco reveal a history of ancient atmospheric carbon dioxide. Philos Trans R Soc Lond B Biol Sci. 2012;367(1588):483–92.

    CAS  Article  PubMed  PubMed Central  Google Scholar 

  54. 54.

    Gharib WH, Robinson-Rechavi M. The branch-site test of positive selection is surprisingly robust but lacks power under synonymous substitution saturation and variation in GC. Mol Biol Evol. 2013;30(7):1675–86.

    CAS  Article  PubMed  PubMed Central  Google Scholar 

  55. 55.

    Berthelot C, Brunet F, Chalopin D, Juanchich A, Bernard M, Noel B, et al. The rainbow trout genome provides novel insights into evolution after whole-genome duplication in vertebrates. Nat Commun. 2014;5:3657.

    Article  PubMed  PubMed Central  Google Scholar 

  56. 56.

    Beckstead WA, Ebbert MT, Rowe MJ, McClellan DA. Evolutionary pressure on mitochondrial cytochrome b is consistent with a role of CytbI7T affecting longevity during caloric restriction. PLoS One. 2009;4(6), e5836.

    Article  PubMed  PubMed Central  Google Scholar 

  57. 57.

    Facon B, Machinle E, Pointier JP, David P. Variation in desiccation tolerance in freshwater snails and its consequences for invasion ability. Biol Invasions. 2004;6:283–93.

    Article  Google Scholar 

  58. 58.

    Flores-Molina MR, Thomas D, Lovazzano C, Núñez A, Zapata J, Kumar M, et al. Desiccation stress in intertidal seaweeds: Effects on morphology, antioxidant responses and photosynthetic performance. Aquat Bot. 2014;113:90–9.

    CAS  Article  Google Scholar 

  59. 59.

    Mittler R. Oxidative stress, antioxidants and stress tolerance. Trends Plant Sci. 2002;7(9):405–10.

    CAS  Article  PubMed  Google Scholar 

  60. 60.

    Shen YY, Liang L, Zhu ZH, Zhou WP, Irwin DM, Zhang YP. Adaptive evolution of energy metabolism genes and the origin of flight in bats. Proc Natl Acad Sci U S A. 2010;107(19):8666–71.

    CAS  Article  PubMed  PubMed Central  Google Scholar 

  61. 61.

    Da Silva C, Tomasco IH, Hoffmann F, Lessa EP. Genes and ecology: accelerated rates of replacement substitutions in the cytochrome b gene of subterranean rodents. The Open Evolution Journal. 2009;3:17–30.

    Google Scholar 

  62. 62.

    Caballero S, Duchene S, Garavito MF, Slikas B, Baker CS. Initial evidence for adaptive selection on the NADH subunit Two of freshwater dolphins by analyses of mitochondrial genomes. PLoS One. 2015;10(5), e0123543.

    Article  PubMed  PubMed Central  Google Scholar 

  63. 63.

    Pfenninger M, Weigand A, Balint M, Klussmann-Kolb A. Misperceived invasion: the Lusitanian slug (Arion lusitanicus auct. non-Mabille or Arion vulgaris Moquin-Tandon 1855) is native to Central Europe. Evol Appl. 2014;7(6):702–13.

    Article  PubMed  PubMed Central  Google Scholar 

  64. 64.

    Weigand AM, Jochum A, Pfenninger M, Steinke D, Klussmann-Kolb A. A new approach to an old conundrum--DNA barcoding sheds new light on phenotypic plasticity and morphological stasis in microsnails (Gastropoda, Pulmonata, Carychiidae). Mol Ecol Resour. 2011;11(2):255–65.

    Article  PubMed  Google Scholar 

  65. 65.

    Steinke D, Albrecht C, Pfenninger M. Molecular phylogeny and character evolution in the Western Palaearctic Helicidae s.l. (Gastropoda: Stylommatophora). Mol Phylogenet Evol. 2004;32(3):724–34.

    CAS  Article  PubMed  Google Scholar 

  66. 66.

    Feldmeyer B, Hoffmeier K, Pfenninger M. The complete mitochondrial genome of Radix balthica (Pulmonata, Basommatophora), obtained by low coverage shot gun next generation sequencing. Mol Phylogenet Evol. 2010;57(3):1329–33.

    CAS  Article  PubMed  Google Scholar 

  67. 67.

    Kearse M, Moir R, Wilson A, Stones-Havas S, Cheung M, Sturrock S, et al. Geneious Basic: an integrated and extendable desktop software platform for the organization and analysis of sequence data. Bioinformatics. 2012;28(12):1647–9.

    Article  PubMed  PubMed Central  Google Scholar 

  68. 68.

    Bernt M, Donath A, Juhling F, Externbrink F, Florentz C, Fritzsch G, et al. MITOS: improved de novo metazoan mitochondrial genome annotation. Mol Phylogenet Evol. 2013;69(2):313–9.

    Article  PubMed  Google Scholar 

  69. 69.

    Katoh K, Standley DM. MAFFT multiple sequence alignment software version 7: improvements in performance and usability. Mol Biol Evol. 2013;30(4):772–80.

    CAS  Article  PubMed  PubMed Central  Google Scholar 

  70. 70.

    Talavera G, Castresana J. Improvement of phylogenies after removing divergent and ambiguously aligned blocks from protein sequence alignments. Syst Biol. 2007;56(4):564–77.

    CAS  Article  PubMed  Google Scholar 

  71. 71.

    Abascal F, Zardoya R, Telford MJ. TranslatorX: multiple alignment of nucleotide sequences guided by amino acid translations. Nucleic Acids Res. 2010;38:W7–13.

    CAS  Article  PubMed  PubMed Central  Google Scholar 

  72. 72.

    Williams ST, Foster PG, Littlewood DT. The complete mitochondrial genome of a turbinid vetigastropod from MiSeq Illumina sequencing of genomic DNA and steps towards a resolved gastropod phylogeny. Gene. 2014;533(1):38–47.

    CAS  Article  PubMed  Google Scholar 

  73. 73.

    Lanfear R, Calcott B, Kainer D, Mayer C, Stamatakis A. Selecting optimal partitioning schemes for phylogenomic datasets. BMC Evol Biol. 2014;14:82.

    Article  PubMed  PubMed Central  Google Scholar 

  74. 74.

    Stamatakis A. RAxML-VI-HPC: maximum likelihood-based phylogenetic analyses with thousands of taxa and mixed models. Bioinformatics. 2006;22(21):2688–90.

    CAS  Article  PubMed  Google Scholar 

  75. 75.

    Stamatakis A, Hoover P, Rougemont J. A rapid bootstrap algorithm for the RAxML Web servers. Syst Biol. 2008;57(5):758–71.

    Article  PubMed  Google Scholar 

  76. 76.

    Miller MA, Pfeiffer W, Schwartz T. Creating the CIPRES Science Gateway for inference of large phylogenetic trees. In: 2010 Gateway Computing Environments Workshop (GCE). New Orleans: IEEE; 2010. p. 1–8.

    Chapter  Google Scholar 

  77. 77.

    Ronquist F, Teslenko M, van der Mark P, Ayres DL, Darling A, Hohna S, et al. MrBayes 3.2: efficient Bayesian phylogenetic inference and model choice across a large model space. Syst Biol. 2012;61(3):539–42.

    Article  PubMed  PubMed Central  Google Scholar 

  78. 78.

    Rambaut A, Suchard MA, Xie D, Drummond AJ. Tracer v1.6. 2014.

    Google Scholar 

  79. 79.

    Douady CJ, Delsuc F, Boucher Y, Doolittle WF, Douzery EJ. Comparison of Bayesian and maximum likelihood bootstrap measures of phylogenetic reliability. Mol Biol Evol. 2003;20(2):248–54.

    CAS  Article  PubMed  Google Scholar 

  80. 80.

    Leaché AD, Reeder TW. Molecular systematics of the Eastern Fence Lizard (Sceloporus undulatus): a comparison of Parsimony, Likelihood, and Bayesian approaches. Syst Biol. 2002;51(1):44–68.

    Article  PubMed  Google Scholar 

  81. 81.

    Tamura K, Stecher G, Peterson D, Filipski A, Kumar S. MEGA6: molecular evolutionary genetics analysis version 6.0. Mol Biol Evol. 2013;30(12):2725–9.

    CAS  Article  PubMed  PubMed Central  Google Scholar 

  82. 82.

    Merker S, Thomas S, Volker E, Perwitasari-Farajallah D, Feldmeyer B, Streit B, et al. Control region length dynamics potentially drives amino acid evolution in tarsier mitochondrial genomes. J Mol Evol. 2014;79(1–2):40–51.

    CAS  Article  PubMed  Google Scholar 

  83. 83.

    Zeileis A, Grothendiek G. zoo: S3 infraestructure for regular and irregular time series. J Stat Softw. 2005;14(6):1–17.

    Article  Google Scholar 

  84. 84.

    Yang Z. PAML 4: phylogenetic analysis by maximum likelihood. Mol Biol Evol. 2007;24(8):1586–91.

    CAS  Article  PubMed  Google Scholar 

  85. 85.

    Yang Z, Nielsen R, Goldman N, Pedersen AMK. Codon-substitution models for heterogeneous selection pressure at amino acid sites. Genetics. 2000;155(1):431–49.

    CAS  PubMed  PubMed Central  Google Scholar 

  86. 86.

    Woolley S, Johnson J, Smith MJ, Crandall KA, McClellan DA. TreeSAAP: selection on amino acid properties using phylogenetic trees. Bioinformatics. 2003;19(5):671–2.

    CAS  Article  PubMed  Google Scholar 

  87. 87.

    McClellan DA, McCracken KG. Estimating the influence of selection on the variable aminoacid sites of the cytochrome b protein functional domains. Mol Biol Evol. 2001;18(6):917–25.

    CAS  Article  PubMed  Google Scholar 

  88. 88.

    Romero PE, Weigand AM, Pfenninger M. Positive selection on panpulmonate mitogenomes provide new clues on adaptations to terrestrial life. FigShare. 2016.

Download references


We would like to thank Dr. Barbara Feldmeyer for providing suggestions to the manuscript and Claudia Nesselhauf for assistance in the laboratory.


The project was supported by the German funding program “LOEWE – Landes-Offensive zur Entwicklung Wissenschaftlich-ökonomischer Exzellenz” of the Hessen State Ministry of Higher Education, Research and the Arts. PER also received a scholarship from “CONCYTEC/CIENCIACTIVA: Programa de becas de doctorado en el extranjero del Gobierno del Perú” (291-2014-FONDECYT).

Availability of data and materials

Mitogenome sequences are available at GenBank, accession numbers KT626607, KT696545, KT696546. Sequence information, multiple sequence alignments, partition schemes, and phylogenetic trees supporting the conclusions of this article are available in the FigShare repository. [88].

Consent for publication

Not applicable.

Ethics approval and consent to participate

Not applicable.

Author’s contributions

PER carried out the phylogenetic and molecular evolution analyses, participated in the mitogenome assembly, conceived the study and wrote the manuscript. AMW carried out the mitogenome assemblies, participated in the sequence alignment, and helped to draft the manuscript. MP participated in the design and coordination of the study and helped to draft the manuscript. All authors read and approved the final manuscript.

Authors’ information

PER is a doctoral student the Faculty of Biosciences - Goethe University Frankfurt am Main and a member of the Molecular Ecology at the Senckenberg Biodiversity and Climate Research Centre. AMW is a post-graduate scientist and project coordinator in the Aquatic Ecosystem Research, University of Duisburg-Essen and a member of the Centre for Water and Environmental Research (ZWU) Essen, University of Duisburg-Essen. MP is a professor at the Faculty of Biosciences - Goethe University Frankfurt am Main and the leader of the Molecular Ecology group at the Senckenberg Biodiversity and Climate Research Centre.

Competing interests

The authors declare that they have no competing interests.

Author information



Corresponding author

Correspondence to Pedro E. Romero.

Additional files

Additional file 1:

Evolutionary rates for nucleotides (NT; blue line) and amino acids (AA; red line) in the euthyneuran mitogenomes. Rates are scaled such that the average evolutionary rate across all sites is 1 (red line). The x-axis shows amino acid positions in the final concatenated alignment. (PDF 481 kb)

Additional file 2:

Amino acid divergence versus nucleotide divergence in mitochondrial genomes of euthyneuran gastropods. Clades are differentiated by colors and symbols as shown in the legend. (PDF 205 kb)

Additional file 3:

Comparison of models M0 (one-ratio), M1a (nearly neutral), and M2a (positive selection), and. Values in bold represent highly significant differences (p < 0.01) from the null model. LRT: Likelihood ratio test, df: degrees of freedom, lnL: log likelihood, np: number of parameters. (XLSX 14 kb)

Additional file 4:

Comparison of models M7, M8 and M8a. Values in bold represent highly significant differences (p < 0.01) from the null model. LRT: Likelihood ratio test, df: degrees of freedom, lnL: log likelihood, np: number of parameters. (XLSX 14 kb)

Additional file 5:

Branch-site test of positive selection on the mitochondrial genes (Model A vs. Null model). Values in bold represent highly significant differences (p < 0.01) from the null model. LRT: Likelihood ratio test, df: degrees of freedom, lnL: log likelihood, np: number of parameters. (XLSX 10 kb)

Additional file 6:

GenBank accession numbers of the 50 sequences used in this study. (XLSX 12 kb)

Additional file 7:

Primer sequences used to close the gaps and complete the mitochondrial genomes of Arion and Carychium. (XLSX 9 kb)

Rights and permissions

Open Access This article is distributed under the terms of the Creative Commons Attribution 4.0 International License (, which permits unrestricted use, distribution, and reproduction in any medium, provided you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons license, and indicate if changes were made. The Creative Commons Public Domain Dedication waiver ( applies to the data made available in this article, unless otherwise stated.

Reprints and Permissions

About this article

Verify currency and authenticity via CrossMark

Cite this article

Romero, P.E., Weigand, A.M. & Pfenninger, M. Positive selection on panpulmonate mitogenomes provide new clues on adaptations to terrestrial life. BMC Evol Biol 16, 164 (2016).

Download citation


  • Codon models
  • Land invasion
  • Mitogenomics
  • Panpulmonata
  • Positive selection