Skip to main content
  • Research article
  • Open access
  • Published:

Prey preference follows phylogeny: evolutionary dietary patterns within the marine gastropod group Cladobranchia (Gastropoda: Heterobranchia: Nudibranchia)



The impact of predator-prey interactions on the evolution of many marine invertebrates is poorly understood. Since barriers to genetic exchange are less obvious in the marine realm than in terrestrial or freshwater systems, non-allopatric divergence may play a fundamental role in the generation of biodiversity. In this context, shifts between major prey types could constitute important factors explaining the biodiversity of marine taxa, particularly in groups with highly specialized diets. However, the scarcity of marine specialized consumers for which reliable phylogenies exist hampers attempts to test the role of trophic specialization in evolution. In this study, RNA-Seq data is used to produce a phylogeny of Cladobranchia, a group of marine invertebrates that feed on a diverse array of prey taxa but mostly specialize on cnidarians. The broad range of prey type preferences allegedly present in two major groups within Cladobranchia suggest that prey type shifts are relatively common over evolutionary timescales.


In the present study, we generated a well-supported phylogeny of the major lineages within Cladobranchia using RNA-Seq data, and used ancestral state reconstruction analyses to better understand the evolution of prey preference. These analyses answered several fundamental questions regarding the evolutionary relationships within Cladobranchia, including support for a clade of species from Arminidae as sister to Tritoniidae (which both preferentially prey on Octocorallia). Ancestral state reconstruction analyses supported a cladobranchian ancestor with a preference for Hydrozoa and show that the few transitions identified only occur from lineages that prey on Hydrozoa to those that feed on other types of prey.


There is strong phylogenetic correlation with prey preference within Cladobranchia, suggesting that prey type specialization within this group has inertia. Shifts between different types of prey have occurred rarely throughout the evolution of Cladobranchia, indicating that this may not have been an important driver of the diversity within this group.


Predator-prey interactions are among the most fundamental processes in ecology and constitute the fabric of community structure and ecosystem function [1, 2]. However, the role of those interactions in evolution, and their impacts on biodiversity, are less well understood in marine systems [3, 4]. The most widely accepted hypothesis to explain the origin of biological diversity traces its origins to Mayr [5, 6], who proposed that the ranges of organisms are fragmented by the formation of physical barriers, resulting in isolation and divergence in allopatry. However, in the marine realm, where barriers to genetic exchange are less obvious than in terrestrial or freshwater systems [7], non-allopatric divergence and speciation may play a fundamental role in the generation of biodiversity (e.g., [8, 9]). In this context, shifts between major prey types (e.g., different cnidarian classes) could constitute important factors explaining the biodiversity of marine taxa, particularly in groups with highly specialized diets. The greatest obstacle to testing these ideas is the lack of well-supported phylogenies for groups of specialized consumers.

In this study, we generate RNA-Seq data to test the role of prey preference shifting in the evolution of Cladobranchia (Mollusca: Gastropoda: Hererobranchia: Nudibranchia), a group of marine invertebrates with at least 1000 species [10]. Cladobranch sea slugs occupy various marine environments, from coastal reefs, where diversity is highest, to the deep sea, as well as highly specialized pelagic and neustonic niches [11,12,13,14]. Species of Cladobranchia are exclusively carnivorous, and exhibit diverse dietary specializations, preying on a variety of animal taxa, including bryozoans and crustaceans, eggs of fishes and molluscs, and cnidarians (Fig. 1) [15,16,17]. However, the vast majority of cladobranchs prey on species of the two most diverse clades within Cnidaria, Anthozoa (e.g., anemones, stony corals, and octocorals) and Hydrozoa (hydroids, siphonophores, and hydromedusae) [15, 18,19,20]. This preference for cnidarian prey is hypothesized to have facilitated the evolution of the ability to sequester cnidarian nematocysts in Cladobranchia [21], which is believed to have evolved only once within this group [16].

Fig. 1
figure 1

Select photographs of cladobranch taxa on their food source, including: (a) Dondice parguerensis on the scyphozoan jellyfish Cassiopea sp., (b) Doto chica on the hydroid Eudendrium sp.; (c) Tritonia hamnerorum on the octocoral Gorgonia ventalina; and (d) Favorinus tsuruganus on an opisthobranch egg mass (Photo credits: Ángel Valdés)

Based on recent classifications, two of the three main groups in Cladobranchia (Aeolidida and Dendronotida [22]; the third being Arminida, as defined in [23]) contain taxa that prey on animals distributed across the list given above. These classifications suggest that shifts in prey type preference are relatively common throughout Cladobranchia over evolutionary timescales. There also exist many well-documented cases where cladobranch species are tightly associated with specific prey types or species (e.g., [24,25,26,27,28,29]). In many of these cases the prey species might even be considered a host as defined by Coyne & Orr [93], due to similarities that many sea slugs share with herbivorous insects, including their small size relative to their hosts and the use of hosts for both food and shelter [3].

Two main hypotheses have emerged regarding the roles that dietary specialization and prey shifts have played in the evolution of heterobranch sea slugs (formerly called opisthobranchs). The first of these hypotheses suggests that increased speciation occurs due to species-specific prey switching in groups where specialization is prevalent [30]. This leads to clades consisting of many taxa that specialize on individual prey species. In many metazoan groups studied, mainly involving terrestrial symbiotic and parasitic systems [31,32,33], host shifting (shifting between prey species) has been implicated as a driver of diversification, with colonization of new hosts often leading to bursts of cladogenesis. Speciation of taxa by host or prey shifting may also be important in a handful of specialized marine consumers such as some bivalves [34], amphipods [35], barnacles [36], gobies [37], and gastropods [24]. We do not assess this hypothesis here, as it requires broad taxon sampling across Cladobranchia and solid evidence of dietary specialization, both of which are lacking in many cases.

The second hypothesis relating to dietary specialization is that major radiations within heterobranch sea slugs may be related to the evolution of particular morphological structures necessary for feeding on different types of prey [16, 17, 38], such as the distinct radular morphology present in members of Aeolidida. This hypothesis suggests that shifting to new prey items leads to an increase in niche availability, similar to the effect of habitat preference shifts in some groups [39]. This hypothesis is broader than the first, in that it refers to switches between prey types at higher levels of organization. The consequences of this type of switching relate to species that, following a switch, are able to prey on multiple taxa within a general prey type, rather than explicitly focusing on those that specialize on certain prey species. This pattern has been found in only a few taxa [34, 40, 41]. Diversification in this context relates more to the expansion of possible prey types rather than specialization.

Given the variety of prey type preferences exhibited by its members, Cladobranchia constitutes an excellent system to explore the relationship between prey shifting and cladogenesis. Until now, it has not been possible to test how prey choice has evolved through the history of this group, because existing cladobranch phylogenies are notoriously poorly resolved. Support for Cladobranchia as a monophyletic group is high [23, 42, 43], but the relationships among major lineages within Cladobranchia have long been problematic [22, 23, 42, 44,45,46,47]. However, a recent phylogenomic study provided evidence that these lineages could be resolved with RNA-Seq data [43]. In addition to the growing availability of RNA-Seq data for this group, the prey preferences of the majority of species within this group have been published (e.g., [15, 16, 18,19,20, 48]).

To address the role of dietary specialization and host shifts in the evolution of this group and resolve outstanding systematics issues, we reconstruct the phylogeny of Cladobranchia using RNA-Seq data. In this study, we increase the taxon sampling compared to previous phylogenomic work on Cladobranchia [43] by incorporating additional diversity from both previously sampled clades (Aeolidida and Dendronotida) and the previously unsampled Arminida (as defined in [23], though we only include members of Arminidae). In addition, we seek to address patterns of prey type switching among and within the major lineages of cladobranchs by assessing prey type preference for each taxon included in the phylogenetic analyses, and using these data to reconstruct the most likely ancestral prey type preference for each node in the tree. These analyses provide the means to examine the prevalence of prey type switching within Cladobranchia in order to provide a framework for studying how dietary preferences may have affected evolution within this group.


Organismal sampling

One or two specimens of each of 16 representative species were collected in tide pools or via snorkeling or SCUBA (self-contained underwater breathing apparatus; under AAUS certification) using a variety of methods including direct collection, substrate collection, and non-destructive collecting under rocks. A visual examination was used for confirmation of identity using field guides for the Caribbean [13] and the Indo-Pacific [14]. Barcode sequences and expert opinions were used when the identity of specimens was still uncertain. Images of select specimens are in Additional file 1: Figure S1 and Additional file 2: Figure S2. One of the two specimens was placed in RNAlater solution (Qiagen, Hilden, Germany) for RNA preservation and frozen at -80 ºC within one week of collection to prevent RNA degradation. Some specimens in RNAlater were instead stored at -20 ºC within 24 h and remained there for up to a month. A second specimen of each species, when available, was fixed as a voucher for morphological analysis, first in 10% formalin and subsequently preserved in 70% ethanol for long-term storage. Voucher specimens were deposited in the Smithsonian National Museum of Natural History (NMNH) and are available for study under the catalog numbers provided in Additional file 3: Table S5.

For Bulbaeolidia alba, Hancockia uncinata, Unidentia angelvaldesi, Bornella anguilla, Dermatobranchus sp., Phestilla sp., and Eubranchus rustyus, we were unable to obtain morphological vouchers. Cella et al. [49] proposed to use the genus name Tenellia for species previously assigned to Catriona, Cuthona and Phestilla, but recognized that the assignment of species to Tenellia is problematic due to the absence of morphological synapomorphies. Thus, we chose to temporarily maintain species in the genera Catriona, Cuthona and Phestilla until the additional studies suggested by Cella et al. [49] are carried out. We generated RNA-Seq data for 16 Cladobranchia species, downloaded data for 19 additional Cladobranchia species from the NCBI Sequence Read Archive (SRA) and obtained two RNA-Seq datasets from colleagues at Georgia State University. Three outgroup RNA-Seq datasets were also obtained from the SRA: two representatives of Anthobranchia (the sister taxon of Cladobranchia; [22, 46]), and one of Pleurobranchoidea (the sister taxon to Nudibranchia [50]). Specimen data, SRA numbers and barcode GenBank numbers are listed in Additional file 3: Table S5.

RNA extraction and sequencing

A 20–100 mg tissue sample was taken from the anterior of each animal and homogenized using a motorized pestle. In some cases, the specimen was so small the entire animal was used. After homogenizing for 1–2 min the tissue was flash-frozen in liquid nitrogen for subsequent homogenizing until tissue mixture was fully uniform. 500 μL of TRIzol Reagent (Life Technologies, Carlsbad, CA, USA) was then added and the mixture was homogenized again. This procedure was repeated until the solution was deemed fully homogenized. Once this process was complete, an additional 500 μL of TRIzol Reagent was added to the solution and the mixture was left at room temperature for five min.

Following the five min incubation, 100 μL of 1-Bromo-3-chloropropane was added to the solution, which was subsequently mixed thoroughly. The mixture was then left at room temperature for five min, and then centrifuged at 16,000 g for 20 min at 8 ºC. The top aqueous phase was then removed and placed in another tube where 500 μL of 100% isopropanol was added, and stored overnight at -20 ºC for RNA precipitation.

After precipitation, the samples were centrifuged at 17,200 g for 10 min at 4 ºC. The supernatant was then removed and the pellet washed with freshly prepared 75% ethanol. The sample was then centrifuged at 7,500 g for 5 min at 4 ºC. The supernatant was removed and the pellet air-dried for 1 to 2 min (or until it looked slightly gelatinous and translucent). The total RNA was then re-suspended in 10–30 μL of Ambion Storage Solution (Life Technologies, Carlsbad, CA, USA), and 1 μL of SUPERase•In (Thermo Fisher Scientific, Waltham, Massachusetts, USA) was added to prevent degradation.

Total RNA samples were submitted to the DNA Sequencing Facility at University of Maryland Institute for Bioscience and Biotechnology Research, where quality assessment, library preparation, and sequencing were performed. RNA quality assessment was done with a Bioanalyzer 2100 (Agilent Technologies, Santa Clara, CA, USA), and samples with a concentration higher than 20 ng/μL were used for library construction. Library preparation used the Illumina TruSeq RNA Library Preparation Kit v2 (Illumina, San Diego, CA, USA) and 200 bp inserts; 100 bp, paired-end reads were sequenced with an Illumina HiSeq1000 (Illumina, San Diego, CA, USA).

Quality control and assembly of reads

Reads that failed to pass the Illumina “Chastity” quality filter were excluded from our analyses. Reads passing the quality filter were assembled using Trinity (version 2.1.1; [51]) with default settings, which required assembled transcript fragments to be at least 200 bp in length.

Orthology assignment

Translated transcript fragments were organized into orthologous groups corresponding to a custom gastropod-specific core-ortholog set of 3,854 protein models [43] using HaMStR (version 13.2.2; [52]), which in turn used FASTA (version 36.3.6d; [53]), GeneWise (version 2.2.0; [54]), and HMMER (version 3.1b2; [55]). In the first step of the HaMStR procedure, substrings of assembled transcript fragments (translated nucleotide sequences) that matched one of the gastropod protein models were provisionally assigned to that orthologous group. To reduce the number of highly divergent, potentially paralogous sequences returned by this search, we set the E-value cutoff defining an HMM hit to 1e-05 (the HaMStR default is 1.0), and retained only the top-scoring quartile of hits. In the second HaMStR step, the provisional hits from the HMM search were compared to the reference taxon, Aplysia californica, and retained only if they survived a reciprocal best BLAST hit test with the reference taxon using an E-value cutoff of 1e-05 (the HaMStR default was 10.0). In our implementation, we substituted FASTA [53] for BLAST [56] because FASTA programs readily accepted our custom amino acid substitution matrix (GASTRO50; [43]).

Construction of data matrix and paralogy filtering

Protein sequences in each orthologous group were aligned using MAFFT (version 7.187; [57]). We used the --auto and --addfragments options of MAFFT to align transcript fragments to the Aplysia californica reference sequence, which was considered the existing alignment. We converted the protein alignments to corresponding nucleotide alignments using a custom Perl script. A maximum likelihood tree was inferred using GARLI (Genetic Algorithm for Rapid Likelihood Inference version 2.1; [58]) for each orthologous group where at least 75% of the taxa were present (716 orthologous groups), and was given as input to PhyloTreePruner (version 1.0; [59]). Orthologous groups that showed evidence of out-paralogs for any taxa (352 orthologous groups out of 716) were pruned according to the default PhyloTreePruner protocol, which removes all additional sequences outside of a maximally inclusive sub-tree. For orthologous groups containing in-paralogs, multiple sequences were combined into a single consensus sequence for each taxon, and orthologous groups for which fewer than 75% of taxa remained were discarded. This process left 406 orthologous groups eligible for inclusion in data matrices. Individual orthologous group alignments were concatenated (nt123 matrix) (Table 1). Codons not represented by sequence data in at least four taxa were then removed (nt123sitesremoved matrix).

Table 1 Data matrix statistics for each of the two data matrices

Phylogenetic analyses

Four separate phylogenetic analyses were completed in this study: (i) an analysis with the nt123 data matrix partitioned by codon position (nt123partitioned) by assigning different model parameters and rates to the three types of codon positions, (ii) an analysis with the nt123sitesremoved data matrix partitioned by codon position (nt123sitesremoved_partitioned), (iii) an analysis of the nt123 matrix partitioned by codon position, but excluding the third position (nt12partitioned), and (iv) an analysis of the unpartitioned nt123sitesremoved data matrix (nt123sitesremoved_unpartitioned). To conduct all four phylogenetic analyses we used GARLI (version 2.1; [58]) through the GARLI web service hosted at [60]. We used the default settings in GARLI, including a general time reversible substitution model (GTR; [61]) with a rate heterogeneity model with a proportion of invariant sites estimated (+I; [62]) and the remainder with a gamma distribution (+G; [63]), along with stepwise-addition starting trees. Post-processing of the phylogenetic inference results was performed by the GARLI web service at using DendroPy [64] and the R system for statistical computing [65]. For all analyses, 1000 bootstrap replicates were generated and a best tree search was performed with 10 search replicates.

Ancestral state reconstruction

We conducted a literature search to collect prey preference data for all nudibranch taxa in our phylogeny and coded each species as Anthozoa: Octocorallia, Anthozoa: Hexacorallia, Hydrozoa, Scyphozoa, Bryozoa, Crustacea, Gastropoda eggs, or generalist, for a total of eight states (Additional file 4: Table S4). In the cases where more than one type of prey is fed upon by an individual species, we provide that information and run additional analyses to test the effect of these alternatives on the final results. The final analysis incorporates the prey type for each species that that species has been observed to feed on more than 50% of the time. Data was compiled primarily from review papers on feeding and defense in nudibranchs [15, 16, 18,19,20], field guides [13, 66], one additional paper [67], and web sources where necessary [68,69,70,71]. Though limited, the taxon selection in this study represents a large portion of the morphological and ecological diversity of Cladobranchia, including the diversity of prey type preferences. Using these character states, we compared the fit of three discrete trait models using the AICcmodavg 2.0-4 package [72] in R 3.3.1 [65]. We assessed fit for models where: (i) all transition rates were equal (ER); (ii) forward and reverse transitions were equal between states (i.e. symmetrical, SYM); and (iii) all transition rates were different (ARD) using the corrected Akaike information criterion (AICc). The ER model (AICc = 100.61) was a better fit to the data than either the SYM model (AICc = 118.53) or the ARD model (AICc = 165.16). The final ancestral state reconstruction analysis was completed using the ace function, in the APE package [73], under the ER model using default parameters and a joint reconstruction approach. The ace function uses a Markov model employing a maximum likelihood approach. In this analysis, the reconstructed ancestral states that are returned are given as the proportion of the total likelihood calculated for each state for each node.


Assembly and data matrix properties

The raw number of RNA-Seq reads for each species ranged from 25,756,442 to 133,156,930 (x̅ ≈ 49M reads; Additional file 5: Table S1). Once assembled, the number of transcript fragments per sample ranged from 71,967 to 295,127 (x̅ = 146,403; Additional file 6: Table S2). N50 ranged from 395 to 1,058 bp (x̅ = 716 bp). The transcript fragments from each assembly that matched the HaMStR database ranged from 615 to 2,013 (x̅ = 1,282). However, the number of matches to unique orthologous groups ranged from 512 to 1,198 (x̅ = 935). The mean length of transcript fragment matches to the HaMStR database was 282 amino acids. HaMStR results are presented in Additional file 7: Table S3.

Phylogenetic results

Results from all analyses supported Cladobranchia as a monophyletic group with a bootstrap (BS) value of 100% (Fig. 2). Arminidae (Arminida) is also supported as monophyletic (BS = 100%), and is sister (BS = 100%) to Tritoniidae (BS = 100%).

Fig. 2
figure 2

The maximum likelihood topology from the nt123_partitioned analysis, with bootstrap support values from each analysis labeled on some nodes (nt123_partitioned/nt123sitesremoved_partitioned/nt123sitesremoved_ unpartitioned/nt12partitioned). All unlabeled nodes have 100% bootstrap support in all analyses

Dendronotida is non-monophyletic across all topologies, with dendronotid taxa comprising two major clades. The first of these clades is sister to all other cladobranchs (BS = 100%) and contains Doto (Dotidae), Bornella (Bornellidae), Hancockia (Hancockiidae), Scyllaea (Scyllaeidae), Melibe (Tethyidae), Dendronotus (Dendronotidae), and Lomanotus (Lomanotidae).

Aeolidida is supported as monophyletic (BS = 100%) across all topologies, containing Flabellina (Flabellinidae), Berghia, Spurilla, Bulbaeolidia, Anteaeolidiella, and Limenandra (Aeolidiidae), Hermissenda, Dondice, Noumeaella, Favorinus, Palisa, Austraeolis, Learchis, and Phidiana (Facelinidae), Fiona, Cuthona, Catriona, Phestilla, and Eubranchus (Fionidae), and Unidentia (Unidentiidae). All families within Aeolidida where multiple taxa from the same family are included are supported as monophyletic, with the exception of Facelinidae, which is polyphyletic and forms two separate clades.

Two taxa previously unassigned to any of the three major clades, Dirona (Dironidae) and Janolus (Proctonotidae), are supported as sister taxa (BS = 100%) and form a clade that is sister to Aeolidida (BS = 100%).

Ancestral state reconstruction analysis

The ancestral state reconstruction results support the hypothesis that the most recent common ancestor (MRCA) of Cladobranchia preyed upon species of Hydrozoa (Fig. 3). A cladobranch that preyed upon Hydrozoa also appears to be the MRCA for Aeolidida and the clade composed of most of the taxa assigned to Dendronotida, as well as the rest of the MRCAs along the backbone of the tree. However, a taxon that preyed upon Octocorallia species is the likely MRCA for the Arminidae + Tritoniidae clade (88.78% of the scaled likelihood), a taxon that preyed upon bryozoans or hydrozoans is the most likely MRCA for the Dirona + Janolus clade (74.44% of the scaled likelihood; Fig. 3, Table 2), and the MRCA for Aeolidiidae most likely fed upon species within Hexacorallia (98.03% of the scaled likelihood). Additional ancestral state reconstruction analyses were completed to evaluate the effects of alternative prey types for certain taxa on the overall reconstruction of ancestral states (Additional file 8: Table S6, Additional file 9: Table S7, Additional file 10: Table S8, Additional file 11: Table S9, Additional file 12: Table S10). With the exception of the ancestral node of the Dirona + Janolus clade, which changes to >97% of the scaled likelihood supporting a Hydrozoa feeding ancestor in three of the alternative analyses, the results are robust to these changes. The scaled likelihoods across all other nodes within each of the alternative analyses remain within 5% of the value in the original analysis.

Fig. 3
figure 3

Ancestral state reconstruction results for the evolution of diet preference in Cladobranchia. Pie charts on the nodes are scaled likelihoods calculated using the ace function in APE. Alternative states and results are indicated in parentheses with an asterisk at the tips of the tree and nodes, and only alternative node states with greater than or equal to 5% difference from the original reconstruction are shown. Nodes are also labeled with numbers consistent with Table 2

Table 2 Ancestral state reconstruction results for the evolution of diet preference in Cladobranchia. This table provides the percentage (%) of the total likelihood assigned to each state for each node. The node numbers correspond to those provided in Fig. 3


In this study we significantly increased the breadth of RNA-Seq sampling in Cladobranchia in order to generate a robust phylogenetic hypothesis, and provide a framework for the evolution of prey type preference within this group.

Prey preference evolution in Cladobranchia

Well-supported clades recovered within Cladobranchia appear to be strongly associated with prey groups. Most of the larger clades recovered in the phylogenetic tree prey almost exclusively on particular types of organisms (Fig. 3), such as Aeolidiidae on Hexacorallia, Arminida + Tritoniidae on Octocorallia, and multiple clades that prey on Hydrozoa. This result is in opposition to previous studies [22, 45], which indicated that groupings within Cladobranchia contained taxa that fed on a broad range of prey types. The results here support the idea that prey preference within Cladobranchia may be a taxonomically useful trait for placing taxa into some groups. Past taxonomic work on Cladobranchia has focused on different anatomical features to diagnose groups, such as the presence of rhinophoral sheaths (Dendronotida) or oral veils (Arminida) [22]. In the future, incorporating feeding adaptations to particular prey taxa may accelerate taxonomic progress in the group.

Our results indicate that prey preference shifts from one major taxon to another are relatively rare in Cladobranchia. The ancestral state reconstruction unambiguously supports an ancestor for Cladobranchia that preyed upon Hydrozoa (Fig. 3). This analysis also suggests at least five transitions from hydrozoan to other prey taxa, such as Hexacorallia (Anthozoa), Octocorallia (Anthozoa), and Scyphozoa. Interestingly, the clade containing Dirona and Janolus has expanded to feeding on Bryozoa in addition to Hydrozoa, rather than shifting to Bryozoa exclusively. This expansion is mirrored in Phidiana, which is also able to feed on members of both taxa. Overall, expansions to feeding on multiple types of prey have occurred at least six times in Cladobranchia, leading to multiple generalist taxa (Melibe and Hermissenda) and those that can feed on both Hydrozoa and either Bryozoa or Crustacea. These are cases in which diversification might be related more to an increase in options rather than specialization.

The mechanism by which the evolution of prey preference is constrained is unknown, but it could result from relative difficulty in evolving specific traits for protection against nematocysts (or other defenses) from various cnidarian prey groups. Although it is possible for some species (e.g., Phidiana hiltoni) to prey upon different cnidarian species [74], there are few examples of cladobranchs preying on multiple, taxonomically distant cnidarians (only H. crassicornis in this study). Cladobranchs require a series of adaptations to prevent cnidarian nematocysts from firing or minimizing the damage in the case of firing [75], including mucous secretions. These secretions appear to be specific to the prey species in one studied case [76], and may be why switching between types of prey is much more challenging and occurs much less frequently than previously thought. Species of Cladobranchia that do not prey on cnidarians, such as Favorinus, can easily switch between egg masses of distantly related gastropods, including aplysiids, sacoglossans, and other nudibranchs [48], lending support to this hypothesis.

Previous work has suggested that dietary specialization on particular prey types was crucial in the evolution of Euthyneura [16, 17, 77], and has been proposed as a “driving force” in heterobranch sea slug evolution [38, 78]. Dietary specialization has also been considered a contributing factor in the species richness of Nudibranchia as a whole [17], and especially cladobranchs [16, 18] in conjunction with the evolution of nematocyst sequestration. This hypothesis is entirely plausible when looking at the numbers of species in prey groups and how those correlate with species diversity in the cladobranch predators. Although Bryozoa has nearly 6,000 species [79], there are fewer than 50 species within Cladobranchia that prey on members of this group [48]. Conversely, more than 700 species of cladobranchs prey on Hydrozoa, a clade of cnidarians with ~3,500 species [79]. This drastic difference is primarily due to Aeolidida, which contains a large proportion of the taxa that prey on Hydrozoa, and which appears to have diversified primarily while preying on hydrozoans. Our results do not support the hypothesis that prey type shifts lead to morphological adaptations that increase diversity, as the distinct radular morphology found in Aeolidida is not associated with a prey type shift according to our ancestral state reconstruction. However, Aeolidida is also one of two lineages where nematocyst sequestration is known to have evolved within Cladobranchia [80]. Given that hydrozoans are known to have the highest diversity of nematocyst types [81], the ability to sequester nematocysts may have had an impact on diversification.

The hypothesis that the diversity of larger clades within Cladobranchia is related to the frequency of major prey type shifts is not supported by these results. Instead, we suspect that if shifts in diversification associated with diet in Cladobranchia are going to be found, these may occur within groups that prefer a major prey type (e.g., at the family or genus level), where more species-specific prey shifting is likely to occur [16, 24]. The literature indicates that in groups where taxa are specialized on particular prey species, shifting to a new prey (or host) species often leads to speciation and diversification [82,83,84]. This pattern is found in many metazoan taxa, including flies [85], amphipods [35, 86, 87], alpheid shrimp [88, 89], barnacles [36], whelks [90], gobies [37], and sacoglossan gastropods [91, 92], and has been extensively investigated in phytophagous insects (reviewed in [93]). Within Nudibranchia, a large subset of taxa exhibit specialization on a single species, with many others preferring only two or three prey species [94]. We suspect that in the case of Cladobranchia, this specialization and prey shifting at the species level may be the primary impact that prey preference has on the diversification rate across lineages, rather than shifts to new prey types, as is true in many herbivorous insect lineages (reviewed in [93]). Tests of this hypothesis require a broader sampling of members of Cladobranchia for both the phylogenetic inference and species-specific prey preference data.

Systematics of Cladobranchia and prey preference within individual clades

Based on the phylogenetic hypothesis presented here, the monophyly of Cladobranchia is reinforced with full bootstrap support across all analyses. Though monophyly was indicated in previous morphological [22] and molecular [42,43,44,45,46, 95] analyses, there has also been a study suggesting paraphyly [23], though the authors of that paper contended that this might be due to a deletion of a string of nucleotides within one lineage (Melibe) that was biasing the results.


The most significant systematics results from this study involve Arminida, a group not included in the one previous phylogenomic study of Cladobranchia [43]. Arminida, when first described, comprised the genera Janolus and Dirona, among other taxa, including Arminidae [96]. The inclusion of Janolus and Dirona within Arminida renders this group paraphyletic in both morphological and molecular analyses [22, 97], and they, along with others (Charcotiidae and Pinufiidae), have since been removed from Arminida and are considered unassigned members of Cladobranchia [22]. The analyses presented here support this exclusion of Janolus and Dirona from Arminida, consistent with recent studies [22, 42, 98]. Both of these genera primarily prefer bryozoan prey, but also feed on members of Hydrozoa.

There is strong support for Arminidae (one of two families within Arminida) as the sister group to Tritoniidae, which is a novel result. This result is in agreement with only one previous phylogenetic hypothesis, which was generated using 18S rDNA data [42]. In all other previous studies, taxa from Arminida had been either unplaced within the Cladobranchia phylogeny [46, 95, 98] or supported as sister to various other combinations of taxa from Dendronotida and Aeolidida [42]. The position of Tritoniidae in the previous phylogenomics study of Cladobranchia was uncertain [43]. It appears that prey preference is particularly relevant for the evolution of Tritoniidae + Arminidae as species within this group prey exclusively on Octocorallia. Species within Octocorallia are known for their noxious chemical defenses in addition to the nematocysts present in their tissues [99, 100], and these defenses could help explain why a switch to octocorals has occurred rarely within Cladobranchia.

A caveat does exist, however, in regard to the classification of these clades. Both of the taxa from Arminida included within the present analyses (Armina and Dermatobranchus) are members of the family Arminidae; thus, the monophyly of Arminida (containing Arminidae and Doridomorphidae) as a whole has not yet been rigorously tested. That said, the un-sampled family within Arminida, Doridomorphidae, is monotypic and its sole species lives on the blue coral Heliopora coerulea, an octocoral with a massive calcium carbonate skeleton [14, 101]. This is congruent with the dietary evolution results offered here.


With regard to Dendronotida, the analyses presented here strongly contradict monophyly. The majority of families within Dendronotida (Lomanotidae, Hancockiidae, Dotidae, Bornellidae, Scyllaeidae, Tethyidae, and Dendronotidae) form a single clade that is sister to all other species within Cladobranchia. This clade and the relationships within it are fully supported (BS = 100%) by all analyses presented here. Lomanotidae as sister to the rest of the species within this clade is a result novel to this study, with most previous morphological and molecular analyses [23, 46] supporting alternative topologies, though support was mostly low for these hypotheses.

The rest of the group contains two clades, the first of which is one where Bornellidae is sister to Dotidae and Hancockiidae is sister to the Dotidae + Bornellidae assemblage. This result is also novel as compared to most previous studies [46, 98, 102,103,104,105]. In morphological analyses in particular, both Hancockia and Doto have been considered “problematic” genera [22], and as such have mostly been unplaced (Hancockia) [22] or unassigned (Doto) [47] within the Cladobranchia phylogeny. The analyses presented here, however, very strongly support the position of these genera in the tree, and therefore provide a much stronger hypothesis for their relationships. The second clade within this grouping contains sister groups Scyllaeidae and Tethyidae, as well as Dendronotidae, which is sister to the Scyllaeidae + Tethyidae assemblage. These relationships are consistent with most previous studies [23, 43, 46, 103], but similar to the Dotidae + Bornellidae + Hancockiidae clade, in other cases Scyllaeidae has been previously supported as sister to Dendronotidae, with Tethyidae (usually Melibe specifically) as an early branching lineage [98, 102, 106].

This clade containing the members of “Dendronotida” appears to be almost exclusively composed of taxa that prey on hydrozoans, with the exception of Melibe, which prefers crustaceans that it catches with a remarkable oral hood [107]. Species within Tritoniidae (originally assigned to Dendronotida) form a separate monophyletic group in all analyses as sister to Arminida, as discussed above. In addition, this topology supports the hypothesis that nematocyst sequestration evolved at least twice, because the genus Hancockia (the only non-aeolid genus to sequester nematocysts; Martin et al. 2009) does not form a clade with Aeolidida.


Aeolidida is fully supported as monophyletic, consistent with previous studies [22, 42, 43, 45]. The first of two clades within Aeolidida is made up of taxa from Facelinidae, Aeolidiidae and Flabellinidae. The family Facelinidae forms two separate clades within this group, while Aeolidiidae is monophyletic. The relationships within this clade are consistent with most previous studies using molecular data [43, 46, 98, 108,109,110,111]. Based on these results, Facelinidae should likely be split into two separate families, with one clade retaining the name Facelinidae and the other assigned a more appropriate identifier. However, until a member of the genus Facelina (the type genus for this family) is included in the analyses (ideally the type taxon Facelina auriculata), it is impossible to say which clade should receive the Facelinidae designation. These results also include support for Aeolidiidae as a monophyletic group, at the base of which is one of two shifts to Hexacorallia prey within Aeolidida. The other shift occurs within the family Fionidae.

The second clade within Aeolidida is fully supported across all analyses (BS = 100%), and contains taxa from two families (Fionidae and Unidentiidae) [49]. The relationships between these taxa are also fully supported across all analyses, with the exception of the relationships within Fionidae (though the family itself is monophyletic with full support). Sister to the Fionidae is Unidentiidae. This position for Unidentiidae is a novel result. Only one study previously addressed the phylogenetic position of this family, using morphological data, and in that case Unidentiidae was found to be more closely related to members of Flabellinidae, Piseinotecidae, and Babakinidae [67]. This Fionidae + Unidentiidae clade in particular has multiple shifts to different prey types, including shifts to Crustacea (Fiona) and Hexacorallia (Phestilla).


RNA-Seq data have recovered a well-supported phylogeny for Cladobranchia. The results of this study include a robust hypothesis of relationships between the major cladobranch clades, and indicate that some taxonomically diverse groups, such Dendronotida and Facelinidae, are not monophyletic. The ancestral state reconstruction indicates a strong phylogenetic correlation with prey preference within this group, indicating that host shifts are much more infrequent than previously thought. The mechanism causing evolution of prey preference to be constrained remains unknown, but it could result from difficulties in evolving specific traits for protection against nematocysts from various cnidarian prey groups and chemical compounds from Octocorallia. Future research of Cladobranchia would benefit from combined analyses of prey specialization and prey switching, nematocyst sequestration evolution, and diversification using broader sample coverage. The present study provides a framework for understanding major evolutionary trends in Cladobranchia and indicates that prey type specialization within this group has phylogenetic inertia.


  1. Holt RD, Polis GA. A theoretical framework for intraguild predation. Am Nat. 1997;149:745–64.

  2. Sih A, Englund G, Wooster D. Emergent impacts of multiple predators on prey. Trends Ecol Evol. 1998;13:350–5.

    Article  CAS  PubMed  Google Scholar 

  3. Sotka EE. Local adaptation in host use among marine invertebrates. Ecol Lett. 2005;8:448–59.

    Article  Google Scholar 

  4. Poore AGB, Hill NA, Sotka EE. Phylogenetic and geographic variation in host breadth and composition by herbivorous amphipods in the family Ampithoidae. Evolution. 2008;62:21–38.

    PubMed  Google Scholar 

  5. Mayr E. Animal species and evolution. Eugen Rev. 1963;55:226–8.

    Google Scholar 

  6. Mayr E. Systematics and the origin of species: from the viewpoint of a Zoologist. Nature. 1942;151:347–8.

    Google Scholar 

  7. Mayr E. Geographic speciation in tropical echinoids. Evolution. 1954;8:1–18.

    Article  Google Scholar 

  8. Rocha LA, Robertson DR, Roman J, Bowen BW. Ecological speciation in tropical reef fishes. Proc R Soc London B. 2005;272:573.

    Article  Google Scholar 

  9. Churchill CKC, Alejandrino A, Valdés A, Ó Foighil D. Parallel changes in genital morphology delineate cryptic diversification of planktonic nudibranchs. Proc R Soc B. 2013;280:20131224.

  10. Gosliner TM, Valdés Á, Behrens DW. Nudibranch & Sea Slug Identification. Jacksonville, FL: New World Publications; 2015.

    Google Scholar 

  11. Lalli CM, Gilmer RW. Pelagic snails. The biology of holoplanktonic gastropod mollusks. Stanford, California: Stanford University Press; 1989.

    Google Scholar 

  12. Churchill CKC, Valdés Á, Ó Foighil D. Afro-Eurasia and the Americas present barriers to gene flow for the cosmopolitan neustonic nudibranch Glaucus atlanticus. Mar Biol. 2014;161:899–910.

  13. Valdés Á, Hamann J, Behrens DW, DuPont A. Caribbean Sea Slugs. 1st ed. Gig Harbor, Washington, USA: Sea Challengers Natural History Books; 2006.

    Google Scholar 

  14. Gosliner TM, Behrens DW, Valdés Á. Indo-Pacific Nudibranchs and Sea Slugs: A field guide to the World’s most diverse fauna. 1st ed. Gig Harbor, Washington, USA: Sea Challengers Natural History Books; 2008.

    Google Scholar 

  15. McDonald G, Nybakken JA. preliminary report on a world-wide review of the food of nudibranchs. J Molluscan Stud. 1991;57:61–3.

    Article  Google Scholar 

  16. Wägele H. Potential key characters in Opisthobranchia (Gastropoda, Mollusca) enhancing adaptive radiation. Org Divers Evol. 2004;4:175–88.

    Article  Google Scholar 

  17. Göbbeler K, Klussmann-Kolb A. Molecular phylogeny of the Euthyneura (Mollusca, Gastropoda) with special focus on Opisthobranchia as a framework for reconstruction of evolution of diet. Thalassas. 2011;27:121–54.

    Google Scholar 

  18. Putz A, König GM, Wägele H. Defensive strategies of Cladobranchia (Gastropoda, Opisthobranchia). Nat Prod Rep. 2010;27:1386–402.

    Article  CAS  PubMed  Google Scholar 

  19. McDonald G, Nybakken J. A worldwide review of the food of nudibranch mollusks. Part I. Introduction and the suborder Arminacea. The Veliger. 1997;40:157–9.

    Google Scholar 

  20. McDonald G, Nybakken J. A worldwide review of the food of nudibranch mollusks. Part II. The suborder Dendronotacea. The Veliger. 1999;42:62–6.

  21. Goodheart JA, Bely AE. Sequestration of nematocysts by divergent cnidarian predators: mechanism, function, and evolution. Invertebr Biol. 2017;136:75–91.

    Article  Google Scholar 

  22. Wägele H, Willan RC. Phylogeny of the Nudibranchia. Zool J Linn Soc. 2000;130:83–181.

    Article  Google Scholar 

  23. Pola M, Gosliner TM. The first molecular phylogeny of cladobranchian opisthobranchs (Mollusca, Gastropoda, Nudibranchia). Mol Phylogenet Evol. 2010;56:931–41.

    Article  CAS  PubMed  Google Scholar 

  24. Faucci A, Toonen RJ, Hadfield MG. Host shift and speciation in a coral-feeding nudibranch. Proc R Soc B. 2007;274:111–9.

    Article  CAS  PubMed  Google Scholar 

  25. Wagner D, Kahng SE, Toonen RJ. Observations on the life history and feeding ecology of a specialized nudibranch predator (Phyllodesmium poindimiei), with implications for biocontrol of an invasive octocoral (Carijoa riisei) in Hawaii. J Exp Mar Bio Ecol. 2009;372:64–74.

    Article  Google Scholar 

  26. Cronin G, Hay ME, Fenical W, Lindquist N. Distribution, density, and sequestration of host chemical defenses by the specialist nudibranch Tritonia hamnerorum found at high densities on the sea fan Gorgonia ventalina. Mar Ecol Prog Ser. 1995;119:177–90.

    Article  CAS  Google Scholar 

  27. Ritson-Williams AR, Shjegstad S, Paul V, Ritson-williams R, Shjegstad S, Paul V. Host specificity of four corallivorous Phestilla nudibranchs (Gastropoda: Opisthobranchia). Mar Ecol Prog Ser. 2003;255:207–18.

    Article  Google Scholar 

  28. Affeld S, Kehraus S, Wägele H, König GM. Dietary derived sesquiterpenes from Phyllodesmium lizardensis. J Nat Prod. 2009;72:298–300.

    Article  CAS  PubMed  Google Scholar 

  29. Morrow C, Thorpe J, Picton B. Genetic divergence and cryptic speciation in two morphs of the common subtidal nudibranch Doto coronata (Opisthobranchia: Dendronotacea: Dotoidae) from the northern Irish Sea. Mar Ecol Prog Ser. 1992;84:53–61.

    Article  Google Scholar 

  30. Krug PJ. Patterns of speciation in marine gastropods: A review of the phylogenetic evidence for localized radiations in the sea. Am Malacol Bull. 2011;29:169–86.

    Article  Google Scholar 

  31. Ricklefs RE, Fallon SM. Diversification and host switching in avian malaria parasites. Proc R Soc B. 2002;269:885–92.

    Article  PubMed  PubMed Central  Google Scholar 

  32. Hoberg EP, Brooks DR. A macroevolutionary mosaic: Episodic host-switching, geographical colonization and diversification in complex host-parasite systems. J Biogeogr. 2008;35:1533–50.

    Article  Google Scholar 

  33. McLeish MJ, Noort S. van, Tolley KA. African parasitoid fig wasp diversification is a function of Ficus species ranges. Mol Phylogenet Evol. 2010;57:122–34.

    Article  PubMed  Google Scholar 

  34. Goto R, Kawakita A, Ishikawa H, Hamamura Y, Kato M. Molecular phylogeny of the bivalve superfamily Galeommatoidea (Heterodonta, Veneroida) reveals dynamic evolution of symbiotic lifestyle and interphylum host switching. BMC Evol Biol. 2012;12:172.

    Article  PubMed  PubMed Central  Google Scholar 

  35. Sotka EE, Wares JP, Hay ME. Geographic and genetic variation in feeding preference for chemically defended seaweeds. Evolution. 2003;57:2262–76.

    Article  PubMed  Google Scholar 

  36. Mokady O, Brickner I. Host-Associated Speciation in a Coral-Inhabiting Barnacle. Mol Biol Evol. 2001;18:975–81.

    Article  CAS  PubMed  Google Scholar 

  37. Munday PL, Van Herwerden L, Dudgeon CL. Evidence for sympatric speciation by host shift in the sea. Curr Biol. 2004;14:1498–504.

    Article  CAS  PubMed  Google Scholar 

  38. Mikkelsen PM. Shelled opisthobranchs. Adv Mar Biol. 2002;42:67–136.

    Article  PubMed  Google Scholar 

  39. Mckenna DD, Farrell BD, Caterino MS, Farnum CW, Hawks DC, Maddison DR, et al. Phylogeny and evolution of Staphyliniformia and Scarabaeiformia: forest litter as a stepping stone for diversification of nonphytophagous beetles. Syst Entomol. 2015;40:35–60.

    Article  Google Scholar 

  40. Janz N, Nylin S, Wahlberg N. Diversity begets diversity: host expansions and the diversification of plant-feeding insects. BMC Evol Biol. 2006;6:4.

    Article  PubMed  PubMed Central  Google Scholar 

  41. Johnson KP, Weckstein JD, Meyer MJ, Clayton DH. There and back again: Switching between host orders by avian body lice (Ischnocera: Goniodidae). Biol J Linn Soc. 2011;102:614–25.

    Article  Google Scholar 

  42. Wollscheid-Lengeling E, Boore J, Brown W, Wägele H. The phylogeny of Nudibranchia (Opisthobranchia, Gastropoda, Mollusca) reconstructed by three molecular markers. Org Divers Evol 2001;1:241–56.

    Article  Google Scholar 

  43. Goodheart JA, Bazinet AL, Collins AG, Cummings MP. Relationships within Cladobranchia (Gastropoda: Nudibranchia) based on RNA-seq data: An initial investigation. R Soc Open Sci. 2015;2:150196.

    Article  PubMed  PubMed Central  Google Scholar 

  44. Goodheart J, Bazinet A, Collins A, Cummings M. Phylogeny of Cladobranchia (Gastropoda: Nudibranchia): a total evidence analysis using DNA sequence data from public databases. Maryl: Digit. Repos. Univ; 2015.

    Google Scholar 

  45. Wollscheid E, Wägele H. Initial results on the molecular phylogeny of the Nudibranchia (Gastropoda, Opisthobranchia) based on 18S rDNA data. Mol Phylogenet Evol. 1999;13:215–26.

    Article  CAS  PubMed  Google Scholar 

  46. Mahguib J, Valdés Á. Molecular investigation of the phylogenetic position of the polar nudibranch Doridoxa (Mollusca, Gastropoda, Heterobranchia). Polar Biol. 2015;38:1369.

    Article  Google Scholar 

  47. Valdés Á, Bouchet P. Cephalaspidea, Thecosomata, Gymnosomata, Aplysiomorpha, Umbraculida, Acochlidiacea, Sacoglossa, Cylindobullida, Nudipleura. Classif Nomencl Gastropod Fam Malacol. 2005;47:1–397.

  48. McDonald GR, Nybakken JW. List of the worldwide food habits of nudibranchs. Available from:

  49. Cella K, Carmona L, Ekimova I, Chichvarkhin A, Schepetov D, Gosliner TM. A radical solution: The phylogeny of the nudibranch family Fionidae.PLoS One. 2016;11:e0167800.

    Article  PubMed  PubMed Central  Google Scholar 

  50. Zapata F, Wilson NG, Howison M, Andrade SCS, Jörger KM, Schrödl M, et al. Phylogenomic analyses of deep gastropod relationships reject Orthogastropoda. Proc R Soc B. 2014;20141739:281.

    Google Scholar 

  51. Grabherr MG, Haas BJ, Yassour M, Levin JZ, Thompson DA, Amit I, et al. Full-length transcriptome assembly from RNA-Seq data without a reference genome. Nat Biotechnol. 2011;29:644–52.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  52. Ebersberger I, Strauss S, von Haeseler A. HaMStR: profile hidden markov model based search for orthologs in ESTs. BMC Evol Biol. 2009;9:157.

    Article  PubMed  PubMed Central  Google Scholar 

  53. Pearson WR, Lipman DJ. Improved tools for biological sequence comparison. Proc Natl Acad Sci. 1988;85:2444–8.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  54. Birney E, Clamp M, Durbin R. GeneWise and Genomewise. Genome Res. 2004;14:988–95.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  55. Eddy SR. Accelerated profile HMM searches. PLoS Comput Biol. 2011;7:e1002195.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  56. Altschul SF, Gish W, Miller W, Myers EW, Lipman DJ. Basic local alignment search tool. J Mol Biol. 1990;215:403–10.

    Article  CAS  PubMed  Google Scholar 

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

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  58. Zwickl DJ. Genetic algorithm approaches for the phylogenetic analysis of large biological sequence datasets: University of Texas at Austin; 2006.

  59. Kocot KM, Citarella MR, Moroz LL, Halanych KM. PhyloTreePruner: A phylogenetic tree-based approach for selection of orthologous sequences for phylogenomics. Evol Bioinforma. 2013;2013:429–35.

    Article  Google Scholar 

  60. Bazinet AL, Zwickl DJ, Cummings MP. A gateway for phylogenetic analysis powered by grid computing featuring GARLI 2.0. Syst Biol. 2014;63:812–8.

    Article  PubMed  PubMed Central  Google Scholar 

  61. Tavaré S. Some probabilistic and statistical problems in the analysis of DNA sequences. Lect Math Life Sci. 1986;17:57–86.

  62. Hasegawa M, Kishino H, Yano T. Dating of the human-ape splitting by a molecular clock of mitochondrial DNA. J Mol Evol. 1985;22:160–74.

    Article  CAS  PubMed  Google Scholar 

  63. Yang Z. Maximum-likelihood estimation of phylogeny from DNA sequences when substitution rates differ over sites. Mol Biol Evol. 1993;10:1396–401.

    CAS  PubMed  Google Scholar 

  64. Sukumaran J, Holder MT. DendroPy: a Python library for phylogenetic computing. Bioinformatics. 2010;26:1569–71.

    Article  CAS  PubMed  Google Scholar 

  65. R Development Core Team. R: A Language and Environment for Statistical Computing. R Found. Stat. Comput. Vienna Austria. 2016. Available from:

  66. Picton BE, Morrow CC. A field guide to the nudibranchs of the British Isles. London: Immel Publishing Limited; 1994.

    Google Scholar 

  67. Millen S, Hermosillo A. Three new species of aeolid nudibranchs (Opisthobranchia) from the Pacific Coast of Mexico, Panama, and the Indopacific, with a redescription and redesignation of a fourth species. Veliger. 2012;51:145–64.

    Google Scholar 

  68. Rudman WB. Spurilla neapolitana (Delle Chiaje, 1823). Sea Slug Forum. 1999; Accessed 23 Nov 2016.

  69. Rudman WB. Lomanotus vermiformis Eliot, 1908. Sea Slug Forum. 1999; Accessed 23 Nov 2016.

  70. Rudman WB. Tritoniopsis frydis Marcus & Marcus 1970. Sea Slug. Forum. 2001; Accessed 23 Nov 2016.

  71. Rudman WB. Predators of Sea Anemones. Sea Slug Forum. 2002; Accessed 23 November 2016.

  72. Mazerolle MJ, AICmodavg. Model selection and multimodel inference based on (Q)AIC(c). R package version. 2016;2:1–0.

    Google Scholar 

  73. Paradis E, Claude J, Strimmer KAPE. Analyses of phylogenetics and evolution in R language. Bioinformatics. 2004;20:289–90.

    Article  CAS  PubMed  Google Scholar 

  74. Goddard JHR, Gosliner TM, Pearse JS. Impacts associated with the recent range shift of the aeolid nudibranch Phidiana hiltoni (Mollusca, Opisthobranchia) in California. Mar Biol. 2011;158:1095–109.

    Article  PubMed  PubMed Central  Google Scholar 

  75. Greenwood PG. Acquisition and use of nematocysts by cnidarian predators. Toxicon. 2009;54:1065–70.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  76. Mauch S, Elliott J. Protection of the nudibranch Aeolidia papillosa from nematocyst discharge of the sea anemone Anthopleura elegantissima. Veliger. 1997;40:148–51.

    Google Scholar 

  77. Rudman W, Willan RC. Opisthobranchia, Introduction. In: Beesley PL, GJB R, Wells A, editors. Introduction. Mollusca South. Synth: Melbourne: CSIRO Publishing; 1998.p. 915–42.

  78. Thompson TE. Biology of Opisthobranch Molluscs, vol. Volume 1. London: Ray Society; 1976.

    Google Scholar 

  79. Appeltans W, Ahyong ST, Anderson G, Angel MV, Artois T, Bailly N, et al. The magnitude of global marine species diversity. Curr Biol. 2012;22:2189–202.

    Article  CAS  PubMed  Google Scholar 

  80. Martin R, Heß M, Schrödl M, Tomaschko K-H. Cnidosac morphology in dendronotacean and aeolidacean nudibranch molluscs: From expulsion of nematocysts to use in defense? Mar Biol. 2009;156:261–8.

    Article  Google Scholar 

  81. Fautin DG. Structural diversity, systematics, and evolution of cnidae. Toxicon. 2009;54:1054–64.

    Article  CAS  PubMed  Google Scholar 

  82. Bowen BW, Rocha LA, Toonen RJ, Karl SA. The origins of tropical marine biodiversity. Trends Ecol Evol. 2013;28:359–66.

    Article  PubMed  Google Scholar 

  83. Dijkstra K-DB, Monaghan MT, Pauls SU. Freshwater biodiversity and aquatic insect diversification. Annu Rev Entomol. 2014;59:143–63.

    Article  CAS  PubMed  Google Scholar 

  84. Via S. Sympatric speciation in animals: The ugly duckling grows up. Trends Ecol Evol. 2001;16:381–90.

    Article  CAS  PubMed  Google Scholar 

  85. Bush GL. Sympatric Host Race Formation and Speciation in Frugivorous Flies of the Genus Rhagoletis (Diptera, Tephritidae). Evolution. 1969;23:237–51.

    Article  PubMed  Google Scholar 

  86. Stanhope MJ, Connelly MM, Hartwick B. Evolution of a crustacean chemical communication channel: Behavioral and ecological genetic evidence for a habitat-modified, race-specific pheromone. J Chem Ecol. 1992;18:1871–87.

    Article  CAS  PubMed  Google Scholar 

  87. Schiaparelli S, Alvaro MC, Kilgallen N, Scinto A, Lörz AN. Host-shift speciation in Antarctic symbiotic invertebrates: further evidence from the new amphipod species Lepidepecreella debroyeri from the Ross Sea? Hydrobiologia. 2015;761:143–59.

  88. Duffy JE. Resource-associated population subdivision in a symbiotic coral-reef shrimp. Evolution. 1996;50:360–73.

    Article  PubMed  Google Scholar 

  89. Hurt C, Silliman K, Anker A, Knowlton N. Ecological speciation in anemone-associated snapping shrimps (Alpheus armatus species complex). Mol Ecol. 2013;22:4532–48.

    Article  CAS  PubMed  Google Scholar 

  90. Sanford E, Roth MS, Johns GC, Wares JP, Somero GN. Local selection and latitudinal variation in a marine predator-prey interaction. Science. 2003;300:1135–7.

  91. Jensen KR. Evolution of the Sacoglossa (Mollusca, Opisthobranchia) and the ecological associations with their food plants. Evol Ecol. 1997;11:301–35.

    Article  Google Scholar 

  92. Trowbridge CD, Todd CD. Host-plant change in marine specialist herbivores: Ascoglossan sea slugs on introduced macroalgae. Ecol Monogr. 2001;71:219–43.

    Article  Google Scholar 

  93. Coyne JA, Orr HA. Speciation. Sunderland, MA: Sinauer Associates, Inc.; 2004.

    Google Scholar 

  94. Todd CD, Lambert WJ, Davies J. Some perspectives on the biology and ecology of nudibranch mollusks: generalisations and variations on the theme prove the rule. Boll Malacol. 2001;37:105–20.

    Google Scholar 

  95. Thollesson M. Phylogenetic analysis of Euthyneura (Gastropoda) by means of the 16S rRNA gene: use of a “fast” gene for “higher-level” phylogenies. Proc R Soc B. 1999;266:75.

    Article  PubMed Central  Google Scholar 

  96. Odhner N. The Nudibranchiata. British Antarctic (Terra Nova) Expedition, 1910. Br Museum Nat Hist Rep. 1934;7:229–310.

    Google Scholar 

  97. Gosliner TM, Fahey SJ. Previously undocumented diversity and abundance of cryptic species: a phylogenetic analysis of Indo-Pacific Arminidae Rafinesque, 1814 (Mollusca: Nudibranchia) with descriptions of 20 new species of Dermatobranchus. Zool J Linn Soc. 2011;161:245–356.

    Article  PubMed  PubMed Central  Google Scholar 

  98. Hulett RE, Mahguib J, Gosliner TM, Valdes A. Molecular evaluation of the phylogenetic position of the enigmatic species Trivettea papalotla (Bertsch et al. 2009) (Mollusca: Nudibranchia). Invertebr Syst. 2015;29:215–22.

    Article  Google Scholar 

  99. Almeida MTR, Moritz MIG, Capel KCC, Pérez CD, Schenkel EP. Chemical and biological aspects of octocorals from the Brazilian coast. Rev Bras Farmacogn. 2014;24:446–67.

    Article  CAS  Google Scholar 

  100. Coll JC. The chemistry and chemical ecology of octocorals (Coelenterata, Anthozoa, Octocorallia). Chem Rev. 1992;92:613–31.

    Article  CAS  Google Scholar 

  101. Hoon TC, Ong R. New record of the nudibranch Doridomorpha gardineri in Singapore. Singapore Biodivers Rec. 2015;2015:151–3.

    Google Scholar 

  102. Shipman C, Gosliner T. Molecular and morphological systematics of Doto Oken, 1851 (Gastropoda: Heterobranchia), with descriptions of five new species and a new genus. Zootaxa. 2015;3973:057–101.

    Article  Google Scholar 

  103. Pola M, Camacho-García YE, Gosliner TM. Molecular data illuminate cryptic nudibranch species: the evolution of the Scyllaeidae (Nudibranchia: Dendronotina) with a revision of Notobryon. Zool J Linn Soc. 2012;165:311–36.

    Article  Google Scholar 

  104. Ekimova I, Korshunova T, Schepetov D, Neretina T, Sanamyan N, Martynov A. Integrative systematics of northern and Arctic nudibranchs of the genus Dendronotus (Mollusca, Gastropoda), with descriptions of three new species. Zool J Linn Soc. 2015;173:841–86.

  105. Pola M, Rudman WB, Gosliner TM. Systematics and preliminary phylogeny of Bornellidae (Mollusca: Nudibranchia: Dendronotina) based on morphological characters with description of four new species. Zootaxa. 2009;1975:1–57.

    Google Scholar 

  106. Stout CC, Pola M, Valdes A. Phylogenetic analysis of Dendronotus nudibranchs with emphasis on northeastern Pacific species. J Molluscan Stud. 2010;76:367–75.

  107. Espinoza E, Dupont A, Valdés Á. A tropical Atlantic species of Melibe Rang, 1829 (Mollusca, Nudibranchia, Tethyiidae). Zookeys. 2013;66:55–66.

    Google Scholar 

  108. Carmona L, Gosliner TM, Pola M, Cervera JL. A molecular approach to the phylogenetic status of the aeolid genus Babakina Roller, 1973 (Nudibranchia). J Molluscan Stud. 2011;77:417–22.

    Article  Google Scholar 

  109. Carmona L, Pola M, Gosliner TM, Cervera JL. A tale that morphology fails to tell: A molecular phylogeny of Aeolidiidae (Aeolidida, Nudibranchia, Gastropoda). PLoS One. 2013;8:e63000.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  110. Carmona L, Pola M, Gosliner TM, Cervera JL. The Atlantic-Mediterranean genus Berghia Trinchese, 1877 (Nudibranchia: Aeolidiidae): Taxonomic review and phylogenetic analysis. J Molluscan Stud. 2014;80:1–17.

    Article  Google Scholar 

  111. Carmona L, Pola M, Gosliner TM, Cervera JL. Burnaia Miller, 2001 (Gastropoda, Heterobranchia, Nudibranchia): a facelinid genus with an Aeolidiidae’s outward appearance. Helgol Mar Res. 2015;69:285–91.

    Article  Google Scholar 

Download references


We are grateful to Craig Hoover (California State Polytechnic University, Pomona), Hans Bertsch (Universidad Autónoma de Baja California), Karen Cheney (University of Queensland), Allison Fritts-Penniman (California Academy of Sciences), David Fenwick III, and Ariane Dimitris for providing specimens and/or collecting assistance; and Paul Katz, Jonathan Boykin, Aastha Vashist and Amirah Hurst (Georgia State University) for providing two of the RNA-Seq datasets for analysis. We would also like to thank Vanessa Gonzalez at the Smithsonian National Museum of Natural History for help with RNA extractions. Finally, we are grateful to the Laboratories of Analytical Biology of the National Museum of Natural History for use of the laboratory facilities, and the staff of the Smithsonian Tropical Research Institute in Panama and the Richard B. Gump South Pacific Research Station in French Polynesia for use of their facilities and their help in acquiring the proper permits. Lastly, we want to thank four anonymous reviewers for their detailed and constructive criticism on previous versions of this manuscript.


This work was supported by a Smithsonian Institution Peter Buck Pre-doctoral Fellowship, Conchologists of America, Society of Systematic Biologists, University of Maryland Graduate School Dean’s Fellowship and Summer Research Fellowship to JAG, Smithsonian Institution Small Grant award to AGC, funding from University of Maryland to MPC, and NSF Partnerships for International Research and Education program Award 1243541.

Availability of data and materials

RNA-Seq sequence data can be accessed at the NCBI Sequence Read Archive (SRA) with the following accessions: SRA accession numbers SRR1505104, SRR1505108, SRR1505109, SRR1505130, SRR1719366, SRR1721590, SRR1950939–SRR1950954, SRR3726692–SRR3726707, SRR4124996 and SRR4190242. Aligned data matrices and tree files can be accessed at the Dryad Digital Repository (DOI:10.5061/dryad.7kh2n).

Author information

Authors and Affiliations



JAG, AGC, AV and MPC conceived of the study; JAG and AV collected samples and field data; JAG carried out the molecular lab work; JAG and ALB performed the bioinformatics analyses; all authors participated in study design and data analysis, helped draft the manuscript, and gave final approval for publication.

Corresponding author

Correspondence to Jessica A. Goodheart.

Ethics declarations

Ethics approval and consent to participate

Specimen collection was conducted under permits issued by the California Department of Fish and Wildlife (Scientific Collecting Permit #SC-13009), Florida Fish and Wildlife Conservation Commission (Special Activity License #SAL-14-1565-SR), Délégation Régionale à la Recherché et à la Technologie de la Polynesie Francaise (no number for permit), Republica de Panama Ministerio de Ambiente (Permiso Científico #SE/A-64-15 and #SEX/A-65-15), Queensland Government (General fisheries permit #161624), and Oficio Dirección General de Ordenamiento Pesquera y Acuicola (Permit #01020/130214).

Consent for publication

Not applicable.

Competing interests

The authors declare that they have no competing interests.

Publisher’s Note

Springer Nature remains neutral with regard to jurisdictional claims in published maps and institutional affiliations.

Additional files

Additional file 1: Figure S1.

Select photographs of dendronotid and unassigned taxa used in this project, including: A) Scyllaea fulva (SRR3726701), B) Dermatobranchus sp. (SRR3726698; Photo credit: Karen Cheney), C) Lomanotus vermiformis (SRR3726706) and D) Hancockia uncinata (Photo credit: David Fenwick III). (TIFF 19274 kb)

Additional file 2: Figure S2.

Select photographs of aeolid taxa used in this project, including: A) Phidiana lynceus, B) Eubranchus rustyus (SRR3726692; Photo credit: Craig Hoover), C) Learchis evelinae (SRR3726693), and D) Spurilla braziliana. (TIFF 14185 kb)

Additional file 3: Table S5.

List of specimens examined in this study, including species name, locality, and morphological and molecular tissue voucher and barcode information. Sequence Read Archive accession numbers are also provided for each RNA-Seq dataset. (XLSX 41 kb)

Additional file 4: Table S4.

Prey preference data used for the ancestral state reconstruction. (XLSX 48 kb)

Additional file 5: Table S1.

Table of sequence read information for each sample. (XLSX 55 kb)

Additional file 6: Table S2.

Trinity-assembly details, including number of transcript fragments and total number of bases assembled, as well as N50 and L50 statistics for each transcriptome. (XLSX 42 kb)

Additional file 7: Table S3.

HaMStR statistics for each RNA-Seq dataset. (XLSX 43 kb)

Additional file 8: Table S6.

Ancestral state reconstruction results for the evolution of diet preference in Cladobranchia with the alternative prey type states. (XLSX 51 kb)

Additional file 9: Table S7.

Ancestral state reconstruction results for the evolution of diet preference in Cladobranchia with the alternative prey type state for Dirona picta. (XLSX 54 kb)

Additional file 10: Table S8.

Ancestral state reconstruction results for the evolution of diet preference in Cladobranchia with the alternative prey type state for Fiona pinnata. (XLSX 54 kb)

Additional file 11: Table S9.

Ancestral state reconstruction results for the evolution of diet preference in Cladobranchia with the alternative prey type state for Janolus barbarensis. (XLSX 54 kb)

Additional file 12: Table S10.

Ancestral state reconstruction results for the evolution of diet preference in Cladobranchia with the alternative prey type state for Phidiana lynceus. (XLSX 54 kb)

Rights and permissions

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

Reprints and permissions

About this article

Check for updates. Verify currency and authenticity via CrossMark

Cite this article

Goodheart, J.A., Bazinet, A.L., Valdés, Á. et al. Prey preference follows phylogeny: evolutionary dietary patterns within the marine gastropod group Cladobranchia (Gastropoda: Heterobranchia: Nudibranchia). BMC Evol Biol 17, 221 (2017).

Download citation

  • Received:

  • Accepted:

  • Published:

  • DOI: