Phylogenetic relationships within the Phyllidiidae (Opisthobranchia, Nudibranchia).

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 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 is considered an “Informal Group” among the (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 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).

The present study aims to clarify the phylogenetic relationships within the Rafinesque, 1814, belonging to the (Bouchet and Rocroi 2005). This family consists of more than 100 species divided over five genera: Eliot, 1903, Cuvier, 1797, Bergh, 1869, Bergh, 1875, and Brunckhorst, 1990 (Bouchet 2015). The genera JE Gray, 1853, and Yonow, 1986, have been synonymised with (Valdés and Gosliner 1999).

Most nudibranchs of the family 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 into a phylogenetic tree (e.g. Wollscheid-Lengeling et al. 2001) and even fewer deal with phylogenetic relationships among . 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 are the retractile lamellate rhinophores, the compact digestive gland mass, and the triaulic reproductive system (Brunckhorst 1993). Another important character diagnosing the is the possession of numerous subdermal calcareous spicules of different microstructures (Chang et al. 2013). The 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 , 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.

Materials and methods

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 PageBreak2009 PageBreakoff 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 were sampled in 2007.

Figure 3.

Halmahera and Ternate sites where were sampled in 2009.

Table 1.

Information on analysed species: RMNH.MOL catalogue number or field code number in case voucher specimen became lost; Genbank number if available; collection site, station number (RAJ, TER), coordinates.

= Raja Ampat

= Ternate, Halmahera

RMNH.MOL or Field nr.	Genbank accession number	Species	Locality	Station	Coordinates
336464	KX235918	Phyllidia babai	Tanjung Ebamadu	TER08	N0°45'23.4", E127°24'26.5"
336575	KX235920	Phyllidia cf. babai	South Gam, shoal near mangroves	RAJ37	S0°31'08.2", E130°38'28.0"
336614	KX235919	Phyllidia cf. babai	Tanjung Ratemu (South of river)	TER27	N0°54'44.5", E127°29'09.9"
336573	KX235921	Phyllidia coelestis	Eastern entrance of passage	RAJ44	S0°25'44.3", E130°33'56.8"
336574	KX235922	Phyllidia coelestis	Wallace Lake	RAJ13	S0°26'31.1", E130°41'08.0"
58		Phyllidia elegans	Pulau Maka	TER13	N0°54'42.7", E127°18'32.9"
137		Phyllidia elegans	Pulau Pilongga, North	TER34	N0°42'49.8", E127°28'45.4"
156		Phyllidia elegans	Teluk Dodinga; Karang Ngeli West	TER40	N0°46'25.3", E127°32'22.0"
336475	KX073972	Phyllidia elegans	Tanjung Tabam	TER12	N0°50'05.1", E127°23'10.0"
336478	KX073973	Phyllidia elegans	Pulau Maka	TER13	N0°54'42.7", E127°18'32.9"
336488	KX073974	Phyllidia elegans	Tanjung Pasir Putih	TER16	N0°51'50.4", E127°20'36.7"
336514	KX073975	Phyllidia elegans	Dufadufa / Benteng Toloko	TER24	N0°48'49.1", E127°23'21.6"
336515	KX073976	Phyllidia elegans	Idem	TER24	N0°48'49.1", E127°23'21.6"
336554	KX073985	Phyllidia elegans	Passage	RAJ43	S0°25'45.2", E130°33'37.3"
336555	KX073990	Phyllidia elegans	Akber Reef	RAJ14	S0°34'15.2", E130°39'33.7"
336556	KX073988	Phyllidia elegans	Passage	RAJ43	S0°25'45.2", E130°33'37.3"
336557	KX073987	Phyllidia elegans	Idem	RAJ43	S0°25'45.2", E130°33'37.3"
336558	KX073984	Phyllidia elegans	Southwest Pulau Kri	RAJ40	S0°33'58.1", E130°39'46.2"
336559	KX073991	Phyllidia elegans	South Gam, shoal near mangroves	RAJ37	S0°31'08.2", E130°38'28.0"
336560	KX073983	Phyllidia elegans	Southwest Pulau Kri	RAJ40	S0°33'58.1", E130°39'46.2"
336561	KX073986	Phyllidia elegans	Passage	RAJ43	S0°25'45.2", E130°33'37.3"
336562	KX073989	Phyllidia elegans	Akber Reef	RAJ14	S0°34'15.2", E130°39'33.7"
336628	KX073977	Phyllidia elegans	Pulau Gura Ici, East	TER32	S0°01'17.3", E127°14'17.2"
336629	KX073978	Phyllidia elegans	Idem	TER32	S0°01'17.3", E127°14'17.2"
336631	KX073979	Phyllidia elegans	Pulau Pilongga, North	TER34	N0°42'49.8", E127°28'45.4"
336632	KX073980	Phyllidia elegans	Idem	TER34	N0°42'49.8", E127°28'45.4"
336633	KX073981	Phyllidia elegans	Idem	TER34	N0°42'49.8", E127°28'45.4"
336649	KX073982	Phyllidia elegans	Teluk Dodinga; Karang Ngeli West	TER40	N0°46'25.3", E127°32'22.0"
336484	KX235923	Phyllidia exquisita	Tanjung Ngafauda	TER14	N0°54'38.3", E127°29'20.7"
336494	KX235924	Phyllidia ocellata	Southwest of Tobala	TER19	N0°44'56.6", E127°23'13.5"
336563	KX235926	Phyllidia ocellata	Southeast Gam, Friwen Wonda	RAJ11	S0°28'29.9", E130°41'54.8"
336564	KX235925	Phyllidia ocellata	Idem	RAJ11	S0°28'29.9", E130°41'54.8"
336565	KX235927	Phyllidia picta	South Gam, Shoal near mangroves	RAJ37	S0°31'08.2", E130°38'28.0"
336566	KX235929	Phyllidia picta	Passage	RAJ43	S0°25'45.2", E130°33'37.3"
336567	KX235928	Phyllidia picta	North Batanta, West Telok Gegenlol	RAJ29	S0°49'42.5", E130°42'42.0"
336619	KX235930	Phyllidia sp.	Pulau Popaco, East	TER28	S0°01'51.9", E127°14'01.8"
74		Phyllidia varicosa	Tanjung Pasir Putih	TER16	N0°51'50.4", E127°20'36.7"
336489	KX235931	Phyllidia varicosa	Idem	TER16	N0°51'50.4", E127°20'36.7"
336568	KX235942	Phyllidia varicosa	Northeast Pulau Mansuar	RAJ38	S0°34'05.0", E130°38'31.5"
336569	KX235941	Phyllidia varicosa	Idem	RAJ38	S0°34'05.0", E130°38'31.5"
336570	KX235943	Phyllidia varicosa	North Batanta, West Telok Gegenlol	RAJ29	S0°49'42.5", E130°42'42.0"
336571	KX235938	Phyllidia varicosa	South Gam, Eastern entrance Besir Bay, Cape Besir	RAJ25	S0°30'51.5", E130°34'11.5"
336572	KX235940	Phyllidia varicosa	Idem	RAJ25	S0°30'51.5", E130°34'11.5"
336604	KX235932	Phyllidia varicosa	East side Ternate Harbour (outside)	TER25	N0°46'55.3", E127°23'19.9"
336609	KX235933	Phyllidia varicosa	Pasir Lamo (West side)	TER26	N0°53'20.5", E127°27'34.2"
336612	KX235934	Phyllidia varicosa	Idem	TER26	N0°53'20.5", E127°27'34.2"
336617	KX235935	Phyllidia varicosa	Tanjung Ratemu (South of river)	TER27	N0°54'44.5", E127°29'09.9"
336621	KX235936	Phyllidia varicosa	Pulau Popaco E	TER28	S0°01'51.9", E127°14'01.8"
336637	KX235937	Phyllidia varicosa	Teluk Dodinga East; North of Pulau Jere	TER36	N0°50'47.8", E127°37'48.7"
336647	KX235939	Phyllidia varicosa	Teluk Dodinga, Karang Galiasa Kecil West	TER39	N0°51'09.1", E127°35'19.5"
336590	KX235944	Phyllidiopsis fissuratus	Yenweres Bay	RAJ46	S0°29'13.0", E130°40'23.6"
336589	KX235945	Phyllidiella rudmani	Southeast Gam, Friwen Wonda	RAJ11	S0°28'29.9", E130°41'54.8"
336434	KX235946	Phyllidiella nigra	Off Danau Laguna	TER02	N0°45'29.7", E127°20'59.2"
336471	KX235947	Phyllidiella nigra	Maitara Northwest	TER10	N0°44'32.0", E127°21'50.9"
336472	KX235948	Phyllidiella nigra	Idem	TER10	N0°44'32.0", E127°21'50.9"
336501	KX235949	Phyllidiella nigra	Sulamadaha I	TER22	N0°52'03.6", E127°19'33.1"
336505	KX235950	Phyllidiella nigra	Sulamadaha II	TER23	N0°52'02.0", E127°19'45.8"
336576	KX235952	Phyllidiella nigra	South Gam, Eastern entrance Besir Bay, Pulau Bun	RAJ26	S0°30'59.3", E130°33'48.7"
336577	KX235951	Phyllidiella nigra	South Gam, Southeast Besir Bay	RAJ32	S0°30'45.2", E130°35'00.1"
75F		Phyllidiella pustulosa	North Batanta, West Telok Gegenlol	RAJ29	S0°49'42.5", E130°42'42.0"
336436	KX235953	Phyllidiella pustulosa	Off Danau Laguna	TER02	N0°45'29.7", E127°20'59.2"
336460	KX235954	Phyllidiella pustulosa	Desa Tahua	TER07	N0°45'09.1", E127°23'31.3"
336461	KX235955	Phyllidiella pustulosa	Idem	TER07	N0°45'09.1", E127°23'31.3"
336470	KX235956	Phyllidiella pustulosa	Northwest side of Maitara	TER10	N0°44'32.0", E127°21'50.9"
336474	KX235957	Phyllidiella pustulosa	Tanjung Tabam	TER12	N0°50'05.1", E127°23'10.0"
336495	KX235958	Phyllidiella pustulosa	Tanjung Ratemu (South of river)	TER21	N0°54'24.7", E127°29'17.7"
336508	KX235959	Phyllidiella pustulosa	Dufadufa / Benteng Toloko	TER24	N0°48'49.1", E127°23'21.6"
336510	KX235960	Phyllidiella pustulosa	Idem	TER24	N0°48'49.1", E127°23'21.6"
336578	KX235965	Phyllidiella pustulosa	South Gam, Southeast Besir Bay	RAJ32	S0°30'45.2", E130°35'00.1"
336579	KX235971	Phyllidiella pustulosa	South Gam, Besir Bay	RAJ35	S0°48'58.3", E130°59'16.6"
336580	KX235967	Phyllidiella pustulosa	Southwest Pulau Kri	RAJ40	S0°33'58.1", E130°39'46.2"
336581	KX235963	Phyllidiella pustulosa	South Gam, Besir Bay	RAJ35	S0°48'58.3", E130°59'16.6"
336582	KX235968	Phyllidiella pustulosa	Southwest Pulau Kri	RAJ40	S0°33'58.1", E130°39'46.2"
336583	KX235964	Phyllidiella pustulosa	South Gam, East entrance Besir Bay, Cape Besir	RAJ25	S0°30'51.5", E130°34'11.5"
336584	KX235961	Phyllidiella pustulosa	West Pulau Yeben Kecil	RAJ48	S0°29'20.6", E130°30'04.9"
336585	KX235969	Phyllidiella pustulosa	Southeast Gam, Desa Besir	RAJ41	S0°27'48.1", E130°41'14.6"
336586	KX235966	Phyllidiella pustulosa	Idem	RAJ41	S0°27'48.1", E130°41'14.6"
336587	KX235962	Phyllidiella pustulosa	South Gam, Eastern entrance Besir Bay, Cape Besir	RAJ25	S0°30'51.5", E130°34'11.5"
336588	KX235970	Phyllidiella pustulosa	West Pulau Yeben Kecil	RAJ48	S0°29'20.6", E130°30'04.9"
336453	KX235972	Phyllidiopsis krempfi	Kampung Cina / Tapak 2	TER06	N0°47'15.0", E127°23'25.0"
336462	KX235973	Phyllidiopsis krempfi	Tanjung Ebamadu	TER08	N0°45'23.4", E127°24'26.5"
336466	KX235974	Phyllidiopsis krempfi	Idem	TER08	N0°45'23.4", E127°24'26.5"
336469	KX235975	Phyllidiopsis krempfi	West Maitara	TER09	N0°43'47.6", E127°21'44.7"
336512	KX235976	Phyllidiopsis krempfi	Dufadufa / Benteng Toloko	TER24	N0°48'49.1", E127°23'21.6"
336594	KX235979	Phyllidiopsis krempfi	Southwest Pulau Kri, Kuburan	RAJ15	S0°33'42.8", E130°39'40.4"
336595	KX235984	Phyllidiopsis krempfi	Southwest Pulau Kri	RAJ40	S0°33'58.1", E130°39'46.2"
336596	KX235983	Phyllidiopsis krempfi	Northwest Pulau Mansuar, Lalosi reef	RAJ49	S0°32'53.5", E130°29'51.1"
336597	KX235978	Phyllidiopsis krempfi	Southwest Pulau Kri, Kuburan	RAJ15	S0°33'42.8", E130°39'40.4"
336598	KX235980	Phyllidiopsis krempfi	North Batanta, North Pulau Yarifi	RAJ28	S0°46'46.7", E130°42'42.7"
336599	KX235982	Phyllidiopsis krempfi	East Kri, Sorido Wall	RAJ12	S0°33'13.2", E130°41'16.9"
336600	KX235981	Phyllidiopsis krempfi	Northeast Mansuar	RAJ38	S0°34'05.0", E130°38'31.5"
336650	KX235977	Phyllidiopsis krempfi	Teluk Dodinga; West Karang Ngeli	TER40	N0°46'25.3", E127°32'22.0"
336451	KX235985	Phyllidiopsis shireenae	Kampung Cina / Tapak 2	TER06	N0°47'15.0", E127°23'25.0"
336652	KX235986	Phyllidiopsis shireenae	Teluk Dodinga; East Karang Luelue	TER41	N0°46'32.8", E127°33'43.4"
336591	KX235987	Phyllidiopsis xishaensis	Southeast Gam, Pulau Kerupiar, Mike’s Point	RAJ05	S0°30'57.1", E130°40'22.1"
336592	KX235988	Phyllidiopsis xishaensis	East Pulau Kri, Cape Kri	RAJ07	S0°33'22.2", E130°41'28.7"
336593	KX235989	Phyllidiopsis xishaensis	Eastern entrance of passage	RAJ44	S0°25'44.3", E130°33'56.8"
336640	KX235990	Reticulidia fungia	East Teluk Dodinga; North of Pulau Jere	TER36	N0°50'47.8", E127°37'48.7"
336455	KX235991	Reticulidia halgerda	Kampung Cina / Tapak 2	TER06	N0°47'15.0", E127°23'25.0"

Figure 5.

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

Figure 15.

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

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 PageBreakfor 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 based on COI gene sequence data (Figure 4).

Figure 4.

Phylogeny reconstruction of the 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.

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

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

External morphology and colouration of specimens used for COI phylogeny reconstruction: (a–f), (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).

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

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

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

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

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

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

Polymerase Chain Reaction

Sequence alignment and phylogenetic analyses

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 (Cheeseman, 1881) and Bergh, 1880 were collected from GenBank and used as outgroups. Additional COI sequences of Bergh, 1905, Bergh, 1869, Cuvier, 1804, Pruvot-Fol, 1957, Lamarck, 1801, Brunckhorst, 1993, (Cuvier, 1804), Bergh, 1875 were obtained from GenBank (Table 2).

Table 2.

Mitochondrial COI sequences of (and outgroups) obtained from GenBank.

Species	Accession number	Reference	Collection locality
Dendrodoris citrina	GQ292043	Shields et al. (2009 unpubl.)	Ross Sea, Antarctica?
Doriopsilla areolata	AJ223262	Thollesson (2000)	Cadiz, Andalusia, Spain
Phyllidia coelestis	KJ001305	Cheney et al. (2014)	Lizard I., Queensland Australia
Phyllidia elegans	AJ223276	Thollesson (2000)	Tab I., Papua New Guinea
Phyllidia ocellata	KJ001307	Cheney et al. (2014)	Mooloolaba, Queensland, Australia
Phyllidia picta	KJ001304	Cheney et al. (2014)	Lizard I., Queensland Australia
Phyllidia varicosa	KJ001306	Cheney et al. (2014)	Lizard I., Queensland Australia
Phyllidiella lizae	KJ001309	Cheney et al. (2014)	Lizard I., Queensland Australia
Phyllidiella pustulosa	KJ001310	Cheney et al. (2014)	Lizard I., Queensland Australia
Phyllidiopsis cardinalis	KJ001308	Cheney et al. (2014)	Mooloolaba, Queensland, Australia

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 calculated with jModeltest (Posada 2008), resulting in TVM+I+G as the most suitable model. Phylogenetic reconstructions were PageBreakcarried 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 ; 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).

Akaike Information Criterion

Markov Chain Monte Carlo techniques

Initial phylogenetic analyses showed high intraspecific variation on the COI region between specimens identified as . Tests to estimate the average evolutionary divergence over sequence pairs between and within groups were carried out in MEGA 6.06. , , (van Hasselt, 1824), , and 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 sequence from GenBank (KJ001310) was excluded from this analysis: based on its position in the phylogeny reconstruction the identification of this specimen as is doubtful. The web version of ABGD, Puillandre et al. 2012) was used to estimate the genetic distance corresponding to the difference between a speciation process versus intra-specific variation in . 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.

(Automatic Barcode Gap Discovery

All available mitochondrial 16S sequences of 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 ( (Bergh, 1869) and Brunckhorst, 1993) for which no COI data were available. (JG Cooper, 1863) was used as outgroup (Table 3). The sequences were aligned using the Guidance server PageBreak(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 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).

Table 3.

16S sequences of obtained from GenBank.

Species	Accession number	Reference	Collection locality
Doropsilla albopunctata	AF430354	Valdés (2003)	Baja California, Mexico
Phyllidia coelestis	AF430361	Valdés (2003)	Lifou I., New Caledonia
Phyllidia coelestis	KJ018917	Cheney et al. (2014)	Lizard I., Queensland Australia
Phyllidia elegans	AF430362	Valdés (2003)	Lifou I., New Caledonia
Phyllidia elegans	AJ225201	Thollesson (2000)	Tab I., Papua New Guinea
Phyllidia ocellata	AF430363	Valdés (2003)	Lifou I., New Caledonia
Phyllidia picta	KJ018916	Cheney et al. (2014)	Lizard I., Queensland Australia
Phyllidia rueppelii	AF430358	Valdés (2003)	Hurghada, Egypt
Phyllidiella lizae	AF430365	Valdés (2003)	Lifou I., New Caledonia
Phyllidiella lizae	KJ018918	Cheney et al. (2014)	Lizard I., Queensland Australia
Phyllidiella pustulosa	AF249232	Wollscheid-Lengeling et al. (2001)	Great Barrier Reef, Australia
Phyllidiella pustulosa	AF430366	Valdés (2003)	Lifou I., New Caledonia
Phyllidia varicosa	AF430364	Valdés (2003)	Lifou I., New Caledonia
Phyllidiopsis cardinalis	AF430367	Valdés (2003)	Lifou I., New Caledonia
Phyllidiopsis sphingis	AF430368	Valdés (2003)	Lifou I., New Caledonia
Phyllidiopsis xishaensis*	AF430369	Valdés (2003)	Lifou I., New Caledonia
Reticulidia fungia	AF430370	Valdés (2003)	Lifou I., New Caledonia
Reticulidia halgerda	AF430371	Valdés (2003)	Lifou I., New Caledonia

* Re-identification according to Yonow (pers. comm.)

Akaike Information Criterion

Results and discussion

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 PageBreakthe 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 , , , and are retrieved in distinct clades, with as a sister clade to . Brunckhorst, 1993 formed a separate lineage basal to species (albeit without support). does not cluster with its congeners, but instead forms a separate lineage in the .

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 is separated by a long branch. Within the overall clade four main groupings can be distinguished: , , and , and a mixed clade of and . Based on this analysis only the genus is monophyletic. 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).

Figure 17.

Phylogeny reconstruction of the 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 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 . Brunckhorst (1993) distinguished from based on the position of the anus and other anatomical features. (with its synonyms (Pruvot-Fol, 1957), Brunckhorst, 1993, (Brunckhorst, 1993)) was included in our analyses which, according to Brunckhorst, should belong to the genus . 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 in the genus .

Figure 16.

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

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 and . 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 a sister group to all the other genera (Figure 6b). Valdés (2002; Figure 16d) distinguished two larger groupings within the ; and as one group and , , and as the other group. and in turn formed a sister group of (Figure 16d). The cladogram by Brunckhorst (1993) and our cladogram based on COI (Figure 4) both show that is a sister clade of and . In contrast, is not a sister group of but to all the other genera grouped together in the cladogram of Valdés (2003).

Figure 6.

External morphology and colouration of specimens used for COI phylogeny reconstruction: . 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).

Unfortunately no specimens were available to complete our analysis at genus level. Up to this point the phylogenetic position of the genus 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 . formed a highly supported clade. In the clade containing much variation is visible indicating larger genetic differences among individuals. The ABGD analysis shows that four Molecular Operational Taxonomic Units (MOTUs) are present in , suggesting the presence of cryptic species or, alternatively, high intraspecific variation. The of Cheney et al. (2014) falls in between the group consisting of and on one side and Brunckhorst, 1993 on the other and probably represents another species. Our specimen of clustered with the specimen identified as in Cheney et al. (2014). and resemble each other (Brunckhorst 1993) and hence it is possible that the species identified as in Cheney et al. (2014) is in fact .

Specimens of seven nominal species were sequenced. Sequences of 25 individuals of (including one from GenBank) formed a highly supported clade, just like the clades containing , , and was also retrieved as a highly supported clade. An individual identified as by Cheney et al. (2014) was part of this group suggesting that it should probably be identified as . 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 Brunckhorst, 1993 and 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 was retrieved as a sister group of .

Material of four nominal species in the genus was sequenced, with additional data of one species from GenBank (). clusters basal to , without support. Brunckhorst, 1990 and (Lin, 1983) cluster as sister species, in highly supported clades. also formed a clear group. does not cluster with any of the phyllidiid genera based on either the 16S or the COI analysis. This result suggests that should be separated from the other species, but further morphological analyses are needed to confirm this outcome. Brunckhorst (1993) noted that is the type species of the genus , 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 is not monophyletic when molecular characters are used, because is at the base of the clade, and not nested with the other members of ”. Surprisingly, in the analysis of Cheney et al. (2014), based on a concatenated dataset of 16S and COI mtDNA, was retrieved in a highly supported clade with several species of and .

Variation within

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 the base differences were compared per site for all grouped sequences of the species (n = 24), (n = 15), (n = 7), (n = 20), and (n = 13) (Tables 4–5).

Table 4.

Estimates of average evolutionary divergence (p-distance) over sequence pairs within groups, in percentages.

Species	Distance (%)
Phyllidia elegans	0.7
Phyllidia varicosa	0.7
Phyllidiella nigra	0.6
Phyllidiella pustulosa	3.9
Phyllidiopsis krempfi	1.2

Table 5.

Estimates of average evolutionary divergence (p-distance) over sequence pairs between groups, in percentages.

	Distance (%)
Species	Phyllidia elegans	Phyllidia varicosa	Phyllidiella nigra	Phyllidiella pustulosa	Phyllidiopsis krempfi
Phyllidia elegans
Phyllidia varicosa	12.1
Phyllidiella nigra	15.8	15.5
Phyllidiella pustulosa	18.3	18.9	10.5
Phyllidiopsis krempfi	15.8	16.4	14.6	17.2

The genetic variation on the barcoding marker COI is much higher within (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 PageBreak and differ by 10.5%, and the congeners and differ by 12.1%. The observed levels of genetic variation within (Table 4) and between the five species (Table 5) call for additional studies on possible cryptic speciation in .

Conclusions

The barcoding marker COI works well to separate the different species in the , and confirms that the species boundaries in highly variable species, such as , , and , 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 (, ) the studied genera (, , , and ) were retrieved as separate genera within the family. Additional representatives of are needed to indicate the position of this genus within the . The observed groupings within 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 that could be separated slightly according PageBreakto 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.