TY - JOUR AB - Abstract The sister taxon of one of the most diverse and pharmacologically important heterobranch gastropod groups, Nudibranchia, is unclear. Pleurobranchomorpha have been recovered sister to Nudibranchia in most molecular phylogenetic studies to date, forming a clade called Nudipleura. However, this relationship has not been consistently recovered; some phylogenomic studies have recovered Nudipleura as a paraphyletic group, but transcriptome data have been available from only one representative of Pleurobranchomorpha to date. Here, we supplemented the paucity of data available for Pleurobranchomorpha by sequencing the transcriptome of Bathyberthella antarctica Willan & Bertsch, 1987 and analysing transcriptomes and genomes from a total of 27 gastropod species, inferring the phylogeny of Heterobranchia based on 993 genes. Consistent with analyses of datasets dominated by 18 S rDNA and recent phylogenomic studies, all of our analyses unequivocally support Nudipleura, recovering Nudibranchia as sister to Pleurobranchomorpha with maximal support. Our study supports previous hypotheses that nudipleurans evolved from a common ancestor with an androdiaulic reproductive system, a blood gland and a secondarily lost osphradium. INTRODUCTION Heterobranchia are a diverse and important clade of marine, freshwater and terrestrial snails and slugs (Wägele et al., 2014). Some species are model organisms in neurobiology (e.g. Moroz et al., 2006; Senatore, Edirisinghe & Katz, 2015), while others produce or sequester chemicals that have medical applications (e.g. Newman & Cragg, 2014) and many terrestrial snails and slugs are of great conservation concern (e.g. Schilthuizen et al., 2005). Despite their diversity and importance, many questions about the evolutionary relationships among these fascinating and important molluscs remain. Numerous phylogenetic analyses based on morphology and molecular data have attempted to resolve heterobranch phylogeny, but this has proved challenging. Phylogenetic studies on Heterobranchia based on morphological data are made difficult by the parallel evolution of numerous organ systems and other structures in different groups, resulting in homoplasies in the analysed datasets (reviewed by Wägele et al., 2014). Analyses of mitochondrial protein-coding genes have recovered Opisthobranchia, a traditionally recognized group of sea slugs that includes groups such as the familiar and charismatic nudibranchs (Grande, Templado & Zardoya, 2008; Medina et al., 2011; White et al., 2011), as monophyletic, but with Pulmonata, the group that includes most terrestrial snails and slugs, as paraphyletic. In contrast, analyses of nuclear and mitochondrial rRNA genes plus the cytochrome c oxidase I (COI) gene have consistently recovered Opisthobranchia as paraphyletic with respect to a monophyletic Pulmonata (e.g. Klussmann-Kolb et al., 2008; Dinapoli & Klussmann-Kolb, 2010; Jörger et al., 2010; Dayrat et al., 2011, Wilson et al., 2017). These conflicting results have led some to view Opisthobranchia as a valid taxon (Medina et al., 2011; White et al., 2011), while others have explicitly said “bye-bye” to Opisthobranchia (Jörger et al., 2010; Schrödl et al., 2011). More recently, phylogenomics—an approach to the study of evolutionary history using dozens to hundreds of genes derived from genomes or transcriptomes (Telford, 2008)—has been applied to heterobranch phylogeny. Kocot, Halanych, and Krug (2013a) conducted a phylogenomic analysis and recovered Opisthobranchia as paraphyletic with respect to Panpulmonata, a clade in which Sacoglossa was recovered as sister to Pulmonata. Interestingly, Nudipleura, a hypothesized grouping in which Nudibranchia are sister to Pleurobranchomorpha (reviewed by Wägele et al., 2014), was recovered as paraphyletic with strong support. Notably, this result was also observed in some analyses in an earlier phylogenomic analysis that focused on class-level molluscan phylogeny (Kocot et al., 2011). However, most studies addressing heterobranch phylogeny have recovered Nudipleura as monophyletic (Wägele & Willan, 2000; Göbbeler & Klussmann-Kolb, 2010; Dinapoli & Klussmann-Kolb, 2010; Jörger et al., 2010; Kocot et al., 2011 [in part]; Schrödl et al., 2011; Zapata et al., 2014; Oskars, Bouchet & Malaquias, 2015; Kano et al., 2016; Wilson et al., 2017). Zapata et al. (2014) conducted the most comprehensive phylogenomic analysis of Gastropoda to date. In addition to shedding light on some long-standing questions on higher-level gastropod phylogeny and providing further evidence for the paraphyly of Opisthobranchia with respect to Pulmonata, they recovered Nudipleura as monophyletic with maximal support. This result is consistent with most previous molecular studies, but inconsistent with other phylogenomic analyses (Kocot et al., 2011 [in part], 2013a). Despite the diversity, importance and popularity of the Nudibranchia, the sister taxon of this clade remains uncertain. Phylogenomic studies to date may have recovered differing relationships because of the limited amount of data available for Pleurobranchomorpha, a group sometimes called the side-gilled slugs. This clade of sea slugs is characterized by a reduced internal shell, a protruding mantle and, in at least some species, secretion of defensive acidic mucus (Göbbeler & Klussmann-Kolb, 2010; Wägele, Knezevic & Moustafa, 2017). Pleurobranchomorpha are of great interest because they have been recovered as the sister taxon of Nudibranchia in most studies so far, but also because a species of this group is an emerging model organism in the field of neurobiology (Gillette & Brown, 2015). Transcriptome data have only been available from one species, hindering understanding of the phylogenetic position of this clade. Here, we re-evaluated the Nudipleura hypothesis in light of publicly available genome and transcriptome data, with improved sampling for Pleurobranchomorpha using new transcriptome data from the Antarctic pleurobranch Bathyberthella antarcticaWillan & Bertsch, 1987 and a molecular marker selection strategy that includes careful manual screening to exclude exogenous contamination. MATERIAL AND METHODS A single specimen of Bathyberthella antarctica was collected by Agassiz trawl during R/V Polarstern cruise PS96 in the Weddell Sea (74° 56.49′S, 32° 27.27′W at 621.3 m depth) by K. Kocot on 4 January 2016. Foot tissue was stored in RNAlater and the rest of the specimen was preserved in 10% formalin overnight and stored in 75% ethanol. This specimen is deposited in the Alabama Museum of Natural History Invertebrate Zoology collection (acc. no. 21268). Total RNA was extracted from foot tissue using the E.Z.N.A. Mollusc RNA kit (OMEGA Bio-Tek) and 1 μg was sent to Macrogen Inc. (Seoul, South Korea) for Illumina TruSeq RNA library preparation and sequencing on an Illumina HiSeq 2500 with 2 × 101 bp paired-end reads. Raw data are available for download from NCBI SRA (acc. no. SRR5725019) and the assembled transcriptome is available via NCBI TSA (BioProject no. PRJNA391256). Publicly available gene models and transcriptomes were downloaded for 26 gastropods (Table 1). For unassembled transcriptomes, raw reads were assembled using Trinity v. 2.2.0 (Grabherr et al., 2011) with the ‘trimmomatic’ and ‘normalize_reads’ flags. Paired-end and single-end libraries were treated as appropriate (–left and –right or –single, respectively). We then used Transdecoder v. 2.1.0 (https://sourceforge.net/p/transdecoder/) to translate nucleotide sequences to amino acids (see Supplementary Material for example script with commands used). Table 1. Taxon sampling. Species Clade Abbreviation Genes Source Acc. nos Aplysia californica Anaspidea ACAL 1796 NCBI SRA SRR016568 Architectonica perspectiva Architectonicoidea APER 1540 NCBI SRA SRR1505102 Bathyberthella antarctica Pleurobranchomorpha BATH 1055 NCBI SRA SRR5725019 Bathydoris clavigera Nudibranchia BCLA 2049 NCBI SRA SRR1505104 Clio pyramidata Thecosomata CPYR 1232 NCBI SRA SRR1300146 Clione antarctica Pteropoda CANT 1946 NCBI SRA SRR1505107 Conus miliaris Caenogastropoda CMIL 2079 NCBI SRA PRJNA257931 Crepidula fornicata Caenogastropoda CFOR 1888 NCBI SRA PRJNA285294 Doris kerguelenensis Nudibranchia DKER 2040 NCBI SRA SRR1505108 Elysia timida Sacoglossa ETIM 2192 NCBI SRA SRR1582570 Haminoea antillarum Cephalaspidea HANT 2168 NCBI SRA SRR1505111 Hermissenda crassicornis Nudibranchia HCRA 1921 NCBI SRA SRR1950939 Hydatina physis Cephalaspidea HPHY 2107 NCBI SRA SRR1505113 Lottia gigantea Patellogastropoda LGIG 2257 InParanoid See Kocot et al. (2017) Melibe leonina Nudibranchia MLEO 1957 NCBI SRA SRR1950947 Nucella lapillus Caenogastropoda NLAP 1680 NCBI SRA PRJNA217409 Onchidella floridana Systellommatophora OFLO 2072 NCBI SRA SRR1505123 Oxynoe viridis Sacoglossa OVIR 1439 NCBI SRA SRR1505125 Phallomedusa solida Amphiboloidea PSOL 2115 NCBI SRA SRR1505127 Philine angasi Cephalaspidea PANG 2121 NCBI SRA SRR1505129 Pleurobranchaea californica Pleurobranchomorpha PCAL 2181 NCBI SRA SRR1505130 Pomacea canaliculata Caenogastropoda PCAN 2167 NCBI SRA SRR125493 Rissoella caribaea Rissoelloidea RCAR 2212 NCBI SRA SRR1505135 Strubellia wawrai Acochlidoidea SWAW 2163 NCBI SRA SRR1505137 Tritonia festiva Nudibranchia TFES 1430 NCBI SRA SRR1950941 Turbonilla sp. Pyramidellidae TUSP 1972 NCBI SRA SRR1505139 Tylodina fungina Umbraculida TFUN 1787 NCBI SRA SRR1505140 Species Clade Abbreviation Genes Source Acc. nos Aplysia californica Anaspidea ACAL 1796 NCBI SRA SRR016568 Architectonica perspectiva Architectonicoidea APER 1540 NCBI SRA SRR1505102 Bathyberthella antarctica Pleurobranchomorpha BATH 1055 NCBI SRA SRR5725019 Bathydoris clavigera Nudibranchia BCLA 2049 NCBI SRA SRR1505104 Clio pyramidata Thecosomata CPYR 1232 NCBI SRA SRR1300146 Clione antarctica Pteropoda CANT 1946 NCBI SRA SRR1505107 Conus miliaris Caenogastropoda CMIL 2079 NCBI SRA PRJNA257931 Crepidula fornicata Caenogastropoda CFOR 1888 NCBI SRA PRJNA285294 Doris kerguelenensis Nudibranchia DKER 2040 NCBI SRA SRR1505108 Elysia timida Sacoglossa ETIM 2192 NCBI SRA SRR1582570 Haminoea antillarum Cephalaspidea HANT 2168 NCBI SRA SRR1505111 Hermissenda crassicornis Nudibranchia HCRA 1921 NCBI SRA SRR1950939 Hydatina physis Cephalaspidea HPHY 2107 NCBI SRA SRR1505113 Lottia gigantea Patellogastropoda LGIG 2257 InParanoid See Kocot et al. (2017) Melibe leonina Nudibranchia MLEO 1957 NCBI SRA SRR1950947 Nucella lapillus Caenogastropoda NLAP 1680 NCBI SRA PRJNA217409 Onchidella floridana Systellommatophora OFLO 2072 NCBI SRA SRR1505123 Oxynoe viridis Sacoglossa OVIR 1439 NCBI SRA SRR1505125 Phallomedusa solida Amphiboloidea PSOL 2115 NCBI SRA SRR1505127 Philine angasi Cephalaspidea PANG 2121 NCBI SRA SRR1505129 Pleurobranchaea californica Pleurobranchomorpha PCAL 2181 NCBI SRA SRR1505130 Pomacea canaliculata Caenogastropoda PCAN 2167 NCBI SRA SRR125493 Rissoella caribaea Rissoelloidea RCAR 2212 NCBI SRA SRR1505135 Strubellia wawrai Acochlidoidea SWAW 2163 NCBI SRA SRR1505137 Tritonia festiva Nudibranchia TFES 1430 NCBI SRA SRR1950941 Turbonilla sp. Pyramidellidae TUSP 1972 NCBI SRA SRR1505139 Tylodina fungina Umbraculida TFUN 1787 NCBI SRA SRR1505140 ‘Genes’ refers to the number of orthologous groups analysed. When BioProject accession numbers are provided, all runs for the project were combined in one assembly. Table 1. Taxon sampling. Species Clade Abbreviation Genes Source Acc. nos Aplysia californica Anaspidea ACAL 1796 NCBI SRA SRR016568 Architectonica perspectiva Architectonicoidea APER 1540 NCBI SRA SRR1505102 Bathyberthella antarctica Pleurobranchomorpha BATH 1055 NCBI SRA SRR5725019 Bathydoris clavigera Nudibranchia BCLA 2049 NCBI SRA SRR1505104 Clio pyramidata Thecosomata CPYR 1232 NCBI SRA SRR1300146 Clione antarctica Pteropoda CANT 1946 NCBI SRA SRR1505107 Conus miliaris Caenogastropoda CMIL 2079 NCBI SRA PRJNA257931 Crepidula fornicata Caenogastropoda CFOR 1888 NCBI SRA PRJNA285294 Doris kerguelenensis Nudibranchia DKER 2040 NCBI SRA SRR1505108 Elysia timida Sacoglossa ETIM 2192 NCBI SRA SRR1582570 Haminoea antillarum Cephalaspidea HANT 2168 NCBI SRA SRR1505111 Hermissenda crassicornis Nudibranchia HCRA 1921 NCBI SRA SRR1950939 Hydatina physis Cephalaspidea HPHY 2107 NCBI SRA SRR1505113 Lottia gigantea Patellogastropoda LGIG 2257 InParanoid See Kocot et al. (2017) Melibe leonina Nudibranchia MLEO 1957 NCBI SRA SRR1950947 Nucella lapillus Caenogastropoda NLAP 1680 NCBI SRA PRJNA217409 Onchidella floridana Systellommatophora OFLO 2072 NCBI SRA SRR1505123 Oxynoe viridis Sacoglossa OVIR 1439 NCBI SRA SRR1505125 Phallomedusa solida Amphiboloidea PSOL 2115 NCBI SRA SRR1505127 Philine angasi Cephalaspidea PANG 2121 NCBI SRA SRR1505129 Pleurobranchaea californica Pleurobranchomorpha PCAL 2181 NCBI SRA SRR1505130 Pomacea canaliculata Caenogastropoda PCAN 2167 NCBI SRA SRR125493 Rissoella caribaea Rissoelloidea RCAR 2212 NCBI SRA SRR1505135 Strubellia wawrai Acochlidoidea SWAW 2163 NCBI SRA SRR1505137 Tritonia festiva Nudibranchia TFES 1430 NCBI SRA SRR1950941 Turbonilla sp. Pyramidellidae TUSP 1972 NCBI SRA SRR1505139 Tylodina fungina Umbraculida TFUN 1787 NCBI SRA SRR1505140 Species Clade Abbreviation Genes Source Acc. nos Aplysia californica Anaspidea ACAL 1796 NCBI SRA SRR016568 Architectonica perspectiva Architectonicoidea APER 1540 NCBI SRA SRR1505102 Bathyberthella antarctica Pleurobranchomorpha BATH 1055 NCBI SRA SRR5725019 Bathydoris clavigera Nudibranchia BCLA 2049 NCBI SRA SRR1505104 Clio pyramidata Thecosomata CPYR 1232 NCBI SRA SRR1300146 Clione antarctica Pteropoda CANT 1946 NCBI SRA SRR1505107 Conus miliaris Caenogastropoda CMIL 2079 NCBI SRA PRJNA257931 Crepidula fornicata Caenogastropoda CFOR 1888 NCBI SRA PRJNA285294 Doris kerguelenensis Nudibranchia DKER 2040 NCBI SRA SRR1505108 Elysia timida Sacoglossa ETIM 2192 NCBI SRA SRR1582570 Haminoea antillarum Cephalaspidea HANT 2168 NCBI SRA SRR1505111 Hermissenda crassicornis Nudibranchia HCRA 1921 NCBI SRA SRR1950939 Hydatina physis Cephalaspidea HPHY 2107 NCBI SRA SRR1505113 Lottia gigantea Patellogastropoda LGIG 2257 InParanoid See Kocot et al. (2017) Melibe leonina Nudibranchia MLEO 1957 NCBI SRA SRR1950947 Nucella lapillus Caenogastropoda NLAP 1680 NCBI SRA PRJNA217409 Onchidella floridana Systellommatophora OFLO 2072 NCBI SRA SRR1505123 Oxynoe viridis Sacoglossa OVIR 1439 NCBI SRA SRR1505125 Phallomedusa solida Amphiboloidea PSOL 2115 NCBI SRA SRR1505127 Philine angasi Cephalaspidea PANG 2121 NCBI SRA SRR1505129 Pleurobranchaea californica Pleurobranchomorpha PCAL 2181 NCBI SRA SRR1505130 Pomacea canaliculata Caenogastropoda PCAN 2167 NCBI SRA SRR125493 Rissoella caribaea Rissoelloidea RCAR 2212 NCBI SRA SRR1505135 Strubellia wawrai Acochlidoidea SWAW 2163 NCBI SRA SRR1505137 Tritonia festiva Nudibranchia TFES 1430 NCBI SRA SRR1950941 Turbonilla sp. Pyramidellidae TUSP 1972 NCBI SRA SRR1505139 Tylodina fungina Umbraculida TFUN 1787 NCBI SRA SRR1505140 ‘Genes’ refers to the number of orthologous groups analysed. When BioProject accession numbers are provided, all runs for the project were combined in one assembly. Orthologous sequences were identified among the translated genes (which we simply call ‘genes’ henceforth for brevity) using HaMStR v. 13.2.6 (Ebersberger, Strauss & Von Haeseler, 2009; see HaMStR directory in Supplementary Material for example commands used and input/output files). We used the ‘Trochozoa’ core orthology set (2,259 orthologous groups; Kocot et al., 2017) and a modified version of the bioinformatic pipeline used by Kocot et al. (2017; see orthology_script.sh in Supplementary Material). This pipeline first removes contigs that are too short (<50 bp) and groups with few taxa (<15 taxa). We then removed redundant identical sequences using the programme UniqHaplo (http://raven.iab.alaska.edu/~ntakebay/) and aligned sequences with Mafft v. 7.273 (Katoh et al., 2005). Subsequently, we identified ambiguously aligned positions in alignments with ZORRO (Wu, Chatterji & Eisen, 2012) and positions with a ZORRO score less than 0.5 were deleted. Next, a consensus sequence was inferred for each alignment using the EMBOSS programme infoalign (Rice, Longden & Bleasby, 2000). For each sequence in each single-gene amino acid alignment, the percentage of positions of that sequence that differed from the consensus of the alignment were calculated using Infoalign’s ‘change’ calculation. Any sequence that deviated from the consensus sequence by more than 75% was deleted. Subsequently, a custom script was used to delete any mistranslated sequence regions of 20 or fewer amino acids in length surrounded by ten or more gaps on either side. This step was important, as sequence ends were occasionally mistranslated or misaligned. Alignment columns with fewer than four non-gap characters were subsequently deleted. Alignments shorter than 50 bp were excluded and sequences that overlapped by fewer than 20 amino acids with other sequences were removed from the remaining alignments using AlignmentCompare (github.com/kmkocot). At this stage, groups with fewer than 20 taxa were removed. Single gene trees were then inferred using FastTree v. 2.1.7 (Price, Dehal & Arkin, 2010) with the ‘–slow’ and ‘–gamma’ options, and PhyloTreePruner (Kocot et al., 2013b) was used to help remove undetected paralogs and contamination, keeping only the longest ortholog for each taxon. Nodes with support values below 0.95 were collapsed by PhyloTreePruner prior to application of the paralogy pruning algorithm. To further screen for and reduce contamination, post-PhyloTreePruner alignments were visualized using MEGA v. 7 (Kumar, Stecher & Tamura, 2016) and sequences appearing as highly divergent were searched against the NCBI nr database using BLASTP (Altschul et al., 1990) on the NCBI web interface in March 2017. Any sequences with a hit to a non-mollusc in the top 3 hits sorted by e-value were eliminated (see Supplementary Material). Sequences that had top hits to mollusc sequences, but were very divergent from all other sequences in an alignment, were also excluded due to suspicion that the ‘molluscan’ sequence in the database also represented exogenous contamination. Finally, the remaining sequences were concatenated with FASconCAT-G v. 1.02 (Kück & Meusemann, 2010). Although this manual evaluation approach may seem subjective, it was generally very easy to tell by eye when a contaminant sequence was present. Using the resulting matrix, we inferred heterobranch phylogeny by reconstructing a maximum-likelihood (ML) tree with raxmlHPC-PTHREADS-AVX-8.2.4 (Stamatakis, 2014) using the best-fitting model for each partition with a GAMMA model of rate heterogeneity with ML estimation of the alpha-parameter (–m PROTGAMMAAUTO). Nonparametric bootstrapping was conducted using the optimal number of replicates as determined using the RAxML autoMRE criterion (100 replicates for the complete dataset and 50 for the MARE-reduced dataset; see below). We also assembled a reduced dataset in which only the most phylogenetically informative genes were retained, according to the weighted geometry quartet mapping approach implemented in MARE v. 0.1.2-rc (Misof et al., 2013). The default options were used except that all taxa were constrained to be kept with the ‘–c’ option. The resulting data matrix was analysed in RAxML as described above. We additionally explored the consistency of our data using a supertree approach in ASTRAL v. 4.10.12 (Mirarab et al., 2014) with 100 bootstrap replicates. This analysis was based on RAxML single-gene trees reconstructed as described above for the complete dataset, except that these single-gene alignments were not partitioned. All alignments, gene partition files, commands used and resulting trees are available in Supplementary Material. RESULTS AND DISCUSSION Our bioinformatic pipeline resulted in a final concatenated data matrix (‘complete dataset’) of 993 genes that was 163,616 amino acids in length and had an overall completeness of 76% (i.e. 24% missing data). The percentage of missing sequences per taxon ranged from 0.6% (Lottia gigantea) to 55% (Bathyberthella antarctica), with a median percentage per taxon of 13.6% and a mean of 18.7%. Unfortunately, with 1,055 genes, B. antarctica is our ‘worst’ taxon in terms of gene sampling, but we note that all taxa are sampled for roughly half or more of the genes and overall matrix occupancy is high, suggesting that missing data are unlikely to be a problem with this dataset (Roure, Baurain & Philippe, 2012). Out of 24 total nodes in the ML tree inferred from the complete dataset, only three had bootstrap values below 100 (Fig. 1). Results of the analysis of the MARE-reduced dataset and the supertree approach implemented in ASTRAL differed only slightly from the analysis of the complete dataset, and thus we focus on the analysis of the complete dataset here (but see below for discussion of differences in results among analyses). Figure 1. View largeDownload slide Maximum-likelihood tree inferred based on 993 genes. Only bootstrap support values less than 100 are shown. Scale bar indicates 0.05 substitutions per site. Figure 1. View largeDownload slide Maximum-likelihood tree inferred based on 993 genes. Only bootstrap support values less than 100 are shown. Scale bar indicates 0.05 substitutions per site. In the analysis of all 993 genes, we recovered Nudipleura as monophyletic with maximal support, confirming this widely supported hypothesis (e.g. Wägele & Willan, 2000; Göbbeler & Klussmann-Kolb, 2010; Dinapoli & Klussmann-Kolb, 2010; Jörger et al., 2010; Zapata et al., 2014; Kano et al., 2016; reviews by Schrödl et al., 2011; Wägele et al., 2014; Oskars et al., 2015; Wilson et al., 2017). The clade Nudipleura has been supported in most studies to date, but this hypothesis was called into question by two phylogenomic analyses (Kocot et al., 2011, 2013a) that recovered the group as paraphyletic. However, in these studies the sampling for nudipleurans (and Heterobranchia in general) was rather limited. In the study by Kocot et al. (2013a), Pleurobranchomorpha were only represented by Pleurobranchaea californica, and Nudibranchia was only represented by Tritonia diomedea and Hermissenda crassicornis, and these were relatively low-coverage 454 transcriptomes. Our strongly supported results are based on increased taxon sampling and a much more complete dataset in terms of matrix occupancy, and consist entirely of deeply sequenced (e.g. Illumina) transcriptomes and genomes. Our results suggest that recovery of Nudipleura as paraphyletic in earlier phylogenomic studies may have been due to limited taxon sampling, which is well-known to be problematic in phylogenetic analyses (review by Wiens & Tiu, 2012), and/or limited phylogenetic signal in the relatively small (18,101 amino acid positions) and sparse (61% matrix occupancy) data matrix of Kocot et al. (2013a). Morphological characters shared by nudibranchs and pleurobranchs supporting this relationship include an androdiaulic reproductive system, the presence of a blood gland and secondary loss of the osphradium (Wägele & Willan, 2000). Within Nudipleura, our analysis of the complete dataset recovered both Pleurobranchomorpha and Nudibranchia as monophyletic with maximal support. Consistent with recent phylogenomic analyses focused on evolutionary history of cladobranch nudibranchs (Goodheart et al., 2015, 2017), our results indicate that H. crassicornis (Facelinidae) has a more recent common ancestor with Tritonia festiva (Tritoniidae) than Melibe leonina (Tethyidae), although the Hermissenda + Tritonia clade was only moderately well-supported (bootstrap BS = 74) in the analysis of the complete dataset. This result is at odds with morphology-based hypotheses of nudibranch phylogeny, which place Tritonia and Melibe in Dendronotoidea and Hermissenda in Aeolidoidea (review by Wägele & Willan, 2000 and references therein), but our analyses and other recent phylogenomic studies (Goodheart et al., 2015, 2017) found strong support for this relationship. We recovered Anthobranchia (Bathydoris clavigera + Doris kerguelenensis) in agreement with Wägele & Willan (2000), but in contrast to Goodheart et al. (2015, 2017), who recovered this hypothesized grouping as paraphyletic. However, all phylogenomic analyses to date have included just two representatives of this putative clade and these differing results may be due to limited taxon sampling. Considering higher-level heterobranch phylogeny, our ML analysis of the complete dataset was generally consistent with the results of Zapata et al. (2014). We recovered Architectonicoidea (Architectonica) as sister to all other sampled heterobranchs, consistent with previous studies (e.g. Haszprunar, 1985; Dinapoli & Klussmann-Kolb, 2010; Zapata et al., 2014). Nudipleura were recovered in a group that also included Acteonacea (Acteonoidea + Rissoelloidea) and Tectipleura (Euopisthobranchia + Panpulmonata; Schrödl et al., 2011), but the node placing Acteonacea sister to Tectipleura received no support (BS = 48). Notably, placement of Acteonacea was also weakly supported in some of the analyses by Zapata et al. (2014). Within Tectipleura, Euopisthobranchia [represented by Umbraculoidea (Tylodina), Cephalaspidea (Philine and Haminoea), Anaspidea (Aplysia) and Pteropoda (Clione and Clio)] were recovered with maximal support in agreement with previous studies (Jörger et al., 2010; Schrödl et al., 2011; Zapata et al., 2014; Kano et al., 2016; review by Wägele et al., 2014). We recovered Panpulmonata, represented here by Sacoglossa (Elysia and Oxynoe), Acochlidia (Strubellia), Onchidiidae (Onchidella), Amphiboloidea (Phallomedusa) and Pyramidellidae (Turbonilla), with maximal support in congruence with previous studies (review by Schrödl et al., 2011; Schrödl, 2014), but we had limited taxon sampling for this group, whose internal relationships are beyond the scope of the present study. Our ML analysis of the MARE-reduced matrix (Fig. 2A) and the ASTRAL analysis (Fig. 2B) recovered the same higher-level branching order as our analysis of the complete dataset. The only differences between the ML analysis of the MARE-reduced matrix and the analysis of the complete dataset concerned internal relationships within Panpulmonata, which received weaker support in the analysis of the MARE-reduced matrix. Interestingly, support for Acteonacea as sister to Tectipleura was much stronger in the analysis of the MARE-reduced dataset (BS = 91). Relationships recovered in the ASTRAL analysis were also generally consistent with the ML tree based on the complete matrix. However, the ASTRAL tree differed from the analysis of the complete dataset and most studies to date in the placement of Tylodina, which was weakly supported as sister to Aplysia + Pteropoda. The ASTRAL tree also differed with respect to relationships within Panpulmonta. Figure 2. View largeDownload slide A. Maximum-likelihood tree constructed from dataset generated by MARE. Only bootstrap values less than 100 are shown. Scale bar indicates 0.05 substitutions per site. B. ASTRAL tree constructed from individual maximum-likelihood gene trees. Only bootstrap values less than 100 are shown. Figure 2. View largeDownload slide A. Maximum-likelihood tree constructed from dataset generated by MARE. Only bootstrap values less than 100 are shown. Scale bar indicates 0.05 substitutions per site. B. ASTRAL tree constructed from individual maximum-likelihood gene trees. Only bootstrap values less than 100 are shown. Manual examination of single-gene alignments revealed a small number of non-molluscan contaminant sequences in several different taxa, but proved to be particularly important for P. californica (NCBI SRR1505130) which, based on our BLAST searches, was apparently infected with an apicomplexan parasite. Contamination in this transcriptome was initially suspected because of an abnormally long branch length for this taxon in a preliminary analysis (data not shown). Notably, the apicomplexan contamination in this transcriptome would not have affected the results of Kocot et al. (2011, 2013a), which recovered Panpulmonata as paraphyletic, because a different (smaller, 454 pyrosequencing-based) Pleurobranchaea transcriptome (NCBI SRR108976) was used. After detecting this contamination, we opted to screen all alignments manually and exclude this contamination—instead of using the Kocot et al. (2011) 454 transcriptome, because the latter dataset is much smaller. Based on data available in the NCBI nucleotide database from medically important apicomplexans (e.g. Plasmodium falciparum), contaminated sequences were easily detected and excluded using the approaches described above (see Methods). Despite apparent apicomplexan contamination in the Pleurobranchaea transcriptome of Zapata et al. (2014), they observed unremarkable branch lengths for this taxon and recovered Nudipleura with maximal support in all analyses without employing any special contamination filtering approaches. This suggests that the AGALMA (Dunn, Howison & Zapata, 2013) pipeline used for gene selection by Zapata et al. (2014) may be more efficient than the approach used herein at excluding such contaminant sequences. In summary, consistent with analyses of datasets based on PCR-amplified gene fragments (e.g. Göbbeler & Klussmann-Kolb, 2010; Dinapoli & Klussmann-Kolb, 2010; Jörger et al., 2010; Schrödl et al., 2011; Zapata et al., 2014; Kano et al., 2016) and some phylogenomic analyses conducted to date, our results support the clade Nudipleura, recovering Nudibranchia as sister to Pleurobranchomorpha with maximal support. Relationships within Nudipleura and Euopisthobranchia were consistent with previous studies (Jörger et al., 2010; Schrödl et al., 2011; Zapata et al., 2014; Kano et al., 2016). However, our results for relationships among the relatively few sampled members of Panpulmonata differed among analyses. Ongoing research in multiple laboratories to obtain transcriptome data from previously unsampled heterobranch lineages will undoubtedly continue to improve our understanding of this diverse and important gastropod clade. SUPPLEMENTARY MATERIAL Supplementary material is available at Journal of Molluscan Studies online. ACKNOWLEDGEMENTS This work was supported by start-up funds from the University of Alabama to KMK. We thank C. Held and the scientists and crew of R/V Polarstern cruise PS96 for facilitating collection of Bathyberthella antartica. We thank C. Laumer for sharing a script to trim alignment files based on ZORRO output. We thank two anonymous reviewers, Associate Editor H. Wägele, and R. Varney for helpful feedback on earlier drafts of this manuscript. We also thank the University of Alabama Randall Research Scholars Program. AUTHORS’ CONTRIBUTIONS EAP conducted the laboratory work and led the bioinformatic analyses. KMK conceived the study and assisted with bioinformatic analyses. Both authors contributed to writing the manuscript. REFERENCES Altschul , S. , Gish , W. , Miller , W. , Myers , E. & Lipman , D. 1990 . Basic local alignment search tool . Journal of Molecular Biology , 125 : 403 – 410 . Google Scholar CrossRef Search ADS Dayrat , B. , Conrad , M. , Balayan , S. , White , T.R. , Albrecht , C. , Golding , R. , Gomez , S. , Harasewych , M.G. & de Frias Martins , A.M. 2011 . Phylogenetic relationships and evolution of pulmonate gastropods (Mollusca): new insights from increased taxon sampling . Molecular Phylogenetics and Evolution , 59 : 425 – 437 . Google Scholar CrossRef Search ADS PubMed Dinapoli , A. & Klussmann-Kolb , A. 2010 . The long way to diversity—phylogeny and evolution of the Heterobranchia (Mollusca: Gastropoda) . Molecular Phylogenetics and Evolution , 55 : 60 – 76 . Google Scholar CrossRef Search ADS PubMed Dunn , C.W. , Howison , M. & Zapata , F. 2013 . Agalma: an automated phylogenomics workflow . BMC Bioinformatics , 14 : 330 . Google Scholar CrossRef Search ADS PubMed Ebersberger , I. , Strauss , S. & Von Haeseler , A. 2009 . HaMStR: Profile hidden Markov model based search for orthologs in ESTs . BMC Evolutionary Biology , 9 : 157 . Google Scholar CrossRef Search ADS PubMed Gillette , R. & Brown , J.W. 2015 . The sea slug, Pleurobranchaea californica: a signpost species in the evolution of complex nervous systems and behavior . Integrative and Comparative Biology , 55 : 1058 – 1069 . Google Scholar PubMed Goodheart , J.A. , Bazinet , A.L. , Collins , A.G. & Cummings , M.P. 2015 . Relationships within Cladobranchia (Gastropoda: Nudibranchia) based on RNA-Seq data: an initial investigation . Royal Society Open Science , 2 ( 9 ): 150196 . Google Scholar CrossRef Search ADS PubMed Goodheart , J.A. , Bazinet , A.L. , Valdés , Á. , Collins , A.G. & Cummings , M.P. 2017 . Prey preference follows phylogeny: evolutionary dietary patterns within the marine gastropod group Cladobranchia (Gastropoda: Heterobranchia: Nudibranchia) . BMC Evolutionary Biology , 17 : 221 . Google Scholar CrossRef Search ADS PubMed Grabherr , M. , Haas , B. , Yassour , M. , Levin , J. , Thompson , D. , Amit , I. , Adiconis , X. , Fan , L. , Raychowdhury , R. , Zeng , Q. , Chen , Z. , Mauceli , E. , Hacohen , N. , Gnirke , A. , Rhind , N. , Palma , F. , Birren , B. , Nusbaum , C. , Lindblad-Toh , K. , Friedman , N. & Regev , A. 2011 . Full-length transcriptome assembly from RNA-Seq data without a reference genome . Nature Biotechnology , 29 : 644 – 652 . Google Scholar CrossRef Search ADS PubMed Grande , C. , Templado , J. & Zardoya , R. 2008 . Evolution of gastropod mitochondrial genome arrangements . BMC Evolutionary Biology , 8 : 61 . Google Scholar CrossRef Search ADS PubMed Göbbeler , K. & Klussmann-Kolb , A. 2010 . Out of Antarctica?–New insights into the phylogeny and biogeography of the Pleurobranchomorpha (Mollusca, Gastropoda) . Molecular Phylogenetics and Evolution , 55 : 996 – 1007 . Google Scholar CrossRef Search ADS PubMed Haszprunar , G. 1985 . The Heterobranchia—a new concept of the phylogeny of the higher Gastropoda . Journal of Zoological Systematics and Evolutionary Research , 23 : 15 – 37 . Google Scholar CrossRef Search ADS Jörger , K. , Stoger , I. , Kano , Y. , Fukada , H. , Knebelsberger , T. & Schrödl , M. 2010 . On the origin of Acochlidia and other enigmatic euthyneuran gastropods, with implications for the systematics of Heterobranchia . BMC Evolutionary Biology , 10 : 323 . Google Scholar CrossRef Search ADS PubMed Kano , Y. , Brenzinger , B. , Nützel , A. , Wilson , N.G. & Schrödl , M. 2016 . Ringiculid bubble snails recovered as the sister group to sea slugs (Nudipleura) . Scientific Reports , 6 : 30908 . Google Scholar CrossRef Search ADS PubMed KATOH , K. , KUMA , K. I. , TOH , H. & MIYATA , T. 2005 . MAFFT version 5: improvement in accuracy of multiple sequence alignment . Nucleic Acids Research , 33 : 511 – 518 . Google Scholar CrossRef Search ADS PubMed Klussmann-Kolb , A. , Dinapoli , A. , Kuhn , K. , Streit , B. & Albrecht , C. 2008 . From sea to land and beyond–new insights into the evolution of euthyneuran Gastropoda (Mollusca) . BMC Evolutionary Biology , 8 : 57 . Google Scholar CrossRef Search ADS PubMed Kocot , K. , Cannon , J. , Todt , C. , Citarella , M. , Kohn , A. , Meyer , A. , Santos , S. , Schander , C. , Moroz , L. , Lieb , B. & Halanych , K. 2011 . Phylogenomics reveals deep molluscan relationships . Nature , 477 : 452 – 456 . Google Scholar CrossRef Search ADS PubMed Kocot , K. , Citarella , M. , Moroz , L. & Halanych , K. 2013 b. PhyloTreePruner: a phylogenetic tree-based approach for selection of orthologous sequences for phylogenomics . Evolutionary Bioinformatics , 2013 : 429 – 435 . Kocot , K. , Halanych , K. & Krug , J. 2013 a. Phylogenomics supports Panpulmonata: opisthobranch paraphyly and key evolutionary steps in a major radiation of gastropod molluscs . Molecular Phylogenetics and Evolution , 69 : 764 – 771 . Google Scholar CrossRef Search ADS PubMed Kocot , K. , Struck , T. , Merkel , J. , Waits , D. , Todt , C. , Brannock , P. , Weese , D. , Cannon , J. , Moroz , L. , Lieb , B. & Halanych , B. 2017 . Phylogenomics of Lophotrochozoa with consideration of systematic error . Systematic Biology , 66 : 256 – 282 . Google Scholar PubMed Kumar , S. , Stecher , G. & Tamura , K. 2016 . MEGA7: molecular evolutionary genetics analysis version 7.0 for bigger datasets . Molecular Biology and Evolution , 33 : 1870 – 1874 . Google Scholar CrossRef Search ADS PubMed Kück , P. & Meusemann , K. 2010 . FASconCAT: convenient handling of data matrices . Molecular Phylogenetics and Evolution , 56 : 1115 – 1118 . Google Scholar CrossRef Search ADS PubMed Medina , M. , Lal , S. , Vallès , Y. , Takaoka , T.L. , Dayrat , B.A. , Boore , J.L. & Gosliner , T. 2011 . Crawling through time: transition of snails to slugs dating back to the Paleozoic, based on mitochondrial phylogenomics . Marine Genomics , 4 : 51 – 59 . Google Scholar CrossRef Search ADS PubMed Mirarab , S. , Reaz , R. , Bayzid , M.S. , Zimmermann , T. , Swenson , M.S. & Warnow , T. 2014 . ASTRAL: genome scale coalescent based species tree estimation . Bioinformatics , 30 : i541 – i548 . Google Scholar CrossRef Search ADS PubMed Misof , B. , Meyer , B. , von Reumont , B. , Kück , P. , Misof , K. & Meusemann , K. 2013 . Selecting informative subsets of sparse supermatrices increases the chance to find correct trees . BMC Bioinformatics , 14 : 348 . Google Scholar CrossRef Search ADS PubMed Moroz , L. , Edwards , J. , Puthanveettil , S. , Kohn , A. , Ha , T. , Heyland , A. , Knudsen , B. , Sahni , A. , Yu , F. , Liu , L. , Jezzini , S. , Lovell , P. , Iannucculli , W. , Chen , M. , Nguyen , T. , Sheng , H. , Shaw , R. , Kalachikov , S. , Panchin , Y. , Farmerie , W. , Russo , J. , Ju , J. & Kandel , E. 2006 . Neuronal transcriptome of Aplysia: neuronal compartments and circuitry . Cell , 127 : 1453 – 1467 . Google Scholar CrossRef Search ADS PubMed Newman , D. & Cragg , G. 2014 . Marine-sourced anti-cancer and cancer pain control agents in clinical and late preclinical development . Marine Drugs , 12 : 255 – 278 . Google Scholar CrossRef Search ADS PubMed Oskars , T.R. , Bouchet , P. & Malaquias , M.A.E. 2015 . A new phylogeny of the Cephalaspidea (Gastropoda: Heterobranchia) based on expanded taxon sampling and gene markers . Molecular Phylogenetics and Evolution , 89 : 130 – 150 . Google Scholar CrossRef Search ADS PubMed Price , M.N. , Dehal , P.S. & Arkin , A.P. 2010 . FastTree 2—approximately maximum-likelihood trees for large alignments . PLoS One , 5 : e9490 . Google Scholar CrossRef Search ADS PubMed Rice , P. , Longden , I. & Bleasby , A. 2000 . EMBOSS: the European molecular biology open software suite . Trends in Genetics , 16 : 276 – 277 . Google Scholar CrossRef Search ADS PubMed Roure , B. , Baurain , D. & Philippe , H. 2012 . Impact of missing data on phylogenies inferred from empirical phylogenomic data sets . Molecular Biology and Evolution , 30 : 197 – 214 . Google Scholar CrossRef Search ADS PubMed Schilthuizen , M. , Liew , T. , Elhan , B. & Lackman-Ancrenaz , I. 2005 . Effects of karst forest degradation on pulmonate and prosobranch land snail communities in Sabah, Malaysian Borneo . Conservation Biology , 19 : 949 – 954 . Schrödl , M. 2014 . Time to say “bye-bye Pulmonata” . Spixiana , 37 : 161 – 164 . Schrödl , M. , Jörger , K. , Klussmann-Kolb , A. & Wilson , N. 2011 . Bye bye “Opisthobranchia”! A review on the contribution of mesopsammic sea slugs to euthyneuran systematics . Thalassas , 27 : 101 – 112 . Senatore , A. , Edirisinghe , N. & Katz , P. 2015 . Deep mRNA sequencing of the Tritonia diomedea brain transcriptome provides access to gene homologues for neuronal excitability, synaptic transmission and peptidergic signalling . PLoS One , 10 : e0118321 . Google Scholar CrossRef Search ADS PubMed Stamatakis , A. 2014 . RAxML Version 8: a tool for phylogenetic analysis and post-analysis of large phylogenies . Bioinformatics , 30 : 1312 – 1313 . Google Scholar CrossRef Search ADS PubMed Telford , M.J. 2008 . Resolving animal phylogeny: a sledgehammer for a tough nut? Developmental Cell , 14 : 457 – 459 . Google Scholar CrossRef Search ADS PubMed White , T.R. , Conrad , M.M. , Tseng , R. , Balayan , S. , Golding , R. , de Frias Martins , A.M. & Dayrat , B.A. 2011 . Ten new complete mitochondrial genomes of pulmonates (Mollusca: Gastropoda) and their impact on phylogenetic relationships . BMC Evolutionary Biology , 11 : 295 . Google Scholar CrossRef Search ADS PubMed Wiens , J.J. & Tiu , J. 2012 . Highly incomplete taxa can rescue phylogenetic analyses from the negative impacts of limited taxon sampling . PLoS One , 7 : 42925 . Google Scholar CrossRef Search ADS WILLAN , R. C. & BERTSCH , H. 1987 . Description of a new pleurobranch (Opisthobranchia: Notaspidea) from Antarctic waters, with a review of notaspideans from southern Polar seas . Veliger , 29 : 292 – 302 . Wilson , N.G. , Jörger , K.M. , Brenzinger , B. & Schrödl , M. 2017 . Phylogenetic placement of the enigmatic worm-like Rhodopemorpha slugs as basal Heterobranchia . Journal of Molluscan Studies , 83 : 399 – 408 . Google Scholar CrossRef Search ADS Wu , M. , Chatterji , S. & Eisen , J.A. 2012 . Accounting for alignment uncertainty in phylogenomics . PLoS One , 7 : e30288 . Google Scholar CrossRef Search ADS PubMed Wägele , H. , Klussmann-Kolb , A. , Verbeek , E. & Schrödl , M. 2014 . Flash back and foreshadowing—a review of the taxon Opisthobranchia . Organisms, Diversity & Evolution , 14 : 133 – 149 . Google Scholar CrossRef Search ADS Wägele , H. , Knezevic , K. & Moustafa , A.Y. 2017 . Distribution and morphology of defensive acid-secreting glands in Nudipleura (Gastropoda: Heterobranchia), with an emphasis on Pleurobranchomorpha . Journal of Molluscan Studies , 83 : 422 – 433 . Google Scholar CrossRef Search ADS Wägele , H. & Willan , R. 2000 . Phylogeny of the Nudibranchia . Zoological Journal of the Linnean Society , 130 : 83 – 181 . Google Scholar CrossRef Search ADS Zapata , F. , Wilson , N. , Howison , M. , Andrade , S. , Jörger , K. , Schrödl , M. , Goetz , F. , Giribet , G. & Dunn , C. 2014 . Phylogenomic analyses of deep gastropod relationships reject Orthogastropoda . Proceedings of the Royal Society of London B , 281 : 20141739 . Google Scholar CrossRef Search ADS © The Author(s) 2018. Published by Oxford University Press on behalf of The Malacological Society of London, all rights reserved. For Permissions, please email: journals.permissions@oup.com This article is published and distributed under the terms of the Oxford University Press, Standard Journals Publication Model (https://academic.oup.com/journals/pages/about_us/legal/notices) TI - Phylogenomics confirms monophyly of Nudipleura (Gastropoda: Heterobranchia) JF - Journal of Molluscan Studies DO - 10.1093/mollus/eyy013 DA - 2018-04-30 UR - https://www.deepdyve.com/lp/oxford-university-press/phylogenomics-confirms-monophyly-of-nudipleura-gastropoda-DkHX0knyuJ SP - 1 EP - 265 VL - Advance Article IS - 3 DP - DeepDyve ER -