Specimens were collected by SCUBA diving in West Papua by Gerard van der Velde in 2007, mostly in the coastal areas of Gam, Kri, Mansuar, and Batanta (Figures 1–2; see Hoeksema and van der Meij 2008). Additional specimens were mainly collected by Joris van Alphen and Nicole de Voogd, and also by Bert Hoeksema, Sancia van der Meij, and other expedition members (Hoeksema and van der Meij 2010) in 2009 off Halmahera (northern Moluccas), especially around Ternate (Figures 1, 3). A locality list of the sampling stations is provided in Table 1. Collected slugs were first photographed and subsequently preserved in 96% ethanol (West Papua 2007). Halmahera specimens were transferred into fresh 96% ethanol and labelled in order to prepare them for DNA analysis. These have been deposited in the mollusc collection of Naturalis Biodiversity Center, Leiden (coded as RMNH.Mol.), with the exception of some specimens that dried out after sequencing (Table 1; Figures 5–15; Suppl. material 1: COI sequences).

Figure 1. 

Location of field areas: Halmahera (including Ternate) and West Papua (including Raja Ampat).

Figure 2. 

Raja Ampat sites where Phyllidiidae were sampled in 2007.

Figure 3. 

Halmahera and Ternate sites where Phyllidiidae were sampled in 2009.

Collected specimens were identified according to their external morphology using Brunckhorst (1993), Yonow et al. (2002), and Yonow (2011). In addition, field guides showing in situ photographs were used (Gosliner et al. 2008). All individuals except for three could be identified to species level. All specimens were photographed alive or in the preserved state (Figures 5–15); these photos can be linked to the phylogeny reconstruction of the Phyllidiidae based on COI gene sequence data (Figure 4).

Figure 4. 

Phylogeny reconstruction of the Phyllidiidae based on COI gene sequence data of 109 specimens (including outgroups). Topology derived from Bayesian inference 50% majority rule, significance values are posterior probabilities / bootstrap values. Numbers refer to GenBank accession numbers / RMNH.Moll catalogue numbers.

Figure 5. 

External morphology and colouration of Phyllidiidae specimens used for COI phylogeny reconstruction: Phyllidia elegans. Order of specimens (a–h) according to Figure 4 (f, h dorsal and ventral sides). Numbers refer to RMNH. Moll catalogue numbers.

Figure 6. 

External morphology and colouration of Phyllidiidae specimens used for COI phylogeny reconstruction: Phyllidia elegans. Order of specimens (a–i) according to Figure 4 (d dorsal and ventral sides). Numbers refer to RMNH.Moll catalogue numbers and locality codes (137 and 156, dried-out).

Figure 7. 

External morphology and colouration of Phyllidiidae specimens used for COI phylogeny reconstruction: Phyllidia elegans (a–f), Phyllidia sp. (g dorsal and ventral sides), P. exquisita (h), P. coelestis (i). Order of specimens (a–i) according to Figure 4. Numbers refer to RMNH.Moll catalogue numbers or locality code (058, dried-out).

Figure 8. 

External morphology and colouration of Phyllidiidae specimens used for COI phylogeny reconstruction: Phyllidia coelestis (a), P. varicosa (b–i). Order of specimens (a–i) according to Figure 4 (d dorsal and ventral sides). Numbers refer to RMNH.Moll catalogue numbers.

Figure 9. 

External morphology and colouration of Phyllidiidae specimens used for COI phylogeny reconstruction: Phyllidia varicosa (a–f), P. ocellata (g–i). Order of specimens (a–i) according to Figure 4 (c dorsal and ventral sides). Numbers refer to RMNH.Moll catalogue numbers or locality code (074, dried-out).

Figure 10. 

External morphology and colouration of Phyllidiidae specimens used for COI phylogeny reconstruction: Phyllidia picta (a–c), Phyllidia babai (d), Phyllidia cf. babai (e–f), Reticulidia fungia (g), Reticulidia halgerda (h), Phyllidiopsis fissuratus (i). Order of specimens (a–i) according to Figure 4 (e dorsal and ventral sides). Numbers refer to RMNH.Moll catalogue numbers.

Figure 11. 

External morphology and colouration of Phyllidiidae specimens used for COI phylogeny reconstruction: Phyllidiella rudmani (a), Phyllidiella nigra (b–h), Phyllidiella pustulosa (i–j). Order of specimens (a–j) according to Figure 4. Numbers refer to RMNH.Moll catalogue numbers.

Figure 12. 

External morphology and colouration of Phyllidiidae specimens used for COI phylogeny reconstruction: Phyllidiella pustulosa. Order of specimens (a–j) according to Figure 4. Numbers refer to RMNH.Moll catalogue numbers or locality code (75F, dried-out).

Figure 13. 

External morphology and colouration of Phyllidiidae specimens used for COI phylogeny reconstruction: Phyllidiella pustulosa (a–h), Phyllidiopsis xishaensis (i–j). Order of specimens (a–j) according to Figure 4. Numbers refer to RMNH.Moll catalogue numbers.

Figure 14. 

External morphology and colouration of Phyllidiidae specimens used for COI phylogeny reconstruction: Phyllidiopsis xishaensis (a), Phyllidiopsis shireenae (b–c), Phyllidiopsis krempfi (d–i). Order of specimens (a–i) according to Figure 4 (c dorsal and ventral sides). Numbers refer to RMNH.Moll catalogue numbers.

Figure 15. 

External morphology and colouration of Phyllidiidae specimens used for COI phylogeny reconstruction: Phyllidiopsis krempfi. Order of specimens (a–g) according to Figure 4 (f, g dorsal and ventral sides). Numbers refer to RMNH.Moll catalogue numbers.

For each species encountered in the field surveys one or more individuals were chosen for DNA analysis as well as from the morphologically distinct unidentified specimens, resulting in a total of 99 samples (Table 1). DNA was extracted from tissue of small foot fragments with the DNeasy Blood & Tissue Kit (Qiagen, Germany) according to the manufacturer’s protocol. DNA was eluted in DEPC treated water. The quality of the extracted DNA was tested by agarose gel (0.7%) electrophoresis.

Extracted DNA was used for Polymerase Chain Reaction (PCR) to amplify fragments of the mitochondrial gene COI (cytochrome c oxidase subunit 1). The primers used for the amplification of the COI gene were: LCO1490 (5’GGT CAA CAA ATC ATA AAG ATA TTG G 3’) and HCO2198 (5’TAA ACT TCA GGG TGA CCA AAA AAT CA 3’) (Folmer et al. 1994). Thermal cycling conditions used for the amplification of the COI gene were: initial denaturing at 94 °C for 3 min followed by 38 amplification cycles of denaturation at 94 °C for 15 sec, primer annealing at 50 °C for 30 sec, and elongation at 72 °C for 1 min. A final elongation step at 72 °C for 5 min was performed. After checking by agarose (1%) electrophoresis if the PCR resulted the unique PCR fragments of the expected size (approximately 658 bp), the fragments were purified using the GeneJET PCR Purification Kit (Thermo Scientific, Landsmeer, NL). Purified PCR products were sequenced with corresponding primers.

The quality of the sequences was checked using Chromas Lite (Technelysium Pty Ltd.). Subsequently the sequences were edited in MEGA 6 (Tamura et al. 2013) and analysed by BLAST searches (http://www.ncbi.nlm.nih.gov). COI sequences of Dendrodoris citrina (Cheeseman, 1881) and Doriopsilla areolata Bergh, 1880 were collected from GenBank and used as outgroups. Additional COI sequences of Phyllidia coelestis Bergh, 1905, Phyllidia elegans Bergh, 1869, Phyllidia ocellata Cuvier, 1804, Phyllidia picta Pruvot-Fol, 1957, Phyllidia varicosa Lamarck, 1801, Phyllidiella lizae Brunckhorst, 1993, Phyllidiella pustulosa (Cuvier, 1804), Phyllidiopsis cardinalis Bergh, 1875 were obtained from GenBank (Table 2).

The newly obtained COI sequences and the sequences from GenBank were aligned using the Guidance server (Clustal W; Penn et al. 2010), resulting in an alignment score of 1.000. There were no unreliable columns. Prior to the model-based phylogenetic analysis, the best-fit model of nucleotide substitution was identified by means of the Akaike Information Criterion (AIC) calculated with jModeltest (Posada 2008), resulting in TVM+I+G as the most suitable model. Phylogenetic reconstructions were carried out with Bayesian inference in MrBayes 3.1.2 (Ronquist and Huelsenbeck 2003) using the most complex GTR+I+G model of nucleotide substitution. Bayesian inference coupled with Markov Chain Monte Carlo techniques (MCMC; six chains) were run for 5,000,000 generations with a sample tree saved every 1000 generations. The burnin was set to 25%. Likelihood scores stabilized at 0.007476. Consensus trees were visualized in FigTree v.1.3.1 (Rambaut 2009). A maximum likelihood analysis (GTR+I+G; 1000 bootstraps) was carried out with Phyml 3.1 (Guindon et al. 2010) using the Seaview platform (Gouy et al. 2010).

Initial phylogenetic analyses showed high intraspecific variation on the COI region between specimens identified as Phyllidiella pustulosa. Tests to estimate the average evolutionary divergence over sequence pairs between and within groups were carried out in MEGA 6.06. Phyllidia elegans, P. varicosa, Phyllidiella nigra (van Hasselt, 1824), P. pustulosa, and Phyllidiopsis krempfi Pruvot-Fol, 1957 were used as representatives for each of the species groups, because of the larger number of available sequences for these species. The Phyllidiella pustulosa sequence from GenBank (KJ001310) was excluded from this analysis: based on its position in the phylogeny reconstruction the identification of this specimen as P. pustulosa is doubtful. The web version of ABGD (Automatic Barcode Gap Discovery, Puillandre et al. 2012) was used to estimate the genetic distance corresponding to the difference between a speciation process versus intra-specific variation in Phyllidiella pustulosa. Runs were performed using the default range of priors (pmin = 0.001, pmax = 0.10) using the JC69 Jukes-Cantor measure of distance. The analysis involved 20 nucleotide sequences with a total of 588 positions in the final dataset.

All available mitochondrial 16S sequences of Phyllidiidae on GenBank (Tholesson 2000, Wolfscheid-Lengeling et al. 2001, Valdés 2003, Cheney et al. 2014, Shields et al. unpublished) were used for a phylogeny reconstruction based on this marker, which allowed us to study the phylogenetic position of 17 phyllidiid species including two species (Phyllidia rueppelii (Bergh, 1869) and Phyllidiopsis sphingis Brunckhorst, 1993) for which no COI data were available. Doriopsilla albopunctata (JG Cooper, 1863) was used as outgroup (Table 3). The sequences were aligned using the Guidance server (ClustalW; Penn et al. 2010), resulting in an alignment score of 0.996281. All unreliable columns (confidence score below 0.93) were removed. Prior to the model-based phylogenetic analysis, the best-fit model of nucleotide substitution was identified by means of the Akaike Information Criterion (AIC) calculated with jModeltest (Posada 2008), resulting in TVM+I+G. Because of the unavailability of TVM in MrBayes 3.1.2 (Ronquist and Huelsenbeck 2003), we used the most complex GTR+I+G model of nucleotide substitution. Bayesian inferences coupled with MCMC techniques (six chains) were run for 3,000,000 generations, with a sample tree saved every 1000 generations and the burnin set to 25%. Likelihood scores stabilized at a value of 0.005654. Consensus trees were visualized in FigTree v.1.3.1 (Rambaut 2009). A maximum likelihood analysis (GTR+I+G; 1000 bootstraps) was carried out with Phyml 3.1 (Guindon et al. 2010) using the Seaview platform (Gouy et al. 2010).

The reconstruction based on COI (Figure 4) is derived from the Bayesian inference 50% majority rule consensus. This topology is congruent with the one resulting from the maximum likelihood analysis. Three large groupings can be discerned (indicated as A, B, and C in Figure 4), albeit with low support for the higher taxonomic levels. The support values in the distal branches are high. The genera Phyllidia, Phyllidiella, Phyllidiopsis, and Reticulidia are retrieved in distinct clades, with Reticulidia as a sister clade to Phyllidia. Phyllidiopsis fissuratus Brunckhorst, 1993 formed a separate lineage basal to Phyllidiella species (albeit without support). Phyllidiopsis cardinalis does not cluster with its congeners, but instead forms a separate lineage in the Phyllidiidae.

The 16S phylogeny reconstruction is also derived from the Bayesian inference 50% majority rule consensus of the trees remaining after the burnin. There are low support values in the basal part of the tree and high support values in the distal phylogenetic branches (Figure 17). The Bayesian inference topology is congruent with the topology resulting from the maximum likelihood analysis. The outgroup Doriopsilla albopunctata is separated by a long branch. Within the overall clade four main groupings can be distinguished: Phyllidiella, Phyllidiopsis, and Reticulidia, and a mixed clade of Phyllidiella and Phyllidia. Based on this analysis only the genus Reticulidia is monophyletic. Phyllidiopsis cardinalis does not cluster with any of the other analysed taxa, and holds a separate position in the phylogeny reconstruction. The latter is in accordance with the COI reconstruction (Figure 4).

The arrangement of the four phyllidiid genera based on the molecular data (Figures 4, 16a) is similar to that of Brunckhorst (1993) that was based on morphological and anatomical data (Figure 16b). The only exception is the position of the genus Fryeria. Brunckhorst (1993) distinguished Fryeria from Phyllidia based on the position of the anus and other anatomical features. Phyllidia picta (with its synonyms Fryeria picta (Pruvot-Fol, 1957), Fryeria menindie Brunckhorst, 1993, Phyllidia menindie (Brunckhorst, 1993)) was included in our analyses which, according to Brunckhorst, should belong to the genus Fryeria. Valdés and Gosliner (1999) synonymized both genera, which was later followed by Valdés (2003) and Cheney et al. (2014). The present reconstruction based on COI (Figure 16a) reconfirms the inclusion of Fryeria in the genus Phyllidia.

Figure 16. 

a Cladogram based on COI gene sequence data showing topology of four genera of Phyllidiidaeb Cladogram according to Brunckhorst (1993) based on morphological data showing topology of six genera of Phyllidiidaec Cladogram based on 16S mtDNA sequence data showing topology of four genera of Phyllidiidae (Valdés 2003) d Cladogram based on morphological data (Valdés 2002) showing topology of five genera of Phyllidiidae.

Figure 17. 

Phylogeny reconstruction of the Phyllidiidae based on 16S mtDNA of 17 specimens of 14 species (including outgroup). Topology derived from Bayesian inference 50% majority rule, significance values are posterior probabilities/bootstrap values. Numbers refer to GenBank accession numbers. *Re-identification according to Yonow (pers. comm.)

The cladogram of the genera based on 16S mtDNA sequence data collected by Valdés (2003) (Figure 16c) is roughly similar to the cladogram based on COI, except for the different positions of Phyllidiopsis and Phyllidiella. The cladogram based on morphological and anatomical data as shown by Valdés (2002; Figure 16d) is different from the other proposed classifications (Figures 16a–c). Brunckhorst (1993) considered Ceratophyllidia a sister group to all the other genera (Figure 6b). Valdés (2002; Figure 16d) distinguished two larger groupings within the Phyllidiidae; Ceratophyllidia and Phyllidiopsis as one group and Phyllidia, Phyllidiella, and Reticulidia as the other group. Phyllidia and Phyllidiella in turn formed a sister group of Reticulidia (Figure 16d). The cladogram by Brunckhorst (1993) and our cladogram based on COI (Figure 4) both show that Phyllidiella is a sister clade of Reticulidia and Phyllidia. In contrast, Phyllidiella is not a sister group of Phyllidia but to all the other genera grouped together in the cladogram of Valdés (2003).

Unfortunately no Ceratophyllidia specimens were available to complete our analysis at genus level. Up to this point the phylogenetic position of the genus Ceratophyllidia remains unclear, and additional molecular analyses are necessary to establish its position.

Species level analysis was mainly based on COI (Figure 4). Four nominal species were sequenced in the genus Phyllidiella. Phyllidiella nigra formed a highly supported clade. In the clade containing P. pustulosa much variation is visible indicating larger genetic differences among individuals. The ABGD analysis shows that four Molecular Operational Taxonomic Units (MOTUs) are present in Phyllidiella pustulosa, suggesting the presence of cryptic species or, alternatively, high intraspecific variation. The P. pustulosa of Cheney et al. (2014) falls in between the group consisting of P. nigra and P. pustulosa on one side and P. rudmani Brunckhorst, 1993 on the other and probably represents another species. Our specimen of P. rudmani clustered with the specimen identified as P. lizae in Cheney et al. (2014). Phyllidiella rudmani and P. lizae resemble each other (Brunckhorst 1993) and hence it is possible that the species identified as P. lizae in Cheney et al. (2014) is in fact P. rudmani.

Specimens of seven nominal Phyllidia species were sequenced. Sequences of 25 individuals of Phyllidia elegans (including one from GenBank) formed a highly supported clade, just like the clades containing P. ocellata, P. picta, and P. varicosa. Phyllidia coelestis was also retrieved as a highly supported clade. An individual identified as P. picta by Cheney et al. (2014) was part of this group suggesting that it should probably be identified as P. coelestis. Brunckhorst (1993) already noticed the close similarity between the two species but still confused them (Yonow 1996), and hence identification errors are likely to occur. Individuals identified as Phyllidia babai Brunckhorst, 1993 and P. cf. babai were retrieved in two different clades. Specimens 336464 and 336614 differ in 75 base pairs, 336464 and 336575 by 68 base pairs and 336614 and 336575 by 32 base pairs. Differences based on COI suggest that they represent two, or possibly three, different species. The genus Reticulidia was retrieved as a sister group of Phyllidia.

Material of four nominal species in the genus Phyllidiopsis was sequenced, with additional data of one species from GenBank (P. cardinalis). Phyllidiopsis fissuratus clusters basal to Phyllidiella, without support. Phyllidiopsis shireenae Brunckhorst, 1990 and P. xishaensis (Lin, 1983) cluster as sister species, in highly supported clades. Phyllidiopsis krempfi also formed a clear group. Phyllidiopsis cardinalis does not cluster with any of the phyllidiid genera based on either the 16S or the COI analysis. This result suggests that P. cardinalis should be separated from the other Phyllidiopsis species, but further morphological analyses are needed to confirm this outcome. Brunckhorst (1993) noted that P. cardinalis is the type species of the genus Phyllidiopsis, and that it has a unique and complex coloration totally different from that of any other known phyllidiid species, as well as a different anatomy, especially in the foregut. Valdés (2003) states “Additionally, the genus Phyllidiopsis is not monophyletic when molecular characters are used, because Phyllidiopsis cardinalis is at the base of the Phyllidiidae clade, and not nested with the other members of Phyllidiopsis”. Surprisingly, in the analysis of Cheney et al. (2014), based on a concatenated dataset of 16S and COI mtDNA, P. cardinalis was retrieved in a highly supported clade with several species of Phyllidiella and Phyllidia.

Phyllidiella pustulosa is the only species in the COI cladogram (Figure 4) in which highly supported subclades can be discerned. To estimate the average evolutionary divergence within Phyllidiella pustulosa the base differences were compared per site for all grouped sequences of the species Phyllidia elegans (n = 24), P. varicosa (n = 15), Phyllidiella nigra (n = 7), P. pustulosa (n = 20), and Phyllidiopsis krempfi (n = 13) (Tables 4–5).

The genetic variation on the barcoding marker COI is much higher within Phyllidiella pustulosa (3.9%) than within the other four species, which showed genetic variations between 0.6 and 1.2% (Table 4). The interspecific genetic variation (involving three different genera) ranges between 10.5 and 18.9% (Table 5). The congeners Phyllidiella nigra and P. pustulosa differ by 10.5%, and the congeners Phyllidia elegans and P. varicosa differ by 12.1%. The observed levels of genetic variation within Phyllidiella pustulosa (Table 4) and between the five species (Table 5) call for additional studies on possible cryptic speciation in P. pustulosa.