Phylogenetic relationships within the Phyllidiidae (Opisthobranchia, Nudibranchia)

Abstract The Phyllidiidae (Gastropoda, Heterobranchia, Nudibranchia) is a family of colourful nudibranchs found on Indo-Pacific coral reefs. Despite the abundant and widespread occurrence of many species, their phylogenetic relationships are not well known. The present study is the first contribution to fill the gap in our knowledge on their phylogeny by combining morphological and molecular data. For that purpose 99 specimens belonging to 16 species were collected at two localities in Indonesia. They were photographed and used to make a phylogeny reconstruction based on newly obtained cytochrome oxidase subunit (COI) sequences as well as sequence data from GenBank. All mitochondrial 16S sequence data available from GenBank were used in a separate phylogeny reconstruction to obtain information for species we did not collect. COI data allowed the distinction of the genera and species, whereas the 16S data gave a mixed result with respect to the genera Phyllidia and Phyllidiella. Specimens which could be ascribed to species level based on their external morphology and colour patterns showed low variation in COI sequences, but there were two exceptions: three specimens identified as Phyllidia cf. babai represent two to three different species, while Phyllidiella pustulosa showed highly supported subclades. The barcoding marker COI also confirms that the species boundaries in morphologically highly variable species such as Phyllidia elegans, Phyllidia varicosa, and Phyllidiopsis krempfi, are correct as presently understood. In the COI as well as the 16S cladogram Phyllidiopsis cardinalis was located separately from all other Phyllidiidae, whereas Phyllidiopsis fissuratus was positioned alone from the Phyllidiella species by COI data only. Future studies on phyllidiid systematics should continue to combine morphological information with DNA sequences to obtain a clearer insight in their phylogeny.


Introduction
Nudibranch gastropod molluscs have traditionally been classified with the Infraclass Opisthobranchia Milne Edwards, 1848, which consists of more than 6000 species (Yonow 2008). Although this taxon is not monophyletic and therefore is considered obsolete (Schrödl et al. 2011), taxonomic works still refer to "opisthobranchs" for practical reasons (e.g. Uribe et al. 2013) and Opisthobranchia is considered an "Informal Group" among the Heterobranchia (Wägele et al. 2014). These animals form, ecologically and morphologically, one of the most diverse groups of marine gastropods (Wägele et al. 2014). To avoid use of their misnomer, this well-known group of marine animals can also be referred to as sea slugs (Yonow 2015). Among these, the Nudibranchia Cuvier, 1817 form the largest order with an estimated number of more than 2000 species (Gosliner et al. 2008), although also estimates of nearly 3000 species are known (Vonnemann et al. 2005).
Much work has already been done to elucidate the phylogeny of the opisthobranchs by molecular analyses (e.g., Wollscheid and Wägele 1999, Grande et al. 2004a, 2004b, Vonnemann et al. 2005, Turner and Wilson 2008, Maeda et al. 2010, Pola and Gosliner 2010, but most of the phylogenetic relationships still remain unclear at family, genus, and species level, especially with regards to the nudibranchs. All nudibranch species and many other sea slugs are predators, which usually can be observed together with their prey (Behrens 2005, Pola and Gosliner 2010, van Alphen et al. 2011. Only rarely they are found together with potential predators such as sea anemones, mushroom corals, and pycnogonids (Piel 1991, Behrens 2005, van der Meij and Reijnen 2012, Mehrotra et al. 2015).
Most nudibranchs of the family Phyllidiidae are commonly encountered on coral reefs, where they can easily be noticed because of their aposomatic colouration, which serves to deter possible predators from eating them (Ritson-Williams and Paul 2007). Nevertheless, only eight phyllidiid COI sequences can be found in GenBank, as well as two 18S sequences and 17 16S sequences. There are only a few published studies that incorporate even a single member of Phyllidiidae into a phylogenetic tree (e.g. Wollscheid-Lengeling et al. 2001) and even fewer deal with phylogenetic relationships among Phyllidiidae. Among the latter, most are using anatomical characters (Brunckhorst 1993, Valdés and Gosliner 1999, Valdés 2001, 2002 and only two are known to include a molecular and phylogenetic analysis (Valdés 2003, Cheney et al. 2014. Phyllidiid slugs are characterized by their oval elongate and tough bodies, which generally possess hard notal tubercles on the dorsal side. Although their colouration is a main character used for their identification, many species cannot be identified based on colouration alone owing to their high intra-specific colour variation. Structure and pattern of the notal tubercles are important characters for identification. Other distinctive features of the Phyllidiidae are the retractile lamellate rhinophores, the compact digestive gland mass, and the triaulic reproductive system (Brunckhorst 1993). Another important character diagnosing the Phyllidiidae is the possession of numerous subdermal calcareous spicules of different microstructures (Chang et al. 2013). The Phyllidiidae have no jaws or radula and lack the dorsal, circumanal circlet of gills that is typical of other dorids (Brunckhorst 1993).
To study the phylogenetic relationships within the Phyllidiidae, a molecular analysis was performed based on DNA sequence data of the mitochondrial cytochrome oxidase I (COI) gene, combined with external morphological assessments of material collected in two areas in eastern Indonesia, the Raja Ampat islands (West Papua) and Ternate, off western Halmahera (Moluccas). Both locations are situated in the centre of maximum marine biodiversity, also known as the Coral Triangle (Hoeksema 2007). In earlier studies, high numbers of phyllidiid species were recorded from this area: 13 from the Bismarck Sea, Papua New Guinea (Domínguez et al. 2007), eleven from Ambon (Moluccas, Indonesia) (Yonow 2011), and eleven from the South China Sea (Sachidhanandam et al. 2000). Therefore, both of our areas were expected to show a high number of phyllidiid species that could be used for the present study.

Sampling
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 (

Morphological study
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).

DNA extraction
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.

PCR amplification, purification, and sequencing
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 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  (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 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

Position of genera
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.
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).  Brunckhorst (1993) based on morphological data showing topology of six genera of Phyllidiidae c 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.
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
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.
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 Phyl- lidiella 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.

Conclusions
The barcoding marker COI works well to separate the different species in the Phyllidiidae, and confirms that the species boundaries in highly variable species, such as Phyllidia elegans, P. varicosa, and Phyllidiopsis krempfi, are correct as presently understood. However, a multi-locus approach, preferably including nuclear markers, is needed to improve the resolution for the higher taxonomic levels. With the exception of a few species that are difficult to place (Phyllidiopsis fissuratus, Phyllidiopsis cardinalis) the studied genera (Phyllidia, Phyllidiella, Phyllidiopsis, and Reticulidia) were retrieved as separate genera within the family. Additional representatives of Ceratophyllidia are needed to indicate the position of this genus within the Phyllidiidae. The observed groupings within Phyllidiella pustulosa suggest that multiple (cryptic) species could be present in this species, for which further analyses are needed including morphological data and multiple markers. Chang and Willan (2015) indicated that at least nine clades could be recognized in Phyllidiella pustulosa that could be separated slightly according to morphological characters. We recommend that future studies combine DNA sequences with morphological characters, which can easily be done by adding pictures of the specimens to avoid increasing confusion in the identification of specimens.