Morphological and Molecular Assessments of Bobtail Squids (Cephalopoda: Sepiolidae) Reveal a Hidden History of Biodiversity

Molecular species delimitation assists taxonomic decisions for challenging species, like cryptic species complexes. Bobtail squids (Family Sepiolidae Leach, 1817) are a very diverse group of benthic and nektonic small to medium size cephalopods with many taxonomic questions to solve. In this study we provided new sequence data for 12 out 17 Mediterranean bobtail squid species including all the genera present i n the area. Other relevant species from other parts of the world were used as comparison. The combined use of several molecular species delimitation methods consistently showed a picture of hidden biodiversity within this family which hinders the use of molecular data isolated from morphological characters. On the one hand, those methods provided contrasting results for the number of recognized species of some morphologically well-defined species. We suggest this can be an effect of recent speciation phenomena followed by an intense morphological drift. On the other hand, cryptic biodiversity was detected among members of several monophyletic clades assigned to the same nominal species, pointing to recent speciation phenomena without a parallel morphological evolution. Although Mediterranean bobtail diversity has been extensively studied for more than a century, a new species of Stoloteuthis Verrill (1881) was discovered and described here, both using molecular and morphological methods. This new research stresses the necessity of combined morphological and molecular studies to correctly assess cephalopod diversity. urn:lsid:zoobank.org:act:57AFBB38-18EA-4F80-B1D4-73519C12694F.

The family Sepiolidae Leach, 1817, commonly known as bobtail squids, is a very diverse group of benthic and nektonic small to medium size cephalopods. In recent years, the taxonomy and systematics of the family has been dynamic, including a few molecular phylogenetic studies (e.g., Groenenberg et al., 2009;Sanchez et al., 2019), major systematics reviews (Bello, 2019(Bello, , 2020 and the description of a relatively large number of new species (e.g., de Heij and Goud, 2010; Kubodera and Okutani, 2011;Sanchez et al., 2019). In the last decade, new species have been even discovered in the Mediterranean Sea (Bello, 2013;Bello and Salman, 2015), where bobtail squid diversity has been extensively studied for more than a century (Bello, 2015(Bello, , 2019. The closure of the Strait of Gibraltar is an historical process that may explain the contemporary Mediterranean sepiolid endemism and high species richness, especially in the western area (Mangold and Boletzky, 1988;Bello, 2003;Rosa et al., 2019). Sepiolid systematics and taxonomy rely mostly on the morphology of the light organs or photophores, the male copulatory organ or hectocotylus, and the female sperm storage organ or bursa copulatrix (Reid and Jereb, 2005;Bello, 2020). Three subfamilies are recognized: Heteroteuthinae Appellöf, 1898, Rossiinae Appellöf, 1898and Sepiolinae Leach, 1817. Even though hectocotylus morphology is a reliable character, recently it was discovered that the previously recognized intraspecific variability in one species undercovered pseudocryptic biodiversity (Groenenberg et al., 2009;de Heij and Goud, 2010). Moreover, identifications based on early life stages and females are challenging and misidentifications are abundant on GenBank (Groenenberg et al., 2009), hindering identification based only on DNA barcoding.
Here, we examined most of the Mediterranean biodiversity for the family Sepiolidae, performed several molecular delimitation methods in order to assess the actual diversity of this group, and provide a solid molecular framework for future studies based on DNA identification methods. Additionally, we discovered a new species of Heteroteuthinae, which is described here using both molecular sequences and morphological characters.

Sample Collection
In total, 77 newly collected bobtail squids were analyzed (Table 1), covering 12 of the 17 known Mediterranean species and all the genera (Bello, 2015(Bello, , 2020. Sampled individuals were mostly collected in the Mediterranean coasts of the Iberian Peninsula from Tiñoso Cape in the south to Blanes in the north, including the Balearic Islands, between the years 2006 and 2015. Most of this material was collected during the Spanish research cruises MEDITS 2006 (https://www.sibm.it/ SITO%20MEDITS/principaleprogramme.htm) and FORMED 5 (Demestre et al., 2017). Additional material were collected by commercial fishing trawlers from Tarragona and Vilanova i la Geltrú fishing ports and from the "sonsera" littoral artisanal fishery from Blanes (Lleonart et al., 2014), which is a small scale littoral boat seine fishing method aimed to catch sand eels. Four individuals of Heteroteuthis dispar (Rüppell, 1844) were collected in NE Atlantic waters near the Canary Islands in the spring of 2015, during the MAFIA research cruise (Olivar et al., 2017). Figure 1 summarized the collection localities. After collection, the individuals were frozen at −20 • C until their study in the lab. In a time lapse ranging from a few weeks to 8 years, the material was defrosted, identified, and photographed. From each specimen, a portion of the mantle was removed for DNA extraction and the rest of the body was preserved as a morphological voucher.

Morphological Identification and Vouchering
Individuals were identified following the morphological key of Bello (1995). In the case of the genera Sepietta (Naef, 1912a) and Sepiola Leach, 1817, only males were identified based on the hectocotylus morphology to ensure reliable identifications. After removal the tissue for DNA extraction, the specimens were fixed in 4% buffered formalin for 3-10 days and transferred to ethanol 70%. The specimens are deposited in the Biological Reference Collections (CBR-ICM) at the Institut de Ciències del Mar (ICM-CSIC, Barcelona, Spain) under the accession numbers provided in Table 1 (Guerrero et al., 2020).
The description provided here for a new butterfly squid is based in a careful morphological examination of all the available specimens (n = 4). The measurements and morphometric characters of these specimens followed Roper and Voss (1983) as: dorsal mantle length (DML), ventral mantle length (VML), fin width (FW), fin length (FL), fin base (FB), head width (HW), head length (HL), funnel length (FnL), arm I-IV length (AIL-AIVL), tentacle length (TeL), tentacle club length (CL), and web depth A-E (WDA-E). Two additional measurements were taken: the occipital band length (mantle-head fusion) (OBL), defined as the length of the fusion between the head and the mantle, and the ventral shield length (VSL), as the shield length along its midline. All the morphological measurements from this study were performed on formaline-fixed individuals stored in 70% ethanol.
Beaks and radulae were extracted from selected individuals. For the beaks, the upper and lower rostral lengths (URL and LRL, respectively) were measured according to Clarke (1986). The radula was observed under a Hitachi S3500N scanning electron microscope (SEM). At the beginning of SEM preparation, the radulae were dehydrated in an increasing concentration of ethanol (80, 90, and 96%) until they were saturated in absolute ethanol. Each ethanol bath lasted 10 min. After complete dehydration in the ethanol series, the samples were dried to a critical point using CO 2 as the transition liquid. After the drying stage, samples were mounted on stubs with double-sided conductive sticky tape to place them in the preferred position. The mounted samples were sputter coated with gold-palladium before SEM observations. The spermatophores from selected males were extracted for the assessment of the spermatophore count (SpC), and the spermatophore length (SpL) based on 30 randomly selected spermatophores that were measured according with 1 | List of specimens used in this work, indicating the sampling locality, the number of studied sequences (n), and the GenBank Accession Numbers.

Species
Locality n Morphological voucher GenBank Accesion Number

DNA Extraction, Amplification, and Sequencing
Tissues for molecular analysis were fixed in 96% ethanol. Total genomic DNA was extracted from an ethanol-fixed piece of the mantle using the NZY Tissue gDNA isolation kit (NZYTech, Lisbon, Portugal), following the manufacturers' protocol and resuspended in a final volume of 100 µL. A negative control that contained no sample was included in every isolation round to check for contamination during the experiments. Sequences from the partial mitochondrial cytochrome c oxidase I (COI) gene were amplified, using the primer pair LCO1490 and HCO2198 FIGURE 1 | Sampling localities of the specimens sequenced in this work sorted by subfamilies. Sampling localities of Stoloteuthis cthulhui sp. nov. are also indicated. Modified from Google Earth Pro. (Folmer et al., 1994). Standard PCR reactions were performed using the NZYTaq Green PCR Master Mix (NZYTech, Portugal) following the manufacturer's protocol in a total volume of 25 mL, which included 0.5 µM of each primer, 25 ng of template DNA and PCR-grade water up to 25 µL. PCRs consisted of an initial denaturation at 95 • C for 5 min, followed by 35 cycles of denaturation at 95 • C for 30 s, annealing at 50 • C for 30 s and extension at 72 • C for 45 s, with a final extension of 5 min at 72 • C. The amplified products were sequenced using both forward and reverse PCR primers on an ABI 3730xl. DNA sequence data were edited and aligned with Geneious 8.1.5 (http://www.geneious. com). The GenBank Accession numbers of the sequences used in this work together with the morphological voucher Accession numbers are summarized in Table 1.

Phylogenetic Analyses
The new sequences obtained in this work together with selected sequences available from GenBank were analyzed (Table 1). Idiosepius pygmaeus (Steenstrup, 1881) was selected as outgroup for the phylogenetic analyses. Some sequences coming from GenBank were shorter, so Ns were added to align them with the complete sequences. The final alignment contained 245 sequences and 658 positions. As an initial analysis, a Maximum Likelihood tree was obtained through the IQTree portal (Trifinopoulos et al., 2016) [available at http://iqtree.cibiv. univie.ac.at] using the automatic model selection feature. The selected model was TIM2+F+I+G4 according to both Akaike and Bayesian Informative Criteria. The support of the branches was calculated after 2,000 ultrafast bootstrap generations.
When possible, species identifications of the clades were based on voucher specimens. Several sequences uploaded to GenBank from previous works were originally based on misidentified specimens. The correct identifications, the original identifications, their GenBank Accession numbers and the original references are summarized in Table 4.

Molecular Species Delimitation
Thirty seven different clades of bobtail squids were identified based on molecular and morphological data (Figure 2). The individuals examined in the present study and identified based on morphological characters always clustered together in a single clade, but some inconsistencies were detected on previously published sequences. Some of them are amended according with our identifications over other members of their respective clades ( Table 4). Sepiolinae sp. 1 was not originally identified at the species level (Strugnell et al., 2004) and the only sequence available did not match with any correctly identified clade. Sepiolinae sp. 2 was originally identified as Sepiola affinis (Naef, 1912c) (Lindgren et al., 2004; Table 4), but this sequence did not group with the correctly identified S. affinis. The sequences AY293710 and AY293715 were described as Lusepiola birostrata (Sasaki, 1918) or Euprymna morsei (Verrill, 1881) and their divergence is compatible with an intraspecific distance. Sanchez et al. (2019) assigned those sequences to Sepiola (= Lusepiola) birostrata. Euprymna hyllebergi, E. scolopes, and E. tasmanica showed 2-4 highly divergent clades. The species delimitation methods provided conflicting results. The eight partitions of the ABGD identified from 26 to 41 groups depending of the prior maximal distance. It identified 28 groups with a maximal intragroup distance of 2.1%, 33 groups with maximal distances of 0.04-1.2% and 41 groups for 0.2-0.1 maximal distances. Based on the distribution of the distances of the whole dataset, no discrete barcoding gap was detected. Figure 2 represents the species assemblages based on the results from ABGD with a maximal distance prior of 1.2%, since it is closer to the observed maximal intra-specific distance observed in the dataset. The TCS analysis with 95% of maximum connectivity identified 34 networks. Sepiola affinis, Sepiola atlantica d'Orbigny, 1842 in Férussac and d'Orbigny, 1834-1848, and Sepiola intermedia (Naef, 1912d) formed a single network, while E. hyllebergi, E. scolopes, and E. tasmanica were split in 2-4 networks (Figure 2). The TCS analyses at 98 and 99% of divergence differed in the relations between S. affinis, S. atlantica, S. intermedia and the sequence DQ646730 (identified as E. tasmanica). In the 99% analysis, the three species and DQ646730 formed four independent networks, while in the 98% analysis S. affinis, S. atlantica, and DQ646730 formed a single network. Both 98 and 99% analyses over-split Rossia macrosoma (delle Chiaje, 1830), Heteroteuthis dagamensis (Robson, 1924), Heteroteuthis hawaiinensis (Berry, 1909), and Sepiola robusta (Naef, 1912c) in two networks each. The Maximum Likelihood solution of the bPTP analysis recovered 38 species majorly consistent with the results of the TCS 95% analysis and the morphological and molecular assignations of the Maximum Likelihood tree (Figure 1). However, two species were detected for E. scolopes OTU 1 and three for H. dagamensis (one including all New Zealand specimens and the two others included North Atlantic specimens); and S. affinis, S. atlantica, and S. intermedia were recognized as a single species. The highest Bayesian supported solution recognized 83 species, and S. intermedia was recognized as a different species than S. affinis and S. atlantica (results not shown). The single threshold method of the GMYC identified 28 clusters (confidence interval 20-32) and 41 entities (confidence interval 26-48) with a significant Likelihood Ratio test (LR = 4.312508e-07). Sepiola affinis including the sequence DQ646730, S. atlantica and S. intermedia were recognized as three independent species. Two species were recognized for H. dagamensis: one including all the Atlantic individuals and another formed by the New Zealand specimens. Rossia macrosoma and Rondeletiola minor (Naef, 1912c) were split in two species each. The multi-threshold method revealed 31 clusters (confidence interval 24-31) and 43 entities (confidence interval 31-43) also with a significant result (LR = 1.928034e-08). In this analysis S. affinis including the sequence DQ646730 and S. intermedia were recognized as independent species, but S. atlantica was recognized as two species. Rossia macrosoma, E. hyllebergi OTU 2 and Sepietta oweniana d'Orbigny in Férussac andd'Orbigny, 1834-1848 were also recognized as two species each, while E. tasmanica OTUs 3 and 4, and Stoloteuthis japonica Kubodera and Okutani, 2011 and Sepiolina nipponensis (Berry, 1911) were merged in single species. It also split Atlantic and New Zealand H. dagamensis in two different species. Uncorrected p-distances across species ranged from 0.81 to 18% (Table 2; mean: 12.7%) and from 0 to 2.3% ( Table 3; mean: 0.3%) at an intraspecific level. Regarding interclade distances, it should be noted that most interspecific distances were values above 3% [in accordance to the values for the family Sepiolidae published by Gebhardt and Knebelsberger (2015)], but the distances were lower between Se. nipponensis and Sl. japonica (2.4%); H. dispar, Heteroteuthis sp. KER [from the Kermadec Islands, see Braid and Bolstad (2019)] and H. hawaiinensis (2.5-2.8%); and between S. intermedia, S. atlantica, and S. affinis (1.45%). The identification of the last three species is assured since it was confirmed morphologically (Groenenberg et al., 2009 for S. atlantica; this work for S. affinis and S. intermedia, Figure 3). Regarding intraclade distances, in H. dagamensis distances above 2% were reported ( Table 2).
Among the analyzed Rossinae species [Neorossia caroli (Joubin, 1902), Rossia macrosoma (delle Chiaje, 1830), and Semirossia tenera (Verrill, 1880)] interspecific distances ranged from 4 to 12.2% ( Table 2) and at intraspecific level ranged between 0 and 1.4%. Among R. macrosoma two different clades diverging 1.4% were identified. One of them was formed by four Mediterranean individual while the other one was formed by one Mediterranean and eight North Sea individuals. Among Heteroteuthinae, interspecific p-distances ranged from 2.4 to 17.8% (Table 2) and the intraspecific p-distances range was 0-2.2% (Table 3). Between the Mediterranean Stoloteuthis individuals and Sl. leucoptera the distance levels were typical for interspecific distances. Among Sepiolinae, the interspecific pdistances ranged from 0.81 to 15.8%, but if S. affinis, S. atlantica and S. intermedia are excluded, the lowest intraspecific distance detected is 3.4% and there is no overlapping between intra-and interclade p-distances (Tables 2, 3), as occurs in the other available sequences of the other two subfamilies. In the genus Euprymna (Steenstrup, 1887), cryptic biodiversity was found in E. hyllebergi and E. tasmanica, formed by 2 and 4 highly divergent clades, respectively (Figure 2, Tables 2, 3). Between the two clades of E. hyllebergi, a distance of 3.9% was reported. Distances of 4.3-11.6% are found between the four divergent clades of E. tasmanica. The TCS and bPTP analyses found cryptic biodiversity in E. scolopes, although the ABGD, the GMYC, and the p-distance analyses did not find it.

Diagnosis
Stoloteuthis with a maximum size of 18 mm of mantle length; with a narrow occipital band; with a ventral shield of around 80% of the ventral mantle surface; with a wide head; tentacles representing 251-379% of the mantle length; males with glands in the first two thirds of both dorsal and ventral margins of arms I; males with rows 2-4 of ventral suckers slightly enlarged, ventral, and dorsal rows 5-6 enlarged to the same level in arms II; males with 3-4 series of suckers at the tip of arms IV.

Etymology
The specific epithet cthulhui was erected in honor to the fiction cosmic horror entity "Cthulhu" created by Howard Phillips Lovecraft (1890Lovecraft ( -1937, which holds both cephalopodlike tentacles and wings (Lovecraft, 1928), resembling the pair of fins of this new species of butterfly squid. Cthulhu was described in the literature as male, so a male gender genitive suffix was applied. Pronunciation: /k e 'θu:lu:/.

Description
Butterfly squid up to 18 mm DML (Figures 4A,B). The body is muscular with a rounded posterior end. The mantle dorsally is fused directly with the head forming a narrow commissure of 5-8 mm and is heavily and uniformly pigmented. Ventral mantle is anteriorly bilobed with roughly the same extension of dorsal mantle or slightly larger, with a ventral shield of dark brown pigment and an iridescent greenish hue roughly representing nearly 80% of the ventral mantle surface. This ventral shield is surrounded by a shallow yellowish bright band. Fins are large, oval, and less pigmented than the mantle, attached dorsally in the lateral mantle. The mantle component of the funnel/mantle locking-apparatus is straight with a low ridge at the mantle edge and the funnel component is a straight simple groove. The funnel reach the anterior edge of the eyes and the pigmentation forms a diamond patch from the funnel tip toward the mantle; near the tip of the funnel there is a round patch devoid of chromatophores, and the funnel sections covered by the lobes of the mantle also lack of chromatophores. The head is bulbous, wide and uniformly pigmented, representing from 88 to 145% of the ML. Eyes are large and occupy most part of the head. Olfactory organs are prominent. Arms are short and muscular with a well-developed keel in arms III. Arm formula: II ≈ III > I > IV. Two rows of suckers in both arms, but in males the arms IV might appear to have four rows. Female arm crown lacks of any enlarged suckers. Arm suckers in all arms more developed in males than in females. Arms I of mature males with well-developed glands in the first two thirds of both dorsal and ventral margins. Arms II suckers from rows 2-6 are slightly enlarged according with the following pattern: rows 2-4 of ventral suckers are slightly enlarged and the ventral and dorsal rows 5-6 have the same degree of enlargement ( Figure 4C). A well-developed web unites the arms until approximately the last third of the arm length from arms I-III, arms IV united by a shallower web reaching less than a quarter of the arm length or absent (individual ICMC000165). Tentacles range from 44 to 46 mm length, representing 251-379% of ML. The distal 5-6 mm of the tentacles are occupied by the club. The tentacle organ slightly overlaps the tentacle club, formed by 10-11 transversal rows of suckers, which are larger in both proximal and terminal rows of the club.
The light organ has the typical shape for Heteroteuthinae (Figure 4B), with a round morphology and two pores in the midline. The visceral mass was not dissected. The upper beak has a long and curved rostrum of 0.8-0.9 mm length with a sharp pointed rostral tip (Figure 4D), representing ∼4.9% of the ML. Beak angle of ∼90 • , with a smooth edge; specimen ICMC000165 has an irregular tooth. Hood is small and fragile. The rostrum and the shoulders are darkly pigmented, while the hood and the lateral wings are moderately pigmented. Lower beak has a short and blunt rostrum of 0.6-0.7 mm length with a rounded rostral tip (Figure 4E), representing ∼3.6-3.7% of the ML. The lower beak angle is of ∼130 • , without teeth. The hood is narrow and the wings are wide. The rostrum and the hood are darkly pigmented, while the wings and lateral walls are moderately pigmented and almost transparent at the edges. The male specimens ICMC000165 and ICMC000166 and the female specimen ICMC000165, 17.9, 16.7, and 12.3 mm ML, respectively, were used for measuring the radulae. The radula has seven transverse rows of teeth displaying a typical homodont morphology ( Figure 4F). Rhachidean teeth are 96-126 µm length, with a wide base, sharply pointed and slightly curved shape. First lateral teeth are 88-130 µm length, with a narrower base, sharply pointed and slightly curved shape. The second lateral and marginal teeth are narrow and strongly curved, of 98-164 µm and 176-240 µm length, respectively. All the spermatophores of specimens ICMC000165 and ICMC000166 were extracted and counted, and a random selection of 30 were used for measuring their length. ICMC000165, 17.9 mm ML, had 236 spermatophores of 2.6 ± 0.1 mm length (range 2.5-3.0 mm). Spermatophores were intact and the spermatophoric reaction was not triggered (Figure 4H). Spermatophore threads were short and in most cases they were broken. Spermatophores head measures 0.28 ± 0.01 mm and holds three loops of the ejaculatory ducts. The ejaculatory apparatus is 0.76 ± 0.03 mm long and the spiral filament occupies ∼¾ of its length. The cement body is relatively long (0.64 ± 0.05 mm) and approximately half of the length of spermatophore is occupied by the seminal reservoir (1.32 ± 0.05 mm). There is no posterior empty part, but in those spermatophores fixed during the first steps of the spermatophoric reaction, a posterior empty part is observed. Specimen ICMC000166, 16.7 mm ML, had a strongly enlarged spermatophoric sac with a massive amount of spermatophores (11,645) ranging from 1.4 to 2.4 mm length. Smaller spermatophores usually look similar to the larger ones, but with reduced seminal reservoirs. Other spermatophores of this specimen were empty and/or deformed and should be considered as tentative spermatophores sensu Nigmatullin et al. (2003).

Molecular Data
The COI p-distances analyses (Tables 2, 3) show a 3.5% divergence between Sl. cthulhui sp. nov. and its sister species, Sl. leucoptera. This level of divergence was reported by Gebhardt and Knebelsberger (2015) and Groenenberg et al. (2009) between different closely-related bobtail squid species. Moreover, the same level of divergence was found between different well-established species in our analysis ( Table 2).

Remarks
Morphological differences between Sl. cthulhui sp. nov. and Sl. leucoptera are subtle but still consistent. The morphology of the second pair of arms of the mature male differs in the assemblage of the enlarged suckers. In Sl. cthulhui sp. nov., the pattern is less evident than in Sl. leucoptera: the first row of suckers is unmodified, rows 2-4 of ventral suckers are slightly enlarged and ventral and dorsal rows 5-6 of suckers showed the same level of enlargement, whereas in Sl. leucoptera, the first two rows of suckers are unmodified and ventral suckers from rows 3-6 are progressively enlarged (Figure 4I, see also Vecchione and Young, 2013). In some Sl. leucoptera specimens the fifth ventral sucker is larger than the sixth or dorsal suckers 5-6 might be larger than the ventral ones. Although the range of values overlaps in all cases, some differences also exist in some morphometric indexes ( Table 6): the fin length and shield length are larger in Sl. leucoptera than in Sl. cthulhui sp. nov., while the head is narrower (68-96% of ML) and the tentacles are shorter (115-299% of ML). The rhachidean teeth of the Sl. leucoptera radulae are smaller, with 41-85 µm length, measured in two mature females and two mature males 13.3-15.5 and 11.5-13.3 mm ML, respectively ( Figure 4G). Although the remaining teeth of Sl. leucoptera also tend to be smaller, the length ranges overlap with Sl. cthulhui sp. nov. Other characters, such as other morphometric measures, the beaks morphology and the spermatophores did not show any other relevant morphological differences between the two species. Stoloteuthis leucoptera was described in the Gulf of Maine, 30 miles east Cape Ann (North-western Atlantic Ocean) (Verrill, 1878). Previously to the present work, Sl. leucoptera was thought to be distributed in NW Atlantic Atlantic to Namibian and Mediterranean waters (Reid and Jereb, 2005;Vecchione and Young, 2013). We did not find any morphological differences between Namibian Sl. leucoptera (Villanueva and Sánchez, 1993; this work) and the descriptions from the literature for NW Atlantic Sl. leucoptera (e.g., Vecchione and Young, 2013). The individuals described by Orsi Relini and Massi (1991 : Table 1) and Cuccu et al. (2010 : Table 1) had a HWI of 105-130%, consistent with Sl. cthulhui sp. nov. and not with Sl. leucoptera. We consider previous records of Mediterranean Sl. leucoptera as Sl. cthulhui sp. nov. (Orsi Relini and Massi, 1991;Volpi et al., 1995;Cuccu et al., 2010;Quetglas et al., 2013;Zaragoza et al., 2015).

DISCUSSION
All of the genera and twelve sepiolid species, covering 70% of the known specific biodiversity of the family in the Mediterranean Sea, have been barcoded and vouchered ( Table 1). All the studied individuals were successfully linked with their taxonomic name and no inconsistencies arose among the newly sequenced material. Nevertheless, wrong identifications of bobtail squids are relatively frequent in GenBank, as previously reported by Groenenberg et al. (2009) andSanchez et al. (2019). The work of Groenenberg et al. (2009) on DNA barcoding on vouchered individuals has solved some of the previously misidentified sequences. However, this work was based on animals from NE Atlantic waters and some species not studied by them remained unsolved. Here, we provided a correct barcode for many of these species, solving some of the previous problematic sequences ( Table 4). Solving those problematic sequences based on morphologically identified animals is extremely important to ensure the quality of the identification based only on molecular data.
For most of the species studied here, there is a tendency of intraspecific distances below 2% and interspecific distances above 2.4-3% (Tables 2, 3). However, this pattern is not universal within the family and some exceptions occur. In fact, no clear barcode gap (Meyer and Paulay, 2005) was identified between intra-and interspecific distances, due to the presence of challenging groups, such as the clade formed by S. affinis, S. atlantica and S. intermedia. The interspecific distances of these three species are the smallest among our dataset (0.81-1.45%). In fact, the largest intraspecific p-distances among Sepiolinae (1.8% in S. robusta and 2.3% in E. scolopes, Table 3) are larger than the interspecific distances between S. affinis, S. atlantica and S. intermedia, thus existing an overlapping between intra-and interspecific distances which complicates the use of DNA barcoding methods based on genetic distances for this group of animals. The three species have very distinctive hectocotylus morphology with discrete morphologies (Figure 3): the morphological variation of those species does not overlap (Bello, 1995;Reid and Jereb, 2005;de Heij and Goud, 2010). As far as we know, no hybrids have been described among them. While S. atlantica is allopatric in reference with the other two species, S. affinis and S. intermedia both occur in Mediterranean waters (Reid and Jereb, 2005), pointing out to the presence of effective reproductive isolation mechanisms acting at least between the two Mediterranean species. For this clade, only the single threshold approach of the GMYC among all the tested molecular species delimitation methods provided the same results as the morphology (Figure 3). If molecular identifications were carried out with no further morphological information, this level of interspecific distances might be mistaken as intraspecific variation. This low level of interspecific distances might be due to recent phenomena of speciation with a fast morphological drift of key morphological characters. Recently, Costa et al. (2021) found that two morphologically different species of coastal squids that diverged in recent times were recognized as a single species by molecular species delimitation methods. Their study and ours highlights the importance of combining studies based on molecular identifications with careful morphological examinations.
Another important phenomenon hindering the direct use of bobtail squid DNA sequences for species identifications is the presence of cryptic biodiversity. Although cryptic biodiversity is an increasingly reported phenomenon in cephalopod biodiversity studies (e.g., Anderson et al., 2007, Cheng et al., 2014, this phenomenon is comparatively unknown in bobtail squids. This is especially true for E. hyllebergi and E. tasmanica, with interspecific distances between different OTUs ranging from 3.9 to 11.6%. Remarkably, the two E. hyllebergi OTUs occur allopatrically in both Indian and Pacific coasts of the Thailand Peninsula, while E. tasmanica OTUs might occur sympatrically. It is known that the hatchling size and mode of life has an important effect on the distribution range of cephalopod species, since species with larger benthic hatchlings tend to have smaller distribution areas (Villanueva et al., 2016). Information on the way of life of bobtail squid during their first days of life is scarce (Villanueva et al., 2016: Tables 1, 2). Available information points out that Rossinae and most Sepiolinae tend to have large benthic hatchlings, although some early stages of some species reported as benthic throughout their lives can be also be found in the water column (Olmos-Pérez et al., 2018b). It is particularly remarkable that the species with larger hatchlings from the subfamily Sepiolinae, E. tasmanica (5 mm ML, Villanueva et al., 2016: Table 1), has more cryptic lineages with larger interspecific distances, accounting for four different OTUs with a divergence of 4.3-11.6% (Table 2). Interestingly, R. macrosoma has slightly larger benthic hatchlings (5.5 mm ML) and also has genetic structure, being taken as two species by several species delimitation methods. Comparisons between the sympatric geographic patterns of E. tasmanica lineages with the allopatric pattern found in Euprymna hyllebergi, a species with smaller planktonic hatchlings (2.2 mm ML), suggest that large hatchlings with direct benthic development have a more intense effect on the dispersal capacity and communication between distant cephalopod populations, and it might be one of the triggers that increase the opportunities for speciation. Rossinae are mostly exclusively benthic species, and Sepiolinae usually are usually reported as benthic species (Reid and Jereb, 2005), but they can also be found in the water column (Bello and Biagi, 1995). Absence of ontogenetic migrations have been suggested for some Sepiolinae species (e.g., Villanueva, 1995), which suggests that this restrictions to movement limit dispersal, contributing to slow down the gene flow between distant populations. This putative lower dispersal during both young and adult stages combined by the fact that some species inhabit relatively non-overlapping bathymetric ranges might have helped to trigger speciation in the Mediterranean Sea and increase their endemic bobtail species (Bello, 2019). Bobtail squids are among the cephalopods with higher tolerances to salinity changes (Mangold and Boletzky, 1988) and young stages can be found in estuarine systems (e.g., Olmos-Pérez et al., 2018b), and adults in the intertidal regions (Fernández-Álvarez, pers. obs.). Thus, it seems that salinity it is not a great limitation to bobtail squid dispersal. The Mediterranean Sea Sepiolinae fauna is characterized by the high number of endemic species (Bello, 2019), while many Atlantic species does not distribute also in the Mediterranean. It is possible that the effect of some well-known oceanographic barriers to genetic exchange, such as the Strait of Gibraltar (Pascual et al., 2017), had an important effect on the evolutionary history and current distribution of European Sepiolinae, while the differences in salinity between the Mediterranean and the Atlantic (Mangold and Boletzky, 1988) might have a comparatively smaller effect. Future population genetic studies focused on those species can answer this question.
Members of the subfamily Heteroteuthinae are exclusively nektonic and benthic-pelagic species (e.g., Orsi-Relini, 1995). Although it is known that oceanic currents can isolate populations of some pelagic squids and lead them to speciation (e.g., Fernández-Álvarez et al., 2020), the opportunities to dispersal and population connectivity are comparatively larger in oceanic environments than in shallow benthic ecosystems. Thus, large distribution areas on wide oceanic basins are commonly reported in Heteroteuthinae species (Reid and Jereb, 2005). According with Vecchione and Young (2007), there are no known morphological differences between H. dispar and H. hawaiinensis, while all the molecular species delimitation methods performed here support their treatment as different species. In the absence of known morphological differences, these two species shall be considered as members of a cryptic species complex. A third undescribed species with a not yet described morphology, Heteroteuthis sp. KER (Braid and Bolstad, 2019), form a clade with H. dispar and H. hawaiinensis. Heteroteuthis dagamensis have a large distribution ranging from the Gulf of Mexico and the South Atlantic (Judkins et al., 2016) to New Zealand (Braid and Bolstad, 2019) and North Atlantic (Taite et al., 2020). Four out six of the molecular species delimitation analyses recognized cryptic biodiversity within this species (Figure 3). The large distribution range of the species in combination with the large intraspecific divergence between New Zealand and the specimens coming from other latitudes (>2% , Table 3) suggests that some processes of speciation in its early stage might be taking place. Between those described cryptic Heteroteuthis species (H. dispar and H. hawaiinensis), undescribed new species (Heteroteuthis sp. KER) and cryptic lineages (H. dagamensis) some well-known oceanic and terrestrial barriers exists, such as the Panama Isthmus and the currents that creates the main oceanographic gyres. New combined morphological and molecular studies are necessary for solving this taxonomic problem.
It is remarkable the fact that a closer relationship exists between Sl. japonica and Se. nipponensis (2.4%) rather than between Sl. japonica and other congeneric species (12.7% with Sl. cthulhui sp. nov. and 14.3% with Sl. leucoptera). Divergence between Sepiolina petasus Kubodera and Okutani, 2011 and Se. nipponensis is 13.9%, similar to that reported between species of the genera Sepiolina and Heteroteuthis. These data suggest that the current generic assignations to genera in Heteroteuthinae species are not fully molecularly supported and should be revised in future combined morphological and molecular studies, as already suggested by Allcock et al. (2014).
All molecular species delimitation analyses consistently identified Sl. leucoptera and Sl. cthulhui sp. nov. as different species (Figure 3), also the COI divergence of 3.5% is typical for different species in bobtail squids. This difference is consistent with interspecific levels in many other invertebrates, such as nemerteans (Fernández-Álvarez and Machordom, 2013), land planarians (Lago-Barcia et al., 2015), crustaceans (Robles et al., 2007), and other cephalopods (Gebhardt and Knebelsberger, 2015;Fernández-Álvarez et al., 2020). Besides, the morphological comparisons between Sl. leucoptera and Sl. cthulhui sp. nov. showed a few differences in key morphological characters that until now remained overlooked, such as the length of the tentacle, the width of the head and slight differences in the sexual modifications of arm II in mature males. In general, it seems that the modifications of the arm II in Sl. cthulhui sp. nov. are less pronounced than in Sl. leucoptera. Even so, the differences between both species are few, likely representing a shallow morphological drift and a recent speciation phenomenon. Species involved in those first steps are commonly referred as in the "gray speciation zone" (Roux et al., 2016). As we showed with S. affinis, S. atlantica and S. intermedia, as well as the molecular data provided for the descriptions of Sl. japonica and Sp. nipponensis, these situations are likely to be frequent in bobtail squids. Based on this observation, the specific status of distant populations with mild levels of genetic divergence, which are typically identified as intraspecific levels of divergence, should be taken with care.
We took special care of including as many morphological characters as possible to avoid overlooking possible morphological differences between the congeneric Sl. cthulhui sp. nov. and Sl. leucoptera. That is the case of beak and radulae morphologies, which are rarely used in bobtail squid taxonomy (e.g., Kubodera and Okutani, 2011;Sanchez et al., 2019), and spermatophores. The use of spermatophores in cephalopod taxonomy should be taken with caution, as our study also shows. It is known that spermatophore size depends on the size of the male that produces it and since they can accumulate spermatophores for long periods of time, the same individual can hold spermatophores of a huge range of sizes and morphologies, according with the somatic size of the squid when each spermatophore was formed (Hoving et al., 2010;Cuccu et al., 2014). The number of spermatophores stored by the male would vary according with when it mated the last time (if anything at all) and so will do the size of those structures. In the present study we found 353-358 spermatophores in the spermatophoric sacs of Sl. leucoptera, while 236 (ICMC000165) and 11,645 (ICMC000166) spermatophores were found in the two examined specimens of Sl. cthulhui sp. nov. It is interesting to point the fact that the specimen ICMC000165, with a lower number of spermatophores, was also the same individual with a narrowest range of spermatophore size. Particularly, the specimen ICMC000166 not only showed the largest range of spermatophore size: it also showed empty and aberrant spermatophores. It is known that in the beginning of their reproductive life, cephalopod males produce tentative spermatophores (Nigmatullin et al., 2003) and in this case both functional and tentative spermatophores were present in the spermatophoric sac. It is not known to us if this was related with the absence of mating events by this individual for any circumstance or if it is just an aberrant individual. Pelagic cephalopods might have difficulties to find conspecifics to mate (Hoving et al., 2012), or just the opposite situation (Fernández-Álvarez et al., 2018;Hoving et al., 2019), so species with promiscuous mating systems are expected to hold a less variable morphology in the spermatophores stored in the spermatophoric sac. Special attention should be taken to the characteristics of a cephalopod species sex life if spermatophores are going to be used as a taxonomic character.

CONCLUSIONS
In this study we provided new sequence data for most of the Mediterranean bobtail squid species and added other relevant sequence data from species from other parts of the world. The combined use of several molecular species delimitation methods consistently showed a picture of hidden biodiversity. On the one hand, most of those methods failed to accurately assess the actual biodiversity of some morphologically different species, hindering the use of molecular data for species identification in the absence of morphological data. On the other hand, cryptic biodiversity was detected among members of the same nominal species, pointing to the fact that in some cases, speciation phenomena might be occurring without a parallel morphological evolution. It is also possible that the comparatively low number of taxonomists working on bobtail squids has hindered the discovery of morphological differences between them or that some morphological differences had been overlooked, as it happened with the new Stoloteuthis species described in this work. The Mediterranean Sea is one of the more diverse and the better studied areas for members of the family Sepiolidae, with a literature production spanning through more than a century (Bello, 2015). Despite this intense biodiversity and taxonomic research, two new Mediterranean species were described in the last few years (Bello, 2013;Bello and Salman, 2015). Here, a new species of Mediterranean Stoloteuthis, previously misidentified, was discovered and described, both using molecular and morphological methods. It is also remarkably that Olmos-Pérez et al. (2018b) found a species ("Heteroteuthinidae" sp.) in European waters whose sequences cannot be assigned to any sequenced species of bobtail squid, while Sanchez et al. (2019) discovered two new species of the genus Euprymna in Pacific waters. All these recent new species discoveries stress the need of new taxonomic studies in both benthic and pelagic bobtail squids on a worldwide basis.

DATA AVAILABILITY STATEMENT
The datasets presented in this study can be found in online repositories. The names of the repository/repositories and accession number(s) can be found at: NCBI GenBank (accession: MW260131-MW260134).

ETHICS STATEMENT
Ethical review and approval was not required for the animal study because the cephalopods we worked with for this study were all dead before research began.

AUTHOR CONTRIBUTIONS
FÁF-Á, PS, and RV conceived the study, participate in sample collection and identification, and performed morphological examinations. FÁF-Á performed the phylogenetic analyses, data curation, and prepared the first draft. All authors contributed to the manuscript preparation and revision and read and approved the submitted version.