Specimens of Calmidae Iredale and O’Donoghue, 1923, Coryphellidae Bergh, 1889, Facelinidae Bergh, 1889, Fionidae Gray, 1857, Flabellinidae Bergh, 1889, and Trinchesiidae F. Nordsieck, 1972 from Mediterranean and eastern Atlantic localities were collected by SCUBA diving (Supplementary Table 1). Samples were preserved in 95% EtOH, and DNA was extracted using the salting out procedure (Aljanabi and Martinez, 1997). Amplifications were performed by PCR using universal primers: LCO1490 and HCO2198 (Folmer et al., 1994) for the COI, and 16Sar-L and 16Sbr-H (Palumbi et al., 2001) for the 16S fragment. The newly designed forward primer ITS2 MAT-03 [5′-CGUCGC(A/G)GACGCCUC(U/C)GCGC-3′] and the reverse primer ITS2 MOD Rev: 5′-AGTTCTTTTCCTCCGCTTA-3′ were used for the ITS2. PCR conditions were the same as reported in Furfaro et al. (2016b). Amplicons were sequenced by Macrogen Inc. (Netherlands).

ITS2, COI, and 16S DNA sequences were aligned together with GenBank sequences. We used the Muscle algorithm implemented in MEGA X (Kumar et al., 2018) for COI dataset, while for ITS2 and 16S, sequence-structure alignments were built using combined sequence-structure models in 4Sale (Seibel et al., 2008) and further optimised considering the consensus sequence of each rRNA domain. The software “WebLogo”² was used to generate a graphical representation of conserved domains and stems based on multiple sequence alignment (Crooks et al., 2004). Genetic distance (p-distance and Kimura 2-paramer, K2p, distance) and segregating sites were calculated using the programme Mega X for each marker (Kumar et al., 2018). Species delimitation analyses based on both genetic distances [Automatic Barcode Gap Discovery (ABGD) (Puillandre et al., 2012) and Species Identifier (Meier et al., 2006)] and on phylogenetic trees [Bayesian and maximum-likelihood analyses and Bayesian Poisson Tree Process (bPTP) (Zhang et al., 2013)] were carried out with the ITS2 dataset. Acanthodoris pilosa (Abildgaard in Müller, 1789) was used as the outgroup in BI and ML analyses. The GTR + I + G model was selected as the best substitution model by JModelTest 0.1 (Posada, 2008) according to the Bayesian information criterion (BIC). BI analyses were carried out with MrBayes 3.2.6 (Ronquist et al., 2012) running four Markov chains of five million generations each, sampled every 1000 generations. Consensus trees were calculated on trees sampled after a burnin of 25%. Maximum-likelihood analyses were performed in raxmlGUI 1.5b2 (Silvestro and Michalak, 2012), a graphical front-end for RAxML 8.2.1 (Stamatakis, 2014), with 100 independent ML searches and 1000 bootstrap replicates. Additionally, the recently developed Assemble Species by Automatic Partitioning (ASAP) analysis (Puillandre et al., 2021) was performed using default parameters.

Complete ITS2 rRNA and partial 16S rRNA (region from L7 to L13 stem-loops) secondary structures were obtained using the programme mfold (Zuker and Jacobson, 1998; Zuker, 2003). Best-supported folding models were predicted combining a thermodynamic approach (Mathews et al., 1999) with a close inspection of paired conserved regions. CBCs of rRNA sequence-structure were calculated and visualised in 4Sale (Seibel et al., 2008).

Genetic distance values for ITS2 rRNA and the mitochondrial COI and 16S rRNA markers are reported in Table 1. The lower genetic distance at ITS2 was obtained for Calma spp. with 5.3 and 5.5% p-distance and K2p, respectively, while the higher values were observed for Flabellina spp. with 36.6 and 50.3%. ITS2 is the more variable marker in most of the cases (Table 1). The lowest intrageneric p-distance and K2p values retrieved for the COI are 11.4 and 12.6% for Trinchesia spp. while the highest ones are 19.7 and 23.2% for Flabellina spp. (Table 1). Lowest values for the 16S marker were 3.0 and 3.1% between congeneric Calma spp., while the highest values were 7.1 and 7.6% for Flabellina spp. Intergeneric ITS2 p-distance and K2p values ranged from 16.8 to 19.1% between Coryphella and Fjordia species to 36 and 50.4% between Facelina and Dicata species. Genetic distance values of the two ITS2 conserved domains D1 and D2 are reported in Supplementary Table 2. Results from all the species delimitation analyses were congruent to each other (Figure 1) revealing the ability of the ITS2 to detect species even if closely related. Results from BI and ML phylogenetic inferences produced the same topology (statistical supports at each node are reported in Figure 1) and recovered two distinct monophyletic clades within Fjordia lineata (Lovén, 1846) that will need further taxonomic assessment.

The ITS2 rRNA alignment consisted of 44 sequences, from 14 species in 10 genera and six families (Supplementary Table 1 and Supplementary Figure 1), and 668 positions, among which 446 segregating sites. The ITS2 sequence length ranged from 275 nucleotides in Calma gobioophaga (Calado and Urgorri, 2002) to 391 nucleotides in Flabellina gabinierei (Vicente, 1975) with a mean length of 321 nucleotides. No differences in intra-specific ITS2 rRNA lengths were observed, except for the Atlantic and Mediterranean populations of F. lineata that show 311 and 313 nucleotides, respectively. The common derived ITS2 secondary structure of the eolid nudibranchs is organised in five domains (D1–5) (Figure 2 and Supplementary Figure 1) and conforms to the general folding observed in several other molluscan (Oliverio et al., 2002; Coleman, 2007; Salvi et al., 2010, 2014; Salvi and Mariottini, 2012, 2017).

The ITS2 D1 and D2 helix-loop regions are always identifiable in terms of sequence/structure and position (Figure 3). In fact, the basal double-strand RNA portion of D1 consisting of the very conserved triplet 5′-CGC/GCG-3′ is preceded by the triplet 5′-CGU-3′ and followed by the conserved single-strand sequence that separates D1 and D2 domains 5′-CUUC-3′. In D1 stem, another conserved base pairing is 5′-G/C-3′ positioned next to the triplet 5′-CGC/GCG-3′, present in all the families except for the Flabellinidae (substituted by the CBC 5′-A/U-3′) and the Facelinidae (substituted by the CBC C/G). The basal double-strand helix of D2 displays the consensus base pairing 5′-GG/CC-3′ (Calmidae, Coryphellidae, Flabellinidae, and Trinchesiidae), but in some families, this paired sequence is substituted by 5′-GA/UC-3′ (Fionidae) and 5′-GC/GC-3′ or 5′-GU/GC-3′ (Facelinidae) where CBCs and semi-CBCs maintain the initial stem base-pairing folding of D2. The ITS2 3′ region containing D3–D5 shows a moderate to high divergence structural folding between congeneric species (Figure 2 and Supplementary Figure 1) and the correct multiple sequence alignments were obtained only considering the derived 2D structures.

In all species examined, there is a very high conserved sequence motif located in the apical double helix-loop region of D5 and identified as the Apical STEM, described in other marine molluscs (Salvi et al., 2010; Salvi and Mariottini, 2012). The consensus Apical STEM of the eolid nudibranchs includes up to nine base pairs: there are five invariant nucleotide positions in the 5′-strand (nnnGCUCGn) out of nine and two invariant ones (nnGnGnnnnn) out of 10 in the 3′-strand (Supplementary Figure 2). The base-pairing between the two strands does not perfectly match the entire sequences of the double helix. Within the eolid families analysed, the Apical STEM consists of nine base pairs with a single U conserved in all families (Supplementary Figure 2). Caloria elegans (Facelinidae) shows a quite different sequence due to insertions of nucleotide stretches between the stem regions (Supplementary Figure 2). The Apical STEMs of the families Coryphellidae and Flabellinidae are very similar with semi-CBCs that preserve nine base pairs (Supplementary Figure 2).

The number of ITS2 CBCs observed between species of six eolid families is reported in Supplementary Table 3. Calma spp. does not show any diagnostic CBC, whereas the number of CBCs in Flabellinidae species ranges from one (Flabellina affinis/Flabellina ischitana) to six (F. gabinierei/F. ischitana).

The 16S rRNA dataset consisted of 41 sequences, obtained from the same species (Supplementary Table 1), with 391 positions, among which 152 segregating sites. The 16S rRNA sequence length ranged from 384 to 391 nucleotides. The 16S secondary structures (Supplementary Table 1 and Supplementary Figure 1) correspond to the domain D5 of the 3′ half portion of the gene (Horovitz and Meyer, 1995; Lydeard et al., 2000).

The sequence-structure of 16S (D5) is highly conserved between congeneric species, except for the L7, L10, and L11 stem-loops (Horovitz and Meyer, 1995), which contain diagnostic nucleotides and can be considered as barcoding regions (Furfaro et al., 2016a,b) (Supplementary Figure 1). In the D5 region, a few diagnostic nucleotides distinguish C. gobioophaga/Calma glaucoides when compared to Flabellina spp. and Trinchesia caerulea/Trinchesia morrowae.