In the present study, the datasets of two mitochondrial (cytochrome c oxidase subunit I [COI] and 16S rRNA) as well as one nuclear (28S rRNA) marker were analyzed. The material was obtained from published articles available by March 2021 (Ritchie et al., 2015; Jażdżewska and Mamos, 2019; Iguchi et al., 2020; Mohrbeck et al., 2021; Weston et al., 2021) supplemented by newly produced sequences. Figure 1 summarizes all datasets that have been used in this study, and information on new material collection, DNA extraction, amplification and sequencing is presented below.

All chromatograms of the newly obtained sequences were visually inspected, edited when reading mistakes happened in Geneious 10.1.2, and primer sequences were trimmed. These sequences were uploaded to GenBank under accession numbers: COI: MZ655819–MZ655889, 16S: MZ655914–MZ655967, 28S: MZ655890–MZ655910. Relevant voucher information, taxonomic classifications, and sequences of all studied genes (except for already published 16S sequences) are deposited in the dataset “DS-PCAPTEN” in the Barcode of Life Data System (BOLD)¹ (Ratnasingham and Hebert, 2007)². Additional summary that includes GenBank accession numbers for all sequences is available in the Supplementary Material S1.

Thirty specimens of Paralicella (15 of each of P. caperesca and P. tenuipes) were analyzed from samples collected at the Porcupine Abyssal Plain Sustained Observatory (PAP-SO), situated in the subpolar North-East Atlantic, at 48°50′N 16°30′W and a water depth of 4850 m. Samples were collected by means of baited trap (for trap details see Horton et al., 2020, and sampling details in Hartman, 2019; Ruhl et al., 2019). In the South Atlantic, 41 individuals were collected at one station in the Brazilian Basin during the DIVA-3 expedition. These individuals were collected using a baited trap (its description and sampling procedure is described in Martínez Arbizu et al., 2015).

In the case of North Atlantic amphipods, the total genomic DNA was extracted from one-two pleopods using a mixture of 150 μl pure H₂O with 0.015 g Chelex^® (SIGMA-ALDRICH Co.) and 10 μl proteinase K. The digestion at 55°C lasted for 6 h. The DNA of individuals from the Brazilian Basin were extracted using an AutoGenprep 965 extraction robot after overnight digestion at 55° in the AutoGen buffers with proteinase-K.

For North Atlantic individuals, the barcoding fragment of the COI gene amplification was conducted with the degenerated LCO1490-JJ and HCO2198-JJ primer pair (Table 2) and the reaction conditions described in Hou et al. (2007). Sequences were obtained by Macrogen Inc., the Netherlands on the Applied Biosystems 3730xl capillary sequencer. One-way (forward) sequencing was the standard procedure for all samples, but in addition, at least one individual of each recognized MOTU was sequenced in both directions.

The barcode fragment of the COI gene for the South Atlantic specimens of the studied species was amplified and sequenced with degenerate primers (dgLCO-1490/dgHCO-2198, Table 2) according to the protocol described in Riehl et al. (2014).

Sequences were initially blasted using default parameters on NCBI BLASTn and translated into amino acid sequences to confirm that no stop codons were present.

The 16S marker of 24 individuals (18 from the Central Pacific and six from North Atlantic, preferably one-two per each recognized BIN) was amplified using the primer pair 16SFt_amp/16SRt_amp2 (Table 2) in the conditions presented by Lörz et al. (2018b). In the case of the Central Pacific specimens, polymerase chain reaction was performed with AccuStart II PCR SuperMix (Quantabio), while for North Atlantic ones DreamTaq Green PCR Mastermix (Thermo Scientific) was used. Sequences were obtained from Macrogen Inc., the Netherlands, with the Applied Biosystems 3730xl capillary sequencer. Similarly to the COI gene, one-way (forward) sequencing was the standard procedure for all samples, but in addition, some individuals were sequenced in both directions. Additionally, 30 sequences of the specimens from the South Atlantic were successfully amplified and sequenced using the primer pair 16Sar/16Sbr (Table 2) at the Smithsonian Institution’s Laboratories of Analytical Biology. The protocol was as detailed in Riehl et al. (2014).

The analysis of the nuclear gene fragment was done on a restricted number of individuals from the central Pacific and North Atlantic. Twenty specimens representing 14 BINs (up to three individuals per BIN) were chosen, for which both COI and 16S sequences were already available. The fragment of 28S gene sequence amplification was conducted with the combination of two forward (28F, 28S-700F) and two reverse (28R, 28S-1000R) primers (Table 2) and reaction conditions published by Hou et al. (2007). Polymerase chain reaction was performed with AccuStart II PCR SuperMix (Quantabio) for central Pacific individuals, and DreamTaq Green PCR Mastermix (Thermo Scientific) for North Atlantic specimens. Sequencing in both directions was performed by Macrogen Inc., the Netherlands, while the editing of the sequences was similar to the procedure for COI and 16S genes, resulting in 20 sequences of 1152–1261 bp length.

Separate alignments of the sequences of each gene were performed with MAFFT 7 (Katoh et al., 2002; Katoh and Standley, 2013) using the G-INS-i algorithm, the sequences were trimmed to have all of them of the same length. The RNA was homologous enough to use MAFFT software, the alignment was carefully inspected by eye and no ambiguous columns were found. Four alignments were produced: (1) full dataset of all COI sequences with the alignment length of 594 bp, (2) dataset of all available 16S sequences (“short”: length 267 bp), (3) restricted dataset of 16S sequences (“long”: length 398 bp), (4) dataset of all available 28S sequences (length: 1168 bp). Additionally, COI haplotypes were identified and a COI haplotype dataset generated using DNA SP v6 (Librado and Rozas, 2009).

Five molecular species delimitation methods were applied to reveal the Molecular Operational Taxonomic Units (MOTUs). Two methods were distance-based: Barcode Index Number (BIN) System (Ratnasingham and Hebert, 2013), and the Assemble Species by Automatic Partitioning (ASAP) (Puillandre et al., 2021). The following three were tree-based phylogenetic approaches using Generalized Mixed Yule Coalescent (GMYC) model-based method (Pons et al., 2006), according to Monaghan et al. (2009), the Bayesian implementation of the Poisson Tree Processes (bPTP) (Zhang et al., 2013) and multirate Poisson Tree Process (mPTP) (Kapli et al., 2017). The dataset of all COI sequences was used for BIN, ASAP and mPTP. Short 16S sequences were analyzed with ASAP and GMYC, while both the long ones and the 28S sequences were used for the above two methods and mPTP. The restricted datasets of haplotypes of all three markers were the basis for ASAP, bPTP, mPTP, and GMYC analyses.

The BIN method is implemented as part of the Barcode of Life Data system (BOLD; Ratnasingham and Hebert, 2007). It compares newly submitted sequences with the sequences already available. They are clustered according to their molecular divergence using distance-based algorithms (single linkage clustering followed by Markov clustering) that aim at finding discontinuities between Operational Taxonomic Units (OTUs). Each OTU receives a unique and specific code (aka Barcode Index Number or BIN), either already available or new if the submitted sequences do not cluster with already known BINs. Each BIN is registered in BOLD.

The Assemble Species by Automatic Partitioning (ASAP) (Puillandre et al., 2021) is a method that uses pairwise genetic distances to assemble individuals into groups and proposes species partitioning ranked according to a scoring system.

The GMYC method defines MOTUs through identification of the switch from intraspecific branching patterns (coalescent) to interspecific species branching patterns (Yule process) on a phylogenetic tree. Because for GMYC an ultrametric tree is required, as an input, a Bayesian tree was reconstructed in BEAST 2.6.3 (Bouckaert et al., 2019). The site model was set up with bModelTest (Bouckaert and Drummond, 2017). The tree prior was set to Birth-Death following Bayes factors. Two runs of Markov chain Monte Carlo (MCMC) were performed each 20 M generations-long, sampled every 2,000 generations. Runs were examined for convergence in Tracer 1.7 (Rambaut et al., 2018). All runs reached the effective sample size (ESS) above 200 and were combined using LogCombiner 2.6.3. The final tree was summarized with TreeAnnotator 2.6.3, all being part of BEAST 2.6.3 package. The Bayesian tree was uploaded into the R 4.0.5 (R Core Team, 2021) software package ‘SPLITS’ (Species Limits by Threshold Statistics) (Ezard et al., 2009) and analyzed using the single threshold model.

For the following methods, a phylogenetic tree was reconstructed with Maximum Likelihood (ML) approach using RAxML 8.2.8 (Stamatakis, 2014) through raxmlGUI 2.0 (Edler et al., 2021). The best-scoring ML trees were produced using the GTRGAMMA substitution model. Bipartition information was drawn from the phylogenies obtained with the rapid hill climbing tree search algorithm. Statistical supports were estimated with thorough bootstrap tests set to 1000 repetitions.

The bPTP incorporates the number of substitutions in the model of speciation and assumes that the probability that a substitution gives rise to a speciation event follows a Poisson distribution (Zhang et al., 2013). The analysis was performed on the bPTP webserver³ with 500,000 iterations of MCMC and 10% burn-in.

The mPTP method incorporates different levels of intraspecific genetic diversity deriving from differences in either the evolutionary history or sampling of each species. The method implements MCMC sampling that provides a fast and comprehensive evaluation of the inferred delimitation (Kapli et al., 2017). Five runs of 100 M MCMC generations long chain with a burn in of 10% were performed on a local server.

Pairwise p-distances between all recognized MOTUs and P. tenuipes were calculated on the complete COI dataset and on the restricted dataset of COI haplotypes. Within MOTUs, distances were also calculated using p-distance separately for six datasets (two for COI, three for 16S, and one for 28S). To visualize the MOTUs Neighbor-Joining (NJ) tree of all COI sequences was generated based on p-distances (transition and transversion substitution included and complete deletion) with 1000 bootstrap replicates (Felsenstein, 1985) in MEGA V7.0.18 (Kumar et al., 2016).

The distribution of each recognized MOTU was plotted on the World map using QGIS 3.16 (QGIS.org, 2021). The distributions included all records of both nominal species extracted from the relevant literature (Supplementary Material S2). Furthermore, to present the molecular divergence of haplotypes and their geographical allocation, Median Joining Networks were generated in PopART 1.7 (Bandelt et al., 1999) separately for the two nominal species and independently for each recognized MOTU within P. caperesca.

Up to three individuals of each BIN were used to produce a time-calibrated phylogeny. It was based on the combined COI and 16S dataset, so only the individuals with sequences of both genes were included (except for the two BINs AEG2603 and ACZ4873, for which only COI sequence was available). The analysis was conducted in Beast 2.6.3 (Bouckaert et al., 2014). The molecular clock was set using a strict clock and general for gammarid amphipods the COI rate based on multiple calibration points (including fossils) of 0.01773 substitutions/site Ma-1 (Copilaş-Ciocianu et al., 2019). The value is in congruence with other works on Amphipoda (e.g. Mamos et al., 2016). The substitution model was selected via bModelTest (Bouckaert and Drummond, 2017). Birth-Death process was selected as a tree prior. Four runs of the MCMC, each 20 million generations long and sampled every 2,000 generations, were performed and examined for convergence in Tracer 1.7 (Rambaut et al., 2018). All runs reached the effective sample size (ESS) above 200 and were combined using LogCombiner2.6.3. The final tree was summarized with TreeAnnotator 2.6.3, all being part of BEAST 2.5.2 package.

Sampling of the studied taxa in different areas of the World Ocean differed, but it was possible to separate geographic populations for some of the MOTUs (MOTU 1, 4, and 5 of P. caperesca and P. tenuipes) and they were further studied for molecular diversity and demography. In this case, the minimum number of individuals of a particular population was seven. The COI molecular diversity of each population was calculated as the total number of haplotypes and haplotype diversity. Historical demographic patterns were explored using the COI data set employing two approaches. First, to test for a recent demographic expansion, Tajima’s D (Tajima, 1989) and Fu’s Fs (Fu, 1997), Fu and Li’s D (Fu and Li, 1993), Fu and Li’s F (Fu and Li, 1993; Achaz, 2009) Ramos-Onsins and Rozas’s (Ramos-Onsins and Rozas, 2002) indices were calculated using DNASP 6. Their statistical significance was evaluated using 1000 simulated samples. Second, the extended Bayesian skyline plot (eBSP) (Drummond et al., 2005) in BEAST 2.5.2 was used to visualize the demographic changes through time. The clock model, as well as the rate and priors on the substitution models for each group, were determined in the same way as for the time-calibrated phylogeny population model was set to 0.5. The MCMC chain was run two times to ensure convergence for 20 million generations, sampled every 20,000 generations. One run for each data set was used to plot the eBSP in R script⁴ after a 10% burn-in phase.

In order to assess putative relation between geographic and molecular distances we used linear regression model in R 4.0.5 (see footnote 4). The visualization was performed for Paralicella tenuipes and two MOTUs of P. caperesca (MOTU 1 and 5). Sampling locations were divided into the following groups (putative populations): CCZ (Clarion-Clipperton Zone, central Pacific), APEI 6 (North of CCZ, Central Pacific), NWP (North-West Pacific), BB (Brazilian Basin, South Atlantic), PAP (Porcupine Abyssal Plain, North Atlantic). Geographic distances were estimated using Google Earth Pro. Within an ocean they were measured by drawing straight line between central points of each station group. Considering that in the case of the distances between Pacific and Atlantic populations such method would draw the line across the continent, for calculating the distances between these groups of stations the line was drawn to surround South America through the Drake Passage. Because of the low number of groups to compare (max 10) we decided only to run a simple linear regression visualization instead of correlation and isolation by distance tests.

All species delimitation methods clearly separated the two main lineages leading to the nominal Paralicella tenuipes and P. caperesca (Figure 2). A few sequences identified as P. tenuipes that grouped with P. caperesca seem to be associated with identification mistakes and will be discussed later in the text. The interspecies distances of COI between the two species varied from 0.156 to 0.188 (Table 3).

Delimitation within P. tenuipes recognized from one (majority of methods) to four molecular units (Figure 2). The least conservative were BIN system and ASAP based on COI sequences that have recognized four units. GMYC recognized either three (16S short sequences, 28S sequences) or two (16S long sequences) MOTUs. The mean intraspecific p-distance of COI sequences was 0.018 (0.000-0.040), while for COI haplotypes – 0.021 (0.002–0.040) (Table 4). For the 16S dataset, in the case of short sequences the mean distance was 0.009, in the case of long sequences it was 0.005 (0.007 for haplotypes). All 28S sequences were identical.

Delimitation of MOTUs within P. caperesca varied strongly from the recognition of a single unit (ASAP on COI haplotypes and on all 16S combinations, mPTP on 16S haplotypes as well as on 28S sequences and haplotypes, bPTP on 16S and 28S haplotypes) to as many as 21 groups (bPTP on COI haplotypes) (Figure 2). It confirms the existence of separate lineages. After removing the extreme unifications and divisions, it may be observed that four lineages represented by the following BINs: BOLD:ADD2929 (MOTU 2), BOLD:AEH6662 (MOTU 3), BOLD:ACZ5628 (MOTU 4), BOLD:ADD2497 (MOTU 6) were most commonly separated from the rest of the groups. The lineages forming the remaining two groups did not present a stable pattern of unification/separation and for this study were treated as two separate MOTUs. MOTU 1 consisted of seven BINs (BOLD:AEG0263, BOLD:ACZ5625, BOD:ACZ5671, BOLD:ACZ5631, BOLD:ACZ4905, BOLD:ACZ4904, BOLD:ACZ5630), whereas MOTU 5 grouped six BINs (BOLD:ACZ5627, BOLD:ACZ5629, BOLD:ACZ4489, BOLD:ADP2618, BOLD:ACZ4873, BOLD:ACZ4903). The COI p-distances between discriminated MOTUs of P. caperesca ranged from 0.064 (between MOTU 1 and 2) to 0.128 (observed for the MOTU 2 and 4 pair) (Table 3). The intraspecific p-distances within nominal P. caperesca varied from 0.000 up to 0.133 of COI, from 0.000 to 0.094 of 16S and from 0.000 to 0.014 of 28S (Table 4). These values were clearly reduced when the species were divided into suggested MOTUs. Within MOTU 1 p-distance of COI reached up to 0.072 (mean 0.028 for sequences and 0.035 for haplotypes). Within 16S, the maximum p-distance was 0.027, when the dataset of short (ca. 250 bp) sequences was used, while it slightly raised up to 0.028 when the longer fragment (ca. 400 bp) was analyzed. The 28S sequences within MOTU 1 are identical. MOTUs 2, 3, 4, and 6 showed intraspecific p-distances not higher than 0.003 no matter the dataset explored. In the case of MOTU 5, the maximum p-distance of COI reached 0.059 (mean 0.029 for sequences, 0.031 for haplotypes). For the 16S gene, the dataset of short sequences showed lower values of p-distance (max – 0.019, mean – 0.009), while they reached up to 0.028 (mean 0.015 for sequences, 0.018 for haplotypes) in case of longer sequences. The nuclear gene was diverse with a mean p-distance of 0.011 (maximum – 0.014). MOTU 1 and 3 shared a 28S sequence, while MOTU 5 was represented by three clearly different sequences, additionally separated by the 28S sequence of MOTU 4 (Supplementary Figure S1).

The literature records of both nominal species come from all three oceans; the molecular study confirmed the pan-oceanic distribution of P. tenuipes (Figure 3) and MOTU 1 of P. caperesca (Figure 4). MOTU 5 was widely distributed in both Atlantic and Pacific oceans, whereas MOTUs 4 and 6 were observed at a few localities in the Pacific only. The most restricted geographically appeared to be MOTU 2 and MOTU 3, each recognized from a single station, the first in South-West Pacific and the second in the North Atlantic. Although in some localities only a single MOTU were observed, in the areas where more individuals of P. caperesca were collected and studied, different lineages co-occurred (Figure 4).

Out of the 19 haplotypes identified within P. tenuipes, seven were singletons, the remaining ones were often shared between the studied regions (Figure 5). As many as 171 individuals represented one of the five dominant haplotypes, 34 individuals characterized the remaining 14 haplotypes. The majority of haplotypes of P. caperesca (56 out of 80) were singletons. Within MOTU 1, there were four haplotypes shared between geographical regions, with one almost equally represented in Pacific and Atlantic oceans. One additional haplotype was shared between North and South Atlantic, while another one was found in South Atlantic and Indian Ocean (Figure 6). The structure of the haplotype network of MOTU 1 has a partially star-like topology with one central haplotype present in Central Pacific and several haplotypes, differing from it by a few mutations only, being widely distributed geographically. MOTU 2, MOTU 3, and MOTU 4 were each restricted to a single region, but all of them were represented by a few individuals only. Among the 22 haplotypes recognized within MOTU 5, only two were shared between regions, one identified in Central and South-East Pacific, while the second present in North-West and South-West Pacific. MOTU6 was represented by two haplotypes, one of them present in North-West and South-West Pacific.

The Bayesian phylogenetic reconstruction confirmed the existence of two main lineages leading to the nominal Paralicella caperesca and P. tenuipes (Figure 7). This separation can be estimated at ca. 37 Ma. The differentiation within P. tenuipes can be estimated as happening within the recent 1 Ma. The diversification within P. caperesca appeared between eight to two Ma. Around eight million years ago, the lineages were divided into two groups: combining MOTUs 1–3 and MOTUs 4–6. Soon after this, the MOTU 6 separated from the rest, and six million years ago MOTU 3 diverged. Division between MOTU 4 and 5 happened ca. four Ma, while the most recent separation was of MOTU 2 and MOTU 1, estimated at ca. two million years ago.

All values of haplotype diversity of P. tenuipes were significant with the highest observed in the Central Pacific and the lowest in the North Atlantic (Table 5). Haplotype diversity of MOTU 1 of P. caperesca (represented by 104 individuals in our study) was highest in the North-West Pacific, but this value was not statistically significant (Table 4). In the Central Pacific population this measure was lower than in the North-West Pacific, but in contrast to this region it was significant. The Atlantic population was the least diverse of the three and the value was insignificant. Within MOTU 5, the haplotype diversity of the North Atlantic population was distinctly higher than that of the Central Pacific (both values significant), however, only nine individuals of MOTU 5 were collected in the whole Atlantic. From the results of neutrality tests, some recent population contraction of Paralicella tenuipes may be seen, particularly in South Atlantic. Whereas, in the case of MOTU 1 of P. caperesca, a slight population expansion in the Atlantic appeared. MOTU 5 seemed to have a stable population in the Central Pacific, while in the North Atlantic it was a sign of recent population contraction. The neutrality tests were significant only in a few cases of MOTU 1 of P. caperesca. The results of eBSP confirmed no significant drops or sudden expansion of populations of both species (Supplementary Figures S2, S3), whereas linear regression revealed no separation of the populations of P. tenuipes and MOTU 5 of P. caperesca and a weak signal of differentiation for MOTU 1 of the second species (Supplementary Figure S4).

Our study undoubtedly confirmed the separation of P. tenuipes and P. caperesca. All individuals identified morphologically by the taxonomists (TH, EH) as certain species were molecularly ascribed to the correct group. This confirms the former findings by Mohrbeck et al. (2021) and is in contrast with Ritchie et al. (2015), who indicated that the characters used for morphological species identification are not separating these two taxa correctly. Because the publication of wrongly identified sequences may have a large influence on future studies, it is important to clarify these issues which we have done here (Table 6). Using the sequence information of both COI and 16S, it can be seen that four individuals from Ritchie et al. (2015) identified as P. tenuipes and two presented as “unidentified primitive lysianassoid” appear to be representatives of different clades of P. caperesca. Moreover, one individual cited as Valettietta anacantha is represented by the COI sequence (GenBank accession number: KP713950) belonging to MOTU 1 (BOLD:ACZ5625) of P. caperesca. The 16S of the same individual does not show similarity to any Paralicella species but it groups with sequences of V. anacantha from the study of Ritchie et al. (2015). The relatively short fragment of the 18S gene (591 bp, GenBank accession number KT372893) is the only available sequence of 18S of the genus Valettietta and the whole family Valettiopsidae and presents 93-95% similarity to various species of Alicellidae from the same study. There are no sequences of this gene provided by other researchers to cross-validate the sequence identity. As a result, it is impossible to determine the reason why the COI sequence is almost identical to those of P. caperesca while the 16S resembles another taxon. As such, we recommend that this record be removed from GenBank, or updated with proper comment, to avoid further confusion. In a study by Iguchi et al. (2020), the identifications were of Amphipoda at the order level. The comparison of COI sequences obtained by these researchers with our data revealed 12 individuals to be P. caperesca (three different MOTUs) and one individual of P. tenuipes (Table 6). There are also two sequences of P. caperesca from the region of the Crozet Islands, provided by Corrigan et al. (2014). They are both short (COI – 262 bp, 16S – 301 bp) and do not group with any available sequence of Paralicella stored either in GenBank or BOLD. These records should also be removed from public databases.

The diversity patterns of the two nominal species studied differ considerably. Recent work by Bribiesca-Contreras et al. (2021) and some species delimitation methods used by us (BINs, ASAP on sequences, and GMYC) suggested further division of P. tenuipes, but this is not supported by the majority of methods and should be treated as so-called oversplitting of taxa. The fact that mostly one unit was recognized and the observation of low intraspecific distances confirm that no further division of P. tenuipes is justified. Oversplitting has been reported for other deep-sea species of Amphipoda from the North-West Pacific (Lörz et al., 2018a; Jażdżewska and Mamos, 2019) and could partly be explained by unequal sampling of all populations of certain species. It is worth noting also that the length of the studied fragment of the gene may have an influence on the final delimitation results. In our study, it can be seen that in the case of 16S sequences, when the shorter fragment was considered (ca. 250 bp), both the mean and maximum uncorrected p-distance were higher than when the longer (ca. 400 bp) sequences were analyzed (Table 4). It derives from the higher concentration of variable sites (4.5%) in the shorter fragment of 16S amplified by former authors studying P. tenuipes (Ritchie et al., 2015; Weston et al., 2021).

The opposite situation can be seen within the nominal P. caperesca that was further divided into separate MOTUs by the majority of analytical methods (Figure 2) that suggests the existence of cryptic speciation. The least conservative methods suggested that P. caperesca consists of as many as 18-21 units (e.g., BINs in BOLD). In this case, it appears that the BIN system applied in BOLD may oversplit the taxa, which derives from the fact that the threshold used is set at ca. 2% of COI sequence similarity (Ratnasingham and Hebert, 2013). The threshold values for this gene more successfully separating amphipod species are set at 3–7% of molecular similarity and are suggested to be family specific (Costa et al., 2009; Knox et al., 2012; Lobo et al., 2017; Tempestini et al., 2018; Jażdżewska and Mamos, 2019). As a result, the best justified division suggests the existence of six separate MOTUs within P. caperesca. Apart from two dominant MOTUs (MOTU 1 and MOTU 5), the remaining ones are represented by one to eight individuals. Samples containing Paralicella generally contain very large numbers of individuals (hundreds to thousands of specimens in each sample) and the sorting and identification is therefore generally carried out by the use of rapidly assessed morphological characters, indicating initially that the specimen belongs to the genus Paralicella. This is followed by separation at the species level as either P. tenuipes (with a distinct small red/brown eye and a strong bevel on P7) or P. caperesca (lacking the distinct small red/brown eye and lacking a bevel on P7). It should be noted that P. caperesca specimens possess a white or orange pigmented diffuse non-ommatidial eye in fresh material, but this often fades after preservation, and therefore it cannot be relied on for species determination.

The fine characters which are needed to distinguish between the two nominal species in this study and the other members of the genus (P. fusiformis, P. microps, P. similis, and P. vaporalis) are not routinely checked during the sorting process. Now that there is a clear indication that the P. caperesca complex comprises more than a single entity, extra care will need to be taken when sorting and identifying Paralicella samples, including examination of coxa 1, uropods and possibly some mouthpart characters to determine the species of Paralicella that is being dealt with. Attention will need to be paid to the means by which it will be possible to clearly discriminate between the MOTUs, particularly considering the likelihood that they co-occur throughout their range. This will mean revisiting the already collected samples of many thousands of individuals, to first determine a means to distinguish between the MOTUs morphologically and to describe each in full (whilst in alignment with the known 6 species), only then will we be able to clarify their distribution and depth ranges. Presently, only a small group of specimens representing four of the MOTUs defined here was available for morphological study. Detailed description of the available molecularly defined specimens is now underway.

Until such work can be completed, where studies involve morphological species delimitation, P. caperesca should be treated as a species complex, and cited as such in all publications. It is also important to note that as for P. tenuipes, the length of the 16S sequence fragment influences the uncorrected p-distance results. In the case of MOTU 1, the use of a short fragment suggests only a slightly higher diversity of sequences in comparison to longer ones because the extension of the studied fragment does not change the proportion of the variable sites reaching in both cases ca. 4.5%. For MOTU 5 the opposite pattern may be observed (Table 4) deriving from the considerable increase of variable sites (from two to six percent) in the longer sequence. The results derived from the longer fragment of 16S seem to be more congruent with the results from the other two genes studied, so it is highly recommended to use the longer fragment of 16S in future molecular studies.

Intraspecific diversification has already been observed by Ritchie et al. (2015) who distinguished four clades within the studied Paralicella specimens. ‘Group 1’ recognized in that study corresponds to P. tenuipes, ‘Group 2’ relates to our MOTU 6 of P. caperesca, ‘Group 3’ corresponds to our MOTU 5, whereas ‘Group 4’ combines individuals from our MOTU 1 and 2. In the subsequent paper, presenting the population structure of two Paralicella species collected in five distant Pacific trenches (Ritchie et al., 2017) this division was only partly used – Group 1 (P. tenuipes) was named RFLP sp. 1, while groups 2–4 (P. caperesca) were combined into a single unit – RFLP sp. 2.

One of the important issues emerging from the use of molecular tools to identify species, particularly when cryptic diversity is noticed and the type collection used for species description is not available for genetic analysis, is to decide which of the molecularly recognized units represents the originally described taxon. For P. tenuipes, described from the region of Canary Islands (Chevreux, 1908), the closest area presently studied was the North Atlantic where representatives of three BINs were identified (Figure 5). Nevertheless, they all constitute a single species, and most probably they all inhabit the type locality of this taxon, which is further justified by the presence of the same BINs and haplotypes in the South Atlantic. Since cryptic diversity has been recognized within P. caperesca, deciding which of the lineages constitutes that originally described (Shulenberger and Barnard, 1976) is more difficult. This species was described from a large collection of amphipods from North of Hawaii; molecular analysis was not considered at that time. France and Kocher (1996) attempted to extract and study DNA of certain deep-sea scavenging species including P. caperesca and P. tenuipes collected in 1977 and 1978 in the Central North Pacific, most probably near the locus typicus of the first species. Their study resulted in four sequences of a short fragment of the 16S gene (ca. 160 bp) (GenBank accession numbers: U92692, U92693, U92694, U92695). The single sequence of P. tenuipes groups with sequences of individuals representing BOLD:ACZ6441 confirming its identity. Two sequences of P. caperesca have affinity to the 16S sequences of individuals ascribed to BOLD:ACZ5625, while the third is almost identical to the sequences of two individuals, one representing BOLD:ACZ6571, the second – BOLD:ACZ4905. All of these BINs belong to MOTU 1 of P. caperesca, which may suggest that the species description was based on individuals of this molecular unit. However, since the original collection studied by Shulenberger and Barnard (1976) consisted of several individuals and our study has already revealed that different MOTUs of P. caperesca co-exist in the Central Pacific, any decision on which MOTU can represent P. caperesca sensu stricto must be preceded by detailed morphological examination and, if possible, molecular study of the types.

Based on morphological studies, both studied species have been considered as widely spread or even cosmopolitan (Table 7; Shulenberger and Barnard, 1976; Brandt et al., 2012). Our study is the first that provides molecular evidence for the cosmopolitanism of P. tenuipes, as well as a very wide distribution of two of the MOTUs within P. caperesca. It confirms also the wide bathymetric ranges (2189–3673 vertical metres) of these species, mainly collected at abyssal depths (Table 7). It should be underlined that the restricted number of individuals (five in total) from the Indian Ocean for which the sequences were available might have an influence on our knowledge of these species’ distributions. However, a recent molecular study of scavenging amphipods from hydrothermal vent fields in South-East and Central Indian Ridge did not reveal the presence of either P. tenuipes or any MOTU of P. caperesca there (Kniesz et al., in review⁵). The specific conditions of the areas adjacent to vent fields may influence the composition of the scavenging fauna and may be the reason for not recording the presently studied species. Additionally, the traps used in that study were placed relatively shallow (2500–2900 m) while P. tenuipes and P. caperesca seem to prefer greater depths. P. tenuipes was not recorded south of Antarctic Convergence, whereas individuals resembling P. caperesca have been recorded in the area of South Shetland and South Sandwich Islands as well as in the Lazarev Sea (De Broyer et al., 2004). Those animals however, possessed certain characters that differed from the described species (De Broyer, personal communication), since they were not studied molecularly, we have not included these distribution records in our study.

Low genetic diversity and wide distribution of deep-sea species have been postulated by various authors (Zardus et al., 2006; McClain and Hardy, 2010; Etter et al., 2011) however, only a relatively low number of peracarid species are reported to present large horizontal ranges (Brandt et al., 2012). Recent studies of widely distributed species identified solely by morphology have revealed species complexes of multiple taxa with more restricted ranges when molecular methods are applied (Havermans et al., 2013; Verheye et al., 2016; Jakiel et al., 2019). Although cosmopolitanism of deep-sea species appears to be less common than previously thought, certain scavenging amphipod species have been confirmed by molecular studies as widely distributed. These include Abyssorchomene distinctus (Birstein and Vinogradov, 1960), Eurythenes magellanicus (Milne Edwards, 1848), E. maldoror (d’Udekem d’Acoz and Havermans, 2015), Bathycallisoma schellenbergi (Birstein and Vinogradov, 1958) and Haptocallisoma abyssi (Oldevig, 1959) (Havermans, 2016; Jażdżewska et al., 2020; Mohrbeck et al., 2021; Weston et al., 2021; Kniesz et al., in review⁵). Very few amphipod species have had a very wide (>2 km) vertical distribution molecularly confirmed, as has been observed here for Paralicella. Similar to the situation of wide geographic ranges inferred from morphological identification, several taxa previously listed as recorded from a wide bathymetric range appear to be more restricted when their genetic intraspecific structure is checked. Among amphipod species of confirmed wide bathymetric range are Abyssododecas styx Takeuchi et al., 2016 and Rhachotropis saskia Lörz & Jażdżewska (Takeuchi et al., 2016; Lörz et al., 2018a). These two species inhabit abyssal and hadal depths of North-West Pacific, whereas the Paralicella species from our study are mainly abyssal taxa that are also collected in bathyal depths. There is a report of the collection of P. caperesca in the Indian Ocean that reached the research vessel deck still alive (Treude et al., 2002) suggesting high resistance of this species to decompression, however, this was an incidental case as the majority of collected individuals were dead on arrival at deck. However, other studies have confirmed the resistance of this species to decompression (Macdonald and Gilchrist, 1980; Yayanos, 1981). It is known that these species do cross the bathyal-abyssal border and have often been collected in the water column several metres above the seafloor (Ingram and Hessler, 1983; Thurston, 1990). Amphipods of the genus Paralicella are known to be obligate scavengers (Havermans and Smetacek, 2018; Horton et al., 2020). Carrion (particularly large carcasses of nekton), although providing nutrition for a long time and being more common in the deep sea than previously expected, is believed to be irregularly deposited and unevenly distributed over the seafloor (Smith and Baco, 2003; Havermans and Smetacek, 2018). Near feeding or spawning grounds as well as along species migration corridors, the availability of carrion is high (Smith, 2007) but in other open ocean areas this may not be the case. The irregularity of food availability implies their special adaptations including, their well-developed chemosensory system and good swimming abilities allowing quick access to available food sources (Hessler et al., 1978; Thurston, 1979; Ingram and Hessler, 1983; Klages et al., 2002; Premke et al., 2003). These adaptations may partly explain the wide horizontal ranges of the studied species, and to consider that resistance to changing pressure may also allow these amphipods to profit from food at different depths.

The separation of the two studied nominal Paralicella species from another species from the family Alicellidae (the supergiant amphipod, Alicella gigantea Chevreux, 1899) appeared to be ca. 80 Ma, which is earlier than recorded by Copilaş-Ciocianu et al. (2020). However, when considering the highest posterior density intervals in both cases the time of lineage separation overlap. The separation of the two nominal Paralicella species is reported here at ca. 35–37 Ma, which is similar to the divergence of Paralicella caperesca and Valettietta anacantha recorded by Copilaş-Ciocianu et al. (2020). However, in that study P. tenuipes was not considered. The time of diversification within Alicellidae, occurring from the Late Mesozoic to Middle Cenozoic (Figure 7), is congruent with the global climatic cooling during that period, and the transition from a halothermal ocean circulation into thermohaline one (McClain and Hardy, 2010; Copilaş-Ciocianu et al., 2020). The thermohaline ocean circulation is driven by high-latitude deep-water formation and results in cold, oxygenated waters at the bottom. Because amphipods are considered as well adapted to cold conditions with low tolerance to hypoxia (Copilaş-Ciocianu et al., 2020) the change of environmental conditions that has “opened” vast deep-sea bottom areas for colonization could have promoted the speciation of the studied species. The speciation of the two nominal species of Paralicella occurred at the time of serious reconfiguration of continents, when the Isthmus of Panama was still open, while the Drake Passage and the West Wind Drift was in the process of formation (Brandt et al., 2007; O’Dea et al., 2016; Straume et al., 2020). The fact that both species emerged already at that time may partly explain their present wide distribution – there was enough time to spread and different oceanic gateways were available. It should be considered that these species might have originated from shallow waters and only during the further stages of speciation submerged, so at the initial stage of expansion they might have been less dependent on the water depth both in the area of Isthmus of Panama and Drake Passage. The further speciation within P. caperesca is predicted as happening from ca. 7 Ma. At that time, only small and shallow water basins were present in the area of Isthmus of Panama, preventing the connection of the Central Pacific populations and may be responsible for the speciation of lineages geographically restricted to the Pacific (MOTUs 2, 4, 6) or to the Atlantic (MOTU 3). The full opening of the Drake Passage (reaching abyssal depths) and the formation of the West Wind Drift would have allowed maintenance of intraspecific connectivity. This contemporary Atlantic and Pacific population connectivity of P. tenuipes, as well as between MOTU 1 and MOTU 5 of P. caperesca, is confirmed by the distribution of haplotypes that are often shared between very distant regions (Figures 5, 6) as well as by the lack, or very weak signal, of population separation within each of them (Supplementary Figure S4). On a more local scale, the constant gene flow within P. tenuipes and the two most widely distributed MOTUs of P. caperesca was recently observed by Bribiesca-Contreras et al. (2021) studying scavenging Amphipoda from three Areas of Particular Environmental Interest in the Clarion-Clipperton Zone. A similar pattern was recorded by Ritchie et al. (2017) for Paralicella RFLP sp. 1 that corresponds to P. tenuipes. In that case, although generally bidirectional, North-West direction of migration predominated. At the same time, a very weak signal of migration was observed for Paralicella RFLP sp. 2, but this may be explained by the fact that the taxon combined four different MOTUs of which two have a very restricted distribution. It should not be overlooked that the authors considered the studied species as the “trench” ones, not considering that from morphology-based literature both were reported from Central Pacific abyssal plain (Ritchie et al., 2017). It is also worth noting no evidence for recent population size changes of the P. tenuipes and the two MOTUs of P. caperesca that remains in contrast with other findings for deep-sea species. The significantly negative Tajima’s D values observed for certain vent and non-vent taxa were interpreted as resulting from recent expansion of populations after disturbance events, suggestive of the instability of deep-sea habitats (Vrijenhoek, 2010; Taylor and Roterman, 2017). Although we are aware that some of the populations studied by us consist of a small number of individuals it appears that the populations of studied Paralicella species are stable over time, implying robustness and high plasticity of these taxa. It may be expected that they are less prone to unfavorable (human or non-human induced) events that occur in the abyss.