- Research article
- Open Access
Positive selection on panpulmonate mitogenomes provide new clues on adaptations to terrestrial life
BMC Evolutionary Biology volume 16, Article number: 164 (2016)
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 . 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 .
The vast majority of mollusks are marine, but several lineages within the clade Gastropoda have conquered freshwater and terrestrial habitats , e.g., Neritimorpha, Cyclophoroidea, Littorinoidea, Rissooidea, and Panpulmonata . Most terrestrial mollusks belong to the clade Panpulmonata (Gastropoda: Euthyneura) . The invasion of the non-marine realms in Panpulmonata occurred multiple times independently in the clades Ellobioidea, Hygrophila, Stylommatophora and Systellommatophora [6–8]. 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 . On the other hand, several examples from vertebrates showed various molecular mechanisms of adaptation: duplication and functional diversification of keratin genes , expansion of genes encoding olfactory receptors to detect airborne ligands , and positive selection on either nuclear genes involved in the urea cycle , or mitochondrial genes responding to the increase of oxygen during the Devonian . 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 . 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 . In addition, amino acid substitutions have been shown to improve aerobic capacity and adaptation to new environments .
Several studies have reported non-neutral changes in each of the 13 mitochondrial genes . 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 . 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 . The changes likely impacted the biomechanical apparatus that affect the electrochemical gradient in the mitochondria . 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 .
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 . 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: cox1–nad6–nad5–nad1–nad4L–cob–cox2–atp8–atp6–nad3–nad4–cox3–nad2. The cluster cox2-atp8-atp6 is conserved among other gastropods and cephalopods . However, in A. rufus, we found the small rRNA subunit between cox2 and atp8 (cox2-rrnS–atp8–atp6). Furthermore, clusters trnD-trnC-trnF and trnY-trnW-trnG-trnH-trnQ-trnL2 are typical in Ellobioidea and Systellomatophora . 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).
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 %.
Our reconstructed tree is congruent with previous comprehensive phylogenetic analyses in Euthyneura, using a combination of mitochondrial and nuclear genes , and phylogenomics [3, 24] (Fig. 1). Tectipleura (Euopisthobranchia + Panpulmonata)  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. . In addition, there is high support for the clade Euopisthobranchia. This clade was defined by Jörger et al.  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.
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 . 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. , comprises the clades Amphiboloidea, Ellobioidea, Glacidorboidea, Hygrophila, Siphonarioidea, Stylommatophora, and Systellomatophora plus Acochlidia and Sacoglossa, previously regarded as opisthobranchs . 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 , i.e., evidence from phylogenomics did not recover “Pulmonata” as a monophyletic group . 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) , and rejects the Geophila hypothesis (Stylommatophora + Systellomatophora) . Finally, the association between Stylommatophora and Hygrophila has been also found using phylogenomic analyses , 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) : (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 . 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 . The secondary contact after allopatric divergence of haplotypes (hypothesis 2) has been found in Candidula  and Cepaea . Moreover, introgression of mitochondrial lineages as a result of hybridization has been observed in two species of the land snail genus Trochulus . The effect of natural selection (hypothesis 3) shaping the genetic diversity has been proposed before for land snails  although it was considered to be uncommon . However, a recent study has shown that mitochondrial DNA undergoes substantial amounts of adaptive evolution, especially in mollusks . 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 , and explained the persistence of ancestral polymorphisms and the extreme divergence in Achatinella , Systrophia , and Xerocrassa . In addition, in the case of Hygrophila, some studies found high divergence rates in Physella , and Radix , 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 . 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].
Positive selection tests based on either sites or branches only, are conservative for many genes . 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 . For these reasons, we used the branch-site test of positive selection  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).
The branch-site model has been shown to detect ancient episodes of positive selection . 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 . 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 . 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 . 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) ; while terrestrialization in Stylommatophora appeared to be older (100–150 Ma) .
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 . 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 . 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 . Since desiccation stress (or abiotic stress in general) is linked to metabolic activity and ROS production , this directly affects the invasion success of an evolutionary lineage.
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 . 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 . 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 .
In case of nad5, site 474 is homologous to site 540 from the alignment of Garvin et al.  (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  and in eutherians . Also, site 474 corresponds to site 519 in subterranean rodents . This codon position is positively selected only in lineages that independently colonized the subterranean niche, a habitat suggested to be energetically demanding . Mutations in the NADH complex, especially in the transmembrane domains, may affect the proton pump activity of this complex . These changes could facilitate the proton flow and improve the efficiency of ATP production - characteristics associated to increased energetic requirements in non-marine habitats .
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.
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) , Carychium tridentatum (Risso, 1826) , and Helicella itala (Linnaeus, 1758)  for DNA shotgun sequencing. We followed the protocol described by Feldmeyer et al.  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 . 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 . 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).
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 . Ambiguous aligned regions were removed using Gblocks . 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 . 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.
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.  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 . Maximum likelihood analyses were conducted in RAxML-HPC2 (8.0.9) [74, 75] implemented on XSEDE  (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  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  was used to evaluate effective sample sizes (ESS > 200). We assume that a bootstrap value of >70 %  and a posterior probability of > 0.95  are evidence of significant nodal support.
Patterns of evolutionary rates
Relative evolutionary rates were calculated in the software MEGA6 , using the nucleotide and amino acid information following the procedure described by Merker et al. . 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 . 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  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 .
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 , while in the branch-site models, the test focuses on the so-called foreground branches . 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 . 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 . We followed the suggestions from George and Blieck  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.  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.  and references therein [14, 19, 20, 49] were downloaded and aligned using the global homology algorithm g-insi in MAFFT .
Laurin M. How Vertebrates Left Water. Berkeley: University of California Press; 2010.
Little C. The Terrestrial Invasion: An Ecophysiological Approach to the Origins of Land Animals. Cambridge: Cambridge University Press; 1990.
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.
Kameda Y, Kato M. Terrestrial invasion of pomatiopsid gastropods in the heavy-snow region of the Japanese Archipelago. BMC Evol Biol. 2011;11:118.
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.
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.
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.
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.
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.
Vandebergh W, Bossuyt F. Radiation and functional diversification of alpha keratins during early vertebrate evolution. Mol Biol Evol. 2012;29(3):995–1004.
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.
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.
George D, Blieck A. Rise of the earliest tetrapods: an early Devonian origin from marine environment. PLoS One. 2011;6(7), e22136.
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.
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.
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.
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.
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.
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.
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.
Stoger I, Schrodl M. Mitogenomics does not resolve deep molluscan relationships (yet?). Mol Phylogenet Evol. 2013;69(2):376–92.
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.
Grande C, Templado J, Zardoya R. Evolution of gastropod mitochondrial genome arrangements. BMC Evol Biol. 2008;8:61.
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.
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.
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.
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.
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.
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.
Schrödl M. Time to say “Bye-bye Pulmonata”? Spixiana. 2014;37(2):161–4.
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.
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.
Chiba S. Accelerated evolution of land snails Mandarina in the oceanic Bonin Islands: Evidence from mitochondrial DNA sequences. Evolution. 1999;53(2):460–71.
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.
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.
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.
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.
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.
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.
Parmakelis A, Kotsakiozi P, Rand D. Animal mitochondria, positive selection and cyto-nuclear coevolution: insights from pulmonates. PLoS One. 2013;8(4), e61970.
James JE, Piganeau G, Eyre-Walker A. The rate of adaptive evolution in animal mitochondria. Mol Ecol. 2016;25(1):67–78.
Wakeley J. The effects of subdivision on the genetic divergence of populations and species. Evolution. 2000;54(4):1092–101.
Thacker RW, Hadfield MG. Mitochondrial Phylogeny of Extant Hawaiian Tree Snails (Achatinellinae). Mol Phylogen Evol. 2000;16(2):263–70.
Romero P, Ramirez R. Intraspecific divergence and DNA barcodes in Systrophia helicycloides (Gastropoda, Scolodontidae). Rev Peru Biol. 2011;18(2):201–8.
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.
Nolan JR, Bergthorsson U, Adema CM. Physella acuta: atypical mitochondrial gene order among panpulmonates (Gastropoda). J Molluscan Stud. 2014;80(4):388–99.
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.
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.
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.
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.
Yang Z, dos Reis M. Statistical properties of the branch-site test of positive selection. Mol Biol Evol. 2011;28(3):1217–28.
Yang Z. Molecular Evolution: A Statistical Approach. Oxford: Oxford University Press; 2014.
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.
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.
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.
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.
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.
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.
Mittler R. Oxidative stress, antioxidants and stress tolerance. Trends Plant Sci. 2002;7(9):405–10.
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.
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.
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.
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.
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.
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.
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.
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.
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.
Katoh K, Standley DM. MAFFT multiple sequence alignment software version 7: improvements in performance and usability. Mol Biol Evol. 2013;30(4):772–80.
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.
Abascal F, Zardoya R, Telford MJ. TranslatorX: multiple alignment of nucleotide sequences guided by amino acid translations. Nucleic Acids Res. 2010;38:W7–13.
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.
Lanfear R, Calcott B, Kainer D, Mayer C, Stamatakis A. Selecting optimal partitioning schemes for phylogenomic datasets. BMC Evol Biol. 2014;14:82.
Stamatakis A. RAxML-VI-HPC: maximum likelihood-based phylogenetic analyses with thousands of taxa and mixed models. Bioinformatics. 2006;22(21):2688–90.
Stamatakis A, Hoover P, Rougemont J. A rapid bootstrap algorithm for the RAxML Web servers. Syst Biol. 2008;57(5):758–71.
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.
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.
Rambaut A, Suchard MA, Xie D, Drummond AJ. Tracer v1.6. 2014.
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.
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.
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.
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.
Zeileis A, Grothendiek G. zoo: S3 infraestructure for regular and irregular time series. J Stat Softw. 2005;14(6):1–17.
Yang Z. PAML 4: phylogenetic analysis by maximum likelihood. Mol Biol Evol. 2007;24(8):1586–91.
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.
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.
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.
Romero PE, Weigand AM, Pfenninger M. Positive selection on panpulmonate mitogenomes provide new clues on adaptations to terrestrial life. FigShare. 2016. https://dx.doi.org/10.6084/m9.figshare.c.3291377.
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. https://dx.doi.org/10.6084/m9.figshare.c.3291377 .
Consent for publication
Ethics approval and consent to participate
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.
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.
The authors declare that they have no competing interests.
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)
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)
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)
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)
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)
GenBank accession numbers of the 50 sequences used in this study. (XLSX 12 kb)
Primer sequences used to close the gaps and complete the mitochondrial genomes of Arion and Carychium. (XLSX 9 kb)
About this article
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) doi:10.1186/s12862-016-0735-8
- Codon models
- Land invasion
- Positive selection