Morphology, nuclear SNPs and mate selection reveal that COI barcoding overestimates species diversity in a Mediterranean freshwater amphipod by an order of magnitude

DNA sequence information has revealed many morphologically cryptic species worldwide. For animals, DNA‐based assessments of species diversity usually rely on the mitochondrial cytochrome c oxidase subunit I (COI) gene. However, a growing amount of evidence indicate that mitochondrial markers alone can lead to misleading species diversity estimates due to mito‐nuclear discordance. Therefore, reports of putative species based solely on mitochondrial DNA should be verified by other methods, especially in cases where COI sequences are identical for different morphospecies or where divergence within the same morphospecies is high. Freshwater amphipods are particularly interesting in this context because numerous putative cryptic species have been reported. Here, we investigated the species status of the numerous mitochondrial molecular operational taxonomic units (MOTUs) found within Echinogammarus sicilianus. We used an integrative approach combining DNA barcoding with mate selection observations, detailed morphometrics and genome‐wide double digest restriction site‐associated DNA sequencing (ddRAD‐seq). Within a relatively small sampling area, we detected twelve COI MOTUs (divergence = 1.8–20.3%), co‐occurring in syntopy at two‐thirds of the investigated sites. We found that pair formation was random and there was extensive nuclear gene flow among the ten MOTUs co‐occurring within the same river stretch. The four most common MOTUs were also indistinguishable with respect to functional morphology. Therefore, the evidence best fits the hypothesis of a single, yet genetically diverse, species within the main river system. The only two MOTUs sampled outside the focal area were genetically distinct at the nuclear level and may represent distinct species. Our study reveals that COI‐based species delimitation can significantly overestimate species diversity, highlighting the importance of integrative taxonomy for species validation, especially in hyperdiverse complexes with syntopically occurring mitochondrial MOTUs.


Introduction
In times of biodiversity loss, assessing and monitoring species diversity as well as detecting new species remains of utmost importance. Traditionally, species identification and description were based on morphological characters. However, sometimes speciation does not lead to morphological differentiation (Da€ ınou et al., 2016), making it challenging or impossible to identify diagnostic differences (Mutanen and Pretorius, 2007) and worsening the effects of taxonomic impediments (Coleman, 2015;Rodman and Cody, 2003). One of the solutions proposed is to integrate molecular information in the process of species identification and delineation (DNA taxonomy). In particular, DNA barcoding is a commonly used approach to identify known species based on standard reference gene sequences and to detect potentially new species through deviations from known sequences. The principle is based on the assumption that intraspecific variation is smaller than interspecific variation and that species have unique barcode gene sequences that are not shared with other species (Godfray, 2007;Hebert et al., 2003Hebert et al., , 2010Hebert and Gregory, 2005). For animals, the standard barcoding marker is a fragment of the mitochondrial cytochrome c oxidase subunit I (COI) gene (Hebert et al., 2003). Using COI for DNA taxonomy has shown great promise for speeding up the process of species discovery (Brower, 2010;Meierotto et al., 2019;Sharkey et al., 2021).
Exploratory DNA barcoding studies revealed the presence of a great number of previously overlooked morphologically cryptic species (Fi ser et al., 2018;Johnson et al., 2008;Kane et al., 2008;Katouzian et al., 2016;Nakano and Spencer, 2007;Pfenninger and Schwenk, 2007;Weiss et al., 2014). Furthermore, an increasing number of studies focusing on the ecology of cryptic species have revealed that there can be significant differences in their ecological demands as well as their resilience against stressors (Feckler et al., 2012;Fi ser et al., 2018;Macher et al., 2016), underlining the importance of correct species identification (Bickford et al., 2007;Fi ser et al., 2018;Pfenninger and Schwenk, 2007). For a long time, DNA-based species delineation relied only on single genetic markers, which can mislead species diversity estimates (Flot, 2015;Fontaneto et al., 2015;Padial et al., 2010;Pante et al., 2015). This holds true in particular for mitochondrial markers, such as COI, that have several unique characteristics (e.g. no recombination, haploid genome, usually only maternally inherited, high mutation rates; see Ballard and Whitlock, 2004). Hence, analysing nuclear in addition to mitochondrial markers has increasingly led to the detection of different conflicting patterns between both marker systemsthe so-called mito-nuclear discordance (Petit and Excoffier, 2009;Toews and Brelsford, 2012). One common discordance pattern is that mitochondrial data do not support distinct clusters whereas nuclear data show separately evolving lineages. This pattern can emerge when speciation has happened relatively recently or is still ongoing (Ivanov et al., 2018;Salokannel et al., 2021;Weiss et al., 2018). Two other discordance patterns are the result of differential gene flow for the mitochondrial and nuclear genomes. In one case, distinct mitochondrial lineages are found (i.e. genetic distances between clusters resemble interspecific distances), whereas nuclear data show no differentiation, indicating ongoing gene flow. This pattern can emerge if populations of a species evolved in geographical separation and then came into secondary contact. Ongoing gene flow eliminates nuclear differentiation, but mitochondrial differences remain due to the lack of recombination (Elbrecht et al., 2014;Hinojosa et al., 2019;Lee et al., 2021;Thielsch et al., 2017). In the other case, the presence of distinct species is indicated by both marker types, but the observed divergence patterns are not congruent between them. This pattern can emerge if distinct species exchange genes through hybridization or introgression (Marshall et al., 2021;Papakostas et al., 2016;Weigand et al., 2017).
Even though these issues have been known for some time now, recent reviews indicate that more than 30% of all studies suggesting presence of cryptic species are based on information derived from a single, often mitochondrial, genetic marker with very few studies using genome-wide evidence (Fi ser et al., 2018;Struck et al., 2018). Furthermore, identification of cryptic species using multiple molecular markers rarely is accompanied by high-resolution morphometric analysis. However, when these analyses are conducted, some of the cryptic species turned out to exhibit consistent morphological differences, so they were in fact not morphologically cryptic (Karanovic et al., 2016;King et al., 2012;Lajus et al., 2015). Therefore, detailed morphometric analyses should follow molecular species delimitation.
Amphipod crustaceans are one taxon where cryptic species diversity has been reported numerous times (Fi ser et al., 2018). Overlooked cryptic diversity in amphipods is a widespread global phenomenon observed in communities inhabiting deep-sea habitats, lowland rivers and high-mountain streams (Baird et al., 2011;Copilas ß-Ciocianu and Petrusek, 2015;Grabowski et al., 2017;Havermans et al., 2013;Hupało et al., 2019;Katouzian et al., 2016;Mamos et al., 2014Mamos et al., , 2021Wattier et al., 2020). This includes the Mediterranean freshwater amphipod Echinogammarus sicilianus Karaman & Tibaldi, 1973. E. sicilianus was originally described from the southeastern part of Sicily and later recognized as widely distributed across the island's freshwater habitats. It was also found in isolated insular populations on Sardinia and Malta (Pinkster, 1993). Echinogammarus sicilianus is considered a species adapted to living in a semi-arid climate. It is assumed to actively disperse during the wet period and then retract to the headwaters during drier summer months (Pinkster, 1993). Recently, extraordinary levels of intraspecific genetic diversity have been discovered within E. sicilianus, revealing formerly overlooked potential cryptic diversity. Based mostly on mitochondrial data, up to fifteen potential cryptic species have been identified on Sicily (Hupało, 2019). Four of them co-occurred in a single river system (i.e. Fosso del Tempio) in syntopy. Since some cryptic species are known to coexist (Fi ser et al., 2018), we wanted to verify if the observed diversity represents a case of cryptic species coexistence. In our study, we used the unified species concept according to De Queiroz (2007), where species are defined as separately evolving metapopulation lineages.
The aim of this study was to investigate the species status of the mitochondrial lineages found in the Fosso del Tempio river system. For that, we performed an extensive fine-scale sampling throughout the entire catchment, as well as in the neighbouring catchments. We then used an integrative taxonomic approach combining DNA barcoding of over 1500 specimens with high-resolution morphometric data for the most common mitochondrial MOTUs, as well as analysing selected specimens of nearly all MOTUs with genomewide double digest restriction site-associated DNA sequencing (ddRAD-seq) data. RAD-seq has been applied in population genomics, phylogeography and phylogenomics (Combosch and Vollmer, 2015;Ortiz et al., 2021;Reitzel et al., 2013;Weiss et al., 2022), as well as species delimitation (Ivanov et al., 2018;Weigand et al., 2017;Weiss et al., 2018), where its main advantage is providing information from multiple loci across the entire genome compared to commonly used single-locus nuclear markers (Andrews et al., 2016;Baird et al., 2008). Furthermore, as correct mate recognition is one of the crucial aspects characterising a distinct species, we also investigated the MOTU composition of precopulatory mating pairs. Precopulatory mate guarding preceding the actual copulation is common in numerous amphipod species (Conlan, 1991;Sutcliffe, 1992). This behaviour allows for easy observation of natural mating patterns and thus may indicate the presence of prezygotic barriers and support species status for amphipods (Lagrue et al., 2014;Sutherland et al., 2010). With that, we aimed to assess whether behavioural, morphological and genomic data support species delineation based on COI data.

Sampling
The samples were collected during two seasonal sampling campaigns conducted in September 2019 (autumn) and January 2020 (winter) mainly in the Fosso del Tempio river system in Sicily. The studied river belongs to the basin of Fiume dei Monaci, which covers ca. 590 km 2 and belongs to the Fiume Simeto catchment. The river originates from the slopes of Monte Moliano e Montagna on the border of the territory of the Municipalities of Aidone and Piazza Armerina. The main stretch of Fosso del Tempio is ca. 30 km long and after the first main confluence, it changes its name to Fosso Pietrarossa.
In order to analyse the fine-scale diversity and distributional patterns of E. sicilianus, 25 sites were chosen in the Fosso del Tempio river system, on average 1 km apart (Fig. 1). To more broadly assess regional diversity, an additional five sites from neighbouring catchments were sampled. Although a total of 30 sites were visited, individuals of E. sicilianus were found only at 25 of them, where they co-occurred with E. adipatus (except for sites E9, A27 and D13, where only individuals of E. sicilianus were found). All sites, except for A25 and A27 (sampled only in winter) were visited in both seasons (Table S1). The individuals were collected from all available stream microhabitats (up to five per site) using a hand net, aiming to collect at least 20 individuals per microhabitat per site. Additionally, up to ten precopulatory pairs per site were collected, if available (Table S1). The collected individuals were stored in 96% ethanol until further processing.

DNA barcoding
DNA was extracted using a 10% solution of Chelex 100 resin (Bio-Rad Laboratories, Hercules, CA, USA) and the following protocol: a small piece of the amphipod's muscle tissue from the dorsal side and/or two to three pleopods were placed in one well of a PCRplate with 100 lL of Chelex solution (10%) and incubated at 95°C for 20 min in a thermocycler, vortexed thoroughly every 5 min. Afterwards, the samples were centrifuged at 13 000 g for 1 min and stored at 4°C.
For DNA barcoding, we amplified the standard barcoding fragment of the COI gene using a 20-lL assay with 10 lL of Dream-Taq TM Hot Start Green PCR Master Mix (Thermo Fisher Scientific, Waltham, MA, USA), 1.6 lL (5 lM) of each primer (LCO1490-JJ and HCO2198-JJ; Astrin and St€ uben, 2008), 4.8 lL of nuclease-free water and 2 lL of DNA template per reaction. All PCR plates were prepared using a Biomek FX P liquid handling workstation (Beckmann Coulter Inc., Brea, CA, USA) and subsequently placed in a thermocycler for amplification. PCR settings for COI amplification were as follows: initial denaturation at 94°C for 60 s; five cycles of denaturation at 94°C for 30 s, annealing at 45°C for 90 s and extension at 72°C for 60 s, followed by 35 cycles with 94°C for 30 s, 51°C for 90 s and 72°C for 60 s, with a final extension at 72°C for 5 min. Before sequencing, 5 lL of the PCR products were purified enzymatically with 10 U of Exonuclease I (ExoI) and 1 U of thermosensitive alkaline phosphatase (FastAP) (both Thermo Fisher Scientific) with an incubation step at 37°C for 25 min followed by an inactivation step at 85°C for 15 min. Sanger sequencing of the forward direction (using the LCO1490-JJ primer) was performed at Eurofins Genomics (Cologne, Germany). For individuals where the sequence quality was low, a reamplification and (bidirectional) sequencing was performed using the same protocol (individuals indicated in Table S2). All generated COI sequences were deposited in GenBank (OP684936-OP686437). Additionally, all sequences were compiled in the dataset and deposited in the public repository of the Barcode of Life Data Systems (BOLD; Ratnasingham and Hebert, 2007), where all of the relevant metadata information and sequence trace files are available (project code: SIZES).

ddRAD-seq libraries
Samples for ddRAD-seq were selected depending on COI MOTU and location. Initially, 20 individuals of each of the four most common COI MOTUs (i.e. S01-S04) were selected; additionally to these 80 specimens, 16 individuals from the other MOTUs (i.e. S05-S12, except for S12) were included adding up to 96 specimens in total. After analysing the first batch of individuals, we decided to increase the sampling size especially at sites where MOTUs were found in syntopy as well as for rare MOTUs. Because higher amounts of high-quality DNA were needed for ddRAD-seq in comparison to COI barcoding, extraction from either muscle tissue and all pereiopods or entire specimens (extraction method for each individual indicated in Table S2) was conducted using a magnetic bead or column-based extraction protocol based on the NucleoMagÒ Tissue (Macherey-Nagel; Thermo Fisher Scientific) extraction protocol.
For specimens of the first library (Table S3), tissue was lysed with 195 lL of lysis buffer (T1) and 5 lL of Proteinase K (250 ng/lL working concentration), followed by an incubation at 55°C for 3 h. Afterwards, 5 lL of RNAse A (10 mg/mL; Thermo Fisher Scientific) was added and the samples were incubated for 30 min at 37°C. Then, DNA was extracted in a 96-well plate format using the Biomek FX P (Beckmann Coulter Inc.), following the protocol as described in Buchner et al. (2021). For the full-body extracts for libraries 2-4, the specimens were first homogenized using 1-and 2mm zirconia beating beads in a solution of 477 lL TNES buffer, 20 lL Proteinase K and 3 lL RNAse A (for Library 4 the amounts of TNES, Proteinase K and RNAse A were doubled) and then incubated at 65°C for 20 min, while shaking at 1400 r.p.m. Afterwards, DNA was extracted using the same extraction protocol using selfmade buffers (Appendix S1) and beads (Oberacker et al., 2019) for specimens belonging to the ddRAD-seq libraries 2 and 3 (Table S3). For specimens included in library 4, the same buffers, but a column extraction plate was used to obtain a higher DNA yield. After the extraction, DNA concentration was measured either with a Qubit dsDNA BR Assay Kit (Life Technologies; Thermo Fisher Scientific) or a SpectraMax ID3 plate reader (QuantIT High Sensitivity Kit; Molecular Devices, San Jose, CA, USA).
Depending on the initial DNA concentration (Table S3), up to 600 ng of DNA were used for double digestion with the FastDigest restriction enzymes Csp6I and SdaI (Thermo Fisher Scientific). Reactions were carried out in either a 50-lL or 30-lL reaction. For each reaction, 1 lL of each enzyme and 10% FD-Buffer (5 or 3 lL) were added to a DNA/nuclease-free water-Mix (Table S3) and incubated for 30 min at 37°C. Directly after digestion, reactions were purified using the NucleoMag NGS Clean-up and Size select kit (DNA clean-up protocol; Macherey-Nagel). This and all of the following clean-up and size selection steps were performed on the Biomek FX P (Beckmann Coulter Inc.). To calculate adapter quantities for ligation (Peterson et al., 2012), expected cut frequencies were estimated by in silico digestion using the script genomecut.pl (Rozenberg, https:// github.com/evoeco/radtools/). These frequencies were based on a standard shotgun genomic library of Gammarus fossarum (Macher et al., 2017) and resulted in average cut frequencies of 279 bp for Csp6I and 178 891 bp for SdaI. Adapter sequences are given in Table S4 and are similar to adapters used in Weigand et al. (2017). P7 adapters contain a degenerate base region (DBR) used to detect and remove PCR duplicates (Schweyen et al., 2014). The reaction was carried out using 22 lL of digested DNA, 0.5 lL of T4 DNA Ligase (2000 U/lL; New England Biolabs, Ipswich, MA, USA), 4 lL of T4-Ligase Buffer, calculated amounts of P5 and P7 adapters (Table S3), filled up to 40.5 lL with nuclease-free water. After mixing, the reaction was incubated for 2 h at room temperature and inactivated by heating (10 min at 65°C). Adapter dimers were removed according to the respective protocol from the NucleoMag NGS Clean-up and Size select kit (purification, 1.0 ratio, followed by size selection, 0.9 ratio).
After clean-up, PCR reactions were conducted using 0.02 U/lL of Q5 High-Fidelity DNA Polymerase (New England Biolabs), 19 Q5 Reaction Buffer, 200 lM of dNTPs, 1 lM of each primer (P5 and P7; Schweyen et al., 2014) and 1 lL template DNA, made up to 25 lL with nuclease-free water. The thermocycling conditions used were as follows: initial denaturation at 98°C for 30 s, 20 cycles of 98°C for 10 s, 65°C for 30 s and 72°C for 30 s, and final elongation for 5 min at 72°C. If PCR amplification was not successful, reactions were repeated using more or less DNA template, depending on the initial DNA quantity. PCR products were purified with NucleoMag NGS Beads and a ratio of 1.8, followed by a double-sided size selection with ratios of 0.55 (right side) and 0.85 (left side) and eluted in 30 lL of EB buffer. The success of amplification and size selection as well as the quantity of each sample was measured on a Fragment Analyzer using the High Sensitivity NGS Fragment Analysis Kit (both Advanced Analytical; Agilent, Santa Clara, CA, USA).
After measurement, samples were pooled equimolarly into four libraries (Table S3) and a final size selection targeting a size range of 300-500 bp for sequencing was carried out using NucleoMag NGS Beads with the ratios 0.67 (right) and 0.75 ratio (left). After final quality control (Fragment Analyzer) and DNA concentration measurement (Qubit dsDNA HS Assay Kit; Thermo Fisher Scientific), libraries were sent for sequencing to Macrogen Korea and pairedend sequenced on four HiSeq X Illumina lanes (read length 2 9 150 bp). In addition to the samples of E. sicilianus, a varying number of samples of E. adipatus was included on each lane, resulting in sequencing of 96, 72, 72 and 96 specimens (80, 64, 64 and 18 E. sicilianus) on the different lanes. Data for E. adipatus were not analysed and discussed in this study.

COI data analysis
The quality of the obtained COI sequences was checked in Geneious Prime 2022.0.2 (https://www.geneious.com) and aligned with MAFFT 1.4.0 (Katoh and Standley, 2013) as implemented in Geneious using the automatic algorithm selection and default settings. Sequences which were too short or where quality was insufficient, were excluded from the dataset. The alignment of the remaining sequences was trimmed to the shortest sequence. MOTU delimitation was done using ASAP (assemble species by automatic partitioning; Puillandre et al., 2021), a new implementation of ABGD (automatic barcode gap discovery; Puillandre et al., 2012) considered to be the most conservative single-locus, distance-based species delimitation method. The algorithm is supposed to be less prone to MOTU oversplitting, often observed in tree-based methods (Dellicour and Flot, 2018), and more reproducible than BINs (barcode index numbers; Meier et al., 2022). For amphipods, MOTUs delimited with ASAP had the highest congruence with morphology in comparison to other methods (Copilas ß-Ciocianu et al., 2022). COI groups were defined by using the ASAP web server applying the Kimura (K80 ts/tv 2.0) model for computing distances and otherwise default settings. To visualise phylogenetic relationships between MOTUs, a neighbour-joining tree was generated in MEGA 11 (Tamura et al., 2021). Branch support was estimated using 1000 bootstrap replicates. The Kimura 2-parameter (K2P) substitution model was used to correct distances based on the maximumlikelihood (ML) model test performed in MEGA. For rooting the tree, a single sequence of the closely related E. adipatus co-occurring throughout the studied river system with E. sicilianus, was used (GenBank accession number: OP684936). To also visualise distances between haplotypes within delimited ASAP groups, a minimumspanning network (Bandelt et al., 1999) was created in PopART v.1.7 (http://popart.otago.ac.nz) and coloured according to ASAP groups (MOTUs). Group definitions of ASAP were used to calculate distances between and within groups in MEGA 11 using the K2P model and uniform rates.

Detection of intermixing in precopulatory pairs
For sites where multiple MOTUs were detected, we checked precopulatory pairs for partner MOTU assignment and calculated the frequency of intermixed pairs (i.e. where male and female belong to different MOTUs). To test if the proportion of observed precopulatory pairs and intermixing rate align with expectations of nonassortative mating, chi-square goodness-of-fit tests were performed. First, to test if precopulatory pairs formed randomly (i.e. irrespective of the MOTU affiliation), the MOTU composition of precopulatory pairs was tested against the expected frequencies deriving from relative abundances of MOTUs observed on the studied sites where precopulatory pairs were collected. Secondly, observed MOTU frequencies forming intermixed pairs were tested against the MOTU frequencies observed in all collected precopulatory pairs from sites where intermixing was observed (12 sites). All statistical analyses were calculated in R v.3.3.2 (R Core Team, 2015).

ddRAD-seq data analysis
Pre-processing of the four ddRAD-seq libraries was conducted as described in Weiss et al. (2018). Subsequently, denovo_map.pl of Stacks v.1.34 (Catchen et al., 2013) was used to identify loci and genotypes. Raw data as well as demultiplexed and pre-processed data for all specimens are deposited in the NCBI sequence read Archive (SRA) in the Bioproject PRJNA894373. Individual accession numbers are annotated in Table S3. To be able to identify the optimal parameter settings for Stacks, the Stacks pipeline was run with eight different parameter settings according to the guidelines in Paris et al. (2017). The parameter m (minimum number of raw reads required to form a stack or putative allele) was fixed at m = 3. The parameter M (number of mismatches allowed between stacks for merging into a putative locus) varied between 2 and 5 and N (number of mismatches allowed to align secondary reads to putative loci) was set to M + 2. The number of mismatches allowed between stacks during catalogue constructions (n) was equal to M or M + 1. Results from the Stacks database were exported with export_sql.pl using a minimum stack depth of 3.
Further analyses were conducted using the workflow management tool Snakemake (K€ oster and Rahmann, 2012) combining stacks2fasta.pl (Macher et al., 2015) and several python and R scripts for data reformatting, filtering and analyses similar to the workflow used in Weiss et al. (2018). The following general filtering settings were chosen for the analyses: maximum number of single nucleotide polymorphisms (SNPs) per locus was 1-12, of which only one was further used; minor allele frequency was 1% and minimum percentage of individuals required for a locus to remain in the dataset was varied between 60%, 70%, 80% and 90%. Individuals with >50% missing data were excluded from the dataset for the final analyses. To evaluate differences between the Stacks parameter and locus filtering settings, basic population statistics, such as observed heterozygosity (H O ), observed gene diversity (H S ), overall gene diversity (H T ) and overall F ST and F IS were calculated using hierfstat (Goudet, 2005) in R v.3.3.2 (R Core Team, 2015). Besides, neighbour-net networks (Bryant and Moulton, 2004) were calculated using SplitsTree v.4.14.5 (Huson and Bryant, 2006). Individuals were coloured according to COI group affiliation using Adobe Illustrator CS6. To analyse population structure in more detail, individual ancestry coefficients were estimated based on sparse nonnegative matrix factorization algorithms (sNMF; Frichot et al., 2014) with the R-package LEA (Frichot and Franc ßois, 2015). For sNMF analysis, one to 17 clusters and 30 replicates with 100 000 iterations per replicate were used. To select the best replicate and the most probable number of clusters (K), cross-entropy values were compared.

Morphology
We examined morphological differentiation by measuring 49 traits that reflect diet, locomotion and reproduction (Copilaș-Ciocianu et al., 2021;Fi ser et al., 2013;Kralj-Fi ser et al., 2020;Premate et al., 2021) (Table S5). A complete description of landmarks and function of the examined traits is provided in Copilaș-Ciocianu et al. (2021). In total, we used 80 adult individuals (originating from precopulatory pairs, COI barcoded before the morphological analyses) from the four most common mitochondrial MOTUs (10 males and 10 females per MOTU; Table S2). Ethanol-preserved specimens were soaked overnight in a 2% lactic acid (LegoChem Biosciences Inc., Daejeon, South Korea) solution with subsequent transfer to a 2:1 mixture of 70% ethanol and glycerol (Carl Roth Gmbh & Co. Kg, Karlsruhe, Germany) in order to soften and partially clear the cuticle. Animals were dissected in glycerol, under a Nikon SMZ1000 stereomicroscope, with the help of fine tweezers and microsurgical scissors. Due to asymmetry of some mouthparts, only the right-sided ones were measured. All of the other traits were measured from either the left or right side of the body, depending on specimen completeness. The dissected appendages were temporary mounted on microscope slides in glycerol and photographed with a Nikon DS-Fi2 camera attached to a Nikon Eclipse Ci-L microscope or a Nikon SMZ1000 stereomicroscope, using the manufacturer's software (NIS-Elements; Nikon, Tokyo, Japan). Morphometric measurements were made with the Digimizer software (https://www.digimizer.com/). Statistical analyses were run at the level of COI MOTUs. To account for the effect of body size, all of the measurements were regressed against body length and regression residuals were used in subsequent analyses. Meristic counts were log-transformed before the regression. Differentiation was assessed using a permutational multivariate analysis of variance (PERMANOVA) test with 9999 permutations and a Euclidean similarity index. Bonferroni correction was used for multiple comparisons. We first ran PERMANOVAs with sex as a factor, and subsequently, within each sex, with COI MOTUs as factors. To gain a deeper understanding of morphological differentiation, analyses also were run at the level of five functional complexes of traits related to: (1) sensorial function (antennae and eyes; eight measurements), (2) food gathering and handling (gnathopods; 10 measurements), (3) food processing and digestion (mouthparts and stomach; 13 measurements), (4) crawling (pereopods; five measurements), and (5) swimming, breathing and reproduction (ventral channel; 11 measurements). To visually inspect the patterns of differentiation, a principal component analysis (PCA) was performed on a correlation matrix based on the size-corrected data. A linear discriminant analysis (LDA) also was conducted in order to investigate to what extent the COI MOTUs can be discriminated against with morphometric data. All statistical analyses were performed with PAST 3 (Hammer et al., 2001). Raw morphometric measurements are available in Table S5

Assessment of putative cryptic species diversity -COI analysis
In total, we generated 1501 high-quality COI sequences for E. sicilianus from 25 sampling sites. The final 554-bp alignment contained 182 variable sites (32.9%) of which 28 were nonsynonymous substitutions. Using the distance-based approach ASAP to identify putative cryptic species, 2-140 groups were reported when considering the ten best partitions ( Table S6). The two partitions with the best ASAPscores indicated 12 putative species, which we used as species hypotheses for subsequent analyses and called MOTUs S01-S12 hereafter. The individual assignment is given in Table S2. Pairwise K2P distances between MOTUs ranged between 1.8% (S02 and S11) and 20.3% (S04 and S09) with a mean distance of 11% (Table S7). Maximum K2P distance within MOTUs was 1% in S04. The most common and diverse MOTUs were S01 (535 specimens; 29 haplotypes), S02 (210; 13), S03 (208; 21) and S04 (485; 54) ( Fig. 2a,  Fig. S1; Table S1). High diversity also was found in S07 with 12 haplotypes in 38 individuals. The other MOTUs each were represented only by two to eight individuals and one to three haplotypes. Differences between S02 and S11 (nine substitutions, 1.8% mean K2P) and between S01 and S10 (17 substitutions, 3.4% mean K2P) were relatively small, yet the other MOTUs differed to the next MOTU at ≥23 substitutions (S03-S12; 5% mean K2P). The most divergent MOTU was S09, which differed in ≥78 positions from S03 (mean K2P distances to other MOTUs between 16.8% and 20.3%). This MOTU was exclusively found in a tributary (A27) in the downstream part of the studied river system, where no other MOTU was detected (Fig. 2b). Similarly, S07 was found exclusively in a tributary (D13), which was close-by regarding air distance, but separated by a waterway distance of ca. 30 km from the next site (A25) in the main river system. At D6, also located in an adjacent river, one of the main MOTUs (S04) as well as a private MOTU (S06) was found. In general, co-occurrence of two or more MOTUs was found at 17 sites, whereas only one MOTU occurred at eight sites ( Fig. 2b; Table S1). When considering MOTU distribution in the main catchment, three geographical groups are visible, which are defined by one or two dominant MOTUs. The first group located in the upstream part of the main river consists of individuals sampled at sites E5, E8, E9, A7, A10, A11 and A12, where S01 was the most abundant MOTU. The second group, where S04 was dominant or exclusively found, consists of specimens from sites A15, A17 and A18 (also upstream part) as well as D7, D8 and D11 located in an adjacent tributary. The third group includes the more downstream sites (A27, A25, A1, A2, E1, E2, E3, E4, EF and A6F) where the highest MOTU diversity was found and all four main MOTUs, as well as most of the rare MOTUs, were detected. MOTUs S02 and S03 were the most dominant here, but only rarely detected further upstream. Another relatively abundant MOTU was S01, which was dominant in the upstream part.

Assessment of species status of COI MOTUs
Precopula data. In order to test if precopulatory pairs are formed primarily by specimens of the same MOTU, we analysed 129 pairs collected at sites where at least two MOTUs co-occurred (Table S8). The MOTU composition of precopulatory pairs was not significantly different from the relative abundance of MOTUs observed across the studied sites where precopulatory pairs were collected (v 2 = 1.8038, df = 8, P = 0.9864). In 50 pairs (39%), we observed precopulatory intermixing (i.e. male and female specimen did not belong to the same MOTU) (Fig. 3); the majority of intermixed pairs (80%) included individuals belonging to one of the three most common MOTUs S01, S02 and S03. Four of six rare MOTUs (S05, S10, S11 and S12), which co-occurred with other MOTUs also formed intermixed precopulatory pairs ( Fig. 3; Table S8). Most of the intermixed pairs (92%) were found at sites where three or more MOTUs were observed. In the intermixed pairs, the MOTU affiliation was not significantly different from frequency expectations based on all analysed precopulatory pairs (v 2 = 12.408, df = 7, P = 0.08792).
Morphology. Morphological analyses were performed at the level of the four most abundant COI MOTUs (S01-S04). PERMANOVA testing revealed a significant differentiation between males and females within each functional complex of traits, and when all traits were considered simultaneously (Table 1). However, within each sex, specimens from the four COI MOTUs were not significantly different when considering all traits simultaneously (Table 1; female MOTUs, F = 2.056, P = 0.794; male MOTUs, F = 1.817, P = 0.093). Further testing within each functional group revealed a similar pattern of weak to no differentiation. A significant differentiation was observed among female MOTUs with respect to length of pereopods (F = 2.503, P = 0.042). However, pairwise tests indicated that none of the comparisons were significant, although MOTUs S02 and S04 were close (F = 5.64, P = 0.065). Likewise, among male MOTUs, the only significant differentiation was observed with respect to mouthparts and stomach (F = 2.414, P = 0.0043). Here, only MOTUs S02 and S04 differed significantly (P = 0.027; Table 1). The difference was driven by the on average 13% longer stomach of MOTU S04 than that of MOTU S02 (F = 13.77, P = 0.008). The PCA supported the PERMANOVA results by revealing a large overlap in morphospace among MOTUs (Fig. 4). The first three PCA axes explained 36.9%, 9.9% and 6.04% of the variance, respectively. Likewise, the LDA performed very poorly at classifying specimens according to MOTUs, with an overall accuracy of 45%.  ddRAD-seq data. When comparing Stacks parameter and locus filtering settings, the number of obtained loci varied strongly between Stacks settings and filtering limits (Table S9). When 60% of the individuals needed to have data for the locus, 1477-2229 loci were obtained and global F ST values varied between 0.31 and 0.22. For the 90% limit, only 149-185 loci remained. As a compromise between missing data and number of loci, we decided to use the 80% limit for final analysis. When using the most stringent Stacks setting (m3 M2 N4 n2) resulting in the highest number of loci, especially individuals of MOTU S09 had many missing data (62-67%), which was reduced with more relaxed Stacks settings. As general results did not differ between Stacks settings, we decided to use the setting (m3 M5 N7 n6) that resulted in the lowest numbers of missing data per S09 individual (35-46%). After excluding individuals with >50% missing data, 219 individuals and 495 loci remained, which were shared by ≥80% of the individuals (Table S9).
In order to evaluate if nuclear data support the COI MOTU assignment, we first created a neighbourjoining network (Fig. 5a,b). With the exception of specimens assigned to the mitochondrial MOTU S09 (Fig. 5a), all analysed specimens clustered together in a stretched star-like network pattern, where genetic distances between individuals were much lower than those to specimens of S09. The only other correspondence with regard to COI MOTU assignment were specimens belonging to S07, which formed a distinct cluster exclusively containing all specimens of this MOTU (Fig. 5b). The same pattern is also visible in the sNMF analysis (Fig. 5c). Here, cross-entropy values indicated that six clusters best represent the population structure with two clusters corresponding to S07 (K5) and S09 (K6), respectively. For the other four clusters, no congruence was found with MOTU assignment. Individuals assigned to the same COI MOTU were divided into different nuclear clusters and also individuals clustered into K1-K4 belonged to different MOTUs. However, when sorting individuals according to sampling site (Fig. S2), a geographical clustering after stream sections similar to the described COI groups becomes visible. While specimens from the upstream COI group 1 mostly belong to K1, specimens from the second upstream group (COI group 2) belong mainly to K4. Specimens from the downstream sites had highest membership proportions to K3, whereas specimens from A6F and some further specimens from downstream sites were grouped into K2. However, there were some exceptions, where specimens were assigned to two or more clusters or did not cluster according to the sampling site. Given the strong divergence of S09, we also performed the analysis excluding S09 individuals to evaluate their effect on the observed patterns, but the results remained very similar (Fig. S3).

Validation of putative cryptic species diversity
Based on our comprehensive COI analysis of E. sicilianus, we detected an extraordinary level of intraspecific mitochondrial diversity within this morphospecies, suggesting that it may consist of 12 putative cryptic species. Comparably high intraspecific COI distances have been found in other amphipod species (Desiderato et al., 2019;Wattier et al., 2020;Weiss et al., 2014), leading to the proposition of many cryptic species complexes (Csap o et al., 2020;Mamos et al., 2016;Wattier et al., 2020). Therefore, a higher genetic distance threshold of 5% was proposed for delimitation of amphipod species (Costa et al., 2009). In E. sicilianus, we found a broad range of distances between MOTUs, but only two (S10 and S11) exhibited genetic distances <5%. Hence, ten of the MOTUs could be regarded as putative cryptic species. However, in view of reported discordances between the nuclear and mitochondrial genomes in other taxa (Elbrecht et al., 2014;Hinojosa et al., 2019;Ivanov et al., 2018) as well as conflict between DNA and morphology (Copilas ß-Ciocianu et al., 2022) in other arthropod species, we consider the integration of multiple lines of evidence essential for delineating species. Nevertheless, we agree with recent propositions that such data-intensive approaches are not feasible in tackling "dark taxa" (groups containing thousands of species of which <10% are described; Hartop et al., 2022) and should be applied mainly in cases of unusual COI distances within or between morphospecies. In order to check for evidence of reproductive isolation between co-occurring MOTUs, we analysed the degree of intermixing in precopulatory pairs. For amphipods, it has been shown that closely related (cryptic) species rarely form intermixed precopulatory pairs under laboratory conditions as well as when naturally occurring in sympatry (Byst rick y et al., 2022;Lagrue et al., 2014;Sutherland et al., 2010). By contrast, we found a high rate of intermixing (close to 39%) between E. sicilianus COI MOTUs. The formation of intermixed precopulatory pairs was not significantly different from random, indicating a lack of prezygotic reproductive barrier. Furthermore, detailed morphometric analyses of nearly 50 traits revealed almost no significant morphological differences between the four most common MOTUs. We observed consistent morphological differentiation only between males and females, which is to be expected considering the pronounced sexual dimorphism in gammarid amphipods (Pinkster, 1993). The observed overall morphological similarity either supports true morphological crypsis or indicates that all studied individuals belong to the same species. As has already been shown for other amphipods, morphological traits accurately reflect occupied trophic niche (Copilaș-Ciocianu et al., 2021;Fi ser et al., 2013;Kralj-Fi ser et al., 2020;Premate et al., 2021) and therefore morphological similarity is likely to indicate the same ecological niche and similar use of available resources (Bock and Von Wahlert, 1965;Ricklefs and Miles, 1994). However, when considering individual traits separately, there was one significant difference in the stomach length between the males of MOTU S02 and S04 (mean K2P = 6.5%). Interestingly, we also detected a slight difference in length of pereiopods in females of those two MOTUs. These findings might suggest subtle ecological differences between MOTUs S02 and S04, probably deriving from their distribution patterns within the system (S02 is present mostly downstream while S04 mostly upstream), but to confirm that more research is needed regarding, for example, the feeding habits of the studied populations.
In accordance with the behavioural and morphological, but in discordance with the mitochondrial data, the genome-wide nuclear data indicate ongoing gene flow between most of the COI MOTUs, especially at localities where they co-occur. In general, the number of shared loci differed quite substantially between Stacks settings and in particular individual filtering limits. Even though the overall F ST value always indicated strong differentiation in the overall dataset, the value also varied depending on used settings. The differences that we found could be due to either biological reasons (i.e. high differentiation between populations) or technical issues, or a combination of both factors. To obtain sufficient target DNA for sequencing, we had to run a relatively large number of PCR cycles (20), resulting in many PCR duplicates. As a consequence of the included DBRs (Schweyen et al., 2014), they could be filtered out before loci calling, but this reduced the number of reads per individual quite substantially and could have led to less overlap of loci between individuals. However, the deduced population structure as shown in the network and the sNMF analyses were very similar between datasets (data not shown), indicating the robustness of the inferred population structure.
Ten of twelve of the COI MOTUs were not supported by the ddRAD-seq data. Of these ten, four were previously known (i.e. MOTUs S01-S04; Hupało, 2019) and six were newly detected in the present study (i.e. S05, S06, S08, S10, S11, and S12). The fact that these ten MOTUs do not represent distinct species suggests extraordinary levels of intraspecific variability in COI. This has been observed in other studies, where divergent mitochondrial clades were not supported by nuclear evidence (Elbrecht et al., 2014;Lee et al., 2021;Thielsch et al., 2017). However, this pattern cannot be generalized because lineages of comparable mitochondrial divergence were found to occur in syntopy while maintaining nuclear isolation in the related Gammarus fossarum (Byst rick y et al., 2022). The remaining two MOTUs (S07 and S09) are supported as distinct by both mitochondrial and nuclear data. Both are known from isolated populations occurring outside the main river stretch and do not co-occur with other lineages. They may indeed represent distinct species, but the presence of a physical barrier makes it difficult to infer whether there also would be reproductive isolation once secondary contact with other MOTUs has been re-established. Although individuals of both MOTUs can be ascribed morphologically to E. sicilianus, further detailed morphometric measurements, similar to those used in this study for MOTUs S01-S04, are needed to further validate their species status.

Potential reasons for observed intraspecific diversity in COI
Our study provides another example of local COI hyperdiversity in animals with limited dispersal abilities, which already has been observed in several other amphipod species (Copilas ß-Ciocianu and Petrusek, 2015; Grabowski et al., 2017;Mamos et al., 2021;Weiss et al., 2014;Wysocka et al., 2013) as well as other invertebrates (Elbrecht et al., 2014;T€ anzler et al., 2016;Weiss et al., 2018). In those cases, the speciation events leading to currently observed diversity were often caused by the historical isolation of individual populations. Although the observed COI variability of E. sicilianus might be the result of a single or multiple divergence or speciation events, the nuclear genomic diversity patterns are likely to point to hybridization between formerly isolated populations (Seehausen et al., 2008). The numerous COI substitutions separating MOTUs indicate an ancient divergence that has taken place throughout the late Miocene to the Pleistocene (Hupało, 2019). It is possible that each MOTU diverged allopatrically, in formerly isolated parts of the basin, which subsequently reconnected, allowing for secondary contact and gene flow between the previously isolated populations that had not yet become reproductively isolated. That also would explain why certain parts of the Fiume del Tempio system are dominated by a particular MOTU, possibly indicating the location of these previously isolated river sectors (e.g. S02 dominating on site A6F; S01 as only MOTU on site A12; S04 dominance on sites A15-A18 and D7-D11).
It is hard to evaluate how formerly isolated parts of the streams became connected and when the secondary contact occurred. However, given that so many mitochondrial lineages were retained and not eliminated by genetic drift, one may assume that the secondary contact occurred recently. For example, it could have been a consequence of alterations of the hydrological network of Fiume del Tempio, such as river capture events (Bishop, 1995) that are known to play an important role in the diversification the freshwater fauna (Burridge et al., 2006;Waters et al., 2020). Because no regional historical data are available, one can only speculate that the observed patterns are a result of Sicily's turbulent geological history rather than recent anthropogenic disturbances to the hydrological network (Broquet, 2016;Guglielmo and Marra, 2011). However, the genetic divergence patterns uncovered at small spatial scales suggest the presence of genetically isolated populations inhabiting certain stream sections. This also is underlined by the slight morphological differences we found between S02 (downstream) and S04 (upstream). However, although a physical barrier (concrete weir) between sites E8 and E5 has been shown to partially disconnect the faunal communities (Hupało et al., 2022), this barrier does not seem to fully explain the genetic structure we detected there. Therefore, it would be interesting to analyse the population structure in greater detail to identify potential, also biological, barriers shaping this strong geographical pattern, especially considering that sites were often only around 1 km distant from each other.

Conclusions
In our study, we showed that taxonomic identification based on COI barcoding alone would lead to a great overestimation of species diversity in otherwise morphologically uniform populations of E. sicilianus. Rather than proposing twelve putative cryptic species, we show that ten of them most likely represent the same species since they randomly mate and are morphologically indistinguishable. In view of our findings as well as other studies identifying mito-nuclear discordances, we showed that one should be extremely cautious when describing morphologically cryptic species solely using COI data. Instead, species descriptions should include additional lines of evidence. This holds true especially for hyperdiverse yet morphologically conserved similar species, where secondary contact might have led to gene flow between formerly isolated entities. However, the mechanisms leading to the observed restricted gene flow between populations at small scales still remain unanswered, and future research should focus on identifying further biotic and abiotic factors hindering gene flow between populations.

Acknowledgements
We would like to thank Saskia Schmidt, Marie-Th er ese Werner and Pedro M. G. Gomes for their valuable help during the sampling campaigns. We further want to thank Dominik Buchner for help with laboratory work and as well as the entire Leese Lab working group for helpful comments that further improved our manuscript. We also thank Bernd Sures, Daniel Grabner and Diego Fontaneto for the logistic support, and Florian Altermatt, Michał Grabowski, Alexander Weigand and Hannah Weigand for helpful discussions throughout the project. We also would like to thank two anonymous reviewers and Brent Emerson for providing valuable comments and remarks that helped to substantially improve the manuscript. This research was funded by German Research Foundation (DFG-project LE 2323/10-1). Open Access funding enabled and organized by Projekt DEAL.

Conflict of interest
None declared.

DATA AVAILABILITY STATEMENT
All generated COI sequences were deposited in GenBank (OP684936-OP686437). Additionally, all sequences were compiled in the dataset and deposited in the public repository of the Barcode of Life Data Systems (BOLD, Ratnasingham and Hebert, 2007), where all the relevant metadata information and sequence trace files are available (project code: SIZES, https://dx. doi.org/10.5883/DS-SIZESIC). Raw ddRAD-seq data as well as demultiplexed and pre-processed data for all specimens are deposited in the NCBI sequence read Archive (SRA) in the Bioproject PRJNA894373. Raw morphometric measurements are available on Figshare (doi: https://doi.org/10.6084/m9.figshare.20067404).

Supporting Information
Additional supporting information may be found online in the Supporting Information section at the end of the article.    Table S1. List of sampling localities (site code, coordinates, catchment) with information of number of individuals analysed with the different methods (COI, precopulatory pairs, morphometrics, ddRAD), number of MOTUs per site and the number of individuals per MOTU at each site. Table S2. List of all individuals including sample name, sampling site, sampling season, sequencing direction for COI, notes about the rerun of COI sequencing, MOTU assignment and GenBank accession number (COI), followed by indication of whether an individual was used for precopulatory, morphometric and/or ddRAD-seq analyses, ddRAD-seq library number and information on ddRAD-seq extraction method. Table S3. Information for ddRAD library preparation per sample (original sample name and name in ddRAD library) including initial DNA concentration, amount of DNA (ng) used for digestion, P5 and P7 adapter names and volumes, DNA concentration after PCR (in target size range of 300-500 bp and total DNA concentration), library number, used amount of DNA for library pooling (target and total amount), indication of whether an individual was used in the final analysis after filtering and NCBI Biosample accession number. Table S4. List of adapter and primer oligonucleotides in the 5 0 -3 0 direction. Table S5. List of all raw morphometric measurements (in mm). Table S6. Results for the 10 best partitions found by ASAP using the K80 (K2P) substitution model. Table S7. K2P genetic distances (COI dataset) calculated for all identified MOTUs. Table S8. List of all precopulatory pairs analysed with annotation if cross-lineage intermixing was detected. Table S9. Summary statistics for different Stacksand filtering settings; max. # snps = maximum allowed SNPS per locus, loci limit = percentage of specimens required to have the loci, ma = minor allele frequency, H O = observed heterozygosity, H S = within-population gene diversity, H T = overall gene diversity; F-statistics F ST and F IS calculated according to Weir & Cockerham (1984). Appendix S1. Recipes for self-made buffers used for DNA extraction.