Integrative taxonomy reveals a cryptic species of the nudibranch genus Polycera (Polyceridae) in European waters

Abstract This work aimed to test whether the colour variability featured by the European nudibranch Polycera quadrilineata is consistent with the concept of a single polychromatic species or may hide multiple lineages. Samples from across the geographic range of P. quadrilineata together with representatives from worldwide species with a focus on Atlantic diversity, were gathered and studied using an integrative taxonomic approach. Morpho-anatomical characters were investigated by light and scanning electron microscopy. Bayesian molecular phylogenetics using MrBayes, the Automatic Barcode Gap Discovery species delimitation method, and haplotype network analysis using the PopArt software were employed to help delimit species using the mitochondrial gene cytochrome c oxidase subunit I (COI). The results supported the existence of a second species, here described and named Polycera norvegica sp. nov., only known from Norway where it is sympatric with P. quadrilineata. The COI uncorrected p-genetic distance between the two species was estimated at 9.6–12.4%. Polycera norvegica sp. nov. differs by exhibiting a black dotted or patchy dotted pattern occasionally with more or less defined orange/brown patches, but never black continuous or dashed stripes as in P. quadrilineata. The two species share a common colouration with a whitish base and yellow/orange tubercles. Anatomically, P. norvegica sp. nov. has a weaker labial cuticle, a smaller radula with fewer rows, and only four marginal teeth, a reproductive system with a single lobed bursa copulatrix, shorter reproductive ducts, and a penis armed with two kinds of spines: needle-like and hook-shaped penile spines.


Introduction
The concept of 'cryptic species' has been debated and its definition remains a topic of controversy (Struck et al., 2018;Korshunova et al., 2019). Several authors questioned whether cryptic species truly exist or are transitional taxonomic misinterpretations resulting most of all from overlooked morphological characters (Horsáková et al., 2019;Korshunova et al., 2019). The concept is further complicated by the existence of multiple definitions across disciplines. For example, in behavioural ecology the concept applies to species with colourations and shapes that blend with the environment (Todd, 1981(Todd, , 1983Claridge et al., 2005;Bickford et al., 2007;Korshunova et al., 2019). Additionally, several derivatives of the concept have been proposed such as 'true cryptic species' (when a priori no morphological differences are recognized regardless of the distribution and ecology of species; Horsáková et al., 2019;Korshunova et al., 2019), 'pseudo-cryptic species' and 'quasi-cryptic species' (when subtle morphological difference can be recognized; Horsáková et al., 2019;Korshunova et al., 2019), 'semi-cryptic species' (morphological differences are very difficult to define; Vondrák et al., 2009;Korshunova et al., 2019), or 'false cryptic species' (morphological differences are obvious, but for some reason were missed or not highlighted in previous studies; Korshunova et al., 2019).
Polycera quadrilineata is a common species between the intertidal zone to depths of 30 m where it feeds on bryozoans of the genera Electra Lamouroux, 1816 and Membranipora d'Blainville, 1830 (Thompson & Brown, 1984;Thompson, 1988;Picton & Morrow, 1994), but specimens have been reported between 60-300 m (Bergan & Anthon, 1977). Around the British Isles the species reproduces and lay eggs during late spring and early summer (Miller, 1961;Bruce et al., 1963;Thompson & Brown, 1984), but in Norway we found aggregations of the species from November in the winter to June in late spring, although the exact time of these peaks of abundance can vary from year to year (Erling Svensen, personal observation). The egg-mass is white and crescent-shaped and the eggs are spherical, small and have been documented to range between 0.06-0.08 mm in diameter (Schmekel et al., 1982;Martínez-Pita et al., 2006;Martínez-Pita & García, 2017), which according to Martínez-Pita et al. (2006) is an egg size typically found in species with planktotrophic larval development.
Whether the highly polychromatic pattern of P. quadrilineata is consistent with a case of extreme intraspecific variability in nudibranchs or hides a possible complex of species was never tested before. In this work an integrative approach combining DNA and morpho-anatomical characters was used to investigate the taxonomic status of P. quadrilineata, i.e. the hypothesis whether it is one single biological lineage with extensive chromatic variability or alternatively, comprises a complex of multiple species.

Taxon sampling
Specimens were obtained by intertidal sampling, scuba-diving down to 30 m deep, and dredging using triangular, epibenthic and kelp dredges on board the research vessel 'Hans Brattström' owned by the University of Bergen, Norway (UiB). Most collected specimens were deposited in the scientific collections of the Department of Natural History, University Museum of Bergen, UiB (ZMBN).
Specimens were photographed alive with a digital SLR camera equipped with macro-lens, measured with a ruler (mm) for their maximum length (H), and frozen overnight inside plastic jars with seawater to ensure the body was kept fully extended for later possible anatomical studies. Afterwards, specimens were defrosted and fixed in absolute ethanol (>96%). Additionally, 47 sequences of the cytochrome c oxidase subunit I (COI) of Polycera quadrilineata and other Polyceridae taxa, representing a total of 17 species, together with the outgroup species Jorunna tomentosa (Cuvier, 1804) were obtained from GenBank and BOLD databases ( Table 1). The nudibranch species J. tomentosa (family Discodorididae) was chosen for outgroup taxa in the phylogenetic analyses based on the results of Palomar et al. (2014).

DNA extraction, amplification, purification and sequencing
Tissue samples for DNA extraction were gathered from 69 specimens of Polycera quadrilineata sensu lato by cutting a small part of their foot or mantle using forceps or a scalpel and were kept in 1.5 ml Eppendorf tubes filled with absolute ethanol. In rare cases, when specimens were too small to cut off enough tissue the whole specimen was used. To prevent contamination forceps and scalpel blades were rinsed with ethanol between each sample.

734
Cecilie Gotaas Sørensen et al. DNA was extracted from tissue samples using the 'Qiagen DNeasy Blood and Tissue Kit' (QIAGEN, catalogue no. 69506), following the protocol for 'Purification of Total DNA from Animal Tissues (Spin-Column)'. Amplification of the gene COI was performed through polymerase chain reaction (PCR) using the universal primers by Folmer et al. (1994) and a total reaction volume of 50 μl using 17.5 μl Sigma water (ddH 2 O), 5 μl buffer, 5 μl dNTP, 10 μl Q-solution, 7 μl MgCl, 2 μl of each primer (10 μM), 0.5 μl TAQ and 1 μl DNA. Some amplifications were carried out with only 25 μl volume using the same cocktail mix, but replacing the standard buffer with CoralLoad (CL) buffer from Qiagen, using only half of each quantity. PCR thermal cycles included an initial denaturation at 95°C for 3 min, followed by 39 cycles of 45 s at 94°C (denaturation), 45 s at 45°C (annealing), 2 min at 72°C (extension), and a final extension step at 72°C for 10 min, before cooling down. In order to rule out contamination, a negative and positive control were added to each PCR run. The negative control consisted of distilled water (ddH 2 O), whereas the positive control used DNA extract from a previously successfully tested sea slug species, namely Aplysia punctata (Cuvier, 1803).
The quality of PCR products was assessed using gel electrophoresis, by running 4 ml of PCR mixed with 1 μl Ficoll 5× loading buffer on a 1.0% agarose gel prepared with halfstrength TAE 1× buffer (tris-acetate-EDTA) containing the staining agent GelRed covered in TAE 1× buffer. For the PCR products already containing a loading buffer (i.e. the CL buffer), 5 μl PCR product were added directly into the gel. To quantify and estimate the length of amplified DNA fragments, a 5 μl FastRuller ladder marker was used. The gel was run for 30 min at 80 V and was posteriorly analysed under a Syngene UV-radiation machine. GeneSnap and GeneTools softwares (Syngene, Cambridge, UK) were used for imaging and manual DNA quantification, respectively.
Successful PCR products were purified using EXO-SAP, a combination of the enzymes Exonuclease I (EXO I) and Shrimp Alkaline Phosphatase (SAP), with each purification sample containing 8 μl PCR product and 2 μl EXOSAP (0.1 μl EXO, 1.0 μl SAP and 0.9 μl ddH 2 O). Reactions were run for 30 min at 37°C (incubation), followed by 15 min at 85°C (enzyme inactivation), and 4°C for cooling in the thermal cycler. Sequencing reactions were prepared in 1.5 ml Eppendorf tubes (kept on ice) using 1 μl of purified PCR product mixed with 6 μl of ddH 2 O, 1 μl primer (3.2 μM), 1 μl BigDye (BD) and 1 μl of sequence buffer. This process was repeated independently for each of the two primers Specimens (S) from the same lot were coded sequentially with the acronym S1, S2, S3, etc., in the column 'Sample no'. Novel sequences are marked with an asterisk.
(forward and reverse). The reactions were conducted in a thermal cycler for 5 min at 96°C (initial denaturation), followed by 25 cycles of 10 s at 96°C (denaturation), 5 s at 50°C (annealing), and 4 min at 60°C, before cooling down at 6°C. Afterwards, 10 μl ddH 2 O were added to each sequencing reaction and PCRs were delivered for Sanger sequencing using a capillary-based Applied Biosystem 3730XL DNA Analyzer at the DNA sequencing facility of the Department of Biological Sciences, UiB.

Phylogenetic and species delimitation analyses
The programme Geneious v. 11.0.3 was used to inspect, assemble and edit the chromatograms of the forward and reverse DNA strands. Sequences of each sample were checked by careful examination of the chromatograms and trimmed at both ends to remove parts of low quality and were translated to protein sequences using the invertebrate mitochondrial genetic code to ensure that there were no stop-codons present. Possible contamination of sequences was accessed by comparing with available sequences of molluscs through the BLAST toll in GenBank. Sequences were aligned using MUSCLE (Edgar, 2004) implemented in Geneious. The alignment was trimmed at both ends to a position where at least 50% of the sequences had nucleotide data yielding a final alignment with 642 base pairs (bp).
The MEGA X software (Kumar et al., 2018) was used to estimate uncorrected pairwise ( p) distances within and between species ( Table 2). The jModelTest2 v. 2.1.10 (Darriba et al., 2012) was used to find the best-fit evolutionary model [model selected = GTR + I + G] under the Akaike information criterion (AIC). The Bayesian analysis was performed using MrBayes (Huelsenbeck & Ronquist, 2001), with three parallel runs of five million generations each, sampling every 1000 generations, with a burn-in set to 25%. MrBayes was run through the portal CIPRES (Miller et al., 2010) and the consensus phylogram was converted into a graphical tree in FigTree v.1.4.3 (Rambaut, 2016).
The Automatic Barcode Gap Discovery (ABGD) species delimitation analysis (Puillandre et al., 2012) was used to help delimiting species using the three evolutionary models available (Simple Distance, Kimura (K80) TS/TV = 2 and Jukes-Cantor (JC69)). Each analysis was run independently using default settings.

Haplotype network analysis
Haplotype network analysis was conducted separately for the two lineages recognized within Polycera quadrilineata sensu lato (P. quadrilineata proper [63 seq.] and Polycera norvegica sp. nov. [17 seq.]), using the programme PopArt v. 1.7 (Population Analysis with Reticulate Trees; Leigh & Bryant, 2015). Prior to PopArt the alignment file combining the two lineages was edited using the text editor programme Notepad++ v. npp.7.6.6, and all empty positions at both ends were removed, yielding a final alignment with 594 bp. Sequence P48 (see Table 1) was excluded due to its reduced size (541 bp). Notepad++ was used to generate species-specific alignments and trait files with geographic area codes (i.e. 0 = specimen absent, 1 = specimen present). Alignments were converted into phylip format (phy) using the program Mesquite v.3.51 (Maddison & Maddison, 2018). Alignments and trait files were finally run in PopArt to create a standard TCS Network analysis (Clement et al., 2000) in order to visualize the genetic relationships and distances between the individual genotypes. The single specimen of P. quadrilineata obtained from Sveio (Norway) was for the sake of geographic proximity considered part of the geographic area Haugesund (Norway). The TCS haplotype networks were edited for more satisfying visualization using both PopArt v. 1.7, Adobe Illustrator, CS6 v.16.0.4 and Gravit Designer v.2019-2.1 at https://gravit.io/.

Morpho-anatomical work and scanning electron microscopy
Dissections were done under a stereo microscope Nikon SMZ-1500 equipped with a camera lucida. The animals were opened by dorsal incision, and the reproductive system and buccal mass, with the radula and labial cuticle were removed. The buccal mass was dissolved in a 10% sodium hydroxide solution until the labial cuticle and radula were cleansed from their surrounding tissue. These structures were then rinsed with water, and examined and photographed under a light microscope using the Life Science Imaging software cellSense v.1.18. The reproductive systems were drawn using the camera lucida, and each penis was isolated, opened, examined and photographed using light microscopy. The labial cuticles and penises were critical point dried using hexamethyldisilane. All parts (radula, penises and labial cuticles) were mounted on metallic stubs for scanning electron microscopy (SEM), and sputter coated with gold-palladium. Observations were done with a Hitachi S-3000N SEM-machine.

Phylogenetic analysis and species delimitation analysis
The molecular phylogenetic analysis included a total of 113 COI sequences containing 66 novel ones and 47 downloaded from GenBank and BOLD databases (Table 1). Eighteen lineages (including the outgroup taxon) were recognized and samples Green box refers to P. quadrilineata and purple box to Polycera norvegica sp. nov. Images refer to main morphotypes of these species and specimen depicted across both boxes represents the shared morphotype. preliminarily identified as Polycera quadrilineata split in two different clades with maximum support. One containing specimens from all over Europe including the type locality (Drøbak, Norway) and colour patterns consistent with the original description (O. F. Müller, 1776Müller, , 1779, and a second clade only with specimens from Norway that was rendered sister (PP = 0.95) to a clade containing individuals from an unidentified Polycera species from South Africa and specimens of Polycera faeroensis from Europe ( Figure 2). The genus Polycera as a whole was not rendered monophyletic due to lack of support (PP = 0.55) and inclusion of samples of the genera Polycerella A. E. Verrill, 1880 and Thecacera J. Fleming, 1828. Nevertheless, a clade with maximum support containing six species of Polycera was retrieved, namely with the type species of the genus P. quadrilineata, plus Polycera aurantiomarginata, Polycera capensis, P. faeroensis, Polycera sp.A (from South Africa) and Polycera norvegica sp. nov. (Figure 2).
COI uncorrected pairwise ( p) genetic distances showed a 9.6-12.4% difference between P. quadrilineata and P. norvegica sp. nov., and a range of intraspecific variability in the former species of 0-2.4% and 0.2-2.3% in the latter ( Table 2). The interspecific genetic distance between all included Polycera species was estimated at a maximum of 18.6-19.7% between P. faeroensis from Portugal and Polycera atra from California, USA and a minimum of 4.3-5.8% between P. capensis from South Africa and Australia and P. aurantiomarginata from Spain and Morocco. Intraspecific COI uncorrected p-distance between all studied Polycera species ranged between 0-2.6%.
The ABGD analyses under the three independent models of evolution retrieved the same number of partitions (=18), which were consistent with the 18 lineages suggested by the topology of the COI phylogenetic tree. Only for P values above 0.012915, the suggested number of partitions was lower, namely 17 groups (Kimura 80) and four with the Simple Distance and Juke-Cantor models (see supplementary material).
Class GASTROPODA   External morphology (Figures 3 & 4): Based on studied specimens H = 5-30 mm (but specimens with 45 mm have been reported; Moen & Svensen, 2014). Body surface smooth, covered with scattered tuberculate blotches; number and size of tubercles variable; tubercles either smoothly or sharply edged. Individuals with smaller tubercles often have larger quantity. Head with four to five, rarely six, smooth, digitiform veil processes projecting anteriorly; frontal veil processes yellow or orange in colour; apical tip often whitish. Non-retractile, lamellated rhinophores; stems thick, slightly leaning forwards; lamellated section slightly leaning
backwards; average of 10-12 lamellae present, but number can range from 6-15. Eyespots present behind rhinophores; inconspicuous in some specimens. Gill circlet with 7-9, sometimes 11 pinnate gills. Gills white with yellow or orange apical edges, sometimes black pigmentation present on upper part. One elongate cylindrical papilla presents on each side of gill circlet, projecting backwards, sometimes short and stubby with rounded apical tip, otherwise slender with pointed tip; distal part yellowish, proximal half to 2/3 of length white. Mid-dorsal yellow or orange line extends from behind the gills to tip of tail. Colouration (Figures 3 & 4): Two main colour morphs present. Yellow/orange colour morph ( Figures 3A-F, 4C, D)ground colour white, translucent, with yellow or orange circular to oval scattered tubercles. Rhinophore stems white with upper half of lamellate part yellow or orange; or rhinophore stems with black, dark brownish or grey pigmentation covering entire stem or restricted to frontal half-part, expanding into upper yellowish lamellate part. Gills with yellow or orange apical edges.
Striped colour morph ( Figures 3G-M, 4A, B)ground colour white, translucent, with small yellow or orange often circular scattered tubercles. Black or dark grey continuous or dashed longitudinal lines present. Thickness and number of dark lines varies between individuals; some almost fully covered, giving specimens a melanistic appearance, while others have fewer lines, with a nearly white background. Rhinophore stems often black, dark or light greyish; lamellae often yellow or orange, sometimes with additional black pigmentation in upper part. Gills, in some darker specimens, can depict dark pigmentation on tips and along edges.
Labial cuticle ( Figure 5C-F): Large, robust with two large and elongated lateral wings. Well-developed brownish centre with jaw elements.
Prostate gland massive, narrowing towards distal vas deferens, closely attached to bursa copulatrix. Inside vas deferens a cupshaped structure indicates end of the prostatic section. Vas deferens long, narrow, folded before reaching large penile bulb. Penis armed with numerous elongated, pointed, chitinous spines, of similar size; some spines bifid at apical tip. Vaginal duct long, folded, similar in width to vas deferens, connected to bursa copulatrix. Bursa copulatrix large, elongate, with two different parts; larger elongate proximal part, distal part more oval. Base of bursa copulatrix connected to pyriform, small receptaculum seminis by long, thin duct. Short uterine duct emerging close to receptaculum seminis and entering female gland, behind elongated portion of bursa copulatrix.
Ecology: Intertidal and subtidal species commonly found in shallow waters associated with algae. Feeds on encrusted bryozoans of the genera Electra and Membranipora. Reported from depths up 300 m (Bergan & Anthon, 1977). Lives in cold to temperate waters between 5-25°C (Picton & Morrow, 1994; Betti et al., 2017).
widest at mid-dorsal length close to anus, gills and papillae, ending in elongated and pointy tail. Body surface smooth, covered with scattered tubercles; number and size of tubercles variable, flattened or wart-like. Individuals with dotted morphotype tend to have rounded, slightly elevated, light-yellow tubercles. Head equipped with four to seven digitiform veil processes projecting anteriorly, yellow or orange coloured with whitish tips. Dotted morphotype often with only yellow spots scattered over a whitishtranslucent background. Non-retractile, lamellated rhinophores ending in cylindrical knob; stems slightly leaning forward; 6-10 lamellae present. Eyespots present behind rhinophores; nearly indistinguishable in some specimens. Gill circlet with 7-9, sometimes 11 pinnate gills; individuals may possess both smaller and larger gills; gills can partially retract into pocket. One elongate, narrow, papillae present on each side of gill circlet, projecting backwards or laterally; papillae often shorter and stubby with rounded apical tip on individuals with dotted morphotype; often slender with sharper apical tips on individuals with yellow/orange morphotype. Colouration (Figures 7 & 8): Two main colour morphs present. Yellow/orange colour morph ( Figures 7A-D, 8B, C)as in yellow/orange morphotype described for P. quadrilineata, but in the studied specimens the rhinophore stems were always whitish lacking black pigmentation. Faint greyish tiny spots may be present on head region and gills.
Dotted colour morph ( Figures 7E-H, 8A)ground colour white, translucent, with black, dark brownish or grey dots scattered over entire surface. Brown/orange patches present on head region or sometimes randomly along body surface. Rhinophores white translucent; lamellae sometimes with a weak hint of lightyellow pigmentation; stem sometimes dark brown or grey Frontal veil processes translucent white, with few scattered yellow patches randomly distributed. Posterior papillae white with yellow tips; sometimes white basal part dotted with yellow or black patches. Gills whitish with yellow, orange, or black spots along rachis and tips.
Labial cuticle ( Figure 9C-F): Small, weak, with two short lateral wings. Weakly developed brownish centre with jaw elements.
Reproductive system (Figure 10): Triaulic; hermaphroditic duct elongated, thin. Ampulla small, kidney-shaped; post-ampullary duct bifurcating into short oviduct leading to large female gland mass and short vas deferens through prostrate portion. Prostate Fig. 11. TCS haplotype network analysis based on the COI gene generated in the programme PopArt, including sequences from 63 specimens of P. quadrilineata. Lines between black dots represent one mutation, while black dots represent hypothetical haplotypes. Each coloured circle represents a unique haplotype, and the size of each circle indicates how many specimens share that haplotype. Different colours represent geographic locations. Fig. 12. TCS haplotype network analysis based on the COI gene generated in the program PopArt, including sequences from 17 specimens of Polycera norvegica sp. nov. Lines between black dots represent one mutation, while black dots represent hypothetical haplotypes. Each coloured circle represents a unique haplotype, and the size of each circle indicates how many specimens share that haplotype. Different colours represent geographic locations.
gland large, massive, narrowing towards distal vas deferens, surrounding bursa copulatrix. Inside vas deferens a cup-shaped structure indicates end of prostatic section. Vas deferens short, narrow, folded before reaching genital pore. Penile bulb reduced. Penis armed with two types of chitinous spines; spines closest to prostate more elongate; spines closest to genital opening hookshaped. Vaginal duct elongated, slightly bent, shorter than vas deferens, but of similar width. Vagina ends in large oval bursa copulatrix. Bursa copulatrix rounded; base of bursa copulatrix connected to pyriform, small receptaculum seminis by short, thin duct; short uterine duct emerging close to receptaculum seminis, entering female gland. Ecology: Sublittoral species occurring between 2-15 m depth, often on kelp (Laminaria spp.).
Distribution: Only known from Norway where it has been reported between Trondheim in the central west coast, along the south-west in Bergen, Haugesund, Stavanger and Egersund, and in Larvik in the southern coast.

Haplotype network analysis
The TCS haplotype network of P. quadrilineata with 63 specimens from Scandinavia, Mediterranean Sea and the Azores revealed the presence of 38 distinct haplotypes of which only nine were shared between individuals (Figure 11). From these, four haplotypes were each shared by two individuals (Norway: Bergen, Trondheim and Egersund), two were shared by three individuals (Sweden, Norway: Bergen, Egersund and Drøbak), one by four individuals (Norway: Egersund and Kristiansund), one by five (Sweden, Norway: Bergen, Drøbak and Egersund), and the most common haplotype by 11 individuals (Sweden, Norway: Bergen, Egersund, Kristiansund and Stavanger). Apart from the Azorean specimens which clustered together, the network showed a general lack of geographic structure. The haplotypes from the Azores differed between themselves by two to three base pairs (0.3-0.5%), whereas the difference to the closest haplotype outside the Azores (two haplotypes from Norway) was four base pairs (0.7%).
For Polycera norvegica sp. nov. (Figure 12) with 17 specimens (three from northern Norway and 14 from southern-western Norway), the TCS haplotype network revealed a lack of shared haplotypes and geographic structure.

Discussion
Molecular taxonomy has contributed in recent years to unravel a large number of cryptic lineages in sea slug molluscs, otherwise difficult to distinguish using traditional morphological and anatomical characters (see Introduction for an overview). Nevertheless, the concept of 'cryptic species' has been debated and remains controversial and subject to several definitions and interpretations (see for reviews Struck et al., 2018;Korshunova et al., 2019). Korshunova et al. (2019) considered cryptic species to be transitional taxonomic problems resulting often from incorrect weighting of morphological characters.
The phylogenetic and ABGD analyses (Figure 2 and supplementary material) retrieved one group containing specimens consistent with the original description by O. F. Müller (1776Müller ( , 1779  and collected in Drøbak in the Oslofjord, Norway (Figure 1), which is the type locality of Polycera quadrilineata. These specimens were therefore attributed to the latter species. Whereas, specimens preliminarily ascribed to P. quadrilineata, but which clustered elsewhere in the phylogenetic analysis and grouped separately on the ABGD analyses were attributed to a new species here described, namely Polycera norvegica sp. nov. In the latter species, one of the recognized colour morphs (the dotted colour morph) is in fact fairly distinct from the colour morphs found in P. quadrilineata, but surprisingly conspecificity was never questioned. However, in specimens of both species with the yellow/orange colour morph, separation based only on external chromatic and morphological features is difficult, virtually impossible. In addition, the variability found in the number of lamellae of the rhinophores and gills and maximum length of both species (Results section; Table 3) do not seem to be taxonomically informative and may in fact be ontogenetic or ultimately an artefact of sampling bias. The species P. quadrilineata and P. norvegica sp. nov. are genetically distinct (Table 2; COI uncorrected p-distance = 9.6-12.4%) and bear several discrete anatomic differences in the radula, labial cuticle and reproductive system. The radula of P. quadrilineata is larger, thicker, more elongate, with a greater number of rows, and has an additional marginal tooth on each side. The labial cuticle in P. quadrilineata is comparatively thicker, more robust and with a stronger central brownish area with jaw elements. The reproductive system is substantially distinct between both species. In P. quadrilineata the ampulla is larger and more robust and the bursa copulatrix depicts two different parts with a larger and elongate proximal part and an oval-shaped distal part, whereas P. norvegica sp. nov. has a bursa copulatrix entirely oval. Polycera norvegica sp. nov. has a larger prostate gland and a less developed penial bulb. Both species have armed penis but we have observed in P. norvegica sp. nov. the occurrence of two types of chitinous spines, whereas only one type was observed in P. quadrilineata. Moreover, the tubercles of P. quadrilineata are frequently more smoothly rounded than those found in P. norvegica sp. nov. (see Table 3 for a synopsis of diagnostic characters).
Interestingly these two species were not rendered sister to each other. The COI gene tree phylogeny suggested a closer phylogenetic relationship of P. norvegica sp. nov. with the European species Polycera faeroensis and an apparently undescribed species from South Africa (Polycera sp.A), whereas P. quadrilineata seems to be closer related to the north-eastern Atlantic Polycera aurantiomarginata and Indo-West Pacific Polycera capensis.
All but one sister pair of Polycera species studied in this work have interspecific genetic distances ranging between 8.8-19.7% (Table 2) and intraspecific variabilities between 0-2.6% showing the existence of a clear DNA barcode gap between species. The single exception is between P. capensis and P. aurantiomarginata with a genetic distance ranging from 4.3-5.8%, yet this does not affect the existence of a molecular barcode gap between sister species among the Polycera species studied in this work.
The haplotype network analyses showed genetic diversity and a lack of geographic structure in both species. The fact that the Azorean specimens showed a structural group in the haplotype network of P. quadrilineata could potentially indicate some degree of genetic isolation (Figure 11), which is not surprising given the geography of the Azorean islands located far off in the middle of the North Atlantic Ocean. Shared haplotypes between species from Norway and Sweden hint at a possible gene flow across populations from these two neighbouring countries. Interestingly, no shared haplotypes between sampling localities in Norway were found for P. norvegica sp. nov., but this may simply be because of the low number of specimens included in the analysis (N = 17; Figure 12).
Polycera quadrilineata and P. norvegica sp. nov. are sympatric in Norway, dwelling together on shallow water Laminaria kelp forests. The species P. quadrilineata has been reported to feed upon the encrusted bryozoans Membranipora membranacea (Linnaeus, 1767) and Electra pilosa (Linnaeus, 1767), but the differences found on the digestive system of P. norvegica sp. nov. with its weaker labial cuticle and distinct radula suggest a possible different diet.
This study showed that even in historically well-studied faunas there is a need to use integrative taxonomic approaches combining DNA and morphological characters together with complementary methodological approaches to truly understand biological diversity. This becomes even more critical when dealing with species that are highly similar morphologicallythe so-called 'cryptic' or 'pseudo-cryptic' species. In those cases, even when subtle morphological differences can be recognized the problem is obviously to access and decide whether the observable differences are part of the natural intraspecific variability of the species or indicate distinct taxa. Drawing this line is often difficult and this is when a barcoding approach can be useful to unravel the putative existence of 'hidden' lineages. However, this does not preclude the need to eventually study several gene markers with different evolutionary rates, and for detailed morphological studies in order to establish helpful discrete taxonomic characters and formalize the description of new recognized lineages. It is now recognized that cryptic species are an important component of biodiversity and the only way to uncover this fraction is by using integrative approaches combining and comparing morphological and genetic characters under a multitude of methodological frameworks.
Supplementary material. The supplementary material for this article can be found at https://doi.org/10.1017/S0025315420000612.