An integrative approach challenges species hypotheses and provides hints for evolutionary history of two Mediterranean freshwater palaemonid shrimps (Decapoda: Caridea)

Abstract The Mediterranean Region is a biodiversity/endemism hotspot whose freshwater fauna remains largely unexplored. Our integrative study challenges the taxonomic status of two freshwater palaemonid shrimps, Palaemon antennarius and Palaemon minos. Three molecular operational taxonomic units (MOTUs) were defined based on 352 cytochrome oxidase subunit I (COI) sequences and 88 haplotypes. Two belonged to P. antennarius: one inhabiting the Apennine Peninsula and Sicily, and the other from the Balkan Peninsula. Palaemon minos was the third MOTU, found on Crete. The Balkan MOTU of P. antennarius was genetically closer to P. minos than to the other conspecific MOTU. Data from a nuclear marker (Histone 3) is congruent with such a pattern. The carapace shape variation (based on 180 individuals) was mainly explained by the geographical distribution. Balkan and Cretan groups were clearly recovered, while other samples clustered along a shape gradient from Sicily, through the Apennine Peninsula to the Balkans. Our results show that, for taxonomic consistency, the MOTU inhabiting the Balkan Peninsula should be either described as a new separate species or synonymised with P. minos. The third possible option would be treating all the populations as part of P. antennarius. Geometric morphometrics supports the first option, phylogenetic reconstructions point to the second one, yet the low genetic divergence favours the third one, illustrating that even emblematic taxa such as shrimps require an in-depth integrative approach.


Introduction
The species is considered a fundamental unit in biological systematics and a basic and convenient unit to measure biodiversity. In addition, speciation itself is a key topic in evolutionary biology (Mayr 1976;De Queiroz 2007). However, there is no consensus concerning the definition of species, as it may refer to various properties of organismal biology, but also to evolution and phylogeny (Mayr 1976;Ghiselin 2001;De Queiroz 2005Wiens 2007). Thus, despite the important need for delimitation of species, the task encounters a lot of difficulties. To tackle such difficulties, combining methods from various fields of studies, as advocated by the integrative taxonomy approach, is considered to be efficient in eliminating failure in the delimitation process (Dayrat 2005;Padial et al. 2010;Schlick-Steiner et al. 2010;Rajaei 2015).
Anatomical quantification used in taxonomic studies commonly employs either the so-called multivariate traditional morphometrics (e.g. Anastasiadou et al. 2009) or geometric morphometrics (e.g. Torres et al. 2014). A number of works have examined the usefulness of the two methods in tackling the same questions and found them to be useful in establishing boundaries between species (e.g. Parsons et al. 2003;Navarro et al. 2004;Fruciano et al. 2011;Schmieder et al. 2015;Ramírez-Sánchez et al. 2016;Lovrenčić et al. 2020). Nonetheless, taxonomic studies supported by geometric morphometrics (i.e. shape variation) remain scarce (e.g. Fruciano et al. 2011;Schmieder et al. 2015;Ramírez-Sánchez et al. 2016;Navarro et al. 2018). In recent years, there has been a tremendous increase in the number of taxonomic studies combining morphological and DNA data (e.g. Grabowski et al. 2017b;Hupało et al. 2018;Rudolph et al. 2018). However, combining traditional morphology (including multivariate morphometrics), geometric morphometrics and molecular data seems rare (e.g. Arnoux et al. 2014;Fruciano et al. 2016;Celik et al. 2019), even though it is proven to be very efficient, not only in helping to formally describe nominal species, but also for better understanding of speciation mechanisms, as shown by e.g. Sangster (2018) or Zheng et al. (2020). Nevertheless, in some cases, e.g. Young et al. (2019) on stoneflies, the integration of delimitation methods can actually lead to disproving species hypotheses.
Known for a complex combination of geological and climatic histories, the Mediterranean region is considered to be an ideal area to conduct studies of speciation and biogeography (Myers et al. 2000;Tierno de Figueroa et al. 2013). Past processes, among other things, significantly influenced the evolution and composition of local freshwater fauna. First, the Alpine orogeny, which started in the Mesozoic and reached its greatest intensity in the Paleogene, is considered a factor of main impact on the landform of that area (Skoulikidis et al. 2009). Second, several eustatic fluctuations of sea level might have had a significant impact. Such events included the regression of the Tethys Ocean that severed the connection between the Indian Ocean and the Mediterranean Basin (Bialik et al. 2019); the major evaporation of the proto-Mediterranean Sea during the so-called Messinian Salinity Crisis, which lasted from 5.96 to 5.33 Ma; and the transgression of the neighbouring epicontinental Paratethys Sea, known as the Lago Mare episode (Hsü et al. 1977;Popov et al. 2004;Krijgsman et al. 2018;Bialik et al. 2019). These events caused changes in hydrological conditions, expressed mainly by recurrent salinity alternations and fragmentations/reconnections of inland aquatic ecosystems (Bianco 1990;Skoulikidis et al. 2009).
According to Mittermeier et al. (2011), the Mediterranean Basin itself could be considered the second largest known hotspot for biodiversity and endemism among the 35 most important hotspots recognised worldwide. Furthermore, some more restricted areas within the basin might be of particular interest and considered local hotspots or key biodiversity areas (KBAs), such as lakes (e.g. Lake Ohrid, Lake Skadar) and sites located on Mediterranean islands (e.g. Sicily or Crete) (Eken et al. 2004;Darwall et al. 2014).
The epigean inland waters of the Mediterranean region are inhabited by a very diverse invertebrate fauna, including many taxa of crustaceans (Balian et al. 2008;Tierno de Figueroa et al. 2013). Among them, two families of freshwater shrimps (Malacostraca: Decapoda: Caridea), Atyidae and Palaemonidae, have attracted a lot of attention during recent years and have undergone intensive taxonomic research. While studies concerning atyids were conducted mostly with integrative methods combining traditional morphology and genetics (Christodoulou et al. 2012;García Muñoz et al. 2014;, palaemonids were investigated mainly based on morphological features (Gottstein Matočec et al. 2006;Anastasiadou et al. 2009;Tzomos & Koukouras 2015), although single specimens with mitochondrial 16S rDNA and/or nuclear Histone H3 markers applied were reported by Ashelby et al. (2012), Cuesta et al. (2012) and Carvalho et al. (2017).
Until recently, it was supposed that there were six species of Palaemon Weber, 1795 inhabiting fresh waters of the Mediterranean, among which P. antennarius H. Milne Edwards, 1837 was believed to have the widest distribution in the region. In their recent study, Tzomos and Koukouras (2015), based on a sampling covering the species range, reexamined morphological variation of P. antennarius and P. migratorius (Heller, 1862), using mostly qualitative features, resulting in the description of two new species: one (P. minos Tzomos & Koukouras 2015) believed to be endemic to Crete, and the other (P. colossus Tzomos & Koukouras 2015) endemic to Rhodes Island and Anatolia. Therefore, the most recent review lists nine Palaemon representatives in inland waters of the Mediterranean region (Christodoulou et al. 2016). Palaemon antennarius is thought to occur in Italy, Slovenia, Croatia, Albania, Montenegro and Greece (Christodoulou et al. 2016). Its range covers areas of the Apennine and Balkan peninsulas (including fresh and brackish waters), surrounding the Adriatic part of the Mediterranean, and consequently belongs to the so-called peri-Adriatic region. The region is of particular concern to taxonomists and evolutionary biologists, due to its geological history, which, among other things, promotes diversification of local freshwater fauna (Jelić et al. 2016;Grabowski et al. 2017a).
The present study focuses on P. antennarius and P. minos. Taking into account that: (1) recurrent regressions/transgressions of the Adriatic and the neighbouring seas from Miocene to Holocene may promote speciation in local fresh waters; (2) P. antennarius is widely distributed in the peri-Adriatic including the Lake Skadar basin which is known as a local hotspot of biodiversity and endemism , as well as in Sicily, an island known to be associated with endemism for many other taxa (Signorello et al. 2018;Hupało et al. 2021); (3) the species is known to show some phenotypic variation related to sex and habitat (Anastasiadou et al. 2009); and (4) P. minos was erected mostly based on qualitative anatomical features (Tzomos & Koukouras 2015), we decided to challenge the species hypotheses for these two taxa using an integrative taxonomy approach, consisting of morphological identification, geometric morphometrics and genetic investigation. We wished to test whether such combined taxonomic methods would be an informative tool to provide consistent results in species recognition of palaemonid shrimps.

Sample collection and identification
A total of 389 individuals of freshwater palaemonid shrimp were examined in this study. All the material was collected in the years 2004-2016 in the Apennine Peninsula (126 individuals from 16 sampling sites), Sicily Island (30 ind., 2 sites), Balkan Peninsula including the Lake Skadar basin (156 ind., 26 sites) and Crete Island (60 ind., 1 site) ( Figure 1; Table I). The samples were gathered with the use of a benthic hand-net or dredge. They were sorted on site and immediately fixed in 96% ethanol. Then the shrimp individuals were Figure 1. Map of the research area. Sampling sites are indicated with colour dots, for the different MOTUs/haplogroups: yellow -northern APS (Apennine Peninsula and Sicily); red -southern APS; green -BP (Balkan Peninsula); black -CI (Crete Island). The black square is indicating the location of Skadar Lake, an enlarged map of which is inset into the upper right corner of the figure. Small grey dots indicate the sampling sites where palaemonid shrimps were not found.

902
A. Jabłońska et al.   morphologically examined under NIKON SMZ 800 stereomicroscope and identified to the species level on the basis of descriptions included in the paper by Tzomos and Koukouras (2015), supported by information from González-Ortegón and Cuesta (2006). In particular, Tzomos and Koukouras (2015) presented the following two features as obvious differences between P. minos and P. antennarius: (i) the rounded, not pointed, fifth pleuron end of P. minos vs. the strongly pointed fifth pleuron end in P. antennarius; and (ii) the plumose setae overreaching the inner spines on the distal part of the telson vs. the shorter or equal-length setae in P. antennarius. These features were quantified on all specimens in the present study to verify whether they are indeed conserved features. The material has been stored in the permanent collection of the Department of Invertebrate Zoology and Hydrobiology (University of Lodz, Poland).

DNA isolation, sequencing, alignment
The total DNA was isolated from pleopod muscle tissue of 352 individuals by a standard proteinase K and phenol/chloroform extraction (Hillis et al. 1996) or by the Chelex procedure (Casquet et al. 2012).
The cytochrome oxidase subunit I (COI) fragment was amplified in Polymerase Chain Reaction (PCR) with the use of the LCO JJ and HCO JJ primer pair (Astrin & Stüben 2008), according to the protocol provided by Hou et al. (2007). A subset of six specimens for each molecular operational taxonomic unit (MOTU, see below) was used for amplification of Histone H3 nuclear DNA fragments with the use of the primer pair HisH3f and HisH3r (Corrigan et al. 2014). The PCR product was amplified with the initial denaturation at 94°C for 1 min 50 s, followed by 30 cycles of 30 s at 94°C, 30 s at 45°C and 1 min at 72°C. The final extension was conducted for 5 min at 72°C. For both markers, all PCR products were purified with exonuclease I (Exo I) and alkaline phosphatase (FastAP) and then sequenced with BigDye terminator technology by Macrogen Inc., Europe.
The amplified markers were verified as belonging to the genus Palaemon against the GenBank resources using BLAST (Altschul et al. 1990). The obtained sequences were aligned and trimmed to the same length of 564 and 310 nucleotides for COI and H3, respectively, using Geneious 10.0.2 software (Kearse et al. 2012). The COI haplotypes, haplotype diversity (Hd) and nucleotide diversity (pi) were defined using DnaSPv5 software (Librado & Rozas 2009). All of the obtained sequences (COI and H3) were deposited in GenBank (Benson et al. 2005) with accession numbers MT517444-MT517795 and MT517796-MT517813 for COI and H3, respectively (Table I), and deposited with BOLD Systems (Ratnasingham & Hebert 2007).

MOTU delimitation
The level of cryptic diversity described by the number of MOTUs was assessed using two distance-based methods: (i) the barcode index numbers (BINs) method (Ratnasingham & Hebert 2013) and (ii) the assemble species by automatic partitioning (ASAP) method based on the Kimura 2-parameter (K2p) model (Puillandre et al. 2020). In addition, between-MOTUs average pairwise K2p genetic distances were calculated in MEGA X software (Kumar et al. 2018).

Phylogenetic relationships among haplotypes
The phylogenetic relationships among the mitochondrial COI haplotypes were explored using the neighbour-net method (Bryant & Moulton 2004) based on the K2p model (default options). The final phylogenetic network was produced and bootstrapped (1000 replicates) with the splits network algorithm method (Dress & Huson 2004). The whole analysis was performed using SplitsTree4 (Huson 1998;Huson & Bryant 2006).
The evolutionary relationships among the nuclear H3 haplotypes were illustrated with the maximum likelihood (ML) tree, based on the K2p model, where the validity of nodes was estimated with the bootstrap test (1000 replicates) (Felsenstein 1985;Saitou & Nei 1987). The analysis was conducted in MEGA X (Kumar et al. 2018). The following sequences were taken from GenBank and used to supplement the tree: Palaemon antennarius (acc. no. KP179081), and, as outgroups, Palaemon adspersus (acc. no. KP179094) and Palaemon elegans (acc. no. KP179102), all deposited by Carvalho et al. (2017).

Defining dataset for geometric morphometrics
A set of 180 individuals, 60 per MOTU, also fitting three defined geographic units (see Results) was chosen for the geometric morphometric analysis (see Table I). These 60 individuals included an equal number of mature males and females, excluding ovigerous females. It should be noted that within these 60 individuals per geographic unit, at least 45 individuals were also DNA barcoded. Although not all individuals processed for morphometrics were also processed for COI, the allopatric distribution of MOTUs combined with our sampling effort implies that a given individual from a given area would belong to the particular MOTU.

Landmark designation
For the landmark-based analysis, the carapace including rostrum was chosen. This part of the body is commonly used in such studies (Ashelby 2012;Zimmermann et al. 2012;Torres et al. 2014;Sganga et al. 2016;De Melo & Masunari 2017) and not only in shrimps (Silva et al. 2010;Simanjuntak & Eprilurahman 2019;Lovrenčić et al. 2020).
Photographs of the right side of the carapace including rostrum were digitised with the use of a Nikon MM60 measuring stereomicroscope (10× magnification), coupled with a Nikon DS-Fi2 (5MPixels) camera. Based on the protocol proposed by Torres et al. (2014), a total of 14 landmarks were used, encompassing the shape of the carapace including the rostrum (Figure 2), and digitised in photographs with the use of tpsDig 2.32 software (Rohlf 2015).

Generalised Procrustes analysis (GPA)
To eliminate the effect of position, size and orientation and to keep only the shape component from individual landmarks recorded, a partial generalised Procrustes analysis (GPA) was performed (Dryden & Mardia 1998), resulting in a set of aligned coordinates.

Principal component analysis (PCA)
The Procrustes-aligned coordinates were projected onto tangent space to the mean shape, and subsequently submitted to principal component analysis (PCA) to summarise the main shape variations in the dataset and also enable the visual exploration of the possible role of factors such as geographic distribution, sex (F, M) and habitat (lagoon, lake, river, spring) in the structure of those variations. Shape variation along PCs was depicted by means of thinplate spline transformation grids (Bookstein 1991). A single outlier individual (from Crete) was removed from subsequent analysis. To evaluate statistical hypotheses concerning patterns of carapace shape covariation with geography, sex and size, a multivariate analysis of variance (MANOVA) was performed on the complete shape space of 24 dimensions.

Linear discriminant analysis (LDA)
Differences in carapace shape in the three studied geographic units, with the Apennine Peninsula and Sicily additionally subdivided into northern and southern/central parts, following the distribution of haplotype groups within the ASAP-MOTU 1 (see results), were afterwards subjected to linear discriminant analysis (LDA; Xanthopoulos et al. 2013) built on individual PC scores. The classification model quality was evaluated with the use of the "leave-one-out cross-validation" (LOOCV) procedure (Hastie et al. 2009). To test possible biases due to that procedure, we performed the analyses Figure 2. Protocol of the 14 landmarks digitised on the lateral view of the carapace: 1 -dorsal posterior carapace margin; 2-5 -four dorsal rostral teeth, starting from the postorbital one; 6 -the tip of the rostrum; 7-8 -two ventral teeth, starting from the posterior one; 9 -orbital margin of the carapace; 10 -branchiostegal spine; 11 -pterygostomial angle; 12 -concavity in the ventral margin of the carapace; 13ventral posterior carapace margin; 14 -concavity in the posterior margin of the carapace.
described by Evin et al. (2013) (see also Navarro et al. 2018 for another application). First, the effect of the number of variables (i.e. numbers of PCs) used in the LDA was tested by computing several LDAs based on an increasing number of PCs (from 1 to 24). Each time, the LOOCV procedure was performed, and the associated prediction error rate was computed. Then, we compared our LDA results with a null model where the classification was due only to chance, so that each individual would have the same chance of being classified in any of the three groups. The null model was simulated by randomly reshuffling individuals among geographic units (simulating a null model where the classification would be due to chance) 100 times. Finally, the quality of the classification process was checked based on the posterior probabilities of each individual to belong to their chosen class. As in Evin et al. (2013), the balance of sampling was tested.

Non-supervised mclust analysis
Model-based Gaussian mixture modelling was used as a non-supervised classification method and density estimation (mclust analysis) to test whether the a priori classification based on geographic units that we used for LDA is the only one structuring the shape variation of our sample, or if other groups exist. For this, we performed an mclust analysis based on finite Gaussian mixture modelling (Scrucca et al. 2016). Modeling of the covariance structure and the number of clusters was optimised based on the Bayesian information criterion (BIC). The search was also carried out over the number of included PCs.

Geometric morphometrics data processing
All morphometric and associated statistical analyses were performed in R v. 3.6.1 (R Core Team 2019) with the use of the following packages: geomorph v.

Morphological identification
Among the 389 individuals analysed we found 192 males, 132 females, 29 ovigerous females and 36 juvenile specimens with sex features not clearly developed (Table II). Based on the combination of the morphological features proposed in the key and species descriptions, all 329 individuals from the Apennine Peninsula and Sicily as well as from Balkan Peninsula are classified as Palaemon antennarius, while all 60 individuals from the Crete Island are assigned to P. minos (Table I).
Nevertheless, when based on a single feature, assignment is not unambiguous (Table III). For example, no Genetic diversity. Out of the 352 COI sequences, the overall nucleotide diversity (pi) in our dataset is 0.023. In total, 88 haplotypes were identified, among which 79 (307 individuals) represent Palaemon antennarius (mtH1-79) and 9 (45 individuals) belong to P. minos (mtH80-88). The BOLD system recognised only one BIN-MOTU (ADI1458) (dx.doi.org/10.5883/BOLD:ADI1458) (but see discussion). On the other hand, the ASAP delimitation method produced three MOTUs (ASAP-MOTU 1-3), each restricted to a particular geographic unit: the Apennine Peninsula and Sicily (APS), the Balkan Peninsula (BP) and Crete Island (CI), respectively (Figures 1 and 3). The nucleotide diversity (pi) and haplotype diversity (Hd) are, respectively, 0.015 and 0.922 for ASAP-MOTU 1, 0.0011 and 0.488 for ASAP-MOTU 2, and 0.0011 and 0.483 for ASAP-MOTU 3. The average K2p genetic distance is 0.037 (SE (standard error) 0.0068) between individuals from APS (ASAP-MOTU 1) and from BP (ASAP-MOTU 2), 0.037 (SE 0.0070) between APS and CI (ASAP-MOTU 3) and 0.021 (SE 0.0056) between BP and CI (Table IV). Furthermore, within APS two main allopatric haplotype groups were found, the northern group and the southern one (Figure 1(a)), occurring allopatrically in the northern/central part of the AP versus the central/southern AP and Sicily, respectively. Only one haplotype (mtH12) of the southern haplogroup is shared between central AP (site 16) and Sicily. Two Sicilian sampling sites (sites 17 and 18, see Figure 1(a)) also have only one haplotype (mtH43) in common.
The 18 individuals (six per ASAP-MOTU) sequenced for H3 belong to six haplotypes. Specimens from APS (ASAP-MOTU 1) are ascribed to four unique haplotypes (nH3-nH6), while the individuals from BP (ASAP-MOTU 2) are represented by nH1 as well as by nH2, the latter haplotype being shared with CI (ASAP-MOTU 3). The K2p distance between APS and BP + CI is very low 0.008 (SE 0.004).

Phylogenetic relationships among haplotypes
The phylogenetic network constructed for COI haplotypes shows three well-distinguished groups congruent with the three ASAP-MOTUs. In the group corresponding with APS (ASAP-MOTU 1), there is also clear separation between the northern and the southern haplogroups (Figure 3(a)). The southern haplogroup appears to be more diverse, possessing a higher number of haplotypes than the northern one. Interestingly, the haplotypes from Sicily do not form a well-defined subgroup, but are intermingled with those from the central Apennine Peninsula.
The phylogenetic tree constructed with the ML method for the H3 marker does not show any supported topology within P. antennarius. However, the existence of two groups is suggested, one composed of the haplotypes from APS and the other of haplotypes from BP and CI, with one haplotype being shared. In addition, a haplotype taken from GenBank, representing an individual from the island of Rhodes, Greece, is phylogenetically close to the haplotypes from BP and CI (Figure 3(b)).

Principal component analysis (PCA)
The main patterns of shape variation depicted on the first two PC axes computed on tangent coordinates illustrate mainly variations located in the rostrum (Figure 4). Thus, PC1, which explains 39% of the initial shape variance in our dataset, describes variation linked mainly to the relative positions of teeth on the ventral side of the rostrum. The second PC, explaining 23% of the initial shape variance, is linked to the rostrum length and width, and is associated with a shift of the ventral spikes relative to the dorsal ones. The eigenvectors of the PCs are presented in Supplementary Table S1.
Incorporating information about geographic units on the first two PC planes (Figure 4) shows that the carapace shape variation seems to be explained mainly by geographical distribution. Indeed, the individuals from Crete form a cloud that almost does not overlap with the clouds grouping the other individuals, mainly with negative scores on the first two PCs, thus exhibiting a teeth-armed rostrum. The individuals from APS and BP form the elongated cloud on the right side of the first PC plane (Figure 4). The studied individuals are characterised mainly by differences in rostrum shape and serration. Individuals from those two peninsulas seem a little less distinct than the individuals from CI in terms of the first two PCs. In contrast, specimens from the APS geographic unit constitute two subclouds, covering the distribution (with a few exceptions) of the northern and southern haplogroups of the ASAP-MOTU 1 in the PCA plot. Moreover, most of the individuals from the more northern part of the APS show generally greater similarity to BP.
To verify whether the shape variation could be connected to the habitat in which the studied shrimp lived, we added information on the habitats in a duplicated plot of the first PCA, which is included in Figure 4. The BP population, in contrast to two other groups (CI only in Kournas Lake; APS almost only in rivers; see Table I), is represented by shrimp inhabiting various types of waters (rivers, lakes, brackish-water lagoons and sublacustrine springs). Within the BP geographic unit, individuals deriving from different habitats are all mixed together, and no specific distributional gradient is shown.
MANOVA, performed on the individual scores from the first 24 PCs (i.e. the complete set of PCs with non-null eigenvalues) and testing the effects of geographic distribution, sex and centroid size on the carapace shape, shows that all of these effects as well as their interactions have significant effects on the carapace shape, even if the geographic factor has the greater effect (Table V).

Linear discriminant analysis (LDA)
No possible bias associated with the LDA or the LOOCV procedure is found for our dataset. First, the prediction error does not increase with the number of PCs included in the analysis; rather, it exhibits an asymptotic decrease ( Figure 5(a), grey curve), quickly reaching an error rate of 2-3%. The error rate of the null models, whatever the number of included PCs, meets the expected error rate of 0.67 ( Figure 5(a), grey confidence interval). Consequently, the LDA is built on the scores from the first 24 PCs. The projections of individuals on the two discriminant axes allow us to distinguish the CI geographic unit from the others on the negative side of the first axis. The individuals deriving from BP and APS form a cloud with a gradient corresponding to the succession of those geographic units, including the division of APS into northern and southern parts ( Figure 5(b)). The proportions of the studied individuals well classified within their own groups by the LDA, with posterior probabilities Figure 4. Ordination of sampling sites projected on the plane defined by the first two principal components (PCs) summarising the main patterns of the carapace shape variation. Thin plate spline (TPS) grids are reported at the extremities of each axis to depict the shape variation along them. Individual symbol shapes stand for the sex factor (circles for females, squares for males), and their colours represent the geographic factor and the designated MOTUs/haplogroups (indicated on the map): yellow -northern APS (Apennine Peninsula and Sicily); red -southern APS (black circles around them indicate individuals from Sicily Island); green -BP (Balkan Peninsula); black -CI (Crete Island). The ordination in the right upper corner corresponds to the same PC1-PC2 plane, but with symbols coloured depending on the environment inhabited by the individuals : purple -river; pink -lake; blue -lagoon; orange -sublacustrine spring.  (Figure 5(c)).

Non-supervised mclust analysis
Non-supervised clustering analysis (mclust) points to EEV (meaning the ellipsoidal, mixture model of the equal shape, equal volume and variable orientation) as the best covariance model, according to which three clusters are consistently detected in the mixture model, until the cost of complexity related to the dimensionality become too high when more than 11 PCs are included (Figure 6(a,b)).
Nevertheless, those units (distinguished based on posterior probability) are partially different from predefined groups corresponding to three investigated populations. Although BP and CI units appear to be clearly qualified, the APS group is partially included in the BP unit. The distribution of individuals assigned to mclust units is illustrated by pie charts in Figure 6(c), and the numbers of individuals from each population in these units are shown in Table VI.

Discussion
Our study shows that when applying the key features proposed by Tzomos and Koukouras (2015) to the studied material, ascribing a specimen to a given taxa was possible, in many cases, only if the features were combined. In those cases, a single key diagnostic character was not sufficient for indisputable distinction between P. antennarius and P. minos (see Results, Table III). For example, the presumably discriminative features, such as the angle of the 5th pleuron end, appeared to show substantial overlap in variation between the two species. However, the outcomes of geometric morphometric analyses confirmed the distinction of P. minos from P. antennarius on the basis of carapace shape variation. What is more, geometric morphometric methods allowed for further subdivision of the individuals ascribed to P. antennarius into two morpho-groups generally fitting the designated geographic units, APS and BP. Also the genetic investigations identified three MOTUs, although the between-MOTU genetic distance did not exceed 0.037, between ASAP-MOTU 1 associated with APS and ASAP-MOTU 2 associated with BP, as well as between ASAP-MOTU 2 and ASAP-MOTU 3 (P. minos, endemic to CI). The distance between ASAP-MOTU 2 and ASAP-MOTU 3 was only up to 0.021. So, interestingly, P. minos appeared to be genetically closer to P. antennarius ASAP-MOTU 2 than the latter was to ASAP-MOTU 1 of the same morphospecies. The fact that all the ASAP-MOTUs belong to only one BIN designated by BOLD is probably associated with the presence of numerous private (i.e. not publicly available) records that are diffusing the barcoding gaps we observe in our dataset. This points to the interest of confronting both genetic and geometric morphometric approaches.
Whilst crustaceans are used as a model system in many different fields of biology, integrative taxonomic investigations of this group remain scarce and their systematics is far from being established. First, some studies combined morphological (but not morphometrics) and genetic approaches. They were mostly used for confirmation or rejection of the morphologically distinguished species, as well as in reference to phylogeography and existing cryptic diversity within nominal species (e.g. Mamos et al. 2016;Rudolph et al. 2018;Rossi et al. 2020;Wattier et al. 2020). Second, some  works combined genetics with morphometrics, but only in terms of body measurement and not analysing shape variation. They were often conducted for decapods and were used to examine the length and width of body parts (e.g. carapace, pleon). One example concerning Palaemonidae is the study by Cartaxana (2015), who on the basis of measurements of the body parts and genetic data suggested the fusion of two Palaemon species. Similar linear distance-based morphometric studies of Palaemon antennarius in north-western Greece by Anastasiadou et al. (2009Anastasiadou et al. ( , 2014Anastasiadou et al. ( , 2017 showed some phenotypic plasticity connected with sex and the type of habitat, although their studies were not supported by genetic analyses. However widely applied in other arthropods (e.g. Ramírez-Sánchez et al. 2016;Lorenz et al. 2017;Ren et al. 2017), geometric morphometrics is much less applied on crustaceans (e.g. Giri & Collins 2004;Accioly et al. 2013;Bissaro et al. 2013;Lovrenčić et al. 2020), and this approach has mainly been used on brachyurans (e.g. Rufino et al. 2004;Silva & Paula 2008;Alencar et al. 2014). For caridean shrimps, such studies are even more scarce, most of them investigating carapace shape variations predominantly in relation to sex (Ashelby 2012;Sganga et al. 2016;De Melo & Masunari 2017) or to habitat (Zimmermann et al. 2012;Torres et al. 2014).
Although geometric morphometric studies combined with genetics have often been conducted on invertebrates, including arthropods (e.g. Silva et al. 2010;Marrone et al. 2013;Zinetti et al. 2013;Marchiori et al. 2014;Kamimura et al. 2020), such an approach has not previously been used on caridean shrimps. Consequently, our study is the first to combine these methods for this group of crustaceans.
Few studies on palaemonid shrimps have tried to relate the carapace shape to environmental factors. While some studies have reported a dependence (e.g. Zimmermann et al. 2012, on Macrobrachium sp.), many others similar to our study have not categorically pointed out such a relation (e.g. Torres et al. (2014) on Macrobrachium sp., or Ashelby (2012) on Palaemon longirostris). Further, our results show that the connection between carapace shape variation and sex is of lower significance than the geographic factor. Sexrelated carapace shape variation was previously evidenced for some palaemonids (Zimmermann et al. 2012;Bissaro et al. 2013;Torres et al. 2014;De Melo & Masunari 2017) but not found in others (Ashelby 2012). This suggests that the nature of such relationships is not obvious and should be studied in detail; however, this is beyond the scope of the present paper.

Overview of geographic distribution patterns integrating geometric morphometry and MOTUs
The spatial distribution of shape variation of Palaemon antennarius and P. minos was illustrated by PCA, LDA and mclust plots. The analyses were congruent in terms of showing that (1) the studied shrimps representing each of the three defined geographic units (APS, BP, CI) were generally characterised by distinct carapace shapes, and each population represented a separate ASAP-MOTU; (2) there was a gradual north-south-oriented variation of carapace shape along AP, which coincides with our finding that no COI haplotypes were shared between the more northern and more southern parts of AP; and (3) on a molecular level, the Sicilian population is part of the southern haplogroup of ASAP-MOTU 1.
Geographical gradient. A geographical gradient of morphological variation seems to be a natural pattern of widely distributed species. In crustaceans in general, however, this phenomenon varies depending on the taxon, ranging from its absence, as in the freshwater amphipod Gammarus balcanicus (Mamos et al. 2014) to a large-scale latitudinal gradient, as in the brachyuran Perisesarma guttatum along the African coast from Kenya to Mozambique (Silva et al. 2010). Specifically for shrimps, our results are consistent with data provided by Ashelby (2012), who studied carapace shape variation in connection to the geographic range of Palaemon longirostris, and confirmed the existence of the northern and southern/central morpho-groups within that species in European coastal waters. The carapace shape variation of the shrimp Macrobrachium borelli was also explained by geographical distribution (Torres et al. 2014).
Focusing on the peri-Adriatic area and adjacent regions, we cannot find much data reporting distributional gradients in relation to shape, regardless of taxon. To our knowledge, the only study involving comprehensive sampling along the Apennine Peninsula and Sicily is that of Pizzo et al. (2011). Although the subject of their interest was two terrestrial beetles of the genus Ontophagus (Scarabaeidae), the pattern of their distribution was quite consistent with our outcomes, pointing to ongoing speciation. Other studies, e.g. by Zaccara et al. (2019) conducted on the cyprinid fish in the southern Apennine Peninsula, showed that body shape variation may be structured according to the hydrographic division of the area. An extensive study of the white-clawed crayfish (Austropotamobius spp.) species complex in Europe, by Jelić et al. (2016), conducted on the basis of molecular analyses, showed a mosaic pattern of distribution of phylogenetic lineages occurring in the Apennine and Balkan Peninsulas. Similarly, an investigation of Austropotamobius torrentium by Lovrenčić et al. (2020) revealed that the reported lineages presented a mosaic composition in northern Croatia, possibly due to the anthropogenic transportation of crayfish among water bodies, which in turn was confirmed by linear and geometric morphometrics. Obviously, in contrast to crayfish and marine shrimps, P. antennarius is not of commercial importance, so human-related translocations probably have much less -if any -impact on the spatial distribution of its genetic variability. A notable exception may be the Cretan population (see below).

Sicilian population
Palaemon antennarius from Sicily Island did not differ from the peninsular individuals in carapace shape. This result is consistent with the work by Pizzo et al. (2011), who on the basis of the same method combined with genetics, questioned the taxonomic status of the beetle Ontophagus massai, endemic to Sicily, and suggested that it belongs to the widely distributed O. fracticornis. The island underwent complex transformations, and due to tectonic movements its geological structure is heterogeneous and consists of three parts of different origins (Broquet 2016;Di Maggio et al. 2017). The Sicilian material used in our study came from the Apenninic-Maghrebian orogen, which was erected after the collision of Southern Apennines and the African plate. So did the material used in Pizzo et al.'s (2011) study. In contrast, Marrone et al. (2013) observed two lineages of the terrestrial tenebrionid beetle, Phaleria bimaculata, but its distribution in Sicily corresponded to the ancient geological division of the island into the African-and European-derived parts. Taking into account the fact that Sicilian shrimps did not differ morphologically from the peninsular ones, it is worth noting that the genetic outcomes of our study reported only one haplotype (mtH12) in common between shrimps collected in central AP (14%) and Sicily (86%). (Figure 1; Table I). Moreover, there was also only one haplotype, mtH43, shared between two sites in Sicily.
Although islands are considered to promote speciation by the isolation of insular populations from the continental ones, and although we observed morphological and genetic divergence between northern/central and southern/central APS geographic units, we also noted the absence of morphological as well as genetic discrepancies within the southern/central APS, which included Sicily Island. This morphological and genetic uniformity of Sicilian shrimp with those of the southern haplogroup on the Apennine Peninsula points to ongoing speciation and colonisation, which can be probably explained by the very recent past land connections between the Apennine Peninsula (Calabria) and Sicily during Pleistocene glaciations and the related eustatic sea level regressions, e.g. during the Last Glacial Maximum (LGM) (Antonioli et al. 2014;Hupało et al. 2021).
Additionally, high haplotype diversity within the southern APS haplogroup indicates two possible scenarios corresponding to colonisation processes. The first presupposes one initial Apennine-Sicilian population subsequently fragmented as a result of the Holocene sea-level rise and subjected to the founder effect and genetic drift, while the second implies multiple colonisation events between Sicily and the southern Apennine Peninsula after they were separated by the Messinian Strait in the Holocene. An airborne dispersal by birds over short distances was evidenced for freshwater shrimp, crayfish and amphipods (Banha & Anastácio 2012;Rachalewski et al. 2013;Águas et al. 2014).

Morphological distinctness in Crete
The palaemonid shrimp occurring on Crete Island were recently described as a distinct species, P. minos, based on traditional morphology (Tzomos & Koukouras 2015). Our study based on geometric morphometrics confirmed this distinctness. Alongside this, taking into account the genetic data, all P. minos individuals were ascribed to a specific ASAP-MOTU: ASAP-MOTU 3. Nevertheless, the low genetic distance between ASAP-MOTU 3 and two other ASAP-MOTUs as well as the H3 haplotype shared between ASAP-MOTU 2 and ASAP-MOTU 3 would suggest a recent or ongoing speciation. Interestingly, the last connection of Crete with the continent was reported from the Miocene, during the Messinian Salinity Crisis (Poulakakis et al. 2015) -thus much earlier than the land connection between the Apennine Peninsula and Sicily. The morphological distinctness but genetic closeness to Balkan population of P. antennarius, as well as the absence of a recent land connection to Crete, would suggest an alternative explanation. It is possible that crustaceans could travel between the Greek mainland and islands using water birds as vectors, as has been demonstrated over short distances by Banha and Anastácio (2012), Rachalewski et al. (2013) and Águas et al. (2014). Or, more likely, they could be introduced to new sites by humans (e.g. fisherman or sailors), which has 918 A. Jabłońska et al. been suggested by Hupało et al. (2020), who reported Gammarus plaitisi, previously considered to be a Cretan endemic (Hupało et al. 2018), from other Aegean islands of Tinos and Serifos, which are equally far from Crete as the mainland. Jesse et al. (2011), based on genetic methods, suggested distinguishing the Cretan endemic Potamon kretaion from P. potamios (Decapoda, Brachyura). However, they also pointed to the recent divergence of these crabs, which strengthens the theory of possible alternatives to land connections.

Conclusions
Our study provides the first evidence of shape differentiation of freshwater palaemonid shrimps in the Mediterranean region based on geometric morphometrics. Moreover, we joined the taxonomic method with COI DNA sequences, a technique that has not been used in Palaemon antennarius and P. minos studies to date. This integrative approach gave us partly conflicting results, which highlighted the need to combine taxonomic methods in species research and provided an approximately holistic view on the given species. Summarising our results, the geometric morphometrics and the COI genetic analysis show two different patterns of divergence on the two investigated islands. Shrimps from Sicily were neither morphologically nor genetically distinct from those of the Apennine Peninsula, while the shrimps deriving from Crete were morphologically and molecularly different from the Balkan ones. Yet, given that P. minos (ASAP-MOTU 3) is genetically closer to the Balkan ASAP-MOTU 2 than the latter is to the Apennine ASAP-MOTU 1, for taxonomic consistency, the MOTU inhabiting the Balkan Peninsula should be either described as a new separate species or synonymised with P. minos. The third possible option would be treating all the populations as part of P. antennarius. Geometric morphometrics supports the first option, phylogenetic reconstructions point to the second one, while the fact that BOLD recognises only one BIN favours the third one. In our opinion, further studies are required to resolve the issue. These studies should employ more molecular markers and crossing experiments as well as more thorough sampling over the area, and should explore the spatial distribution of morphological traits at finer scales, e.g. to explore the potential issue of spatial autocorrelation.
Our study points to the complexity of the species delimitation process and therefore confirms that taxonomy should be based on an integrative approach.
It also shows that doubts and questions concerning the taxonomy and distribution of species in freshwater faunas remain to be answered, even with shrimps, which are quite emblematic representatives of these faunas. Kłosowska, Paula Krzywoźniak, Grzegorz Michoński, Wanda and Manolis Plaitis, Agnieszka Szlauer-Łukaszewska and Anna Wysocka for their help during fieldwork. We also thank Radek Šanda, who provided part of the material used in this study.

Funding
The fieldwork and molecular part of the study were financially supported by the National Science

Disclosure statement
No potential conflict of interest was reported by the authors.

Supplementary material
Supplemental data for this article can be accessed here.