Burrowers from the Past: Mitochondrial Signatures of Ordovician Bivalve Infaunalization.

Bivalves and gastropods are the two largest classes of extant molluscs. Despite sharing a huge number of features, they do not share a key ecological one: gastropods are essentially epibenthic, although most bivalves are infaunal. However, this is not the ancestral bivalve condition; Cambrian forms were surface crawlers and only during the Ordovician a fundamental infaunalization process took place, leading to bivalves as we currently know them. This major ecological shift is linked to the exposure to a different redox environoments (hypoxic or anoxic) and with the Lower Devonian oxygenation event. We investigated selective signatures on bivalve and gastropod mitochondrial genomes with respect to a time calibrated mitochondrial phylogeny by means of dN/dS ratios. We were able to detect 1) a major signal of directional selection between the Ordovician and the Lower Devonian for bivalve mitochondrial Complex I, and 2) an overall higher directional selective pressure on bivalve Complex V with respect to gastropods. These and other minor dN/dS patterns and timings are discussed, showing that the Ordovician infaunalization event left heavy traces in bivalve mitochondrial genomes.

Introduction

The rise of Ediacara-type organisms, as well as the Cambrian explosion, is traditionally connected with the increasing oxygen levels of late-Proterozoic oceans (see, f.i., Anbar and Knoll 2002; Narbonne 2005; Canfield et al. 2007; Sperling et al. 2013; Li et al. 2015; Tostevin et al. 2016; and reference therein). However, knowledge is growing about eukaryotes inhabiting hypoxic, anoxic, or even euxinic environments: such organisms are scattered across many different branches of the tree of eukaryote life (Levin 2003; Engel 2007; Danovaro et al. 2010; Borgonie et al. 2011; Danovaro et al. 2016). Indeed, the existence of several unrelated eukaryotes that are able to live at low oxygen concentration, or even with no oxygen at all (Danovaro et al. 2010, 2016; but see Bernhard et al. 2015), suggests this to be a plesiomorphy of all eukaryotes, whose common ancestor ought to have been a facultative anaerobe (Müller et al. 2012; Mentel et al. 2016).

In any case, the palaeobiochemistry of early eukaryotes and early metazoans is intimately interwoven with the evolutionary history of mitochondria. No doubt can be currently cast on the endosymbiontic origin of mitochondria from α-proteobacterial-like ancestors (Gray et al. 1999; Sicheritz-Ponten and Andersson 2001; Fitzpatrick et al. 2006; Atteia et al. 2009; Koonin 2010; Abhishek et al. 2011; Thrash et al. 2011; Gray 2012; Degli Esposti et al. 2014; Gray 2015), but the detailed steps of this process are largely unknown and several hypotheses have been proposed (Gray 2015).

Organelles of Mitochondrial Origin (OMOs) have been identified in all eukaryotes (Hjort et al. 2010; Shiflett and Johnson 2010; Müller et al. 2012) and notably, some of them work in anaerobic conditions too (Müller et al. 2012; Mentel et al. 2016). Therefore, it is conceivable that the common ancestor of all mitochondria was at least able to survive under hypoxia (Wang and Wu 2014), if not full anoxia (Mentel et al. 2016). As a consequence, all eukaryotes have inherited OMOs (including typical mitochondria) that are in some way effective in anoxic episodes too. Moreover, significant hypoxic/anoxic events are reported in geological strata (f.i., Li et al. 2015). It is therefore conceivable, recall this scenario that eukaryotic life evolved on pre-existent protomitochondrial features, under aerobic, anaerobic, and all intermediate conditions.

Given the gene content of extant genome-owning OMOs, a large proportion of the endosymbiont genome have been lost or transferred to the nucleus, a process termed Genome Reductive Evolution (GRE; Andersson and Kurland 1998; Khachane et al. 2007; Ghiselli et al. 2013; Kannan et al. 2014). Only few protein coding genes (PCGs) are still retained in the organelle genome: typically they are 13 in metazoans’ mitochondria (see, e.g. Boore 1999; Breton et al. 2014), but down to three in apicomplexans (Feagin 1994; Rehkopf et al. 2000). A possible reason for their retention in the mitochondrial genome is that these genes must be directly regulated by mitochondrial environment, being co-location mandatory for an effective redox regulation (CoRR hypothesis; Race et al. 1999; Allen 2003a, 2003b; Lane 2007).

Bivalves and gastropods are the two largest classes of molluscs, with ∼3,550 and ∼7,900 genera, respectively (Millard 2001). Both classes originated in the Lower Cambrian and were initially composed by epibenthic species. However, at the beginning of the Ordovician, their evolutionary pathways diverged, because bivalves experienced a huge change in their living habits. The few Cambrian forms of bivalves slowly disappered and modern bivalves arose (Cope 1996; Fang 2006; Sánchez 2008; Fang and Sánchez 2012; Cope and Kříž 2013; Polechová 2015; Mondal and Harries 2016a, 2016b). This transition was remarkably driven by the invasion of the infaunal (or endobenthic) zone. During the so-called Cambrian Substrate Revolution (Bottjer et al. 2000; Dornbos et al. 2004), Neoproterozoic coherent matgrounds shifted to bioturbated mixgrounds, and the ancestral, surface-crawling bivalve forms, which fed by sediment grazing (Seilacher 1999; Bottjer et al. 2000; Dornbos et al. 2004; Fang 2006; Mondal and Harries 2016a), evolved into filter-feeding sediment burrowers (Fang and Sánchez 2012; Polechová 2015; Mondal and Harries 2016a). This, in turn, was coupled with several major changes in body shape: pedal palps, palp proboscides, gills, as well as foot, were all affected by the infaunalization process and, consequently, heavily reduced or highly modified (Cope 1996; Fang and Sánchez 2012; Cope and Kříž 2013; Mondal and Harries 2016a, 2016b).

At the same time, the infaunalization of bivalves had another, not minor, outcome, as the infaunal zone is also typically hypoxic, and even anoxic in some cases (Anbar and Knoll 2002; Fang and Sánchez 2012). It is not unconceivable that more efficient and enlarged ctenidia, suitable for gas exchange in deeper and possibly anoxic sediments, was itself the trigger that eventually led to the evolution of the distinctive bivalve feeding gills (Morton 1996). However, the adaptation to less oxygenated environments ought to have triggered more profound biochemical changes. If, in some way, the mitochondrial machinery adapted to hypoxic conditions, directional selective pressures on key genes are expected, as well as the ability to exploit different, hypoxia-optimized respiratory substrates.

Specifically, different quinone (Q) types are known from different organisms and they show variable efficiency with respect to oxygen levels. Menaquinone is the dominant membrane Q in some prokaryotes, like gram-positive bacteria (Degli Esposti 2015); plastoquinone is the typical membrane Q of cyanobacteria (Battchikova et al. 2011). According to Aussel and colleagues (2014), the rise of atmospheric oxygen eventually led to ubiquinone, by allowing the required hydroxylation reactions. Ubiquinone is basically ubiquitous across the tree of life: respiratory Complex I is the common name of the largest, multi-subunit complex of the respiratory chain, NADH:ubiquinone oxidoreductase (Wirth et al. 2016).

Rhodoquinone (RQ) is a rare ubiquinone analog and performs poorly as a substrate for respiratory complexes I and II (Lenaz et al. 1968). However, it is known to be correlated with anaerobiosis in many protists and metazoans (Müller et al. 2012; Degli Esposti 2015); notably, bivalves are among those hypoxia-exposed animals that were found to contain RQ (Van Hellemond et al. 1995). Furthermore, a handful of genetic signatures on some Complex I subunits, as well as on the cytochrome b of the cytochrome bc complex (Complex III), were linked to the presence of RQ in the respiratory chain (Degli Esposti 2015). One of the Complex I subunits that show RQ-related mutations is the mitochondrially-encoded NuoH (NAD1; Degli Esposti 2015); the cytchrome b gene (cytb) also maps onto mitochondrial genomes. All this considered, it is conceivable that the use of RQ had a pivotal role in the infaunalization of bivalves and that the examination of nad1 and cytb may show clues of this biochemical fine-tuning, with special reference to the Q reacting chamber.

Indeed, mitochondrial genomes (mtDNAs) of extant bivalves and gastropods are the outcomes of more than 480 Myr of post-Ordovician evolution and they do enclose several genetic signatures of their legacy. In fact, molluscan mtDNAs often show peculiar features. For example, gastropod mtDNAs show high levels of gene-gene overlapping boundaries (White et al. 2011), and are extremely divergent from other mtDNAs (Thomaz et al. 1996; Chiba 1999; Parmakelis and Mylonas 2004; Pinceel et al. 2005; Parmakelis et al. 2013). Bivalve mtDNAs may be very large molecules, up to the 46,985 bp of Scapharca broughtonii (Liu et al. 2013; Plazzi et al. 2016); often contain many regions that are apparently untranslated (Ghiselli et al. 2013); may present the phenomenon called Doubly Uniparental Inheritance (DUI), where two separate, sex-linked mitochondrial lineages are passed from parents to offspring (Breton et al. 2007; Passamonti and Ghiselli 2009; Zouros 2013; Breton et al. 2014); appear to encode supernumerary open reading frames (Breton et al. 2009; Milani et al. 2013); have a terrific degree of gene rearrangement (Vallès and Boore 2006; Simison and Boore 2008; Plazzi et al. 2013); may follow a very unbalanced strand usage (Plazzi et al. 2016).

Even if mitochondria, as aforementioned, most probably retained some sort of ability to live and work in hypoxia, their efficiency with lower concentrations in oxygen was probably improved with infaunalization. Here, we investigated the adaptation to hypoxia in terms of natural selection trends and quinone pool-related signatures on bivalve mitochondrial genomes; moreover, we compared them to gastropods, which did not underwent a massive infaunalization, with respect to the common ancestor of conchiferans (Brusca et al. 2016).

Materials and Methods

Phylogenetic Analyses

Ninety eight bivalves and 110 gastropod complete mtDNAs were downloaded from GenBank. In order to root the Gastropoda tree, five outgroup conchiferan mtDNAs were also downloaded: Katharina tunicata (Polyplacophora), Nautilus macromphalus (Cephalopoda), Graptacme eborea, Siphonodentalium lobatum (Scaphopoda), and Solemya velum (Bivalvia). Bivalve species are detailed in Plazzi et al. (2016); conversely, GenBank Accession Number, as well as systematics of selected gastropod species, are provided in supplementary Additional file S1, Supplementary Material online, for details. The package masking_package v1.1 (available at https://github.com/mozoo/masking_package; Plazzi et al. 2016) was used to 1) align translated PCGs and rDNA genes with T-Coffee (Notredame et al. 2000), 2) mask out phylogenetically noisy and uninformative sites through a consensus among different masking softwares—Aliscore 2.0 (Misof and Misof 2009), BMGE 1.1 (Criscuolo and Gribaldo 2010), Gblocks 0.91b (Castresana 2000), Noisy (Dress et al. 2008), and Zorro (Wu et al. 2012)—and 3) obtain a concatenated alignment of masked sequences. The partitioning scheme, as well as molecular evolution models, was explored through PartitionFinderProtein and PartitionFinder 1.1.0 (Lanfear et al. 2012), opting for the Bayesian Information Criterion and a greedy approach. Furthermore, the simple indel method of Simmons and Ochoterena (2000), implemented in GapCoder (Young and Healy 2003), was used to code indels in the alignments.

A Maximum Likelihood phylogeny was estimated using RA × ML 8.2.0 (Stamatakis 2014). Within the Bivalvia tree, Solemya velum (Opponobranchia) was forced as sister group of all other OTUs. Conversely, Gastropoda were constrained as a monophylum, so that at least some possible artifacts, due to long-branch attraction (LBA) with the 5 outgroup taxa, are avoided. Moreover, the mtDNA of Lottia digitalis (GenBank Accession Number NC_007782) turned out to be very fast-evolving and very prone to create LBA artifacts, therefore this taxon was excluded from the study, thus lowering the number of ingroups to 109. Being Lottia the only known representative of the clade with a published complete mtDNA, unfortunately this also led to the exclusion of the whole Patellogastropoda from our analysis.

The ML analysis was extensively described elsewhere (see, in particular, see supplementary Additional file S4, Supplementary Material online of Plazzi et al. 2016). Briefly, several preliminary inferences were carried out to explore the best-performing combinations of manual/automatic rearrangement radius and number of rate categories for the CAT model accounting for evolutionary rate heterogeneity (Stamatakis 2006); then, the Best-Known Likelihood (BKL) tree was inferred under the selected parameter set calling 10 runs from 10 randomized MP starting trees; finally, 1000 boostrap replicates were run with the same settings and consensus support values were annotated on the BKL tree.

The same datasets were used for Bayesian Inference (BI), exploiting MrBayes 3.2.1 (Ronquist et al. 2012). Two separate analyses (with four chains each) were run for 10,000,000 generations of MC3, sampling every 100 trees. Convergence was assessed by manual inspection of standard deviation of average split frequencies sampled every 1,000 generation and of Potential Scale Reduction Factor (PSRF; Gelman and Rubin 1992). Stable standard deviation of average split frequencies was also used to manually identify the burn-in point.

Time calibration on the BKL tree was carried out with r8s 1.70 (Sanderson 2003). Selected clades were chosen as calibration point and their first appearance was computed thanks to the Paleobiology Database (http://fossilworks.org): namely, these clades are listed in Plazzi et al. (2016) for bivalves, whereas for gastropods they are the complete Gastropoda class, Conidae, Heterobranchia, Neritidae, Planorbidae, and Stylommatophora. We used the Langley-Fitch (clock) method and the Truncated Newton algorithm; several rounds of cross-validation were run in order to estimate the best-performing smoothing parameter up to the sixth decimal digit; five restarts and five guesses were used each round.

Trees were graphically edited and annotated with PhyloWidget (Jordan and Piel 2008), Dendroscope 3.3.2 (Huson and Scornavacca 2012), and FigTree 1.4.2 (Rambaut 2006-2014).

dN/dS Analyses

The software PAML 4.8a (Yang 1997; 2007) was used for all dN/dS analyses, providing the BKL tree (rooted and with fixed branch lengths) as the user tree. The masked amino acid alignments of bivalves and gastropods were retrotranslated into the original nucleotide sequences, thus retaining only codons relative to phylogenetically informative sites that were previously selected; this retrotranslation was carried out in the R environment (R Development Core Team 2008), loading the package seqinr (Charif and Lobry 2007). Genes were clustered by complex; five different datasets were prepared (table 1).

Table 1

Datasets of Codon Alignments

Dataset	Complex	Genes
CI	Complex I	nad1, nad2, nad3, nad4, nad4L, nad5, nad6
CIII	Complex III	cytb
CIV	Complex IV	cox1, cox2, cox3
CV	Complex V	atp6, atp8a
total	Complexes I+III+IV+V	nad1, nad2, nad3, nad4, nad4L, nad5, nad6, cytb, cox1, cox2, cox3, atp6, atp8a

Given that atp8 was originally excluded from bivalve phylogeny due to many phylogenetic issues (Plazzi et al. 2016), this gene was included in the CV and total gastropod datasets only.

As a first step, a dN/dS analysis on each dataset was requested to PAML constraining a single dN/dS ratio along the entire tree; alternatively, a dN/dS analysis was also run allowing a specific dN/dS ratio for each given branch in the tree. We used the likelihood ratio test (LRT) to determine the best-fitting model; all LRTs were computed by R and, to be conservative, a χ2distribution with N − 1 degrees of freedom was used, where N is the number of branches of the rooted tree (Wong et al. 2004).

The selected dN/dS ratios were then plotted against time. At fixed time intervals, starting from the root age, the dN/dS ratios exhibited by all branches present in the age-calibrated tree at that time lapse were plotted against the age of that time lapse: this resulted in a bivariate distribution, with an increasing number of observation towards the recent times, because the number of tree branches is obviously growing in the same direction.

This distribution was determined both for bivalves and gastropods and for all the five datasets. Within-class bivariate distribution of dN/dS ratios > 1 were tested for significant differences through the Cramér test (Baringhaus and Franz 2004) implemented in the cramer R package, using 1,000 permutation Monte-Carlo-bootstrap replicates and ϕCramer as the kernel function.

On the other side, all points with a dN/dS ratio > 1 were considered for each class (Bivalvia and Gastropoda) and for each dataset and their distribution with respect to time was computed. Among-class significant differences in the distribution of dN/dS ratios (>1) over time were tested using the Kolmogorov–Smirnov approach, which is also implemented in R.

Finally, Phanerozoic concentrations of atmospheric oxgen were taken from Glasspool and Scott (2008) and superimposed over the described distributions.

Quinone Pool-Related Signatures

All available cytb and nad1 sequences from metazoan mtDNAs were downloaded from the OGRE database (Jameson et al. 2003) in February, 2017: in this way 1,240 (cytb) and 1,241 (nad1) sequences were obtained and aligned with PSI-Coffee (Notredame et al. 2000). Regions of cytb and nad1 genes that may show RQ-related signatures were taken from Degli Esposti (2015) and located on the complete gastropod protein alignments from the present work and on the complete bivalve protein alignments from our previous one (Plazzi et al. 2016), as well as on total metazoan alignments. Sitewise composition analysis and χ2 tests were carried out with Microsoft Excel®. Conservation and functional analyses were carried out with the TeXshade package (Beitz 2000); given that the complete analysis of the total metazoans alignments turned out to be too computative expensive, a random subset of 124 sequences (10%) was drawn from the whole datasets using a custom R script.

Results

Phylogenetic Analyses of Class Gastropoda

Amino acid and rDNA alignment lengths after masking phase, as well as masking-surviving site percentages with respect to the original alignment and sequence lengths, are shown in supplementary Additional files S2 and S3, Supplementary Material online, for details. The best partitioning scheme selected by PartitionFinder separated cox3 + nad2 genes, all remaining PCGs genes, and rRNAs (which were taken together; see supplementary Additional file S4, Supplementary Material online for details).

The topology yielded by the RA × ML ML analysis is largely concordant with the topology yielded by the BI MC3 run, therefore only the ML tree is shown in figure 1, after collapsing nodes with BP < 60, along with the phylogenetic tree of the class Bivalvia collapsed and redrawn from Plazzi et al. (2016). The tree shown in figure 1 is collapsed at higher taxonomic levels; original trees are presented as supplementary Additional file S5, Supplementary Material online, for details (BKL tree with BP) and see supplementary Additional file S6, Supplementary Material online, for details (BI tree with PP).

F

Outline of gastropod evolution over ages. Redrawn from the BKL tree computed from 114 complete mtDNAs (taking only PCGs and rDNAs into account). Nodal support was estimated by 1,000 bootstrap replicates: only nodes with BP > 60 are shown. Species were clustered by high-level taxa; time scale in million years. Taxa marked with asterisks are polyphyletic. In particular, Hypsogastropoda were also recovered as polyphyletic; however, in that the bulk of Hypsogastropoda was monophyletic, we labeled it as “Hypsogastropoda sensu stricto”. We refer to supplementary Additional file S5, Supplementary Material online for the complete ML tree with BP, see supplementary Additional file S6, Supplementary Material online for the complete BI tree with PP, and see supplementary Additional file S7, Supplementary Material online for the complete ultrametric tree. The phylogenetic tree of bivalves, redrawn from Plazzi et al. (2016) is shown in the top left insert for the sake of comparison. We refer to Neoheterodontei as in Taylor et al. (2007); “Unionoidea (M)” indicates male-transmitted mtDNA, while “Unionoidea (F and F-like)” indicates female-transmitted mtDNA, as well as mtDNAs of non-DUI species (for the DUI phenomenon, see references in text). a–i, representative gastropods (not to scale): a, Fissurella picta (Vetigastropoda); b, Planorbis sp. (Planorbidae); c, Vermetidae sp. (Vermetidae); d, Cymbiola nobilis (Hypsogastropoda sensu stricto); e, Rapana venosa (Hypsogastropoda sensu stricto); f, Helix lucorum (Stylommatophora); g, Haliotis tuberculata (Vetigastropoda); h, Neverita josephinia (Naticidae); i, Ampullaria sp. (Ampullariidae).

Three clades emerge from the Gastropoda node: Neritimorpha, here represented by four species of the genus Nerita (BP = 100; PP = 1.000); Vetigastropoda (BP = 93; PP = 1.000); and a cluster of Caenogastropoda + Heterobranchia (BP = 67; unsupported by BI). Our phylogenetic dataset was not suitable to resolve the deepest node of Gastropoda, that is, it was not able to identify which of these nodes is sister group to all other gastropods. However, it has to be noted that the cluster Neritimorpha + Vetigastropoda was recovered under BI with PP = 0.995.

See supplementary Additional file S7, Supplementary Material online, for details shows the complete ultrametric tree along with the geological scale taken from Cohen et al. (2013), whereas ages and evolutionary rates of single branches are detailed in supplementary Additional file S8, Supplementary Material online, for details. The origin of Neritimorpha was constrained to 265 Ma; the origin of Vetigastropoda (node N104) was estimated to 461.93 Ma (Middle Ordovician). The common ancestor of Caenogastropoda and Heterobranchia (node N92) would have lived 509.20 Ma, in the middle Cambrian.

Caenogastropoda originated earlier than Heterobranchia, being dated to 478.62 Ma in the Tremadoc (node N91), and cladogenetic events took place slowly and continuously for at least 180 Myr: in the second part of the Carboniferous (Pennsylvanian), most of the biodiversity of Caenogastropoda was already established. Conversely, Heterobranchia originated later (409.10 Ma, in the Lower Devonian), but a quick expansion of the clade was already accomplished in its main lineages within 40 Myr, before the end of the period (around the Frasnian).

In all cases, the null hypothesis that a single dN/dS applies to all the tree branches was rejected by the LRT (P = 0). dN/dS over time is shown for bivalves in supplementary Additional file S9, Supplementary Material online, for details, whereas significance levels of pairwise comparisons are given in table 2. Overall, it is evident that atp6 underwent higher degrees of directional selective pressure, as it is demonstrated by a large number of points with dN/dS > 1 along all the phylogenetic history of the Class. Conversely, for Complex I and IV and for cytb directional selective pressure seems to be concentrated before the staggering increase of atmospheric oxygen in the Lower Devonian (400 Ma; Glasspool and Scott 2008): more specifically, whereas for Complex I it lasted up to 400 Ma, for Complexes IV and cytb it substantially ended with the Cambrian–Ordovician transition (485.4 Ma; Cohen et al. 2013). In fact, all distributions are significantly different from each other, with the only exception of the CI/cytb comparison (table 2).

Table 2

Pairwise Cramér Test

	CI	CIII	CIV	CV	Total
Bivalvia
CI		6.49E-02	2.00E-03***	0.00E+00***	9.41E-01
CIII	178.66		0.00E+00***	0.00E+00***	7.39E-02
CIV	1,455.47	1,996.01		0.00E+00***	3.00E-03***
CV	1,214.03	942.23	2,532.72		9.99E-04***
Total	30.47	170.55	1,454.63	1,095.67
Gastropoda
CI		6.99E-02	9.99E-04***	4.00E-03***	5.89E-02
CIII	353.64		2.89E-01	1.60E-02*	8.41E-01
CIV	599.16	131.78		9.99E-04***	6.78E-01
CV	765.16	546.76	712.94		9.99E-03**
Total	408.77	19.67	54.12	560.43

Only point with dN/dS > 1 are considered; see text for details. Values below the diagonal of the matrix are the test statistics, whereas the P-values are given above the diagonal (* < 0.05; ** < 0.01; *** < 0.005).

Data from gastropods are shown in supplementary Additional file S10, Supplementary Material online, for details. With the exception of some scattered events of directional selective pressure on Complex V, none of the mitochondrial complexes seem to have experienced marked degrees of directional selective pressure after the origin of the Class, 516 Ma (Orlowski 1985), and never in the Palaeozoic.

Figure 2 depicts the distribution of all time lapses with dN/dS > 1 over time, whereas level of significance of Kolmogorov–Smirnov tests are given for all datasets in table 3. Although the pairwise comparison is not significant for cytb and IV, it is evident 1) that Complex V experienced, much more frequently in bivalves than in gastropods, significantly higher levels of directional selective pressure, and 2) that directional selective pressure had its acme before the increase in atmospheric oxygen for bivalves, and after it for gastropods.

F

Distribution of dN/dS > 1 intervals over time. For each time lapse, all branches exhibiting a dN/dS > 1 were counted, assigned to time bins of 40 Myr starting from 520 Ma, and normalized. The blue line shows atmospheric values of pO2 (Glasspool and Scott 2008) and is plotted on the right axis. From left to right, red vertical lines show the Cambrian–Ordovician (485.4 Ma), Silurian–Devonian (419.2 Ma), and Devonian–Carboniferous (358.9 Ma) boundaries. A black horizontal line marks a dN/dS of 1; time lapses of 7 Myr.

Table 3

Kolmogorov–Smirnov Test Among Classes

Dataset	P-Value
CI	5.87E-03**
CIII	4.31E-01
CIV	1.12E-01
CV	0.00E+00***
Total	2.57E-01

Tested are the differences in the distribution of all time lapses with dN/dS > 1 over the entire history (* < 0.05; ** < 0.01; *** < 0.005).

Amino Acid Signatures on cytb and nad1

The region selected for cytb gene was from site 176 to site 234 (using yeast numbering; Degli Esposti et al. 1993). Site 204 is a relatively conserved Thr (963 out of 1240 metazoans sequences; 77.66%); in bivalves and gastropods the conservation is much lower (48.98 and 46.94%, respectively) and this residue is often substituted with Lys (15 times in bivalves; 50 in gastropods). Site 206 is a highly conserved Ser (97.66 and 98.26% among metazoans and gastropods, respectively); notably, only 89.80% of the bivalve sequences have S206.

The region selected for nad1 gene was from site 19 to site 66 (Thermus thermophilus numbering; Baradaran et al. 2013). Site 42 is a relatively conserved Met (934 out of 1241 metazoan sequences; 75.26%), which in bivalves and gastropods is much more variable (14.43 and 25.22%, respectively). Q43 is highly conserved both in metazoans (96.54%) and in gastropods (98.26%), but in bivalves is much more variable (71.13%). Finally, and similarly, G52 is highly conserved in metazoans (95.17%) and in gastropods (86.09%), but in bivalves is much less conserved (55.67%) and often substitued with a Ser (37 times).

Almost all χ2 tests are highly significant with the only exception of the bivalve/gastropod comparison on M42: raw data, as well as χ2 levels of significance, are listed in supplementary Additional file S11, Supplementary Material online, for details, whereas protein conservation is shown in figure 3 along with basic functional features.

F

Conservation of key residues part of the Q reacting chamber. Alignments of selected regions of the bivalve (a, d), gastropod (b, e) and metazoan (c, f) cytb (a–c) and nad1 (d–f) genes. Key residues are boxed in black, and are T204 and S206 for cytb (yeast numbering; Degli Esposti et al. 1993) and M42, Q43, and G52 for nad1 (Thermus numbering; Baradaran et al. 2013); see text for further detail. Residue conservation is shown above and below each alignment: in the latter case, conservation is coded following a black (very conserved)/white (not conserved) scale. Hydropathy profiles are shown at the bottom of each alignment. Amino acids are shaded following rasmol conventions. Note that all features are computed on complete alignments (with the exception of the metazoan dataset, where only 10% of the sequences were drawn), but only 25 random sequences are shown. a, bivalve cytb gene; b, gastropod cytb gene; c, metazoan cytb gene; d, bivalve nad1 gene; b, gastropod nad1 gene; c, metazoan nad1 gene.

Discussion

The phylum Mollusca is the second phylum in the world in terms of biodiversity, after the Arthropoda (Brusca et al. 2016). Like arthropods, molluscs are spread in most of the biotas of the planet, and were also able to stably colonize freshwater and subaerial environments. Probably the key of the stunning evolutionary success of molluscs is mainly due to the versatility of their original Bauplan, with special reference to foot muscle, shell (at least for conchiferans), mantle cavity, and bipectinate gills (ctenidia).

In marine environments, the three largest classes of molluscs colonized the three main faunal zones: cephalopods are preferentially pelagic, gastropods are mostly epibenthic, and bivalves are commonly found in the endofaunal zone, although the lamellibranch condition and, among other features, the evolution of byssus allowed different group to exploit other ecological niches.

Bivalves and gastropods are strictly related conchiferans, therefore they share a long array of features, but a key difference, which probably was the trigger of most of the other bivalve autapomorphies (e.g. the bivalve shell, the loss of a head, lamellibranchiate gills; see Morton 1996), is the shift to the infaunal zone itself.

It is known that bivalve mitochondria can turn off Complexes II, III, IV in anaerobic conditions, the proton gradient being established thanks to Complex I alone (Müller et al. 2012). This is in agreement with the peculiar dN/dS pattern of bivalves (fig. 2 and see supplementary Additional file S9, Supplementary Material online): during the early infaunalization phase, Complex I probably underwent massive selective pressures to improve its efficiency to adapt to the new environment. Conversely, gastropods did not enter the infaunal zone and coherently do not show the Ordovician Complex I signature of bivalves (see supplementary Additional file S10, Supplementary Material online).

Establishing a proton grandient by means of the Complex I alone probably requested heavy selective pressure on the Complex V, too—a condition which is typical of the complete evolutionary history of bivalves (see supplementary Additional file S9, Supplementary Material online). Again, nothing similar was detected for gastropods (see supplementary Additional file S10, Supplementary Material online). However, it has to be noted that the catalytic subunits of mitochondrial ATP synthase are all coded by nuclear DNA (Baker et al. 2012; Walker 2013; Zhou et al. 2015), therefore we cannot reject the hypothesis that the peculiar pattern of bivalve atp6 is connected to some feature of this single gene.

Directional selective pressure on Complex I seems to become milder in correspondence with the increasing atmospheric levels of oxygen (see supplementary Additional file S9, Supplementary Material online). It is conceivable that the overall higher oxygen partial pressure led to an overall increase in dissolved oxygen, which ended up in a higher oxygen availability within sediments, where bivalves then stably lived and speciated. In this scenario, selective pressure may have been relaxed for bivalve mtDNAs; once more, nothing similar happened to epibenthic gastropods (see supplementary Additional file S10, Supplementary Material online).

Dramatic differences between bivalves and gastropods are evidenced by distribution of time lapses with dN/dS > 1 (fig. 2). For bivalve Complex I they suddenly drop after the lower Devonian oxygenation event, whereas, on the contrary, they appear for gastropods Complex I at the same time. Notably, no significant differences were detected between bivalves and gastropods when the same distributions are computed for cytochrome b and Complex IV (table 3).

Given that no clear link with environmental redox conditions is known for Lower Devonian gastropods, this may find a taxonomic explanation in the radiation of Gastropoda, instead of a palaeochemical/palaeoecological one as in bivalves: indeed, we found most signal of gastropod—especially heterobranch— radiation exactly in the Devonian period (fig. 1 and see supplementary Additional file S7, Supplementary Material online), which is in broad agreement with current knowledge on gastropod evolution (see, f.i., Dinapoli and Klussmann-Kolb 2010; and references therein). The sudden increase in oxygen availability (Glasspool and Scott 2008) may well have triggered the evolutionary burst of crawling, motile metazoans like gastropods.

Further adaptation to low-oxygen environments may have been achieved through more efficient respiratory substrates in bivalves. The in vivo biochemical detection of various quinone molecules—mainly RQ—is well beyond the scopes of the present paper; however, we investigated (fig. 3) all the residues that are part of the Q reacting chamber and were previously found to be correlated with the presence of RQ (Degli Esposti 2015).

The amine group increases the number of potential H bonds that RQ can form with interacting enzymes: the substitution of T204 with Lys in cytochrome b sequence is a change that also increases the number of potential H bonds and was found to be typical of RQ-containing species, including bivalves like Mytilus edulis and Ostrea edulis (Degli Esposti 2015). Although T204 is highly conserved in metazoans, it is much more variable in bivalves, and 15 bivalve mtDNAs showed a Lys instead, confirming those findings. Notably, this tendence is even stronger in the gastropod dataset, where as many as 50 species showed a Lys: thus, we suggest that a low conservation of T204 was an earlier synapomorphy of Conchifera, and may well have been exapted by bivalves to exploit the potential of RQ in low-oxygen conditions.

Similarly, a conserved G52 in the nad1 gene was found to be the rule for all metazoans (including gastropods); however, it is substituted in about half the bivalve sequences; often (37 times) it is replaced with a H-bond forming amino acid, Ser, again confirming previous findings on M. edulis and Crassostrea angulata (Degli Esposti 2015). We also confirmed the previously noted variability of bivalve residue 43 of the same gene (Degli Esposti 2015; fig. 3), which is also part of the Q reacting chamber (Baradaran et al. 2013).

Summarizing, we were able to detect amino acid trends that constitute strong clues of the presence of RQ as a possible substrate for mitochondrial Complex I (fig. 3). Most interestingly, some of these trends appear to be shared at least with gastropods (K204 in the cytb gene; nonconservativity of M42 in the nad1 gene), while other appear to be bivalve-specific (S52 and nonconservativity of Q43 in the nad1 gene). Beside gathering direct evidence of the respiratory substrates exploited by bivalves, we note that it may be very interesting in the future to investigate the same sites in other molluscan classes (like cephalopods, or polyplacophorans), to identify plesiomorphies of molluscs/conchiferans and spot out bivalve autapomorphies related to the adaptation to hypoxia, which were gained through the intense selective pressure that was described above.

Concluding, the major driving force in bivalve history (Ordovician infaunalization) seems to be deeply connected with a period of intense directional selection on mitochondrial Complex I, which, on the other side, looks totally absent from gastropod data. There are other metazoans group that share with bivalves the infaunal life and the ability to turn off Complexes II, III, and IV during anoxic events, like some flatworms, annelids, and peanut worms (Müller et al. 2012). It would be very stimulating to confirm or disprove the same pattern of selection on mtDNA by applying the present pipeline to these diverse metazoan taxa, as well as to analyze whether the same pattern of selection is mirrored or not in the genes for the other subunits of the respiratory chain complexes, which are coded by the nucleus and poorly studied for bivalves.

Supplementary Material

Supplementary data are available at Genome Biology and Evolution online.

Click here for additional data file.