Phyllidiidae (Nudibranchia, Heterobranchia, Gastropoda): an integrative taxonomic approach including chemical analyses

Abstract

Members of the widely distributed and common nudibranch family Phyllidiidae are often easily spotted in the marine environment because of their conspicuous colours and obvious presence on the reef. They are interesting with regard to their defensive chemical compounds that may lead to new drug discoveries. Despite their abundance, the family is also well known for its taxonomic problems and the difficulties in species identification due to very similarly coloured species and lack of morphological characters. In this study, phyllidiid species were analysed using an integrative approach. Molecular analysis of the mitochondrial genes 16S and CO1 was utilised, running phylogenetic analyses, species delimitation tests, and haplotype network analyses. Additionally, for the first time, external morphological characters were analysed, museum material was re-analysed, and chemical profiles were applied for characterising species. The analyses are based on sequences of 598 specimens collected in Indonesia by the team, with the addition of published sequences available on GenBank. This study comprises 11 species of Phyllidia, seven species of Phyllidiopsis, and at least 14 species of Phyllidiella. Moreover, 11 species belonging to these three genera are probably new to science, Phyllidiopsis pipeki is synonymised with P. krempfi, and Phyllidiella albonigra is resurrected. Some of the external colouration previously used for species identification is shown to not be valid, but alternative characters are provided for most species. Chemical analyses led to species characterisation in a few examples, indicating that these species use particular sponge species as food; however, many species show a broad array of compounds and are therefore characterised more by their composition or profile than by distinct or unique compounds.

Introduction

Nudibranchia is one of the largest clades within the marine Heterobranchia. It comprises two lineages, the Cladobranchia and the Doridina. Within the Doridina, Phyllidiidae is one of the oldest recognised families, described by Rafinesque (1814). Phyllidiidae and Dendrodorididae were united under the Porostomata Bergh, 1876 (now Phyllidioidea in WoRMS Editorial Board, 2021), based on the pore-like mouth opening and the complete loss of radula and jaws in both families. More than one century later, a third family with similar characteristics was described, Mandeliidae Valdés & Gosliner, 1999, and these three families are diagnosed by the absence of the radula and jaws (Valdés, 2003; Valdés & Gosliner, 1999).

Specimens of Phyllidiidae are very commonly encountered on Indo-Pacific reefs. Because of their abundance, often large size, and usually striking colours of bright yellow or orange in combination with black and white, they have been investigated due to their chemical compounds used for defence (Bogdanov et al., 2020; Böhringer et al., 2017; Fisch et al., 2017; Hagadone et al., 1979; Manzo et al., 2004; White et al., 2015). The cloudy white mucous secretions produced in defence have a bad smell, observed as early as 1963, and comprise volatile compounds lethal to crustaceans and fish (Johannes, 1963). Interestingly, Okino et al. (1996) suggested that the sesquiterpenes in this mucus may serve as antifouling in order to maintain the exposed surfaces free of epibionts. In addition to the compounds that are taken up from the sponge food items, phyllidiids also use a dense layer of calcareous spicules arranged in the mantle and foot that provide structural support (Kasamesiri et al., 2011; Wallraff, 2018) and make them difficult to eat (van der Meij & Reijnen, 2012). This combination of defence systems probably allows them to be conspicuously exposed on reef substrates during daytime.

Despite their commonness and conspicuous appearance, many members within the Phyllidiidae are known to be difficult to identify, and cryptic speciation and species complexes have been described recently (Bogdanov et al., 2020; Stoffels et al., 2016). Species of Phyllidiidae were first documented from the tropical Indo-Pacific, which retains the highest number of species, with Phyllidia varicosa Lamarck, 1801, Phyllidia ocellata Cuvier, 1804, and Phyllidiella pustulosa (Cuvier, 1804) among the first species described; the type localities of these are Réunion Island and “Mer des Indes”. Few species are documented in the Atlantic Ocean, Caribbean Sea, and Mediterranean Sea (Brunckhorst, 1993; Valdés & Ortea, 1996; Wägele, 1985; Phyllidiidae in GBIF.org). Most are found in shallow waters, but Valdés (2001a, 2001b) focused on deep sea species, showing that the majority of phyllidiids living in depths from 100 m to below 500 m could be assigned to the genus Phyllidiopsis. A major distribution area lies in the Coral Triangle, which includes Indonesia. Sulawesi is one of the largest islands of Indonesia and is unique in its position between the Wallace and the Webber lines, as well as close to the Lydekker line (Cockey, 2013). Sulawesi is therefore characterised by a species-overlap zone with members from the Asian continent as well as from the Australian continent. The high number of phyllidiid species and specimens in these areas were recently documented in several studies (Eisenbarth et al., 2018; Kaligis et al., 2018; Ompi et al., 2019; Papu et al., 2020; Undap et al., 2019).

Only few studies have begun to analyse phylogenetic relationships between the five recognised genera Ceratophyllidia, Phyllidia, Phyllidiella, Phyllidiopsis, and Reticulidia, containing 81 species (WoRMS Editorial Board, 2021). Brunckhorst (1993) provided the first extensive morphological analysis and review of the Phyllidiidae, distinguishing and defining six genera, the five listed above and Fryeria. Fryeria was later synonymised with Phyllidia by Valdés and Gosliner (1999) who analysed relationships using morphological data based on 11 species covering all five genera; however, the relationships of Phyllidiella and Phyllidia could not be resolved.

In a preliminary molecular study, Valdés (2003) published a phylogenetic tree including 12 specimens from four genera (without Ceratophyllidia). In that work, Phyllidiella grouped with species of the genus Phyllidiopsis, and Phyllidiopsis therefore appeared paraphyletic. Subsequent molecular studies using only the mitochondrial gene CO1 were performed by Stoffels et al. (2016) including 18 species belonging to four genera, also omitting Ceratophyllidia. That analysis recovered the monophyly of Phyllidia and Phyllidiella; however, Phyllidiopsis still appeared paraphyletic. Interestingly, that study presented evidence of cryptic variation in Phyllidiella ‘pustulosa’, which was most recently confirmed by the more extensive molecular and chemical composition study of Bogdanov et al. (2020). At least seven subclades were identified within the Phyllidiella pustulosa complex that were not assigned to species. Moreover, clade-dependent metabolomes became evident, and several unique compounds were identified in one specific clade, confirming the distinctiveness of the cryptic varieties within Phyllidiella pustulosa.

In contrast to the scarce phylogenetic analyses, the chemistry of phyllidiids has evoked interest for many decades. Following an observation of a Phyllidia varicosa specimen secreting “a light-grey mucus” that caused deadly poisoning of a lobster in an aquarium setting (Johannes, 1963), P. varicosa was among the first phyllidiid nudibranchs attracting the attention of natural product chemists who isolated a sesquiterpene isonitrile (Burreson et al., 1975). Several studies focusing on phyllidiid chemistry followed (see review in Fisch et al., 2017) and the published results hint at the variety of compounds isolated from one and the same species, e.g., Phyllidiella pustulosa, implying either a very high degree of metabolomic variation or that species identifications might not be correct in all cases.

In the present study we run a phylogenetic analysis of the family Phyllidiidae by using new and available GenBank sequences from the partial mitochondrial genes CO1 and 16S and perform species delimitation tests to identify putative cryptic variations/species within the genera. We also assign morphological characters to the identified species/clades. Finally, we performed metabolomic analyses of some specimens used in the molecular analyses as an additional chemotaxonomic evidence of species delimitation and identification. With 716 sequences, our study is the most comprehensive phylogenetic analysis of the Phyllidiidae to date and is the first time metabolome analysis has been used, providing unprecedented insights into clade-specific chemical compositions that largely support many of the recognised and new clades in our phyllidiid study.

Materials and methods

Sample collection

Phyllidiidae specimens were collected during several expeditions from 2015 to 2018 in the area of North Sulawesi by SCUBA diving and snorkelling in depths between 0 and 30 m (Fig. 1). The localities include Bunaken National Park, Bangka Archipelago, and Sangihe Island. Species lists of these expeditions have been published in Kaligis et al. (2018), Eisenbarth et al. (2018), Undap et al. (2019), Ompi et al. (2019), and Papu et al. (2020). In addition, phyllidiids were collected in South and Southeast Sulawesi in 2018 (Fig. 1), totalling 610 specimens. Each collected specimen was (usually) documented underwater and additionally in the laboratory with a Sony and/or Olympus TG4 camera and then provided with a unique identifier (abbreviation of name, year, location, and number of specimen). Collected specimens were preserved in 96% EtOH or in seawater formaldehyde. For those specimens preserved in formaldehyde, a tiny piece of the foot was cut off first and preserved in EtOH 96% for barcoding. The collections are registered at UNSRAT under the numbers SRU2015/1, SRU2016/1, SRU2016/2, SRU2017/1, SRU2017/2, SRU2018/1, and SRU2018/2. Metadata, including images of all animals, will be published in the Diversity Workbench (https://www.gfbio.org/; Diepenbroek et al., 2014).

Fig. 1

Sampling locations in Sulawesi, Indonesia

Morphological identification

Six hundred and ten specimens were preliminarily identified with the help of available literature (Domínguez et al., 2007; Fahrner & Beck, 2000; Gosliner et al., 2018; Yonow, 2011, 2012). Additional websites such as the Sea Slug Forum were used to check for external variability and distribution data. Subsequently, original literature was used for final identifications (see details for each species in “Results”). The validity of names was checked against the World Register of Marine Species (WoRMS Editorial Board, 2021). Several type and voucher specimens from the Australian Museum (AM Sydney), Muséum national d’Histoire naturelle (MNHN Paris), Natural History Museum of Denmark (ZMUC Copenhagen), Naturalis Biodiversity Center (Leiden), and images from the British Museum of Natural History (NHMUK London) were used for morphological confirmation (see Table S1 for registration numbers): Phyllidia nobilis Bergh, 1869 syntype, paralectotype, and non-type materials (Fig. 2a–d); Phyllidia pustulosa lectotype (Fig. 2e); Phyllidia melanocera Yonow, 1986 holotype (Fig. 2f), Phyllidiella lizae Brunckhorst, 1993 holotype and paratypes (Fig. 2g-i); Phyllidia annulata Gray, 1853 non-type (Fig. 3a–e); Phyllidiopsis pipeki Brunckhorst, 1993 holotype (Fig. 3f), Phyllidiopsis burni Brunckhorst, 1993 holotype (Fig. 3g), and Phyllidiopsis krempfi Pruvot-Fol, 1957 non-type (Fig. 3h).

Fig. 2

Type and museum voucher material of Phyllidiella: a Phyllidia nobilis syntype NHMD-91773; b Phyllidia nobilis non-type NHMD-633656; c Phyllidia nobilis non-type NHMD-633657; d Phyllidia nobilis paralectotypes NHMD-633658a and b; e Phyllidia pustulosa lectotype MNHN-IM-2000–35147; f Phyllidia melanocera holotype NHMUK 1985206/1; g Phyllidiella lizae holotype C159493; h Phyllidiella lizae paratype C162784; i Phyllidiella lizae paratype C159492

Fig. 3

Type and museum voucher material of Phyllidiella and Phyllidiopsis: a Phyllidiella annulata non-type C159519; b Phyllidiella annulata non-type C159520; c Phyllidiella annulata non-type C159520b; d Phyllidiella annulata non-type C159522; e Phyllidiella annulata non-type C159526; f Phyllidiopsis pipeki holotype C162772; g Phyllidiopsis burni holotype C159452; h Phyllidiopsis krempfi non-type RMNH.Mol.336469

DNA extraction, amplification, sequencing, and alignment

DNA from 361 of Phyllidiidae specimens was extracted successfully from the foot, notum, or whole body, depending upon the size of animals (Table S2). DNA isolation was carried out by means of QIAgen® DNeasy Blood and Tissue-Kit, following the manufacturer’s instructions. Partial sequences of mitochondrial CO1 (ca. 680 bp) and ribosomal 16S (ca. 650 bp) were amplified by polymerase chain reaction (PCR) using the primers LCO1490-JJ (5′–CHACWAAYCATAAAGATATYGG-3′) and HCO2198-JJ (5′-AWACTTCVGGRTGVCCAAARAATCA-3′) (Astrin & Stüben, 2008) for CO1, and the primers 16Sar-L (5′-CGCCTGTTTATCAAAAACAT-3′) and 16Sbr-H (5′-CCGGTCTGAACTCAGATCACGT-3′) (Palumbi et al., 2002) for 16S. Amplification of CO1 was performed by an initial step (95 °C for 15 min) followed by 40 touch-down cycles of denaturation (94 °C for 35 s), annealing (55 °C for 90 s), and extension (72 °C for 90 s), with a final extension step 72 °C for 10 min. For 16S rRNA, the PCR began with an initial step (95 °C for 15 min), denaturation (94 °C for 45 s), followed by 34 touch-down cycles, annealing (56 °C for 45 s), extension (72 °C for 90 s), and a final extension step at 72 °C for 10 min. PCR products were sequenced by Macrogen Europe Laboratory (Amsterdam, Netherlands). GENEIOUS Pro 7.1.9 was used to extract the consensus sequences between the primer regions, and for editing and checking the quality of sequences. Consensus sequences were blasted against the NCBI database for evaluation of identification (NCBI Resource Coordinators, 2018). Sequences are deposited in NCBI GenBank, and the numbers are listed in Table S2.

Additionally, 350 phyllidiid sequences available in GenBank were extracted to increase coverage of both species and genera (see Table S3), including many published specimens from our projects (Ompi et al. 2019; Undap et al. 2019; Bogdanov et al., 2020; Papu et al., 2020). Three sequences of the sister taxon of the Phyllidiidae, the Dendrodorididae [Doriopsilla albopunctata (Cooper, 1863), Doriopsilla bertschi Hoover, Lindsay, Goddard & Valdés, 2015, and Dendrodoris atromaculata (Alder & Hancock, 1864)] were obtained from GenBank and used as outgroups. The two genera Reticulidia and Ceratophyllidia were included in the molecular analysis by extracting available sequences from NCBI GenBank since no specimens were collected from Sulawesi. Some available sequences assigned to four species not present in our collections were also extracted from GenBank (Table S3): Phyllidia rueppelii (Bergh, 1869), Phyllidia larryi (Brunckhorst, 1993), Phyllidia babai Brunckhorst, 1993, and Phyllidiopsis cardinalis Bergh, 1876.

Sequences were aligned using the web program MAFFT Alignment, by applying the algorithm FFT-NS-2, scoring matrix 200 PAM/k = 2. The single-gene data sets (632 sequences with 635 bp of CO1, 583 sequences with 538 bp of 16S) were analysed separately, but both genes were also concatenated using GENEIOUS Pro 7.1.9, resulting in a total alignment of 721 sequences with a length of 1173 bp.

Phylogenetic analyses and species delimitation tests

The single-gene data sets as well as the concatenated data set were processed in the online version of IQ-TREE in http://iqtree.cibiv.univie.ac.at/ using default settings for determination of the best-fit substitution model. For our tree reconstruction for the concatenated data set, TVM + F + I + G4 was selected according to the Bayesian information criterion (BIC) (Burnham & Anderson, 2009; Hoang et al., 2017; Kalyaanamoorthy et al., 2017; Nguyen et al., 2015). Additionally, we also ran analyses with models using fewer assumptions (GTR and HKY). The models K3Pu + F + I + G4 and GTR + F + I + G4 were determined according to BIC for the CO1 and 16S data sets, respectively. Phylogenetic reconstruction was carried out using the maximum likelihood algorithms implemented in the IQ tree web server. The following settings (default) were applied: Ultrafast Bootstrap analyses with 1000 replicates, maximum number of iterations 1000, and minimum correlation coefficient for UFBoot convergence criterion 0.99 (Minh et al., 2013). For visualisation of the trees, Dendroscope Ver. 3.5.10 was used (Huson et al., 2012). Collapsing of clades and tree editing was subsequently processed using FigTree v1.4.4 (Rambaut, 2009) and edited in Inkscape (Draw Freely | Inkscape).

To provide further evidence of species hypotheses, data sets of CO1 and 16S were analysed separately with the help of the ABGD (Automatic Barcode Gap Discovery) methodology (Puillandre et al., 2012), using the web server https://bioinfo.mnhn.fr/abi/public/abgd/abgdweb.html. This method helps to characterise barcode gaps between hypothetical species, especially when these are difficult to identify based solely on morphology. The following default settings were used: Initial Partition prior maximal distance (P) = 1.29e-02, Distance K80 Kimura Min Slope = 1.000000. A second program, the bPTP delimitation test (https://species.h-its.org/) using maximum likelihood and Bayesian algorithms was used for the concatenated data to identify cryptic variabilities (Zhang et al., 2013). To visualise the relationship of clades with incongruences in morphological and chemical data, or with short distances between the clades, haplotype networks were investigated using statistical parsimony network methodology (TCS Network method by Clement et al., 2002) in PopART http://popart.otago.ac.nz (Leigh & Bryant, 2015).

Distance values were calculated using MEGA version X (Kumar et al., 2018), using the Kimura 2 parameter model with G4.

Chemical analyses

Based on the molecular tree, seventy phyllidiid specimens were chosen for chemical investigation. These covered 19 clades from three genera. The emphasis is placed on the genus Phyllidiella based on its complexity, and specimens are selected to cover all three locations in North Sulawesi, marked in Fig. S1 by stars.

Extraction and general experimental procedures

The ethanol storage solution was decanted and evaporated under reduced pressure. The bodies were cut into small pieces and repeatedly extracted with methanol (3 × 30 mL) and dichloromethane (1 × 30 mL). All fractions were combined to give the crude extract.

Mass spectra were recorded on a micrOTOF-QIII mass spectrometer (Bruker) with ESI-source coupled with a HPLC Dionex Ultimate 3000 (Thermo Scientific) using an EC10/2 Nucleoshell C₁₈ 2.7 µm column (Macherey–Nagel). The column temperature was 25 °C. MS data were acquired over a range 100–3000 m/z in positive mode. Auto MS/MS fragmentation was achieved with rising collision energy (35–50 keV over a gradient of 500–2000 m/z) with a frequency of 4 Hz for all ions over a threshold of 100. HPLC began with 90% H₂O containing 0.1% acetic acid. The gradient started after 1 min to 100% acetonitrile (0.1% acetic acid) in 20 min. A 5 µL amount of a 1 mg/mL sample solution (MeOH) was injected at a flow of 0.3 mL/min. All solvents were of LCMS grade. NMR spectra were recorded in MeOH-d₄ using Bruker Avance 300 DPX or Bruker Ascend 600 spectrometers. Spectra were referenced to residual solvent signals with resonances at δ_H/C 3.35/49.0. Preparative HPLC was performed on a Merck Hitachi HPLC system equipped with a L-6200A pump, a L-4500A PDA detector, a D-6000A interface with D-7000A HSM software, a Rheodyne 7725i injection system and a Nucleodur C₁₈ 5 µm Pyramid 250 mm × 10 mm column (Macherey–Nagel).

MS data acquisition and analysis

HR-LCMS analysis was performed using a 1 mg/mL LCMS grade methanol solution of the crude extracts. The obtained raw.d LCMS/MS files were converted using the DataAnalysis 4.2 software (Bruker Daltonic) to.mzXML format and uploaded to the MassIVE server (https://massive.ucsd.edu) via FileZilla (https://filezilla-project.org) FTP client. LCMS data files are publicly accessible under ftp://massive.ucsd.edu/MSV000087294/.

A classical molecular network was created using the online workflow (https://ccms-ucsd.github.io/GNPSDocumentation/) on the GNPS website (http://gnps.ucsd.edu). The data were filtered by removing all MS/MS fragment ions within + / − 17 Da of the precursor m/z. MS/MS spectra were window filtered by choosing only the top six fragment ions in the + / − 50 Da window throughout the spectrum. The precursor ion mass tolerance was set to 0.05 Da and a MS/MS fragment ion tolerance of 0.02 Da. A network was then created where edges were filtered to have a cosine score above 0.6 and more than three matched peaks. Feature-Based Molecular Network (FBMN) was created with their workflow (Nothias et al., 2020) on GNPS (Wang et al., 2016). The mass spectrometry data were first processed with MZmine 2.53 (Pluskal et al., 2010), and the results were exported to GNPS for FBMN analysis. The precursor ion mass tolerance was set to 0.05 Da and the MS/MS fragment ion tolerance to 0.05 Da. The network was created with edges filtered to have a cosine score above 0.6 and more than four matched peaks. Furthermore, edges between two nodes were retained in the networks if and only if each of the nodes appeared in each other’s respective top ten most similar nodes. Finally, the maximum size of a molecular family was set to 100, and the lowest scoring edges were removed from molecular families until the molecular family size was below this threshold. The spectra in the network were then searched against GNPS spectral libraries (Horai et al., 2010; Wang et al., 2016). The library spectra were filtered in the same manner as the input data. All matches retained between network spectra, and library spectra were required to have a score above 0.7 and at least six matched peaks. The molecular networks were visualised using Cytoscape software (Shannon et al., 2003).

Isolation of pure compounds

The following protocol was used for the isolation of pure compounds: crude extracts were re-suspended in 30 mL H₂O and extracted three times with 30 mL ethyl acetate (EtAc) in a separation funnel. EtAc solubles were collected and solvent-evaporated under reduced pressure. The obtained lipophilic salt-free part underwent fractionation on a reversed phase SPE Bakerbond 2000 mg column using a stepwise gradient (50:50, 70:30, 90:10, 100:0 MeOH:H₂O v/v, 20 mL each) to obtain four fractions. In specimens Phpu15Bu14 and Phpu16Sa66 fractions SPE3 were found to be pure 3-isocyanotheonellin (18 mg) and 7-isocyano-7,8-dihydrobisabolene (14 mg), respectively. Theonellin formamide (1 mg) was isolated from the Phpu15Bu14 SPE2 fraction using HPLC and isocratic MeOH:H₂O (70:30 v/v) mobile phase. Crude salt-free extracts from five Phyllidiella nigra (van Hasselt, 1824) specimens (Phpu15Bu-6, -10, -18, -19, and -43) underwent the same SPE fractionation procedure. SPE fractions were further separated using HPLC to yield 3-isocyanobisabolane-7,10-diene and 3-isocyanobisabolane-7(14),11-diene in mixture (1.0 mg), epimeric 9-thiocyanatopupukeananes (1.0 mg), and 10-isocyano-5-cadinen-4-ol (1.4 mg). Fractions containing the same metabolites from several individually processed specimens were combined to obtain sufficient amounts for unambiguous structure elucidation with NMR.

Obtained LC–MS data were analysed with web-based cheminformatic tools (classical molecular networking and feature-based molecular networking) associated with the GNPS platform (gnps.ucsd.edu) that allowed rapid compound identification via comparison of the experimental MS/MS spectra with an integrated library of more than 220,000 curated spectra.

Chemical structures of known compounds could be determined after their isolation and ¹H and ¹³C NMR comparison with literature values. In addition to the GNPS-aided identification, structures of compounds could be putatively assigned via manual analysis of mass spectrometric data and database search (MarinLit, Dictionary of Natural Products), especially in cases where minute extract amounts, chemical instability, or high volatility hindered isolation.

Results

Morphological identification

External identification of the specimens collected during the four expeditions allowed recognition of three genera, namely Phyllidia with 11 species (Figs. 4–6), Phyllidiopsis with seven species (Figs. 7, 8), and Phyllidiella with at least 11 species (Figs. 9–12).

Fig. 4

Phyllidia species and specimens with identifiers. Scale bars: 10 mm. 1a–d Phyllidia elegans: a Phel18Ba2; b Phel18Bu1; c Phel15Bu4; d Phaly16Bu1. 2a–e Phyllidia coelestis: a Phco15Bu5; b Phco18Po1; c Phco15Bu12; d Phco15Bu16; e Phpic16Sa2. 3a Phyllidia sp. 3: a Phex15Bu1. 4a–f Phyllidia picta: a Phpi18Bl3; b Phsp18Bu1; c Phpi16Sa1; d Phsp18Ba3; e Phpi18Po1; f Phsp616Sa1

Fig. 5

Phyllidia species and specimens with identifiers. Scale bars: 10 mm. 1a, b: Phyllidia exquisita: a Phex17Ba1; b Phex18Bl2. 2a, b Phyllidia cf. babai: a Phoc18Ba7; b Phoc16Sa3. 3a–e Phyllidia haegeli: a Phco15Bu1; b Phpic16Sa4; c Phva16Bu2; d Phva16Sa51; e Phva16Sa19. 4a Phyllidia sp. 9: a Phsp15Bu1. 5a–g Phyllidia ocellata: a Phoc18Ba4; b Phoc15Bu1; c Phoc16Sa2; d Phoc16Sa4; e Phoc16Sa7; f Phoc17Ba1; g Phma16Sa

Fig. 6

Phyllidia species and specimens with identifiers. Scale bars: 10 mm. 1a Phyllidia sp. a: a Phsp17Bu1. 2a–i Phyllidia varicosa: a Phel15Bu1; b Phva15Bu8; c Phel18Bu6; d Phel18Ko3; e Phva18Po1; f Phva18Ba1; g Phva15Bu2; h Phva16Sa37; i Phva16Sa38

Fig. 7

Phyllidiopsis species and specimens with identifiers. Scale bars: 10 mm. 1a–j Phyllidiopsis krempfi: a Phskr18Ba1; b Phfi16Sa2; c Phspi18Ba1; d Phpu16Sa72; e Phsh15Bu1; f Phfi16Sa1; g 18TMphni7; h Phpu15Bu37; i Phpu16Sa58; j Phpu16Sa66. 2a–d Phyllidiopsis burni: a Phssp.a18Ba1; b Phpu16Sa19; c Phpu16Sa47; d Phpu16Sa55. 3a–d Phyllidiopsis sp. a: a Phskr18Ba2; b Phpu16Sa42; c Phpu16Sa45 d Phpu16Sa59

Fig. 8

Phyllidiopsis species and specimens with identifiers. Scale bars: 1–3 1 mm; 4 10 mm. 1a–d Phyllidiopsis xishaensis: a Phst16Bu1; b Phst16Bu2; c Phsxi18Ba1; d Phsxi18Bu1. 2a–c Phyllidiopsis annae: a Phsan18Ba1; b Phsan18Ba2; c Phsan18Ba4. 3a Phyllidiopsis sphingis: a. Phsph15Bu1. 4a–c Phyllidiopsis shireenae: a 18LePhsh1; b Phsh16Sa2; c Phpi15Bu4

Fig. 9

Phyllidiella species and specimens with identifiers. Scale bars: 10 mm. 1a–d Phyllidiella pustulosa: a Phpu15Bu8; b Phpu15Bu33; c Phpu18Ba12a; d Phpu17Ba2. 2a–d Phyllidiella cf. pustulosa: a Phpu18Ko7; b Phpu18Ba3; c Phpu18Ko8; d Phpu18Sm1. 3a–f Phyllidiella nigra: a Phpu15Bu7; b Phni18Po1; c Phpu18Ko9; d Phpu18Ba19; e Phpu17M3; f Phpu17Bu1. 4a–d Phyllidiella sp. a: a Phpu18Po2; b Phpu18Bu4; c Phpu16Sa9; d Phpu16Sa1

Fig. 10

Phyllidiella species and specimens with identifiers. Scale bars: 10 mm. 1a–f Phyllidiella zeylanica auctt.: a Phan16Bu1; b Phli18Bu1; c Phze18Bu1; d Phpu15Bu41; e Phpu18Bl4; f Phpu15Bu25. 2a, b Phyllidiella rudmani: a; Phpi15Bu2; b Phpi15Bu5. 3a–d Phyllidiella sp. b: a Phpu16Sa25; b Phpu16Sa39; c Phpu16Sa27; d Phpu18Bu5. 4a–d Phyllidiella sp. c subclade 1: a Phli16Sa8; b Phpu15Bu35; c Phli18Ba3; d Phpu18Po4

Fig. 11

Phyllidiella species and specimens with identifiers. Scale bars: 10 mm. 1a–e Phyllidiella sp. c subclade 2: a Phpu18Ba12; b Phan15Bu3; c Phli16Sa5; d Phli18Ba1; e Phpu16Sa11. 2a–f Phyllidiella albonigra: a Phpu15Bu5; b Phpu16Bu8; c Phpu16Bu3; d Phpu15Bu21; e Phpu16Sa32; f Phpu17Ba4. 3a–d Phyllidiella sp. d: a Phpu18Po1; b Phpu18Ba1; c Phpu18Po3; d. Phpu18Bl5. 4a–d Phyllidiella species complex e: a Phpu15Bu28; b Phpu16Sa35; c Phpu15Bu39; d Phpu18Bl2

Fig. 12

Phyllidiella species and specimens with identifiers. Scale bars: 10 mm. 1a–d Phyllidiella hageni: a Phllsp15Bu1; b Phru18Ko1; c Phpu16Sa30; d Phha18Ko1

Phylogenetic analyses

The trees reconstructed from the concatenated data set of the two mitochondrial genes 16S and CO1 using three different model assumptions (GTR, HKY, and TVM + F + I + G4) resulted in similar trees at genus level; however, they differed in many ways at the species level. The tree based on the highest number of assumptions (TVM + F + I + G4), allowing different substitutional modes of the gene partitions, best reflected the results obtained from morphological analyses and was mainly used for the subsequent analysis. Figure S1 illustrates this tree in full with specimen identifiers, and Fig. 13 summarises the results of the concatenated data set by collapsing the nodes at species level. Additionally, Figs. 14 and 15 show the summarised results based on single-gene analyses, CO1 (using the model K3Pu + F + I + G3) and 16S (GTR + F + I + G4), respectively. Phylogenetic reconstructions using the mitochondrial genes in single (Figs. 14, 15) or concatenated (Figs. 13, S1) data sets usually resulted in the monophyletic genera Phyllidia, Reticulidia, Phyllidiella, and Phyllidiopsis: Phyllidia (with the exception of Phyllidia larryi Brunckhorst, 1993) is sister taxon to Reticulidia in all analyses with the various settings with high bootstrap support values. Unfortunately, only CO1 genes were available in GenBank for Phyllidia larryi (type locality Guam), which was not present in our collections. This species does not group with Phyllidia, but as a sister taxon to the Phyllidia/Reticulidia clade (Figs. 13, 14).This position is obtained in the tree reconstruction based only on the CO1 gene, as well as in the concatenated data set, and its assignment to the genus Phyllidia needs further investigation. Morphologically, this species also differs considerably from all the others in being almost smooth and having a pale yellow dorsum with random red lines. Phyllidiella and Phyllidiopsis (with exception of Phyllidiopsis cardinalis Bergh, 1875, type locality Tonga) are sister taxa in the concatenated and 16S data sets, whereas in the CO1 analysis, Phyllidiella is polyphyletic, with one group being sister taxon to the Phyllidia/Reticulidia clade, and the other part being sister to Phyllidiopsis. In our analysis, Phyllidiopsis is monophyletic with one exception: the specimens and sequences assigned to P. cardinalis by Valdés (2003; AF430367) and Cheney (2014; KJ001308) branch off first within the Phyllidiidae tree when using only 16S and in the concatenated data set (Figs. 13, 15). However, when applying only the CO1 data set, P. cardinalis groups within its genus in the tree (Fig. 14). Since information on 16S is lacking for KJ001308 and CO1 is lacking for AF430367, the correct affiliation needs to be investigated in future studies with more complete data sets. The single sequence of a Ceratophyllidia from GenBank was always nested within the genus Phyllidiella with high bootstrap support (100).

Fig. 13

Maximum likelihood tree of Phyllidiidae based on the concatenated data set of partial CO1 and 16S sequences. Terminal taxa on species level collapsed. Numbers at nodes are bootstrap values (1000 replicates). The genera are marked with the same colours as in Fig. S1. Numbers behind the species names indicate results of the various species delimitation tests. One number indicates that all species delimitation tests showed the same result on species numbers; four numbers indicate the species calculated according to (1) ABGD test based on CO1, (2) ABGD test based on 16S, (3) ML bPTP, and (4) Bayesian bPTP tests based on the concatenated data sets

Fig. 14

Maximum likelihood tree of Phyllidiidae based on partial CO1 sequences. Terminal taxa on species level collapsed. Numbers at nodes are bootstrap values. The genera are marked with the same colours as in Fig. S1

Fig. 15

Maximum likelihood tree of Phyllidiidae based on partial 16S sequences. Terminal taxa on species level collapsed. Numbers at nodes are bootstrap values. The genera are marked with the same colours as in Fig. S1

Species delimitation tests

Four species delimitation tests were run with different algorithms and different data sets. In general, the ABGD test based on the 16S gene resulted in the most conservative species hypotheses, with 16 clades of Phyllidia, 7 Phyllidiopsis (without Phyllidiopsis sp. a where no CO1 sequences could be retrieved), and 20 Phyllidiella clades. Interestingly, the bPTP tests when applying the concatenated data set (ML and Bayesian approach) resulted in similar species delimitation as the ABGD test using only the CO1 data set. However, the results of both bPTP tests appear to overrate the genetically identified clades as distinct species when compared with morphological results. Both bPTP analyses resulted in more than 100 species, with particularly heavy splitting in Phyllidia ocellata, P. varicosa, and Phyllidiella hageni Fahrner & Beck, 2000. The Bayesian approach led to even more splitting than the ML analysis, especially in P. varicosa. The ABGD test based on CO1 mainly showed a stronger splitting in Phyllidiella pustulosa. We provide species descriptions based on the results of the ABGD test using the 16S gene (Fig. S1) and consulted the haplotype network analyses, which supported our more conservative species hypotheses.

Many phyllidiid species were investigated chemically for the first time in this study: Phyllidia elegans Bergh, 1869, P. cf. babai, P. haegeli Fahrner & Beck, 2000, Phyllidiopsis burni Brunckhorst, 1993, P. shireenae Brunckhorst, 1990a, Phyllidiella nigra (van Hasselt, 1824), P. rudmani Brunckhorst, 1993, and several unnamed species of Phyllidiella. High-resolution LC-ESIMS analysis provided unprecedented insights into the metabolomes of 70 specimens across the 19 species from three phyllidiid genera. The generated molecular networks based on MS/MS spectra similarity visualise the complex chemical composition of these 70 metabolomes (see Figs. 16, S2, S3). A large array of natural products (Fig. S4a, b) is detected in individual extracts which are outlined in more detail in the species descriptions and summarised in Fig. 16. The majority of identified compounds belong to sesquiterpenes functionalised with isonitrile, isocyanate, formamide, thiocyanate, and isothiocyanate moieties (see reviews by Garson & Simpson, 2004 and Emsermann et al., 2016). In ESI mass spectrometry, these molecules ionise favourably under loss of the respective nitrogen-containing functional group (X = CN, CNO, HCONH₂, SCN, CNS) and exhibit characteristic high-intensity in-source fragment ions with m/z of 205.195 [M-X]⁺ with low-intensity molecular peaks with m/z [M + H]⁺ 232.205, 248.200, 250.216, and 264.178 that are diagnostic and allow assignment of the corresponding functional group (Fig. S4c). Interestingly, m/z values 205.195 and 232.205 are often accompanied by m/z 259.220 ions. Mass differences of 27 Da can be explained by the presence of an additional isonitrile moiety in the molecule. Bi- and even tri-functionalised diterpene isonitriles are well documented from marine sponges (Patra et al., 1984); however, there are no reports of naturally occurring sesquiterpenes bearing two isonitrile moieties to date.

Fig. 16

Phylogeny of Phyllidiidae and chemical analysis of 70 preserved specimens from three genera (Phyllidia, Phyllidiopsis, and Phyllidiella) and 19 species. a Coloured clusters were obtained with classical or feature-based molecular networking and depict families of secondary metabolites characteristic for certain clades, e.g., brominated alkaloids for the genus Phyllidia, amphilectene-type diterpenes for P. coelestis and P. elegans, kalihinol-type diterpenes for P. rudmani, as well less specific secondary metabolites (sesquiterpene isonitriles and derivatives) that were detected throughout the Phyllidiidae. The chemical structures were established either with NMR after compound isolation, or were assigned based on HR-LCMS/MS analysis, MarinLit database search, and comparison with reference spectra in the GNPS repository. b Visualisation of the chemical space as present in phyllidiid extracts using classic molecular networking with HR-LCMS/MS data. Each node represents a parent ion (molecule). Based on similarities in fragmentation patterns; the nodes are connected with edges producing clusters/compound families. Nodes are colour-coded according to the legend. When the same compounds are detected in more than one species, the nodes are displayed as pie charts reflecting the sum intensity of the observed ions in the respective species

Fig. 17

Haplotype network of Phyllidia picta including 34 CO1 sequences. Note the separation of the two clades by 28 mutational steps. The two clades are coloured red and green, respectively

Sesquiterpene isonitriles and formamides were the most common metabolites and could be detected in almost every single chemically analysed specimen. Besides isonitrile sesquiterpenes a series of halogenated compounds with m/z values 266.166 [M + H] (containing chlorine) and 310.115 [M + H]⁺ (containing bromine), characteristic isotopic patterns could be detected in specimens across all investigated phyllidiid species and genera. MS/MS spectra indicate that these molecules share the same backbone and differ only in the type of halogen attached. Intriguingly, no halogenated natural products having the above-mentioned masses are known to date. These compounds were encountered in low concentrations, and despite their wide distribution, isolation of a sufficient amount for unambiguous structure elucidation was not possible.

Whereas the results of metabolomic analyses of some investigated phyllidiids were difficult to interpret, chemical extract composition of Phyllidia coelestis, P. elegans, P. ocellata, P. cf. babai, Phyllidiella albonigra (Quoy & Gaimard, 1832), P. nigra (van Hasselt, 1824), P. pustulosa, and P. rudmani Brunckhorst, 1993 was found to be clade specific. Thus, it can provide chemotaxonomic evidence for species delimitation and will be described in detail.

Integrative species descriptions

After a short summary of the molecular analyses of the respective genera with major results on species and chemical compounds, species results are presented in detail, the listing presented in the species order of the tree shown in Figs. 13 and S1. After the detailed morphological and colour descriptions, including discussions of similar species, the molecular results, species delimitations, and/or haplotype results are presented, followed by results on chemical compounds.

Phyllidia Cuvier, 1797

Morphological identification resulted in 11 Phyllidia species (Figs. 4–6), seven of which are recognised species: P. elegans Bergh, 1869, P. coelestis Bergh, 1905, P. picta Pruvot-Fol, 1957, P. exquisita Brunckhorst, 1993, P. haegeli Fahrner & Beck, 2000, P. ocellata Cuvier, 1804, and Phyllidia varicosa Lamarck, 1801. Two species have been illustrated as new undescribed species in the available literature, Phyllidia sp. 3 and Phyllidia sp. 9 (Gosliner et al., 2015, 2018). The third new species, which we call Phyllidia sp. a in this work (Fig. 6.1a), was figured by Gosliner et al. (2008: 300, top right figure only) and Venkataraman et al. (2015: pl. 75), both as Phyllidiopsis monacha (Yonow, 1986), and the same specimen (Phsp17Bu1) was illustrated as Phyllidia sp. in Eisenbarth et al. (2018). A recent corrigendum validated the original generic placement of Phyllidia monacha based on the original description and drawings, in combination with recent photographs of the type specimen (Yonow, 2021), thereby indicating possibly greater similarities with our new species.

Molecular analyses confirmed our morphologically based assignments of species to a large extent. Our study also confirms that species previously assigned to the genus Fryeria based on their ventral anal opening are part of Phyllidia and do not form a separate clade. Although species are clearly marked by long branches, branch lengths between species are rather short with low bootstrap values (Figs. 13, S1). Higher support values are provided for the sister-taxa relationship of P. elegans/P. coelestis (99), P. haegeli/P. sp. 9 (100), and P. ocellata/P. sp. a. However, we would like to emphasise here that statements on relationships between species are preliminary because many recognised worldwide species could not be included in this analysis.

LCMS analyses revealed a range of brominated compounds that occur in several Phyllidia species in addition to commonly encountered sesquiterpene isonitriles. These comprise bromoindole alkaloid 5-bromotryptophan (m/z values 283.006 and 285.005 [M + H]⁺, and 265.983 and 267.978 [M-NH₃]⁺) which was detected in all Phyllidia species (Fig. S4b). Obtained MS data perfectly match values in the original description of Smenospongia sp. (Porifera) metabolites by Tasdemir et al. (2002) and allowed confident structure assignment of bromoindoles. In addition, double-charged ions of unidentified polar mono- (m/z 235.020 and 242.029 [M + 2H]²⁺) and dibrominated (m/z 219.035 [M + 2H]²⁺) compounds with retention times of 3.5–5.5 min were observed in all Phyllidia extracts except those of P. varicosa. Some of these brominated metabolites could be detected in low amounts in the extracts of six other specimens belonging to the genera Phyllidiopsis (Phyllidiopsis krempfi) and Phyllidiella (Phyllidiella sp. c subclade 1, P. zeylanica auctt., and P. rudmani).

Phyllidia elegans Bergh, 1869

The 40 specimens of Phyllidia elegans (Fig. 4.1a–d) are oval and broad, somewhat flattened in shape, and their dorsal background colours are white, grey, or pale green with narrow black lines. The white granulated tubercles form three interrupted ridges along the back and are usually capped with yellow or orange. The rhinophores are golden yellow or orange, the same colour as the tubercles. The foot sole possesses a longitudinal black line; the elongated tapering oral tentacles are pale white to cream. Some of our specimens were green and not pink as described in the literature (e.g., Brunckhorst, 1993), and some had no gold colouration on the dorsum but the rhinophores were gold (e.g., Phaly16Bu1, Fig. 4.1d). Of the 40 specimens collected from Sulawesi, we were able to sequence 34 that group with an additional 22 sequences of P. elegans downloaded from NCBI. Intraspecific variability is up to nearly 4% within the monophyletic clade (Table S4). The metabolome of P. elegans is very similar to P. coelestis and is described in detail below.

Phyllidia coelestis Bergh, 1905

Fifty specimens of Phyllidia coelestis (Fig. 4.2a–e) collected around Sulawesi exhibit the typical colouration as described by Brunckhorst (1993) and Yonow (2012). The animals all have three longitudinal black lines on the notum and a general bluish background that is granulated. The middle line is usually interrupted by yellow-capped tubercles, but these central tubercles are encircled by the median black line in only two of our 50 specimens, rendering the usually interrupted line into a continuous line (Phco18Po1, Fig. 4.2b). The median line is divided by an anterior tubercle, thus forming the species-specific Y-shaped black marking in front of the rhinophores. The lateral black lines do not unite in either the anterior or posterior parts of the body. Other large and usually single rounded tubercles, capped by yellow, lie on blue ridges between the black lines. The tubercles lying laterally are smaller and can be yellow-capped or bluish; additionally, there are tiny black flecks in this lateral blue area. The species always lacks a black line on the foot. One pale specimen has yellow tubercles protruding only from the central black line (Phco15Bu16, Fig. 4.2d). There was only one interesting and different specimen within the P. coelestis clade: one with the identifier Phpic16Sa2 (Fig. 4.2e) differs in its black background colouration, resembling Phyllidia madangensis Brunckhorst, 1993. A dark form of P. coelestis was described and illustrated from Australia in Yonow (2011): externally they differed in having a completely black central area with three rows of isolated, orange-tipped tubercles, and differences in the internal anatomy from the typical P. coelestis. Some individuals of P. coelestis with darker colouration are illustrated on Sea Slug Forum (Cobb & Rudman, 2008) and GBIF (2019); however, it should be emphasised that these refer only to externally identified individuals. Nonetheless, 50 specimens of P. coelestis from Sulawesi are included in our molecular analyses with further eight sequences retrieved from NCBI from Lizard Island, Australia (Cheney et al., 2014), and New Caledonia (Valdés, 2003) which clearly group with our sequences with a bootstrap value of 100 (concatenated tree, Figs. 13, S1). Intraspecific variability lies within 2%. Phyllidia coelestis and P. elegans form a monophyletic clade supported by a bootstrap value of 99.

Metabolomic analysis of P. coelestis and P. elegans specimens collected in Bunaken National Park in 2015 demonstrated very similar characteristic chemotypes predominant in these two species (Fig. S7a, b), thus supporting their sister-taxa relationship (bootstrap in the concatenated data set is 99). HR-LCMS analysis of crude extracts revealed a series of major peaks with m/z values that can be assigned to amphilectene-type diterpene isonitriles (m/z [M + H]⁺ 298.250), diisonitriles (m/z [M + H]⁺ 325.259), and formamides (m/z [M + H]⁺ 316.261) (Fig. S4a) that were recently reported from two Hainan P. coelestis populations (Carbone et al., 2019). The amphilectene diterpenes detected in P. elegans and P. coelestis form a separate cluster in molecular networks (Figs. 16, S2, S3). Metabolomic analysis using GNPS tools also facilitated an intriguing detection of dichloroimidic sesquiterpenes in P. elegans that were characterised by us previously (Bogdanov et al., 2020). Even though present in trace amounts, dichloroimidic sesquiterpenes can provide a chemotaxonomic separation between the two closely related species P. coelestis and P. elegans. Two specimens (Phco15Bu6, Phel15Bu8) among the chemically analysed specimens of both species possess an additional unknown metabolite with an m/z 404.187 [M + H]⁺ and a retention time of 13.9 min. Interestingly, the same metabolite is detected in only one specimen of Phyllidiella sp. c subclade 2 (Phpu15Bu27). Minute extract amounts prevented us from isolating and characterising this compound.

Phyllidia sp. 3

Our two specimens grouped with one sequence extracted from GenBank (AF430363), an animal collected from 35 m depth in New Caledonia and assigned to P. ocellata by Valdés (2003). These three specimens form a distinct clade, confirmed by species delimitation tests, and form a sister clade to P. elegans/P. coelestis. One of our specimens, Phex15Bu1 (Fig. 4.3a), has the same pattern as Phyllidia sp. 3 in Gosliner et al. (2018): it has a bluish tinge with an irregular yellow mantle rim and two longitudinal black lines that are slightly wavey, a black curved U-shaped line in front of the rhinophores, and black marks on the dorsal midline. The tubercles are large, single, and round, not topped with yellow pigment: in fact, all tubercles are white in this species. The animal depicted in Gosliner et al. (2018) only differs in so far as the yellow parts appear more orange. Unfortunately, no photograph is available of our living second specimen, but the black pattern is similar. The clade is supported by a bootstrap value of 100 (Fig. 13) and shows no intraspecific variability (Table S4). No chemical analyses were performed for specimens of this clade.

Phyllidia picta Pruvot-Fol, 1957

Thirty-four sequences of Phyllidia picta were included in our molecular phylogeny, which grouped in two distinct subclades, indicating an intraspecific variation. This is confirmed by the haplotype network analysis in which the two clades are separated by 28 differing nucleotides which are shared only by specimens in one of the respective clades and not with the other (Fig. 17). Species delimitation tests, except the one based on 16S, also recognised two distinct species. Since we follow the more conservative species concept in our study, and in addition to the lack of support for these two clades as separate species in the 16S data set and morphologically, we discuss these two clades as one species. In the first clade, all sequences from GenBank group together with sequences from our collection from North and South Sulawesi (Fig. 4.4a–c). The second clade includes specimens from our expeditions to North and Southeast Sulawesi (Fig. 4.4d–f). The intraspecific genetic variability of both clades within P. picta is 9%. Both clades are very similar in their external morphologies, with the anal opening positioned ventrally. They resemble each other in colouration with a typical black trapezoidal pattern with elongations towards the edge of the mantle. The tubercles are single, the larger ones in three median rows usually capped with yellow to orange. The rhinophores are also yellow to orange with one tubercle anterior to them. The oral tentacles are grey, the tips with or without yellow in both clades; the foot sole is grey to dark grey. The colouration of the living specimens matches the colouration described for Phyllidia menindie (Brunckhorst, 1993), a species that was subsequently synonymised with P. picta Pruvot-Fol, 1957 by Yonow (1996) after her examination of the holotype in the Natural History Museum (BMNH 1854.7.19.91). No chemical analyses were performed for this species.

Phyllidia rueppelii (Bergh, 1869)

The two clades of Phyllidia picta group in an unresolved relationship with a single available sequence (AF430358) assigned to Phyllidia rueppelii from Valdés (2003; 16S gene only) (Figs. 15, S1). This species is only known from the Red Sea and Gulf of Oman (Valdés, 2003; Yonow, 2020), but its identity has been confused with P. schupporum occurring in the same geographical range (see Yonow, 2020). Like P. picta, P. rueppelii has a ventral anal opening but differs in having a yellow-orange margin to the mantle and thus can be easily confused with P. schupporum. Unfortunately, no CO1 sequence was available for this Red Sea specimen, nor can we confirm correct identification by photograph (Á. Valdés pers. comm.).

Phyllidia exquisita Brunckhorst, 1993

Only three specimens of Phyllidia exquisita were present in our collections (Fig. 5.1a, b). The main background colour is grey to white, with two or four distinct longitudinal black lines, forming a reticulate pattern across the midline in one specimen (Phex18B12, Fig. 5.1b). Each lateral black line forms 4–6 scallops to the mantle margin. In some specimens, some single, black, small- to middle-sized, round tubercles are present, highly disguised, on the black lines. The white tubercles are usually single and the larger ones are capped by yellow-gold. The yellow mantle margin is characteristic of this species, mainly along the anterior rim, more seldom along the sides. The anal opening lies on the greyish white part of the notum, behind the last large posterior tubercle. Our three specimens group with the only P. exquisita sequence available from GenBank (KY235923) provided by Stoffels et al. (2016) from a specimen (not illustrated) from Ternate, Indonesia.

Phyllidia cf. babai

Two of our specimens with the identifiers Phoc18Ba7 and Phoc16Sa3 (Fig. 5.2a, b) have a very similar colouration to P. ocellata, but actually group with two sequences from Stoffels et al. (2016), who assigned their sequences to Phyllidia cf. babai. This monophyletic clade is supported by a bootstrap value of 100 with an intraspecific genetic variability of 5.67%. The animals that Stoffels et al. (2016) depicted are similar to specimens originally described as P. babai Brunckhorst, 1993 with a white ground colour and only the central row of tubercles being yellow or orange. Brunckhorst (1993) also mentioned a distinct yellow margin for his newly described P. babai, a feature that seems to be lacking in the specimens analysed by Stoffels et al. (2016). We assume that the second specimen depicted in their figure 10e is from a preserved animal and therefore lost the yellow pigment. Our specimens differ from their P. cf. babai specimens by having a yellow or orange background colour identical to both their P. ocellata and our P. ocellata. Moreover, our specimens of this clade have white lines surrounding the black patches that are always continuous between the black rings and form a central white oval line around the three tubercular rows. In P. ocellata, the coloured rings are always separated by the orange background colour (see P. ocellata below). One of our P. cf. babai specimens (Phoc18Ba7, Fig. 5.2a) had an orange foot and orange oral tentacles, and not the typical cream to white ventral colouration of P. babai, nor the darker grey colouration of foot and oral tentacles of P. ocellata. Therefore, our specimens do not match P. cf. babai in colour, although they group with the two sequences of Stoffel et al. (2016); they also do not match the colour description of P. babai in Brunckhorst (1993). Investigation of more specimens is needed to find better distinguishing characters for this molecularly distinct clade, which we retain as P. cf. babai.

Chemical analysis of one P. cf. babai specimen (Phoc16Sa3, Fig. 5.2b; Undap et al., 2019: fig. 3B) and one P. ocellata specimen (Phoc16Sa7, Fig. 5.5e) demonstrates clearly different metabolomes despite external similarities (see Fig. S7c). A unique feature detected only in the metabolome of P. cf. babai among all chemically analysed phyllidiids is a series of polar brominated compounds with m/z values matching masses of the bromopyrrole alkaloids manzacidin A ([M + H]⁺ 344.024 and 346.022) and manzacidin B ([M + H]⁺ 360.017 and 362.014; see Fig. S4b). These ions also form a separate cluster in the GNPS molecular network. Furthermore, a MS² spectrum of a peak (m/z 403.917 M⁺) with a retention time of 4.9 min and an isotopic pattern of a dibrominated compound matched perfectly with the MS² spectrum of spongiacidin A in the GNPS library (see Fig. S5). Its monobrominated derivative spongiacidin B (m/z 324.009 M⁺, retention time of 3.7 min) was also detected. These alkaloids were described from an Okinawan sponge species of the genus Hymeniacidon (Inaba et al., 1998; Kobayashi et al., 1991). Additionally, the MS/MS spectrum of a dibrominated compound with an m/z 389.938 [M + H]⁺ and retention times of 4.8 and 6.8 min matches the GNPS repository spectrum of dibromophakellin (Fig. S6), an alkaloid from a sponge identified as Acanthella flabellata (Tanita, 1961) (Sharma & Magdoff-Fairchild, 1977 as Phakellia). Common sesquiterpene isonitriles and formamides were also detected in the crude extract of P. cf. babai.

Phyllidia haegeli (Fahrner & Beck, 2000)

Phyllidia haegeli (Fig. 5.3a–e) is represented in our collections by seven specimens. The species possesses an anal opening on the ventral side, but the colouration and tubercular pattern resemble that of P. varicosa. The species is characterised by four black lines running along three longitudinal tubercular ridges, which meet behind the rhino-tubercle and usually meet behind the second median tubercle. Shorter horizontal black lines may originate from the longitudinal stripes and extend to the margin resembling the trapezoidal markings of P. picta or P. elegans. Single rounded tubercles on the ridges are capped by yellow, often with confluent bases, whereas the small tubercles on the blue marginal band may lack this yellow cap. One distinct yellow tubercle lies anteriorly between rhinophores. The marginal band resembles that of P. coelestis in being blue, somewhat granular, with tiny tubercles and flecks of black. The foot sole is grey to blue with a somewhat darker stripe or band down the middle. We assume that there was no such band remaining on the type specimen when it was examined: Fahrner and Beck (2000) clearly state that the specimen was not photographed alive (as they thought it was a P. varicosa, presumably with a ventral black line?) and that the black pigment disappeared after some time. Figure 5.3d depicts one of our P. haegeli specimens (Phva16Sa51) in the field showing the high degree of similarity to P. varicosa and the probable misidentifications that will have occurred when only photographs are examined (e.g., on internet sites such as Sea Slug Forum and iNaturalist) when not specifically searching for the position of the anal opening. In future, each specimen of P. varicosa will have to be carefully examined. Phylogenetic analyses support the monophyletic clade with a bootstrap value of 100 and a rather low intraspecific variability of 0.7% (Table S4).

One P. haegeli specimen (Phco15Bu1, Fig. 5.3a) was analysed chemically using LC-HRMS (Fig. S7a). Its metabolomic profile is dominated by sesquiterpenes functionalised with isonitrile and formamide moieties, whereas amphilectene-type diterpenoids found in P. coelestis and P. elegans are absent.

Phyllidia sp. 9

The sister group to P. haegeli is a clade composed of only two specimens. One is from our collection (Phsp15Bu1, Fig. 5.4a), and one is from the study of Stoffels et al. (2016: fig. 7g, RMNH.Moll.336619), both with a colour pattern similar to the individual depicted by Gosliner et al. (2015: 282) under the name Phyllidia sp. 9. The colour of our animal is bright white, with a distinctive black oval area in the middle of the mantle. This black area is disrupted by single white tubercles, usually capped with yellow, arranged in three rows. The yellow rhinophores arise from a broad white base encompassing the rhinophore and a rhino-tubercle in the black area. The animal depicted by Stoffels et al. (2016: fig. 7g) does not show any yellow traces, presumably because it is a photograph of the preserved specimen, and is completely white and black, but the pattern is identical. This species has a ventral anal opening, like both P. picta and P. haegeli. Sister-taxa relationship to P. haegeli is supported by a bootstrap value of 100. Intraspecific variability is less than 0.35%.

Phyllidia ocellata Cuvier, 1804

Twenty-seven specimens of P. ocellata were collected from North Sulawesi. Phyllidia ocellata is considered a very variable group with many colour morphs (Brunckhorst, 1993; Gosliner et al., 2018; Pruvot-Fol, 1957; Yonow, 1996), and it differs from P. babai in having an orange ground colour while P. babai has a white dorsum. Phyllidia ocellata specimens have a yellow (or orange) background and large white tubercles, the lateral ones placed in the centre of black rings which are ocellated with a white line on both sides of the black ring (Fig. 5.5a–f). One of these black rings with a white tubercle is located anteriorly in front of the rhinophores, and there are two or three pairs present on each side of the midline. The foot and oral tentacles are grey, the latter with yellow tips. One specimen (Phma16Sa1, Fig. 5.5g) was first assigned to P. madangensis Brunckhorst, 1993 based on the overall black colouration with rounded white and yellow tubercles arising from a black background (Undap et al., 2019: fig. 5E); however, our phylogenetic study clearly places this specimen within the clade of P. ocellata composed of another 30 P. ocellata sequences, and which is supported by a bootstrap value of 100. All four species delimitation tests resulted in a single species referable to P. ocellata (Fig. S1), with an intraspecific variability of approximately 5% and up to 10% when including the long branch leading to PhocSa6 (Table S4).

Previous chemical investigations of P. ocellata reported isolation of bicyclic sesquiterpene isonitriles (Fusetani et al., 1992; Okino et al., 1996). More recently, bicyclic and tricyclic sesquiterpene isonitriles and isothiocyanate were isolated from Australian P. ocellata specimens (White et al., 2015). In the crude extract of a P. ocellata specimen analysed in this study (Phoc16Sa7, Fig. 5.5e), several minor peaks could be assigned to sesquiterpene isonitriles, which is in accordance with published reports. Polar brominated compounds were also detected, but these differ from those detected in P. cf. babai (see above). A derivative of the commonly occurring 5-bromotryptophane, 5-bromoabrine (m/z 297.000 and 298.996 [M + H]⁺) was detected only in P. ocellata.

Phyllidia sp. a

The sister taxon of P. ocellata is an undescribed species that we call Phyllidia sp. a (Phsp17Bu1, Fig. 6.1a). It is a unique specimen, collected from deeper areas in Bunaken Island at 25.3 m (for details see Eisenbarth et al., 2018). It differs from the other Phyllidia species by having only smooth and very low papillate ridges on the notum instead of tubercles. Three white bands lie on the central notum covered by these papillae, the median ridge starting behind the rhinophores along the midline. The mantle margin is yellow-gold. A smooth black band on the lateral notum forms an elongated ring. Gosliner et al. (2008: 300, top right photograph) depicts an individual identified as Phyllidiopsis monacha that looks like our specimen, and we assume that this was erroneously identified as Phyllidiopsis monacha due to its synonymisation with Phyllidiopsis dautzenbergi (Vayssière, 1912) by Brunckhorst (1993). Phyllidia monacha Yonow (1986) has distinct if tiny tubercles which are absent in our new Phyllidia sp. a, and the central black ring extends as black radiating lines towards the mantle margin in P. monacha. Phylogenetic analysis of the 16S data set revealed a sister-taxa relationship of this new species with P. ocellata, supported by a bootstrap value of 99. Unfortunately, we were not able to retrieve the CO1 sequence of Phyllidia sp. a; however, the species delimitation tests using 16S information clearly render Phyllidia sp. a as a distinct species (Fig. S1).

Phyllidia varicosa Lamarck, 1801

The most frequently collected species was Phyllidia varicosa with 114 specimens, and they were also the largest specimens (up to 87 mm) in our collections. Smaller animals of P. varicosa (Phel15Bu1, Fig. 6.2a: 7 mm) resemble the whiter P. elegans (Phel18Bu2, Fig. 4.1a) and the spikier P. haegeli (Phco5Bu1, Fig. 5.3a). Larger animals have a greater number of simple but acute ridges, and only the tubercles on the median ridges are capped with yellow. The oral tentacles are grey with a faint or bright yellow tip. The foot exhibits a black stripe, which can be broken up in black sections (PhvaBu2, Fig. 6.2g). Fahrner and Schrödl (2000) have shown that the black pigment on the foot can fade during preservation and thus might lead to misidentification. Adult living specimens of P. varicosa usually curl their lateral mantle rim underneath, and so they look rather elongate on the substrate (Fig. 6.2f, h, i). The 130 sequences from our collections combined with those from GenBank are supported as a monophyletic group by a bootstrap value of 100. However, this species, with the highest number of sequences included in our analysis, also has the second greatest intraspecific genetic variability of up to 16% mainly caused by long branches like Phva18LS3 (Table S4).

Extracts of four P. varicosa specimens collected in 2015 from two locations at Bunaken National Park were analysed with LCMS. The metabolomes varied for each locality (Fig. S7d), but the major compounds detected can be assigned to sesquiterpene isonitriles. The herein studied P. varicosa specimens seem to lack chemotaxonomic markers that would allow the clear chemical distinction of P. varicosa from other phyllidiid nudibranchs.

Phyllidia babai Brunckhorst, 1993

A single sequence is available from the studies of Stoffels et al. (2016), assigned to Phyllidia babai, that groups neither with the similarly coloured P. cf. babai nor with P. ocellata in our tree reconstructions. It actually represents the first offshoot within the clade Phyllidia. Stoffels et al. (2016: fig. 10d, RMNH.Moll.336464) illustrated a single specimen which was white and similar to one of their P. cf. babai specimens on the same plate (2016: fig. 10e, f), but has more and smaller black rings and patches, as well as two black rings on the midline. We did not collect any specimens that could be assigned to this clade, which may be the “true” Phyllidia babai.

Reticulidia Brunckhorst, 1990b

Specimens of the genus Reticulidia were not represented in our collection; however, we included sequences from two species (R. halgerda Brunckhorst & Burn in Brunckhorst, 1990b and R. fungia Brunckhorst & Gosliner in Brunckhorst, 1993) from Stoffels et al. (2016) in our analysis for completeness. This monophyletic genus is sister taxon to Phyllidia with a bootstrap value of 92 (Fig. 15) to 99 (Figs. 13, 14).

Phyllidiopsis Bergh, 1876

The genus Phyllidiopsis is characterised by having fused oral tentacles that form an oblong structure. Thirty-one species are currently recognised (WoRMS Editorial Board, 2021). In our collection, seven species could be distinguished, six of which can be assigned to recognised species (P. krempfi Pruvot-Fol, 1957, P. burni Brunckhorst, 1993, P. xishaensis (Lin, 1983), P. annae Brunckhorst, 1993, P. sphingis Brunckhorst, 1993, and P. shireenae Brunckhorst, 1990a). Affiliation of one clade, Phyllidiopsis sp. a, remains uncertain. All molecular analyses result in two distinct clades within the genus, both supported by highest bootstrap values. One clade comprises species that exhibit a tubercular appearance (P. krempfi, P. burni, P. sp. a), and the second clade comprises species with ridges or a smooth notum (P. shireenae, P. xishaensis, P. annae, P. sphingis). The sister-taxa relationship within this clade is also supported by highest bootstrap values. Lower bootstrap values (96) can only be seen with regards to the position of P. xishaensis, which is a sister taxon of the clade P. annae/P. sphingis (in the concatenated and 16S data sets) or forms a sister taxon to P. shireenae (in the CO1 data set). Interspecific variability was highest between P. krempfi and P. annae (Table S5). Chemical analyses revealed no genus-specific compounds. Interspecific variability was highest between P. krempfi and P. annae (Table S5).

Phyllidiopsis krempfi Pruvot-Fol, 1957

Our 16 specimens of P. krempfi (Fig. 7.1a–j) have a pink notum with irregular black lines that meander between the tubercles. Two main longitudinal black lines merge in front of the rhinophores and run backwards where they can merge behind the anal opening or may stay separate. The rhinophores are black on the tip and posterior part, but pink on the anterior side. The hyponotum and foot sole are grey; the oral tentacles are rectangular in shape, with grey to black along the lateral grooves. Phyllidiopsis krempfi is very similar in its external appearance to P. pipeki Brunckhorst, 1993, which was described as having two longitudinal black lines and large complex pink tubercles forming a crest on the midline. The foot sole, oral tentacles, and hyponotum are pale pink to grey. In general, it is very difficult to distinguish between these two species, and several of our specimens (Phskr18Ba1, Fig. 7.1a; Phspi18Ba1, Fig. 7.1c; Phsh15Bu1, Fig. 7.1e; Phfi16Sa1, Fig. 7.1f) match the original and subsequent descriptions of P. pipeki (see Brunckhorst, 1993; Yonow, 2011). However, these specimens clearly group with typical P. krempfi colour morphs in our molecular analyses (Fig. S1). Furthermore, some specimens of P. krempfi depicted by Stoffels et al. (2016) have oral tentacles without the black-grey colouration along the sides; these specimens of P. krempfi (RMNH.Moll.336453, 336512, and 336469) show a colouration that matches the original description of P. pipeki, but the tuberculate notum matches the description of P. krempfi. We also have four specimens with characters of P. pipeki on the ventral side, but characters of P. krempfi on the dorsal side, further supporting this synonymy. In our molecular investigations, all P. pipeki colour morphs grouped together with P. krempfi from Sulawesi (this study), Ternate, and Papua New Guinea (Stoffels et al., 2016), and Terengganu, Malaysia (Alqudah et al., 2015). Based on these results, P. pipeki is here considered a synonym of P. krempfi. Phyllidiopsis krempfi as a monophyletic group is supported by a bootstrap value of 100 and shows an intraspecific genetic variability of 2.43% (Table S5).

Three P. krempfi specimens were selected for chemical analysis. Specimens Phpu16Sa58 and Phpu16Sa66 (Fig. 7.1i, j) from Sangihe Island had identical chemical profiles with one major metabolite (see Fig. S8a). Fractionation and HPLC separation of the crude extract of Phpu16Sa66 led to the isolation of 7-isocyano-7,8-dihydrobisabolene (Fig. S4a) as the major metabolite. This natural product was first reported from a Ciocalypta sponge (Gulavita et al., 1986). The structure was unambiguously established by comparison of its ¹H and ¹³C NMR spectra with published data. Intriguingly, the MS data of the main compound in the crude extract of Phpu16Sa66 show a series of m/z values, i.e., 205.191, 232.202, and 259.215 [M + H]⁺. The m/z 259.215 can be reasonably explained with an undescribed corresponding sesquiterpene core bearing two isonitrile functions (Fig. S4c), and thus, the isolated compound could be regarded as an artefact.

Interestingly, the metabolomic profile of the third P. krempfi specimen (Phpu15Bu37, Fig. 7.1h) collected in Bunaken Island was significantly different from the other two analysed specimens Phpu16Sa58 (Fig. 7.1i) and Phpu16Sa66 (Fig. 7.1j) from Sangihe Island, which is approximately 200 km north of Bunaken Island. In addition to 7-isocyano-7,8-dihydrobisabolene, its extract contained several other major constituents that could be identified as sesquiterpene isonitriles based on diagnostic ions.

Phyllidiopsis burni Brunckhorst, 1993

Ten specimens are assigned to P. burni (Fig. 7.2a–d), characterised by having an ovate body with a narrow but bright pink band around the mantle margin. The anterior and posterior ends of each animal tend to be slightly pointed. Large pink compound tubercles are present on a black background. Small single pink tubercles lie around the mantle margin. The rhinophores are predominantly black. The black anal opening is generally on the top of a tubercle. The hyponotum is dark grey with pink patches. The foot sole is grey and brighter toward the edge. The colour of our specimens is very similar to the external colouration of Phyllidiella pustulosa; however, P. burni can be distinguished by its fused oral tentacles typical of the genus Phyllidiopsis. The rectangular oral tentacles are pale white to pink with dark grey grooves on the outer sides.

The extracts of three investigated P. burni specimens (Phpu16Sa19, Fig. 7.2b; Phpu16Sa55, Fig. 7.2d; Phpu16Sa93) revealed various isonitrile sesquiterpenes in LCMS analysis (see Fig. S8b). While chemotypes of the first two specimens are rather unspecific, the extract of the third specimen was dominated by a unique major sesquiterpene isonitrile with a retention time of 17.1 min that was not encountered in any other chemically investigated phyllidiid.

Phyllidiopsis sp. a

Sister taxon to Phyllidiopsis burni (supported by a bootstrap value of 100 in the 16S analysis) is a clade composed of six specimens which look very similar to P. krempfi. This clade, Phyllidiopsis sp. a (Fig. 7.3a–d), is characterised by black rhinophores with a pink to pale pink base, whereas P. krempfi (including the P. pipeki colour morphs) has rhinophores with black tips and posterior sides, but pink on the anterior faces. The boundary between the black rhinophore and the pink base runs diagonally in P. krempfi and not horizontally as in Phyllidiopsis sp. a (compare Fig. 7.1i–j vs. Fig. 7.3a–d). The rims of the rhinophores are angled and very low, with pustules on the posterior side. The narrow mantle edge is pink. The anal opening occurs on a white to pink tubercle. All six specimens of this clade were only collected in North Sulawesi (Lembeh Strait, Sangihe Island, and Bangka Island). Interestingly, molecular barcoding resulted only in 16S sequences and no CO1 sequence in all six specimens. Species delimitation tests based on the 16S data set clearly indicate this clade as a distinct species with a support value of 100 in the 16S tree (Fig. S1). No chemical analyses were performed on any specimens of this clade.

Phyllidiopsis clade xishaensis /annae / sphingis / shireenae

One monophyletic clade (bootstrap support value 100) of Phyllidiopsis consists of four species characterised by white colouration and longitudinal black lines. This clade comprises P. xishaensis (Lin, 1983) (Fig. 8.1a–d), P. annae Brunckhorst, 1993 (Fig. 8.2a–c), P. sphingis Brunckhorst, 1993 (Fig. 8.3a), and P. shireenae Brunckhorst, 1990a (Fig. 8.4a–c), and all species are supported by bootstrap values of 100. Whereas P. shireenae can reach a large size (more than 100 mm), the other three species are usually small, less than 30 mm (Gosliner et al., 2018). Besides the species-specific black pattern and ridge(s) on the dorsum, our seven specimens of P. xishaensis and one specimen of P. sphingis are characterised by white to yellow rhinophores. Phyllidiopsis shireenae has entirely pink rhinophores, and P. annae has black rhinophores which are pink at the base. Phyllidiopsis annae and P. sphingis can be distinguished by their blue mantle colouration, while P. xishaensis and P. shireenae have pale white mantle colouration. Valdés (2003) assigned one of his specimens (AF430368) from New Caledonia to P. sphingis, but in our study, it is grouped together with one P. annae sequence published by Hallas et al. (2017: MF958283) as well as with our seven P. annae specimens; the misidentification is understandable as they have similar dorsal patterns. Intraspecific variability of all four species is low, with a maximum of 1.89% in P. annae (Table S5). However, the number of included specimens was also low with a maximum of 11 specimens of P. xishaensis (intraspecific variability 1.52%).

Only one species of this Phyllidiopsis clade was investigated chemically. The crude extracts of two P. shireenae specimens contained sesquiterpenes typical for Phyllidiidae and did not provide any species-specific chemical cues.

Phyllidiella Bergh, 1869

According to our morphology and molecular-based results, we were able to identify six species of Phyllidiella: P. pustulosa (Cuvier, 1804), P. nigra (van Hasselt, 1824), P. zeylanica auctt. (not the original description by Kelaart, 1858, but used by authors in error, see below), P. rudmani Brunckhorst, 1993, P. albonigra (Quoy & Gaimard, 1832), and P. hageni Fahrner & Beck, 2000, and at least ten clades or subclades that could not be assigned to a named species (Figs. 9, 10, 11, 12). One unexpected result in our tree was the position of the single specimen of an unidentified Ceratophyllidia species from the Philippine Islands (Hallas et al., 2017). It is placed within Phyllidiella (bootstrap value 79), which would render the genus Phyllidiella paraphyletic if the specimen was correctly identified to genus (Figs. 13, 14, 15, S1). The preserved animal has pigment spots on the fleshy papillae consistent with being identifiable with Ceratophyllidia sp. 4 in Gosliner et al. (2015; T.M. Gosliner pers. comm.). All clades of Phyllidiella, including the un-named clades, were supported by high bootstrap values. A few species relationships are supported with highest bootstrap values, e.g., the sister-taxa relationship of P. pustulosa/P. nigra (100) (not supported by the CO1 data set), P. sp. c/P. albonigra (99), and P. hageni as sister of P. sp. d/P. sp. complex e (100).

Phyllidiella pustulosa (Cuvier, 1804)

Our results strengthen the need to formulate adequate diagnostic characters for P. pustulosa, the type species of the genus Phyllidiella, as described and illustrated by Cuvier (1804). Only then can we characterise the other clades to distinguish them from the nominal species and solve the P. pustulosa “species complex” problem (see Stoffels et al., 2016; Bogdanov et al., 2020). All eight synonymies by Brunckhorst (1993) must be reconsidered in light of these results. However, of the type material of P. pustulosa, only one syntype is available that was selected as lectotype by Brunckhorst (1993: 49). Pruvot-Fol (1957) re-described the material of Phyllidiella pustulosa deposited by Cuvier in the Muséum nationale d’Histoire naturelle in Paris, as did Brunckhorst (1993), and our re-analysis of this lectotype did not reveal any further information. The single animal is completely white with no traces of any black pigment (Fig. 2e). Brunckhorst (1993) revised the genus and its species based on his material as well as some museum specimens. He described P. pustulosa as having an elliptical to oval shape, with a black dorsal background, rounded tubercles of various numbers grouped into clusters of two or three, and a pink margin. This description refers to the same set of characters previously described and illustrated by Cuvier (1804) and matches the majority of our specimens in the first Phyllidiella clade in our tree (Figs. 13, S1), so we assign the nominal name P. pustulosa to this clade, comprising 74 specimens from our collection. Our animals all have a black background with a narrow white mantle margin (Fig. 9.1a–d). The tubercles are rather small and conical, and they usually cluster in two, three, or more in confined tubercular fields. Their colour can vary from pale pink to dark pink and even green. The rhinophores are black; the black anal opening lies in a black field. The hyponotum is grey, with pale pink patches, and the gills are dark grey-pink. The foot sole is pale grey to white, as are the oral tentacles, which have a dark grey line along the lateral grooves. Regarding the structure of the tubercles, most of our specimens also match the type materials of Phyllidia nobilis Bergh, 1869 synonymised with P. pustulosa. Examination of the three paralectotypes deposited in Copenhagen (NHMD-633656, Fig. 2b; NHMD-633657, Fig. 2c; NHMD-633658, Fig. 2d) revealed that the tubercles cluster in groups of three or more on a black background and the mantle has a white band along the rim. These specimens all appear to have median clusters with a clear separation from the lateral clusters. Examination of the lectotype of P. pustulosa deposited in Paris (MNHN, IM-2000–35147, Fig. 2e) also showed the same tubercular dorsal morphology as P. nobilis, but the black pigmentation is lost. Our material is here identified as the true P. pustulosa and groups with most other sequences retrieved from GenBank under the name P. pustulosa. The 94 concatenated sequences from our collection and GenBank cluster with a support value of 96 and are sister taxon to Phyllidiella cf. pustulosa with a bootstrap support value of 100 (Figs. 13, 14, S1). Intraspecific genetic variability is 9.45%.

Numerous chemical studies of P. pustulosa have been published (Hirota et al., 1998; Jomori et al., 2015; Kassühlke et al., 1991; Lyakhova et al., 2010; Manzo et al., 2004; Okino et al., 1996; Sim et al., 2020; White et al., 2017; Wright, 2003) as well as our recent publication (Bogdanov et al., 2020) highlighting the chemical variability that is found in nudibranchs commonly encountered and identified as P. pustulosa. Eight P. pustulosa specimens were utilised in this chemical analysis: LCMS revealed three main chemotype variants within the “true” P. pustulosa (see Fig. S9a–c). As expected from previous reports, all investigated specimens possess sesquiterpene isonitriles. The first chemotype is found in specimens Phpu15Bu14 and Phpu16Sa5 (from Bunaken National Park and Sangihe Island), both lying within the first branch of the clade. The major constituent 3-isocyanotheonellin was isolated from the crude extract of Phpu15Bu14 with its formamide derivative (Gulavita et al., 1986), and the structures were confirmed with ¹H and ¹³C NMR experiments. The second chemotype is characterised by having unidentified isonitriles in addition to 3-isocyanotheonellin and is shared by three specimens (Phpu16Sa3, Phpu16Sa4, Phpu16Sa6) collected in Sangihe Island and grouping in a subclade of P. pustulosa supported by a bootstrap value of 100. The third chemotype is found in specimens Phpu15Bu8 (Fig. 9.1a), Phpu15Bu9, and Phpu16Sa52 from different localities (Bunaken National Park and Sangihe Island) and also occurs in different subclades. Intriguingly, the observed chemotypes represent the subclades within P. pustulosa rather than the collection locality. Figure 18 illustrates the relationship between these clades within P. pustulosa in a network analysis. The clades are separated by a few shared mutations but also confirm the closer relationship of specimens sharing the same chemotype.

Phyllidiella cf. pustulosa

Twenty specimens group together in a separate clade (bootstrap value 98) that is sister taxon to P. pustulosa (Cuvier, 1804), and this relationship is supported by a bootstrap value of 100 (Figs. 13, 14, S1). Specimens of Phyllidiella cf. pustulosa (Fig. 9.2a–d) are characterised by a black mantle covered in short, conical, compound, pink tubercles which have a somewhat broader base than the tubercles in all other specimens assigned to P. pustulosa. This also renders the specimens paler because the black lines are narrower. These are arranged in longitudinal and diagonal lines, forming an X-shaped pattern in the first cluster on the medial dorsum. A narrow pink margin is present. The rhinophores are black with pale white–pink bases, their openings lying at the transition of a black line and a posterior pink tubercular field, and there are no rhinophoral sheaths. The anal opening lies in the connection of the two longitudinal black lines. Specimens of this clade have white or pale pink oral tentacles and foot sole. The colour pattern somewhat resembles the pattern of Phyllidiella lizae Brunckhorst, 1993, but analysis of the holotype and the paratypes of P. lizae (Australian Museum: C159493, C162784, C159492, Fig. 2g-i) allows the clear differentiation of our material from this species. Examination of the type material of P. lizae results in the following external characters: a white to pale pink background of the notum with a narrow cross-like black pattern, the margin is white, and the foot and oral tentacles are also pale pink. Phyllidiella lizae is also characterised by low pink rhinophoral sheaths that are absent in specimens of P. cf. pustulosa. It also appears that the tubercular field between the rhinophores is continuous in P. lizae but not in P. cf. pustulosa. One of the available P. lizae sequences (AF430365) retrieved from GenBank and identified by Valdés (2003) is clearly misidentified and clusters with P. rudmani (see below); unfortunately, there is no description or photograph of this specimen available (Á. Valdés pers. comm.). Our P. cf. pustulosa seems to match with one syntype and one paralectotype of P. nobilis (NHMD 91773, Fig. 2a; NHMD-633658a, Fig. 2d).

The species delimitation tests of the concatenated data set as well as of CO1 clearly distinguish between Phyllidiella pustulosa and Phyllidiella cf. pustulosa, the specimens of the latter united with a support value of 98. However, the ABGD test using only the 16S data set considers these two clades as one species (bootstrap value 100). The maximum distance between P. pustulosa and P. cf. pustulosa using the CO1 data set is 9%; this genetic distance value is lower in comparison to other interspecific genetic distances of Phyllidiella species (lowest minimum values between P. pustulosa and P. nigra is 9.84%, but in most cases higher than 10%, see Table S6), thus not providing adequate evidence for separate entities. Nevertheless, we consider P. cf. pustulosa as a separate subspecies because of the consistent external colour pattern, despite its broad distribution around Sulawesi. Two chemically analysed specimens (Phpu16Sa76, Phpu16Sa79) each represent a different branch of P. cf. pustulosa. LCMS analysis shows very different extract compositions in these two specimens (Fig. S9d), but both are dominated by sesquiterpene isonitriles and their derivatives, which nevertheless differ from those of true P. pustulosa.

Using the haplotype network analyses of CO1 and plotting chemical results on this network (Fig. 18) highlight discrepancies in the results between the various methods. In the network analysis, P. cf. pustulosa cannot be considered as an independent species from P. pustulosa because several specimens that form a subclade of P. pustulosa in the tree cluster with P. cf. pustulosa. This subclade is still unique as it shares at least 18 nucleotide states, thus distinguishing it from P. cf. pustulosa. However, it does not share unique nucleotide states only with P. pustulosa. Metabolomic analysis distinguishes P. cf. pustulosa from true P. pustulosa but also provides evidence that a particular subclade (Phpu15Bu8 and Phpu16Sa52) of P. cf. pustulosa is more similar to the P. pustulosa clade (Fig. 18). All these data indicate an ongoing process of speciation within Phyllidiella pustulosa and the lack of clear species boundaries at this moment in time.

Fig. 18

Haplotype network of Phyllidiella pustulosa (74 sequences, white circles) and P. cf. pustulosa (20 sequences, grey circles) based on the CO1 gene, with the four described chemotypes (coloured circles). Note the closer relationship of specimens with similar chemotypes, except chemotype 3 (purple), which is shared by distantly related specimens

Fig. 19

Haplotype network of Phyllidiella sp. a (29 sequences, white circles) and Phyllidiella zeylanica auctt. (24 sequences, grey circles) based on the CO1 gene. Note the closer relationship of specimens with similar chemotypes, coloured orange and blue, respectively

Fig. 20

Haplotype network from concatenated sequences of Phyllidiella sp. c with subclade 1 (32 sequences, white circles), subclade 2 (24 sequences, light grey circles), and Phyllidiella albonigra (9 sequences, dark grey circles), based on the concatenated data set. Note that different chemotypes (in different colours) are shared by more closely related specimens

Phyllidiella nigra (van Hasselt, 1824)

All 25 collected specimens of Phyllidiella nigra (Fig. 9.3a–f) exhibit a greater extent of black background compared to P. pustulosa and most clades. Additionally, the mantle margin is black in contrast to the white rim in all subclades of P. pustulosa and P. cf. pustulosa. The specimens have irregular single or coalesced tubercles distributed all over the mantle. The coalesced tubercles are smaller and usually conical, whereas the single tubercles are relatively larger and rounded. Their colour can vary between white to pink and green. The rhinophores are black, and the black anal opening exits from the black notum, sometimes in-between two small tubercles. The foot sole is grey to black. The oral tentacles are black; however, they can also be pale pink with dark grey-black tips or with dark lateral grooves. Originally, P. nigra (type locality Tahiti) was described as having tall pink-red rounded tubercles irregularly dispersed on the dorsal notum; the living specimen is illustrated in Bergh (1887: pl. 6, fig. 6). The oral tentacles were black and the foot sole was grey. We also collected animals that have the rather reticulate black pattern typical for P. pustulosa (e.g., Phpu15Bu7, Fig. 9.3a), but these specimens group together with the eight P. nigra sequences of Stoffels et al. (2016) (see Table S3) with highest bootstrap values (100) in the concatenated and CO1 data set.

LC-HRMS analysis of Phyllidiella nigra (Phpu15Bu7, Fig. 9.3a; Phu15Bu10, Phpu15Bu18, Phpu15Bu19) shows that all four specimens share the same major metabolites present as sesquiterpenoids (Fig. S9e). Relatively large extract amounts allowed fractionation and subsequent isolation and structure elucidation of five known sesquiterpenes (Fig. S4a), i.e., two bisabolene-type isonitriles (Iwashima et al., 2002; Manzo et al., 2004), 10-isocyano-5-cadinen-4-ol (Hirota et al., 1998), and an epimeric mixture of 9-thiocyanatopupukeananes (Yasman et al. 2003). However, the major metabolites were highly volatile or unstable and repeated isolation attempts failed (extraction of chemically identical large P. nigra specimens Phpu15Bu6 and Phpu15Bu43).

Phyllidiella sp. a

The clade Phyllidiella sp. a is composed of 29 specimens and contains 26 specimens from our collections (Fig. 9.4a–d) that are externally very similar to the sister clade P. zeylanica auctt. (not the original description by Kelaart, 1858 but used by authors in error), described in more detail below. Some specimens of Phyllidiella sp. a have tuberculate ridges, and some have small, rounded, compound tubercles on the notum. A white band is also present along the mantle margin, typical of most Phyllidiella species. The rhinophores are black, but the base is red in nearly all specimens (Fig. 9.4a–d). All specimens of Phyllidiella sp. a have two or three black X-shaped lines on the middle of the dorsal notum. The species is further characterised by a grey to black anal opening on a black part of the notum. The oral tentacles, the foot sole, and the hyponotum are grey, and the oral tentacles have dark grey tips or lateral grooves. The clade Phyllidiella sp. a is supported by a bootstrap value of 99–100. However, molecular analyses provide evidence of three subclades within Phyllidiella sp. a. Two subclades are composed of our specimens and include one further sequence retrieved from GenBank (Cheney et al., 2014; KJ001310, Lizard Island; assigned to P. pustulosa). Network analyses (Fig. 19) clearly show the distinctiveness of the two clades and thus confirm species delimitation when using the CO1 data set only.

Two further sequences retrieved from GenBank identified as Phyllidiella pustulosa by Wollscheid-Lengeling et al. (2001; AF249232, GBR, Australia) and Valdés (2003; AF430366, New Caledonia) form a sister clade to the two subclades mentioned above. The separation of this sister clade as a distinct species is not supported by our species delimitation test using the 16S data, and we therefore consider these two specimens as part of our Phyllidiella sp. a. We could not include these two specimens in our haplotype network analyses (Fig. 19) because CO1 data were lacking.

LCMS analysis of the crude extracts obtained from four specimens of Phyllidiella sp. a revealed a high degree of variation. We could not identify any unique features except for presence of several non-polar metabolites (retention time of more than 20 min; Fig. S9f), most probably sesqui- and diterpene isonitriles and sesquiterpene thiocyanates with m/z 264.18 [M + H]⁺. However, Phpu16Sa9 and Phpu16Sa90 showed a higher similarity to each other than to the other two specimens PhpuSa24 and PhpuSa26, thus reflecting the closer relationship of these two specimens, also according to molecular data (Figs. 13, S1) and network analysis (Fig. 19).

Phyllidiella zeylanica auctt.

Twenty-four specimens are grouped in a clade tentatively named Phyllidiella zeylanica auctt. These specimens (Fig. 10.1a–f) vary in their background colour from grey to pink. Many specimens have the tubercles arranged into longitudinal series; however, some specimens display single conical tubercles which can vary in size. A few have small, rounded, compound tubercles. The black lines can be arranged very differently, forming a reticulate pattern, but the pattern is basically a linear one surrounding larger central tubercular fields (Fig. 10.1a–f). The pink or white margin is followed by a black submarginal line. The rhinophores are black with white base. In contrast to that of many other species, the anal opening is large in all specimens, placed on grey to black dorsum, and surrounded by a ring of tubercles; it has a pink outer surface and black inner surface (Fig. 10.1e). By this unique character, it can be easily distinguished from Phyllidiella sp. a, which is characterised by a grey to black anal opening on a black area of the notum. The oral tentacles are white or pale pink with faded grey on the lateral grooves. Their outer appearance matches the description of Phyllidiella zeylanica of Brunckhorst (1993), but not the original description of P. zeylanica by Kelaart (1858: type locality Sri Lanka, figured by Eliot, 1906: pl. 52, fig. 10). The original description clearly states that there are two black rings and one central black stripe. The two rings are narrower than the central stripe with the outermost ring being the narrowest one. Burn (1970) subsequently described and illustrated specimens of P. zeylanica (sensu Kelaart) from Gulf of Kutch (India) with three black circular rings on the notum. Phyllidiella zeylanica (sensu Kelaart) was also described and illustrated by Yonow (2012 and references therein) based on specimens from the Maldives, Seychelles, Thailand, and Chagos with two or three circular black lines. Therefore, P. zeylanica was certainly misidentified in several publications, including Brunckhorst (1993), Domínguez et al. (2007), and Gosliner et al. (2015, 2018).

Unfortunately, several species were synonymised with P. zeylanica by Brunckhorst (1993) based on a pattern of black lines that coalesce in the posterior part of the body and pink compound tubercles which form longitudinal ridges and these will now need careful re-examination in light of current results: P. catena (Pruvot-Fol, 1956), P. seriata (Pruvot-Fol, 1957), and P. empelia (Yonow, 1984). Phyllidiella catena (type locality Mayotte, Mauritius), and P. empelia (type locality Sri Lanka) are species only recorded from the Indian Ocean, while the type locality of P. seriata is not known. The animals in our clade are not like Phyllidiella empelia, varying in details such as the numbers of black rings around the lateral mantle (which may be an ontogenetic variation), whether or not the black lines are united around the anal opening, the position of the black line surrounding the rhinophores, the presence of a central black line or band, and the anterior tubercles directly surrounding the rhinophores which may form an angular rectangle-like line instead of a curved line. They also appear to differ significantly from Pruvot-Fol’s species and P. rosans (Bergh, 1873). Thus, we cannot assign this clade to any of the above species without molecular evidence of these nominal Indian Ocean species. We urgently need more investigations involving these species to clarify the difficult taxonomy of this group. Our molecular analyses render Phyllidiella zeylanica auctt. as sister taxon to Phyllidiella sp. a with a very low support value of 51 (concatenated data set; Figs. 13, S1). This result reflects the contradicting results of the single-gene analyses. A joined network analysis of Phyllidiella sp. a and P. zeylanica auctt. (Fig. 19) also clearly shows the distinctiveness of these two species with 50 mutational steps separating them, despite their similar colourations.

Crude extracts of seven P. zeylanica auctt. specimens were analysed with LCMS. The chemotypes found showed a high degree of variation (see Fig. S9g) with dominant peaks attributable to sesquiterpene isonitriles by mass spectrometry. Careful inspection of the feature-based molecular network pointed out a few features that were detected in all investigated P. zeylanica auctt. specimens, e.g., ions with m/z 262.16 (M + H) and retention times of 14.6 min and 15.8 min match the accurate masses of bicyclic unsaturated isothiocyanate sesquiterpenes axisothiocyanate-4 and 10-isothiocyanato-4,6-amorphadiene (Fig. S4a) reported from various sponges (Adinolfi et al., 1977; Alvi et al., 1991). These features are very rare among the 70 chemically analysed phyllidiids and were detected in trace amounts in only five other specimens belonging to three genera, Phyllidia varicosa, Phyllidiopsis krempfi, Phyllidiopsis burni, and Phyllidiella sp. c subclade 2. Small specimen size and the resulting low extract amounts hindered us from isolation and unambiguous structure elucidation of the extract constituents.

Phyllidiella rudmani Brunckhorst, 1993

Two specimens in our collection could be confidently assigned to P. rudmani based on colouration of a pale white to bright pink dorsum with two characteristic narrow black lines on the lateral mantle (Fig. 10.2a, b). These lines do not unite in the anterior or posterior parts of the body in either specimen. The isolated tubercles in-between these two lines are tall and multi-compound, and a few small single rounded tubercles occur between these tall complex tubercles. The rhinophores are pink on the basal half and black on the upper clavus. A pink anal opening looks like a small tubercle lying between the two black lines. The hyponotum is white to bright pink, whereas the gill lamellae are black, only interrupted by the pink genital opening. The oral tentacles are entirely pink. Two sequences identified as P. rudmani and P. lizae from GenBank grouped with our specimens (bootstrap value 100; Fig. S1): one specimen from Raja Ampat, Indonesia (RMNH.Moll.336589), included in the analyses by Stoffels et al. (2016) looks like our specimens of P. rudmani. No information on the external features is available for the P. lizae specimen from Cheney et al. (2014) collected from Lizard Island, Australia.

The crude extract of one Phyllidiella rudmani specimen (Phpi15Bu2, Fig. 10.2a) was submitted to LCMS analysis. Several unique clusters in the molecular network (see Figs. 16, S2) indicated secondary metabolites that are found exclusively in this specimen. Careful inspection of the mass spectrometry data and marine natural product database search (MarinLit) allowed assignment of the major constituents of this complex extract (see Fig. S9h). Unlike the other chemically analysed Phyllidiidae, P. rudmani possesses kalihinol-type diterpenes (see Figs. 16, S4a, S9h) that are known from the sponges Acanthella cavernosa Dendy, 1922 and A. pulcherrima Ridley & Dendy, 1886 (Braekman et al., 1994; Chang et al., 1984; Hirota et al., 1996; Patra et al., 1984; Rodriguez et al., 1994; Wolf & Schmitz, 1998). Whilst none of our specimens was photographed on Acanthella, there are a number of photographs on iNaturalist illustrating P. rudmani on cf. Acanthella including two from Sulawesi: https://www.inaturalist.org/observations/59862922 (Bitung, North Sulawesi), https://www.inaturalist.org/photos/34319534 (Bangka, Sulawesi), https://www.inaturalist.org/observations/41101119 (Malé, Maldives), and https://www.inaturalist.org/observations/69467738 (Phi Phi islands, Thailand).

Phyllidiella sp. b

Eight specimens are characterised by a range of external appearances (Fig. 10.3a–d), and no single specific or distinguishing character permits an unproblematic identification. The notum of the animals possesses compound tubercles which are rather pointed. In general, the black lines form a reticulate pattern with vague X-shaped black lines in the middle of the tubercular groups on the central notum. The rhinophores are black with a white base, and low, pink-rimmed openings. Most specimens have a pink foot sole, but one specimen has a foot sole that is grey in the middle and pink on the sides (Phpu16Sa39, Fig. 10.3b). The oral tentacles are pink with dark grey apices and/or lateral grooves. Specimen Phpu16Sa25 (Fig. 10.3a) has dark pink oral tentacles. The clade is confirmed by a bootstrap value of 100.

Of the eight specimens in Phyllidiella sp. b, five were analysed with LCMS. Despite a high level of external diversity, the observed chemical profiles were very similar (see Fig. S9i), thus strengthening molecular results uniting these specimens. Many of the major peaks observed in UV-chromatograms could be attributed to sesquiterpene isonitriles and formamides. Additionally, a prominent peak with the retention time of 19.3 min that can be attributed to an oxygenated molecule with five degrees of unsaturation (molecular formula C₂₀H₃₂O calculated from m/z values 289.249 [M + H]⁺ and 271.239 [M-H₂O]⁺) was detected in all analysed Phyllidiella sp. b. This feature was also found in P. zeylanica auctt. and in one specimen of Phyllidiella sp. c subclade 2 (see below).

Phyllidiella sp. c

Phyllidiella sp. c consists of two different subclades (subclades 1 and 2) which are not considered separate species according to the molecular phylogenetic analyses or the ABGD test based on the CO1 and 16S data sets. Nevertheless, these two clades are considered as two different species when applying bPTP tests (both maximum likelihood and Bayesian analyses). Moreover, they are morphologically different and therefore described below as two separate entities.

Phyllidiella sp. c subclade 1 comprises 32 specimens. The black pattern of the specimens consists of two longitudinal lines and one outer ring (Fig. 10.4a–d). The obvious margin is pink. The inner median lines unite with the outer one anteriorly in front of the rhinophores and posteriorly and are connected to each other by several black X-shaped marks. There is no black band between the rhinophores, and the tubercular clusters in front of and behind the rhinophores are continuous. Several pink compound tubercles lie between these black lines forming short ridges arranged in clusters in the central part of the notum and in the lateral parts. The long pointed rhinophores are black with white to pink or red bases. The anal opening is on the black part of the notum. The hyponotum is grey with pink-white patches corresponding to the dorsal tubercles and ridges, the foot sole is white to grey and often with two distinct darker longitudinal lines (e.g., Phpu18Po4, Fig. 10.4d), and the oral tentacles are whitish pink with black along the lateral grooves. Some of these characters are similar to parts of the original description of P. granulata Brunckhorst, 1993 (type locality Guam), but future molecular analyses including specimens from Guam are needed.

All 24 specimens of Phyllidiella sp. c subclade 2 have one outer black ring and two lateral black lines (Fig. 11.1a–e), similar to subclade 1. The outer ring and the inner longitudinal lines unite anteriorly and posteriorly, similar to subclade 1. The inner lines also connect to each other transversely, forming two central tubercular areas, but these are H-shaped in subclade 2, and there is a black mark in the middle of these regions. The pink anal opening lies in a pink field, not in a black area as in subclade 1. There is a very obvious white–pink margin. In the central part of the notum, the compound pink tubercles are usually arranged in a circular form, with only four specimens exhibiting tubercles arranged in elongate, but short ridges interrupted by smooth black areas (Fig. 11.1e). The compound lateral tubercles are arranged in curved or straight series and separated by black lines. The rhinophores are predominantly black with a white base, but in Phpu18Ba12 (Fig. 11.1a), the rhinophores have a red base. The hyponotum is pink with grey patches caused by the dorsal colouration showing through, the gills are grey to black, the foot sole is white to pinkish grey, and the outer sides of the oral tentacles are black.

The specimens somewhat match the original description of P. melanocera Yonow, 1986 from the Red Sea, but our specimens have a distinct black submarginal line. There are also vague similarities with P. annulata (Gray, 1853; type locality French Polynesia), but this species was described with the absence of a white margin and an overall black colouration with a high number of pink rings (4–14) on the notum (Brunckhorst, 1993). Re-investigated type material of P. annulata (Fig. 3a-e) indicates similarities to the description of the holotype and paratype (only) of P. melanocera (Fig. 2f) in all other respects. In Phyllidiella sp. c, all specimens from both subclades 1 and 2 lack the black stripe between the rhinophores that is characteristic for most Phyllidiella species, and therefore, the pale pink tubercular field is continuous. The anal opening is not in a black area as usual, but lies between two compound tubercles in a pink to grey field.

Molecular evidence, including haplotype network analyses (Fig. 20), does not strongly separate P. sp. c subclade 2 from P. sp. c subclade 1. Both clades are united with a bootstrap support of 100. Species delimitation tests also reject subclade 2 as a separate species from subclade 1, except when applying bPTP algorithms. Within the network analysis based on CO1, only seven mutational steps are present between the two clades.

Four specimens of Phyllidiella sp. c subclade 1 (Phpu15Bu35, Fig. 10.4b; Phpu16Sa14, Phu16Sa37, Phpu16Sa92) and two members of subclade 2 (Phli16Sa5, Fig. 11.1c; Phpu15Bu27) were chemically investigated. Metabolomes found in both subclades exhibited a high degree of variation rendering P. sp. c “difficult” from a chemical point of view, thus reflecting the discrepancy when comparing morphological and molecular data (see Fig. S9j, k). The crude extracts contained, like many other phyllidiids, various sesquiterpenes with usual functionalities, i.e., isonitriles and formamides. We consider two major sesquiterpene isonitriles with retention times of 18.2 min and 18.4 min as a unifying chemical feature of subclade 1. Even though these metabolites are found in varying amounts in all four subclade 1 specimens, they could not be detected in the two specimens of subclade 2. Crude extracts of both specimens of P. sp. c subclade 2 had different chemical profiles (see Fig. S9k) which also differed from those of subclade 1 (Fig. S9j). Interestingly, the chemical profile of Phpu15Bu27 (subclade 2) was identical with those found in one specimen of Phyllidia coelestis (Phco15Bu6) and one specimen of Phyllidia elegans (Phel15Bu8) (see Fig. S10). All three specimens did not contain the usually observed sesquiterpenes, but several amphilectene-type diterpenes with an unknown major metabolite with an m/z 404.187 [M + H]⁺. In contrast, the extract of the other Phyllidiella sp. c subclade 2 specimen (Phli16Sa5) is composed of sesquiterpene isonitriles and corresponding formamides with a major compound eluting at 18.1 min in LCMS analysis. This compound is probably 7-isocyano-7,8-dihydrobisabolene, also isolated from a P. krempfi specimen (Phpu16Sa66). Both Phyllidiella sp. c subclades seem to lack obvious chemotaxonomic markers other than individual isonitrile compositions, and thus, they are difficult to clearly separate from other phyllidiids. This high degree of variability makes chemical characterisation of Phyllidiella sp. c on subclade level, as well as on species level, inconclusive and does not contribute to a taxonomic solution.

Phyllidiella albonigra Quoy & Gaimard, 1832

All nine specimens of Phyllidiella albonigra (Fig. 11.2a–f) are more elongated in shape than the other Phyllidiella clades, with rather pointed ends and a greater proportion of black background to the mantle than all other species except P. nigra, and a narrow pink, green, or white band along the mantle margin. The dense complex tubercles are small but compound, arranged in small clusters on the notum, and their colour varies from pink to green or grey. The rhinophores are black with a white base. The black anal opening arises on black background. This pattern most resembles the external appearance of the figured type material of P. albonigra (pl. 21, figs. 26, 27; type locality Tonga islands): the original description and drawing depict a specimen that is also rather black with irregular clusters of small tubercles. We therefore remove P. albonigra from synonymy with P. pustulosa (see Brunckhorst, 1993) and reinstate it here. Phyllidiella albonigra is clearly separate from other Phyllidiella species (bootstrap value of the clade is 100). Minimum genetic differences varied from 10.52 (from Phyllidiella sp. c, subclade 1) up to nearly 18% (from P. rudmani; Table S6). Its sister-taxa relationship with Phyllidiella sp. c is also supported in the tree reconstruction based on the concatenated data set and CO1 (bootstrap value 99), but is not resolved in the 16S analysis. Haplotype analysis of this species, together with Phyllidiella sp. c, provides further evidence of its distinctiveness (Fig. 20).

The chemical analysis of P. albonigra was described in detail in a recent publication (see Bogdanov et al., 2020 as Phyllidiella pustulosa clade 6). We were able to isolate the very rare dichloroimidic sesquiterpenes (Figs. 16, S4a) and their derivatives from one specimen (Phpu15Bu21, Fig. 11.2d) of this clade. These dichloroimids were identified as the major extract constituents in three out of four analysed specimens of Phyllidiella albonigra. Importantly, all chemically studied specimens were collected at different localities around North Sulawesi during the course of three years. Whereas the chemical profiles (UV chromatograms obtained during LCMS analysis) of two specimens (Phpu15Bu21, Fig. 11.2d; Phpu16Sa32, Fig. 11.2e) were identical, the chromatogram of a third specimen (Phpu17Ba4, Fig. 11.2f) was very different at first glance (see Fig. S9l). However, after obtaining pure compounds from the in-depth analysis, the different metabolomes could be explained by the presence of degradation products of the chemically unstable dichloroimids. Intriguingly, no dichloroimids or derivatives thereof could be detected in the crude extract of the fourth analysed specimen Phpu16Bu8 (Fig. 11.2b), which was smaller (30 mm) compared to the other three specimens (38–80 mm). Chemical analyses thus provide further evidence that P. albonigra is a separate species and can be distinguished from its closest relative Phyllidiella sp. c. (Fig. 20).

Phyllidiella sp. d

The ten specimens united in Phyllidiella sp. d (Fig. 11.3a–d) are externally quite different from each other. They have a notum covered by clusters of pink to grey low, conical tubercles surrounded by a network of narrow black lines. Sometimes, these lines form an X-shape on the central part of the notum. In one specimen (Phpu18Po3, Fig. 11.3c), the tubercles in the central part are arranged in a hexagonal pattern. There is a horizontal black line just anterior to the rhinophores. A white–pink mantle margin is present in all specimens. Each rhinophore arises from a compound pink to white or grey tubercle and is surrounded by a distinct sheath. The oral tentacles are white to pink with black tips. Two of our specimens (Phpu18Po3, Fig. 11.3c; Phpu18Bl5, Fig. 11.3d) resemble specimens of Phyllidiella nobilis illustrated from Thailand by Bergh (1904: pl. 2, fig. 15). This clade is supported by the highest bootstrap values in all analyses (100).

Phyllidiella species complex e

This clade comprises a complex of four different species, represented by one or two specimens each (Fig. 11.4a–d). We could not find any external characters that are unique to this clade, which most resembles Phyllidiella sp. d (Fig. 11.3a–d). Nevertheless, based on our molecular analyses, this clade is a distinct evolutionary lineage supported by a bootstrap value of 100 (Figs. 13, 14, 15), as are all four species, confirmed by the species delimitation tests using the 16S data set (Fig. S1). To assign any of these species to a nominal species is not possible due to the low number of specimens, and lack of intraspecific variability.

One member of this species complex (Phpu16Sa50) was analysed chemically (Fig. S9m), and its profile is dominated by sesquiterpene isonitriles. Its overall composition most closely resembles that found in Phyllidiella nigra (van Hasselt, 1824). However, thiocyanates, typical for many phyllidiids, were not detected in the extract of Phpu16Sa50. The extract also lacks distinctive chemotaxonomic markers, and since it is the only investigated specimen of the Phyllidiella species complex sp. e, further metabolomic studies are required to characterise this evolutionarily distinct clade.

Phyllidiella hageni Fahrner & Beck, 2000

Five specimens of Phyllidiella hageni (Fig. 12.1a–d) collected at several localities in Sulawesi demonstrate its diagnostic external appearance as described and illustrated by Fahrner and Beck (2000; three specimens from Lombok) and subsequently by Domínguez et al. (2007; two specimens from Papua New Guinea). However, the background is green in our material and not pink, and the oral tentacles are white without black tips; however, in all other colouration, our animals match the original and subsequent descriptions and illustrations. Two black longitudinal lines run the length of the mantle. Tubercles are rather flat and single, but sometimes clustered together. Black lines connecting the longitudinal lines meander between these groups of tubercles. A thin black line runs around the mantle margin. The rhinophores are black with bright green bases arising from low rhinophoral sheaths on green notum. The green anal opening is disguised between the green tubercles on green ground. The hyponotum, gills, foot sole, and oral tentacles are pale white (as described by Domínguez et al., 2007) to bright green. Some specimens with simple black lines, typical of P. rudmani Brunckhorst, 1993 can be distinguished from the latter by their tubercles: the tubercles of P. rudmani are tall and sparse, while those of P. hageni are low and densely distributed. One specimen from Raja Ampat (Stoffels et al., 2016: fig. 10i, KX235944) was assigned to Phyllidiopsis fissurata Brunckhorst, 1993, but clusters within our P. hageni sequences, and was probably misidentified. Despite the very distinctive fused oral tentacles in Phyllidiopsis, which are separate in Phyllidiella, these can contract to a smaller unit during fixation, thus leading to generic misidentification. Fahrner and Beck (2000), Domínguez et al. (2007), and Stoffels et al. (2016 as P. fissuratus, RMNH.Moll.336590) described P. hageni as having a pink dorsum. Pink notal tubercles for Phyllidiella hageni should no longer be a defining character due to variability of dorsal colour. The clade is supported by a bootstrap value of 100 (intraspecific variability up to 9%) and shows the highest interspecific distances between the other Phyllidiella species (17–32%).

Ceratophyllidia Eliot, 1903

Ceratophyllidia is only represented by one sequence obtained from GenBank (Hallas et al., 2017, CASIZ181247, Beatrice, Philippine Islands; identifiable with Ceratophyllidia sp. 4 in Gosliner et al., 2015; T.M. Gosliner pers. comm.) and groups within the genus Phyllidiella (Figs. 13, 14, 15, S1). The genus can be easily distinguished from Phyllidiella species by the shape of the swollen tubercles which are stalked (Brunckhorst, 1993). Because of the lack of data for this genus within our data set, we cannot assess this unexpected relationship of Ceratophyllidia with Phyllidiella in detail.

Nomenclatural changes

We synonymise Phyllidiopsis pipeki Brunckhorst, 1993 with Phyllidiopsis krempfi Pruvot-Fol, 1957 based on the results of both the molecular and morphological analyses. We are able to resurrect Phyllidiella albonigra (Quoy & Gaimard, 1832) that had been synonymised with P. pustulosa (Cuvier, 1804) by Brunckhorst (1993). We recognise the species often identified as Phyllidiella zeylanica as distinct and different from the true P. zeylanica Kelaart, 1858 and herein label it P. zeylanica auctt. We were also able to identify several problematic assignments of sequences in GenBank that are listed in Table S3 together with their NCBI accession numbers, sources, and collection localities.

Discussion

Phylogenetic results

Taxonomy and systematics of the Phyllidiidae has been considered a difficult issue due to the high number of very similarly coloured species, the lack of clearly distinct morphological characters, species variability, and/or the lack of detailed species descriptions and illustrations (Brunckhorst, 1989, 1993; Stoffels et al., 2016; Yonow, 1986, 1996). Broader comparative morphological studies with cladistic analyses have led to the synonymisation of the genera Fryeria Gray, 1853 and Reyfria Yonow, 1986 with Phyllidia, leaving the family with five accepted genera (Phyllidia, Phyllidiella, Phyllidiopsis, Ceratophyllidia, and Reticulidia). Few of the early studies on the Phyllidiidae included molecular data (Alqudah et al., 2015; Hallas et al., 2017; Valdés, 2003), and usually they only included a few sequences with a focus on deeper phylogenetic relationships, and it is not clear whether identifications were correct or if the sequence identities in NCBI are reliable, as shown in this work. Stoffels et al. (2016) were the first to provide deeper insights into intrageneric relationships and to provide information about some problematic groups, e.g., Phyllidiella pustulosa, which is now shown to be represented by numerous clades. The provision of images of their specimens made the interpretation of sequences assigned to their clades very valuable, and we therefore follow their example and provide images for many specimens, especially those with aberrant colouration. Stoffels et al. (2016) discussed the relationships of some phyllidiid genera, comparing their CO1 results with the few 16S results obtained from Valdés (2003) and Cheney et al. (2014) as well as on morphologically based cladograms by Brunckhorst (1993) and Valdés (2003).

Our trees are based on both mitochondrial genes and confirm the synonymisation of Fryeria with Phyllidia as per Valdés and Gosliner (1999), and the sister-taxon relationship of Phyllidia with Reticulidia, but our trees show a sister-taxon relationship between Phyllidiella and Phyllidiopsis without bootstrap support (57) in the concatenated data set in contrast to previous studies (Cheney et al., 2014; Stoffels et al., 2016; Valdés, 2003). Furthermore, the only sequences assigned to the type species of the genus Phyllidiopsis, P. cardinalis Bergh, 1876, does not group with all other Phyllidiopsis species, rendering the genus paraphyletic. In the morphological analysis of Valdés (2001b), P. cardinalis is sister taxon to one of the subclades identified in this study. If we interpret the identification of the sequences of P. cardinalis by Valdés (2003; AF430367, New Caledonia) and Cheney et al. (2014; KJ001308, Queensland, Australia) as an error, we do not need to reconsider all other Phyllidiopsis species as members of a different genus. However, if their identifications were correct (the species is in fact very distinctive), then all other Phyllidiopsis species will need to be assigned to a new genus. To clarify the relationships between genera, more data on Ceratophyllidia and Phyllidiopsis, especially correctly identified type species, need to be included and analysed using nuclear genes. Support values for some species groups are also often very low (lowest within Phyllidia with values between 45 and 60), and therefore, relationships within the genera are not reliable at present. Similar results were recently obtained by Layton et al. (2018, 2020) for Chromodoris within the Chromodorididae. These authors suggested that their study group experienced a recent radiation, thus lacking a clear diversification of the analysed gene. Recent radiations are also discussed for the species complex of the doridid Doris kerguelenensis (Bergh, 1884) in the Southern Polar Ocean (Wilson et al., 2009). However, in this case, the morphological plasticity of the species (see Wägele, 1990) could not be related to molecularly defined clades, and again, results on CO1 data clearly differed from 16S data (Wilson et al., 2009). These authors discuss the glaciation event as a driving force of diversification. Wilson and Burghardt (2015) analysed the recent radiation of the facelinid genus Pteraeolidia Bergh, 1875, which was considered monotypic for more than 100 years. Molecular analyses of specimens identified as Pteraeolidia ianthina (Angas, 1864) widely distributed in the Indo-Pacific clearly showed cryptic diversification. In this case, the symbiotic relationship with the dinoflagellate genus Symbiodinium clades and the probable switch to other species are considered as the driving force for rapid radiation. The situation in the Phyllidiidae also indicates a rapid radiation, as shown by the morphological plasticity, repeated colour patterns, and clear molecular species boundaries, despite a broad range of intraspecific variability in several clades. Driving forces could have been food specialisation on certain sponges with specific compounds and subsequent adaptation to different defensive molecules. However, we aim first to delimit species clades both molecularly and morphologically. Our second aim is to attempt to clarify species identities. Both of these aims reveal the diversity of phyllidiids in North Sulawesi, but we do not imply any relationships between the species at this point.

Species identification

For separating species, the barcode gap of a certain percentage of genetic difference between the species is considered a valuable boundary for identification in many organism groups (Hebert et al., 2003; Pentinsaari et al., 2016; Puillandre et al., 2012). In our study, applying the Kimura 2 parameter model for running the evolutionary distances, minimum intrageneric values between the species were similar in the genera Phyllidia, Phyllidiopsis, and Phyllidiella, with minimum values of 7.41% difference between Phyllidia sp. 9 and P. haegeli (Table S4), 8.65% between Phyllidiopsis krempfi and P. burni (Table S5), and 9.07% between Phyllidiella sp. b and Phyllidiella sp. c subclade 1 (Table S6). Maximum interspecific values were more than 30% with the highest value of 33.98% occurring between two species of Phyllidiella, P. rudmani and P. hageni. Intraspecific genetic variability was lower than 6% for most specie/clades, except Phyllidia picta that had a maximum distance value of 9.03%, P. ocellata (10.62%), P. varicosa (16,34%), Phyllidiella pustulosa (11.94%), Phyllidiella sp. a (8.02%), Phyllidiella zeylanica auctt. (6.04%), Phyllidiella albonigra (7.98%), Phyllidiella sp. e (11.84%), and Phyllidiella hageni (8.67%) (Tables S4, S5, S6).

By including similarities in nucleotide composition, which is taken into consideration in the haplotype network analyses, additional information is provided beyond the genetic distance. This can help in distinguishing species, especially when intraspecific genetic variability is broad and maximum values are close to or even overlapping with interspecific genetic distances. This is the case in, for example, the species Phyllidia picta, where intraspecific genetic variability was up to 9.03% in our study (Table S4). The haplotype network analysis (Fig. 17) confirmed the existence of two different clades that appear to be separated by 28 mutational steps. Problematic cases mainly occurred within the genus Phyllidiella. The P. pustulosa clade showed a maximum of 11.94% genetic difference within the clade and a minimum interspecific distance from P. cf. pustulosa of only 5.65% (Table S6). This value might still define a separate species, as indicated by the species delimitation tests when using only the CO1 gene; however, these tests split P. pustulosa into many more species. The haplotype network analyses clearly show that these subclades share nucleotide states and thus contradict the species concept based only on CO1. One group which clusters with P. cf. pustulosa (and therefore differs from our phylogenetic analysis) shows the typical chemistry observed in P. pustulosa. This specific situation might be caused by hybrid introgression, and analysis of nuclear genes is necessary for further insights into a probably ongoing speciation process.

For the closely related Phyllidiella sp. c subclades 1 and 2, the maximum intraspecific genetic variability was 1.28% and 3.33%, respectively; however, the minimum interspecific genetic distance between these two clades was only 2.92%. This value is far below the typical interspecific genetic distances observed between the other Phyllidiella species, which was a minimum of 9.07% between Phyllidiella sp. b and Phyllidiella sp. c subclade 1 (Table S6). This indicates that both subclades 1 and 2 genetically belong to the same species, despite the differences in their external appearances and metabolomes. Looking at the maximum values between the various species, the largest distance was observed between Phyllidiella rudmani and Phyllidiella hageni (34%) (Table S6).

Using molecular data in different methodological approaches (phylogenetic as well as haplotype network analyses), in combination with chemical compounds to characterise clades, we were able to identify variations within species (clades) and show that colour patterns, thought to be species specific, can occur repeatedly in a number of species, and that in many cases molecular investigation is necessary for correct assignment. Nevertheless, most species can be identified based on the combination of external morphological characters, with important distinguishing specific and generic morphological and colour features, such as the position and colour of the anal opening, the shape of the tubercles, the colour(s) of the rhinophores, the presence of a rhinophoral rim, or the form of the oral tentacles. Very often, the overall black colour pattern is also characteristic; however, there may be limitations in identification based only on colour. In contrast to the black patterns, a variation of the background colour (green, pink, white) might result from different food items (Brunckhorst, 1993), or variations in water quality; change in colouration during ontogeny is also possible and could lead to wrong species assignments. Valdés (2001b) discussed similar colour patterns across the genus Phyllidiopsis as a form of Müllerian, or Batesian mimicry, but indicated the lack of information on chemical compounds involved in the system. Similar intraspecific colour variation by mimicking other species has been shown in the family Chromodorididae: colour morphs within the same species, as well as in similarly coloured species, can occur sympatrically (Layton et al., 2018, 2020; Oskars & Malaquias, 2020; Padula et al., 2016), leading to the conclusion that identification may be limited when based only on colour patterns. Layton et al. (2018) discussed Müllerian mimicry in Chromodoris sea slugs with aposematic colouration as one explanation. However, the various colour patterns, especially in the genus Phyllidiella, which might be obvious for scientific identification, are probably not distinct enough for a natural predator. We assume that hybrid introgression might have occurred or is still occurring, resulting in an introgression of colour morphs into various clades and subclades, as recently shown for Chromodoris species when analysing genomes (Layton et al., 2020). This would explain a certain variability within a species or clade, as well as the repeated occurrence of particular colour patterns in other clades. Hybrid introgression is increasingly found in gastropods (e.g., Costa et al., 2020; Proćków et al., 2021), highlighting the complex nature of speciation, and genomic evidence was provided only recently for the first time for a nudibranch genus, Chromodoris (Layton et al., 2020). These questions can only be addressed by nuclear genomic analyses of phyllidiids in the future.

Novel species / clades

We were able to assign many of our specimens to formally described species, and we clearly show that several phyllidiid species are new to science and in need of formal descriptions or available names found. Phyllidia now comprises 11 described species, with the addition of four species recognised as unnamed (including P. cf. babai). Within the genus Phyllidiopsis, our data comprises six described species (excluding the atypical P. cardinalis), and one unnamed species. The species Phyllidiopsis krempfi and Phyllidiopsis pipeki seem to be indistinguishable on a molecular basis, and P. pipeki is here considered a junior synonym of P. krempfi: we have specimens across the full range of external morphologies, including a mixture of characters from both species in one and the same animal. Within Phyllidiella, 14 genetically distinct species can be separated based on slight differences in colouration, on molecular data, and in some instances on the presence of chemical compounds. Chang and Willan (2015) recognised at least nine clades within the P. pustulosa complex based on the mitochondrial genes CO1 and 16S; Stoffels et al. (2016) were also able to identify cryptic varieties of P. pustulosa. Our much larger data set including many more specimens clearly shows that many cryptic species and subspecies exist with colour patterns similar to the nominal P. pustulosa. Consequently, we no longer advocate the term Phyllidiella pustulosa species complex in order to avoid the hypothesis of a monophyletic clade. However, we emphasise here, again, that the interspecific relationship of the species and clades was not the goal of this study and probably cannot be solved with these two mitochondrial genes. We were able to assign members of our collection to five described Phyllidiella species, including P. cf. pustulosa within the nominal species P. pustulosa, and at least nine clades are recognised as currently unnamed species. This includes the species we call P. zeylanica auctt., as well as the species complex forming clade e comprising four species.

Chemical analyses

Chemical analyses of phyllidiids collected from Sulawesi using high-resolution LCMS, metabolomic tools (GNPS platform), and “classical” isolation and structure elucidation of secondary metabolites aimed to provide chemotaxonomic evidence of clade entities and/or species as revealed by molecular analysis combined with morphological characters. Since phyllidiids sequester their secondary metabolites from sponge prey, the extract composition is directly dependent on the respective food source. Specialisation on a certain prey with morphological adaptation is regarded as a major force that drives evolution in the Heterobranchia (Cimino & Ghiselin, 1999; Wägele, 2004). Precise food sources of many phyllidiids are unknown, and metabolomic analyses are an indirect method of providing information on this important ecological factor. A complex chemical profile of a phyllidiid extract can indicate either a broader array of food or a complex metabolome of a specific food source, as in several species discussed above. Extracts composed of either only few or rarely encountered compounds certainly indicate a restricted food array. We often do not know to what extent chemical composition of the same sponge species may vary, or if and how nudibranchs might chemically alter (= metabolise) acquired chemicals. Additionally, collected specimens were stored in ethanol under field conditions with exposure to environmental factors such as light and high temperatures, and unstable secondary metabolites might thus degrade during the first days after fixation. Taking in account these difficulties associated with metabolomic analyses, some investigated phyllidiids nevertheless were found to have highly specific extract composition, demonstrating the feasibility of a chemotaxonomic approach as an additional character in species delimitation.

The most commonly detected secondary metabolites across the three analysed phyllidiid genera are sponge-derived sesquiterpene isonitriles and related compounds, which is in accordance with previous chemical investigations of phyllidiid nudibranchs (Burreson et al., 1975; Fusetani et al., 1992; Hirota et al., 1998; Jomori et al., 2015; Kassühlke et al., 1991; Lyakhova et al., 2010; Manzo et al., 2004; Okino et al., 1996; Sim et al., 2020; White et al., 2017; Wright, 2003). Despite extensive chemical studies on marine heterobranchs (see recent reviews by Avila, 2020 and Avila & Angulo-Preckler, 2020), sesquiterpene isonitriles are known only from phyllidiids and two dorid nudibranchs, Cadlina luteomarginata MacFarland, 1966 (Burgoyne et al., 1993; Thompson et al., 1981) and Hexabranchus sanguineus (Rüppell & Leuckart, 1830) (Zhang et al., 2007); both species are considered generalists, feeding on many different sponges.

In addition to the typical array of sesquiterpenes throughout the Phyllidiidae, we detected specific compounds, or compound compositions, that allow chemical characterisation of some species or genera. A common feature that is detected across all chemically analysed Phyllidia species is the presence of somewhat polar brominated natural products. Whereas some of them are unique and were found only in P. cf. babai, we detected a range of brominated compounds that occur in several Phyllidia species. One of these, bromoindole alkaloid 5-bromotryptophan, was detected in all Phyllidia species; however, its derivative 5-bromoabrine was detected only in P. ocellata. Double-charged ions of unidentified polar mono- and dibrominated compounds were observed in all Phyllidia extracts except those of P. varicosa. Other than in Phyllidia, some of these brominated metabolites were detected in low amounts in the extracts of only six other specimens (all three Phyllidiopsis krempfi specimens and one specimen each of Phyllidiella sp. c subclade 1, P. zeylanica auctt., and P. rudmani).

Within the genus Phyllidia, extracts of P. coelestis and P. elegans showed very similar characteristic chemotypes with major metabolites being diterpene isonitriles, which were rarely found in other phyllidiids, thus uniting these two species and confirming their sister-taxa relationship. Unusual dichloroimidic sesquiterpenes described from Phyllidiella albonigra in our previous work (Bogdanov et al., 2020 as P. pustulosa clade 6) were detected in very low amounts in P. elegans, providing a chemical clue for its separation from P. coelestis. The brominated alkaloids spongiacidin A and B are here reported from a nudibranch, Phyllidia cf. babai, for the first time and differentiate this species from the very similarly coloured P. ocellata. Within the genus Phyllidiella, P. albonigra and P. rudmani were found to have the most characteristic metabolomes. Phyllidiella rudmani was chemically analysed for the first time and contained an array of kalihinol-type diterpenes not detected in any other phyllidiid in this study. Rare dichloroimidic sesquiterpenes were found to be a characteristic feature of P. albonigra and were detected in three of the four specimens collected from three different areas of our collections, Sangihe Island, Bunaken Island, and Bangka Island. However, these special compounds were absent in the smaller fourth animal from Bunaken Island.

We could not identify such obvious chemical features for the other investigated phyllidiids that would allow a clear separation/identification of the species based on chemical analysis only. Rather than possessing compounds of a specific chemical class, the extracts of many phyllidiids are characterised by varying compositions of sesquiterpene isonitriles and related compounds (thiocyanates, isothiocyanates, formamides). Data obtained with high-resolution LC-ESIMS coupled with a UV/Vis diode array detector (DAD) provided insights into the chemical extract composition, but was limited to UV absorbance, molecular mass/formula, fragmentation pattern, and retention times and could be used mainly for identification of known molecules. Currently, MarinLit database has 54 entries for sesquiterpene isonitriles with a molecular formula C₁₆H₂₅N, 33 entries for formamides (C₁₆H₂₇NO), and 60 for thio- and isothiocyanates (C₁₆H₂₅NS). In mass spectrometry, the major fragmentation is due to the loss of the functional group that results in a fragment with exactly the same mass for all above-listed compounds. Thus, laborious fractionation of crude extracts, isolation of pure compounds, and their structural elucidation with NMR are ideally required for metabolome characterisation of the majority of phyllidiid species. Extract amount is a limiting factor for such in-depth chemical analysis, and it could be performed on only three species, Phyllidiella pustulosa, P. nigra, and Phyllidiopsis krempfi. Isolated and unambiguously identified sesquiterpenes were different in each species, providing evidence that an individual assembly of structurally closely related compounds might be very characteristic for certain species.

Interestingly, there is a similarity of compounds within the subclades of P. pustulosa, or in the two subclades of Phyllidiella sp. a, or Phyllidiella sp. c, but with distinct differences between the respective subclades, indicating a specialisation within the subclades that is not associated with location. For example, the four P. sp. c subclade 1 specimens possess sesquiterpene isonitriles that were not found in the two specimens of subclade 2. Unfortunately, we do not know enough about food specificity nor about the metabolomes of the food organisms. Our data indicate that speciation of phyllidiids began as a generalist feeding on a broad array of sponge species, with subsequent specialisations on specific food organisms during the speciation processes, thereby restricting the number of compounds present and exhibiting the more specialised compound composition of the particular food.

We assessed geographical variations of metabolomes by the investigation of (usually) several specimens from different localities for each targeted species. In most clades, no geographic associations were detected, but there were a few exceptions. Phyllidiopsis krempfi was investigated including two specimens from Sangihe Island and one specimen from Bunaken Island: the two specimens from Sangihe, with nearly the same chemical profile, differed considerably from the third specimen collected from Bunaken Island. This “Sangihe” chemotype (see Fig. S11) was also found in four other specimens from Sangihe belonging to P. zeylanica auctt. (Phpu16Sa60, Phpu16Sa70), Phyllidiella sp. a (Phpu16Sa1), and Phyllidiella sp. c clade 2 (Phli16Sa5), all species with high metabolome variability. However, the number of analysed individuals was low, and therefore, the assumption of geographical differences is preliminary.

Difficulties encountered with museum material

Our re-investigation of museum material clearly revealed the problems of imprecise old descriptions, the loss of holotypes, the subsequent designation of syntypes and lectotypes, and synonymisation of available names in previous revisions. Our study demonstrates the necessity of molecular barcoding to verify species identities as well as careful re-examination of the literature and museum material. Unfortunately, our attempts to barcode more recently collected type material did not provide sequences that could be used to assign our unnamed clades to the respective types and thus to an available name. Colour, apart from black melanin, is not preserved in the type materials and often this melanin also fades. Thus, even identification of specific colour features on, e.g., the mantle rim, oral tentacles, or foot was not possible and hindered a putative assignment of our material to available type material. Further problems occurred in syntypes where the individual animals appear to belong to different species (e.g., Phyllidiella nobilis). Future analyses, and hopefully the development of methodologies to analyse old and formalin-preserved specimens, will solve some of these problems.