Phylogenomic analyses of deep gastropod relationships reject Orthogastropoda

(a) Taxon sampling, RNA isolation and sequencing

We collected new transcriptome data for 40 species, including 34 gastropods and six other molluscs. All new datasets are paired-end Illumina reads, except for single-end Illumina datasets for Hinea brasiliana, Philine angasi and Strubellia wawrai. Samples were prepared for sequencing with TruSeq RNA Sample Prep Kit (Illumina) or a previously described custom protocol [27]. We deposited all these new sequence data, along with associated specimen collection information, voucher accession numbers, RNA extraction methods and library preparation details, in NCBI Sequence Read Archive (BioProject PRJNA253054). Vouchers for most specimens were deposited at the Museum of Comparative Zoology, Harvard University (Cambridge, MA, USA) and Scripps Institution of Oceanography (La Jolla, CA, USA). The publicly available data for Siphonaria pectinata are here shown in the figures as Siphonaria naufragum, according to a recent revision [28].

(b) Data analyses

These data were analysed in combination with publicly available data for 16 additional species to generate 56-taxon matrices. All Illumina reads (new and publicly available) were assembled with Agalma (v. 0.3.4–0.3.5) [29], 454 datasets were assembled externally with Newbler (v. 2.3 and 2.5p1) and gene predictions from Lottia gigantea [30] and Pinctada fucata [31] were imported directly into Agalma. Source code for most analysis steps as well as sequence alignments, tree sets, summary trees and voucher information are available in a git repository at https://bitbucket.org/caseywdunn/gastropoda.

Two methods were used to generate the supermatrices within Agalma. In method 1, after assembly, translation and removal of mtDNA loci, the sequences from all taxa were compared to each other using an All-By-All BLASTP search, and a phylogenetic approach to identify orthologous sequences [32]. We refer to this method as ABA. In method 2, the sequences from all taxa were compared using OMA v. 0.99t [33] to directly assign sequences to groups of orthologues using an entirely phenetic approach [34]. We refer to this method as OMA.

For each method (ABA and OMA), we constructed three progressively smaller and denser amino acid supermatrices, creating a total of six matrices (figure 2). Supermatrix 1 was constructed by concatenating all orthologue sequences until the cumulative gene occupancy was 50% (49 752 sites/862 loci for ABA and 190 752 sites/1245 loci for OMA; 425 loci in common between ABA and OMA) and then removing Pyropelta sp. and Paralepetopsis sp., which were poorly sampled. Supermatrix 2 was constructed by removing taxa with less than 20% gene occupancy from Supermatrix 1. The removed taxa include Haliotis kamtschatkana, Perotrochus lucaya, Littorina littorea, Siphonaria naufragum, Chaetoderma sp. and Pomacea diffusa for both OMA and ABA matrices, as well as Amphiplica gordensis for the ABA Supermatrix 2. This taxon was removed from the ABA Supermatrix 2 bootstrap replicates and the ABA Supermatrix 2 posterior probability tree sets prior to summary so that they could be consistently displayed (electronic supplementary material, figures S1 and S2). Supermatrix 3 was constructed by trimming genes from Supermatrix 2 until the cumulative gene occupancy reached 70% (15 735 sites/300 loci for ABA and 45 084 sites/364 loci for OMA; 110 loci in common between ABA and OMA).

• Download figure
• Open in new tab
• Download PowerPoint

We inferred phylogenetic relationships using both maximum-likelihood (ML) and Bayesian approaches, for a total of 12 phylogenetic analyses on the six supermatrices. For ML, we used ExaML v. 1.0.11 [35] with a WAG+Γ model of amino acid evolution. Bootstrap values were estimated with 200 replicates. Bayesian analyses were conducted with PhyloBayes MPI v. 1.4e [36] using the CAT-Poisson model of amino acid evolution. Two independent MCMC chains were run on each matrix, adjusting the number of cycles until convergence was achieved. Convergence was determined with time-series plots of the likelihood scores, time-series plots of the cumulative split frequencies, maximum bipartition discrepancies across chains less than 0.1, and an estimated effective sample size of tree likelihoods of at least 100. Post-burn-in sampled trees were combined and summarized with a majority rule consensus tree.

Tree dating was conducted with MCMCTree v. 4.7 [37] using the approximate likelihood calculation algorithm [38], and the WAG+Γ model of evolution. A birth–death speciation process was specified as tree prior with default parameters (death and growth rate parameters equal 1, and sampling parameter equals 0). Rate heterogeneity among lineages was modelled using an uncorrelated lognormal relaxed molecular clock [39] with a diffuse gamma Γ(1,1) prior for the substitution rate and the rate-drift parameter. We used fossil calibrations to set prior densities on the ages of five nodes (see figure 4) using minimum soft bounds with a left tail probability of 2.5% [40]. Because MCMCTree always needs a calibration point on the root [37], we used 550 Ma (ca Terreneuvian; [41]) to set a prior density on the root age using a maximum soft bound with 2.5% tail probability. We ran MCMCTree twice each time for 1.2 × 10⁷ generations, sampling every 1.0 × 10³ and discarding 20% of the samples as burn-in. Convergence was determined with time-series plots of the likelihood scores and assessing for correlation of divergence times between runs.

(c) Hypothesis testing for Orthogastropoda

We statistically compared the Orthogastropoda hypothesis to our maximum-likelihood tree using the SOWH test [42]. To carry out this analysis, we used SOWHAT [43] specifying a constraint tree and the WAG+Γ model on supermatrix 1 (OMA). We used the automatic stopping criterion implemented in SOWHAT to determine an appropriate sample size for the null distribution.

(a) Deep relationships among major gastropod clades

Our datasets strongly support the monophyly of gastropods. This result is not surprising in itself but has only recently been supported by molecular analyses of large datasets [32,44,45] (see also [14,15,17,18]), or in the total evidence analysis of Aktipis et al. [6]. Our analyses also support the monophyly of all major gastropod clades represented by multiple taxa: Vetigastropoda, Neritimorpha, Caenogastropoda and Heterobranchia (figure 3a). Patellogastropoda is represented by a single species, so its monophyly could not be evaluated. The deep internal relationships of gastropods therefore can be reduced to a 5-taxon problem (figures 1 and 3b). Our 12 analyses (two inferences methods on two types of supermatrices each subsampled in three different ways) all recover the same unrooted ingroup relationships for these five clades (figure 3b; electronic supplementary material, figures S1b and S2b). These ingroup relationships are strongly supported by all methods except the ABA ML analyses, which have lower support than the other methods for a bipartition Vetigastropoda + Patellogastropoda (58, 75 and 56% for Supermatrices 1, 2 and 3) and recover Vetigastropoda + Neritimorpha in a minority of replicates. The lower support in these analyses may be due to the poor sampling of Patellogastropoda. These ingroup relationships allow us to reject the hypotheses for gastropod relationships indicated in figure 1a,b.

• Download figure
• Open in new tab
• Download PowerPoint

Although the ingroup relationships found broad consistent support, the rooting of gastropods is still not well resolved. Our results are congruent with three possible rootings (orange arrows in figure 3b; electronic supplementary material, figures S1b and S2b). This is akin to other recalcitrant animal phylogeny questions, including the root of Metazoa [46,47] and the root of arthropods [48]. Though the hypothesis indicated in figure 1c is compatible with the ingroup relationships supported here, we never recover this rooting and it can be excluded. This reduces the possible alternatives for deep gastropod relationships to the three hypotheses (figure 1d–f). Two of these remaining hypotheses have been proposed before [6,14]; the other (figure 1f) is recovered for the first time here.

The rejection of several widely held hypotheses for deep gastropod phylogeny (figure 1a–c) has major implications for the understanding of gastropod evolution. All our analyses reject the Orthogastropoda hypothesis (a clade comprised Vetigastropoda, Neritimorpha, Caenogastropoda and Heterobranchia) and the placement of Patellogastropoda as the sister group to other gastropods (figure 1a,b). Even in the minority of ABA ML replicates that recover an ingroup partition Vetigastropoda + Neritimorpha, the rooting is inconsistent with Orthogastropoda. The broadly accepted Orthogastropoda hypothesis has been proposed in multiple configurations [6,9,10,49,50]. The placement of Patellogastropoda as the sister group to Orthogastropoda has been driven by considerable anatomical research. One potential character supporting this placement is the ciliary ultrastructure of the cephalic tentacles, which also occurs in Bivalvia and Solenogastres but is lacking from other gastropods [51]. In this scenario, this character is plesiomorphic for Mollusca, retained in Patellogastropoda and was lost a single time in Orthogastropoda. However, because enforcing the monophyly of Orthogastropoda is significantly worse (SOWH test: n = 152, Δ-likelihood = 374.0137, p = 0) than our most likely tree (figure 3), our results indicate that this character may be convergent between Patellogastropoda and outgroup taxa, or was lost more than once within Gastropoda. We also reject another recent hypothesis for gastropod rooting, the sister group relationship of Neritimorpha to other Gastropoda (figure 1c) [15,17,18].

Our reduction of deep gastropod phylogeny to three alternative hypotheses (figure 1d–f) clarifies multiple open questions. All three of these hypotheses include the monophyly of Apogastropoda (i.e. Heterobranchia + Caenogastropoda), reinforcing this widely accepted aspect of gastropod relationships. Other relationships supported here have been found earlier (figure 1d,e; see review by Aktipis et al. [6]; for different views, see e.g. the gastropod classification by Bouchet et al. [52] and the mitogenomic study by Grande et al. [53]; see also [6,14]). However, it is now acknowledged that mitogenomic data are not appropriate for resolving deep gastropod relationships [54]. To the best of our knowledge, no molecular analysis supported the placement of Vetigastropoda as sister group to all other gastropods.

Because we ran 12 phylogenetic analyses, we can explore the differences in support between these three alternative hypotheses for gastropod rooting across inference method (Bayes and ML), matrix construction method (OMA and ABA) and matrix subsampling (Supermatrices 1, 2 and 3; figure 2). Matrix subsampling had little effect on deep relationships. Analyses of Supermatrix 1 (figure 3) and Supermatrix 2 (electronic supplementary material, figure S1) were consistent with all three rooting positions (figure 1d–f). Analyses of Supermatrix 3 (electronic supplementary material, figure S2) found support for only two of these rootings (figure 1d,f). Unlike the other analyses, it did not recover Apogastropoda as the sister group to all other gastropods (figure 1e). This particular hypothesis (figure 1e) is interesting because it includes Archaeogastropoda, which was proposed nearly a century ago by Thiele [19]. Bayesian analyses recovered Neritimorpha as the sister group to Apogastropoda (figure 1d,f) in all analyses, but ML analyses found very low support for this relationship. Analyses of OMA matrices provided strong support for a clade comprising Patellogastropoda and Vetigastropoda (figure 1d,e), but analyses of ABA matrices did not.

These results suggest clear strategies for distinguishing between the remaining hypotheses for deep gastropod relationships. Since these hypotheses differ only in their rooting, improved outgroup sampling will be critical. To maximize gene sampling and matrix density, we limited our sampling of non-gastropod molluscs to those for which Illumina data [32] or genomes are available. Future analyses of gastropod relationships will need to include more outgroups to resolve the remaining open questions. Previous phylogenomic analyses of molluscs that also included extensive 454 and Sanger data [32,44] had much broader non-gastropod sampling but minimal gastropod sampling. The rooting of gastropods was not fully supported in these analyses either, but the strongly supported ingroup relationships are compatible with the three hypotheses supported here. In addition, improved sampling of Patellogastropoda (here represented by a single species with a complete genome) and Neritimorpha, and the addition of the unsampled Neomphalina and Cocculiniformia, will be critical.

(b) Relationships within major gastropod clades

Although our sampling is focused on resolving deep relationships between the major gastropod clades, our results do find strong support for some previously unresolved relationships within Vetigastropoda, Caenogastropoda and Heterobranchia.

A key question within Vetigastropoda is the placement of Pleurotomarioidea (Perotrochus in our analyses), which appears outside Vetigastropoda in some previous studies [14,15]. Here we find strong support for the placement of Pleurotomarioidea as the sister group to all other vetigastropods (figure 3a), resolving this issue. We also resolve the position of Seguenzioidea (Granata) as the next clade to diverge from other vetigastropods (figure 3a; electronic supplementary material, figures S1 and S2). Our analyses also recover a well-supported clade of deep-sea taxa (Pseudococculinidae (Amphiplica) and Lepetodriloidea (Lepetodrilus); figure 3a). We also find strong support (figure 3a; electronic supplementary material, figures S1 and S2) for a clade comprising Phasianellidae (Phasianella) and Trochoidea (Prothalotia and Monodonta). The position of Haliotis is not resolved (figure 3a), perhaps due to relatively poor gene sampling.

Caenogastropoda is a megadiverse clade comprising about 60% of living gastropod species [55], so our limited sampling can address only a small fraction of open questions about internal relationships of this group. The relationships we can test are largely in agreement with prior morphological [10,55–57] and molecular [58–60] analyses. We find a sister group relationship of Ampullarioidea (represented by Pomacea) to Sorbeoconcha, which comprise the remaining sampled caenogastropods (figure 3a). Within Sorbeoconcha, Cerithioidea (Hinea) is the sister group to Hypsogastropoda, the latter dividing into a siphonate (in our case the two Neogastropoda: Urosalpinx and Ilyanassa) and an asiphonate group (Janthina, Littorina, Euspira plus Rubyspira), similar to the findings of Ponder et al. [55].

The basic structure of internal Heterobranchia relationships has only recently gained some agreement [23,24]. Our strong support for the placement of Architectonica as sister group to the other sampled heterobranchs is consistent with most other analyses [24]. Nudibranchia (Fiona + Bathydoris + Doris) and Nudipleura (Pleurobranchaea + Nudibranchia) were monophyletic [61,62], despite some suggestion that Pleurobranchoidea may not be the sister group to the remaining nudibranch lineages [26]. Our results recovered a monophyletic Cephalaspidea (Philine + Haminoea), sister group to Anaspidea (Aplysia) + Pteropoda (Clione). Umbraculoidea (Tylodina) was sister group to Cephalaspidea + Anaspidea + Pteropoda; all four taxa together represent the well-supported Euopisthobranchia (figure 3a; electronic supplementary material, figures S1 and S2). We find support for Panpulmonata (figure 3a; electronic supplementary material, figures S1 and S2), but their internal relationships are mostly unresolved and clearly require future attention. Like previous Sanger sequencing-based studies, our analyses consistently recover a Panpulmonata + Euopisthobranchia clade or Tectipleura [25,62]. The relationship of Tectipleura to other heterobranchs has been of particular interest. We recover two conflicting hypotheses for these relationships, neither of which has been previously proposed. Our likelihood analyses place the unnamed clade Rissoelloidea (Rissoella) + Acteonoidea (Hydatina) as the sister group to Tectipleura (figure 3a; electronic supplementary material, figures S1 and S2). Our Bayesian analyses, however, place Nudipleura and this Rissoelloidea + Acteonoidea clade together with strong support, and place this clade as sister to Tectipleura (electronic supplementary material, figure S3). Previous analyses have instead favoured Euthyneura, a clade comprising Tectipleura and Nudipleura (but excluding Rissoelloidea + Acteonoidea). We do not recover Euthyneura in any of our analyses. Tectipleura is united by a monaulic reproductive system [23], but even Euthyneura is not entirely defined by euthyneury, as there are internal reversals [63] and the condition is known in other groups (Rhodopemorpha [64]). Rissoelloidea + Acteonoidea and Euthyneura share giant neurons in macroscopic animals [63], and if necessary, a simple redefinition of the taxon Euthyneura to include Rissoelloidea + Acteonoidea would maintain stability.

(c) Chronogram

Our dated phylogeny (figure 4) shows a Cambrian origin of stem gastropods with crown diversification into its five main lineages during the Ordovician to the Devonian, as well shown in the gastropod fossil record [72]. From the well-sampled groups, crown Vetigastropoda diversified first, around the Devonian–Carboniferous, followed by Neritimorpha and Heterobranchia at similar periods. Crown Caenogastropoda seem to have diversified later, around the Permian–Triassic, perhaps initiating its explosive diversification after the end-Permian mass extinction ca 254 Ma, responsible for the extinction of 95–99% of marine species and to change the ecosystems and their faunal composition forever. Such drastic post-extinction diversifications have been recently shown for other modern clades of marine organisms (e.g. Crinoidea [73] and Protobranchia [74]). This could also explain other explosive radiations in gastropods, especially within the euthyneuran Heterobranchia clades such as Nudipleura (figure 4). However, denser sampling will be required to derive accurate diversification curves to test these hypotheses.

• Download figure
• Open in new tab
• Download PowerPoint