An integrative insight into the diversity, distribution, and biogeography of the freshwater endemic clade of the Ponticola syrman group (Teleostei: Gobiidae) in the Caucasus biodiversity hotspot

Abstract Freshwater habitats of the Caucasus biodiversity hotspot represent a center of endemism for the gobiid genus Ponticola Iljin, 1927. Hitherto, large‐scale molecular studies, owing to restricted taxon and geographical sampling, have failed to give an elaborate picture of diversity and evolutionary history of these species. Here, to contribute to filling this gap, we assessed taxonomic diversity, phylogeography and evolutionary history for the south Caspian populations of Ponticola presently classified as P. iranicus and P. patimari, using an integrative taxonomic approach comprising an entire geographic range sampling, and analyses of mitochondrial DNA haplotypes, the head lateral line system, otolith shape, and meristic and morphometric variation. All freshwater samples of the P. syrman group belong to a monophyletic clade with two main subclades: a small subclade confined to the upper Sefidroud sub‐basin including the type locality of P. iranicus and a large subclade with three geographically constrained haplogroups (Hg1, Hg2, and Hg3), comprising the rest of the distribution. Hg1 showed an eastern distribution including the type locality of P. patimari, while Hg2 and Hg3 are sister groups with central and western‐central distributions, respectively. The freshwater clade diverged from P. syrman during the Tyurkyanian low stand (~150 m b.s.l. lasting ~0.1 Myr), while the divergence of P. iranicus and P. patimari and radiations within P. patimari took place during the Bakunian high stand (up to 50 m a.s.l. lasting ~378–480 kya). Species delimitation analyses indicated two distinct species, corresponding to each main subclade. Although the otolith shape and lateral line analyses did not reflect with phylogeographic pattern, PCA and DFA plots of meristic and morphometric data showed a clear separation of the two major subclades corresponding to P. iranicus and P. patimari, suggesting the presence of significant morphological variation meriting formal taxonomic recognition. Overall, our findings (i) reveal the presence of two freshwater endemic species in the P. syrman group, and pending further investigation, hypothesize the presence of a third cryptic species; (ii) revise and document a narrow distributional range and low diversity for P. iranicus, in contrast to a wider distributional range and high diversity for P. patimari; (iii) suggest that the climatic oscillations of the Pleistocene were associated with the cladogenesis within the P. syrman group; and (iv) allowed for the recognition of conservation units and proposition of management measures.


| DNA extraction and PCR
Genomic DNA was extracted from fin samples preserved in 96% ethanol using a salt method protocol following Bruford et al. (1992). PCR amplification of the standard vertebrate DNA barcode region, the mitochondrial cytochrome c oxidase I (COI) was performed using the primer pairs, FishF1 and FishR1 (Ward et al., 2005), or FISH-BCL and FISH-BCH (Baldwin et al., 2009). The 25 μl PCR reaction mixes included 12.5 μl of a ready 2X Taq PCR Master Mix (Parstous TM ), 0.5 μl of each primer (10 pmol/μl), 6 μl of the target DNA, and 5.5 μl DNasefree distilled water. The PCR conditions were initial denaturation for 3 min at 94°C followed by 35 cycles of 94°C for 45 s, 52-56°C for 45 s, and 72°C for 45 s, and final extension for 5 min at 72°C. PCR products were purified using ExoSAP-IT TM and sequenced by the Niagene Lab. (Tehran, Iran) using Applied Biosystems TM BigDye TM Terminator v3.1 Cycle Sequencing Kit on an Applied Biosystems TM ABI PRISM 3730xl.
In addition to the newly sequenced material, our molecular investigation also analyzed the following COI datasets (Tables 1 and S1): (i) P. iranicus COI barcodes mined from Zarei, Esmaeili, Schliewen, et al. (2021b) Stepien, 2009a, andZarei, Esmaeili, Schliewen, et al., 2021b] and ZFMK (Table S1). In general, COI does not give enough resolution for deep phylogenetic analyses, but for the taxonomic scope presented here it provides sufficient phylogenetic resolution for a first step for integrative taxonomy.
All newly obtained sequences are deposited in GenBank (Table 1).
The saturation test of Xia et al. (2003) in DAMBE 7 (Xia, 2018) was used to test the nucleotide substitution saturation in the sequence under study. The best-fit nucleotide substitution model for the dataset was estimated based on the Bayesian information criterion (BIC) in jModelTest 2.1.3 (Darriba et al., 2012). A Maximum Likelihood (ML) phylogeny using 5000 bootstrap replicates (fast bootstrap) was generated in RAxML 7.2.5 (Stamatakis, 2006).
Furthermore, cutoff value of 2% K2P sequence divergence for COI (Ward, 2009) was used as indicator of distinct species.

F I G U R E 1
The distribution ranges of P. iranicus according to Vasil'eva et al. (2015) and Zarei, Esmaeili, Schliewen, et al. (2021b) (green area), P. patimari according to Eagderi et al. (2020) (pink area), and P. syrman according to Zarei, Esmaeili, Kovačić, et al. (2022a) and Pinchuk et al. (2003b). The numbered circles refer to the locations of samples used in this study. The map was originally designed using the subbasins layer of HydroBASINS 1.0 (Lehner & Grill, 2013) in DIVA-GIS 7.5 and Surfer 11. TA B L E 1 Studied samples of P. iranicus and P. patimari (see the initial and revised classifications) from the SCB in the molecular (COI), otolith outline, meristic, morphometric and head lateral line analyses  Stepien, 2009a). A secondary calibration of 0.9 Mya was assigned to the node subtending P. syrman + P. iranicus + P. patimari (based on Zarei, Esmaeili, Schliewen, et al., 2021b) and the analysis was run in four independent runs of 100,000,000 generations, with trees sampled every 1000 generations; the first 10% were discarded as burn-in. Finally, convergence and effective sampling sizes (ESS) were checked in Tracer 1.6 (Rambaut et al., 2018), and a maximum clade credibility consensus tree was built in Tree Annotator 1.8.2 (Drummond et al., 2012).

| Otolith outline analysis
The left sagittal otolith of 213 specimens from seven south Caspian localities ( Table 1) were extracted using fine tweezers, cleaned and washed in distilled water, dried at room temperature, then placed on a dark plate and digital images at 4× magnification were taken using a 14MP Industrial Microscope Camera 180× equipped with an S-EYE 1.2.4.128 image processing system. To measure otolith shape variation, we used statistical functions in R 4.0.5 (Ihaka & Gentleman, 1996) using the packages shapeR 0.1-5 , vegan 2.5-7 (Oksanen et al., 2020), ipred 0.9-12 (Peters et al., 2021), and MASS 7.3-54 . The otolith images were read into R. ShapeR analyzes otolith shape by extracting otolith outlines from these images. A matrix of coordinates (x, y) from all otolith outlines was calculated. Evenly spaced radii with length as a univariate shape descriptor were drawn from the otolith centroid to its outline. Using the wavelet transformation on radii, the wavelet coefficients were extracted with wavethresh 4.6.8 (Nason, 2016). In order to remove the influence of allometric growth, wavelet coefficients that showed significant (p < .05) interaction between samples and SL were excluded , and the remaining coefficients were imported into the R statistical packages. The mean otolith shapes for samples were plotted based on wavelet coefficients. To determine areas of otolith shape variation, mean shape coefficients and their standard deviations were plotted against the angle of the otolith outline using wavelet transform with gplots 3.1.1 (Warnes et al., 2020). Because the proportion of variation among groups (intraclass correlation) was more informative to measure between-population differences, intraclass correlation was estimated along the otolith outline.
The radii length was employed to test for the significance of differences between samples using an ANOVA-like permutation test (1000 permutations) in vegan. To investigate the statistical significance of otolith shape difference between males and females, sex-related stability of otolith shape was analyzed in shapeR.

| Meristic and morphometric analysis
Morphometric methods represent a composite of Schliewen and Kovačić (2008) and Miller (2003). Eight meristic and 41 morphometric variables were recorded for 107 specimens from six south Caspian localities ( To study the distribution of the measured specimens in a multivariate morphometric space, discriminant function analysis (DFA) was performed on the size-adjusted morphometric and meristic data to extract the most important characters for differentiating lineages, using the F-value criterion. DF scores 1 and 2 were selected and plotted in the discriminant space. Also, group centroids were produced by DFA to visualize associations between lineages.
Classification success was taken from DFA to correct assignment of individual fish to original lineages. The size-adjusted morphometric and meristic data were used as input for a dendrogram based on the Squared Euclidean distance as a measure of dissimilarity.

| Head canals and sensory papillae pattern
The sensory canals, pores, and sensory papillae patterns of the head lateral line system provide important diagnostic characters to distinguish goby species (Kovačić, 2020). The head lateral line systems of the samples were checked and mapped using a Zeiss TM Stemi SV6 stereomicroscope. When necessary, the specimens were stained in 2% KMnO 4 solution for 5 s, which allowed a better examination of sensory papillae rows. Drawing of the head lateral line system was prepared with CorelDRAW® Graphics Suite 2020 (https://www. corel draw.com/). The terminology of papillae rows, and of canal pores follows Miller (1986) based on Sanzo (1911) and Miller (2003).

| Genealogy, species delimitation, phylogeography, and divergence time estimation
GTR+I+G was determined as the best fitting nucleotide substitution model for the COI dataset (108 sequences, 648 bp length). The nucleotide substitution pattern based on Xia et al. (2003) method showed that the sequences have not reached substitution saturation and are suitable for phylogenetic analysis (Iss < Iss.cS and Iss < Iss. cA; Table S3). In total, the resulting ML tree included COI sequences for 95 samples from the SCB initially classified as Ponticola iranicus (N = 77) and P. patimari (N = 18) based on the geographical origin of the samples and current literature (Table 1; Figure 1), combined with 13 archived sequences (Table S1) from another twelve Ponticola species from the Ponto-Caspian basin ( Figure 3). The deepest split in the tree is between P. ratan and the remainder of Ponticola species with two fully resolved clades. One clade comprising P. syrman and all individuals from the SCB initially tagged with the names P. iranicus and P. patimari, and a second clade comprising the remainder of the species. Within the clade corresponding to the P. syrman group, P. syrman is sister to the clade containing all the south Caspian individuals initially assigned to P. iranicus and P. patimari. This P. iranicus + P. patimari clade divides into two subclades with a 2.32% average K2P divergence (Figures 3 and 4). The SP and ASAP (the partition with the highest ASAP-score) species delimitation analyses found support for two putative species in the freshwater clade of the P. syrman group: one species included all specimens from Sefidroud at IM, and another species including the rest of the distribution (Hg1 + Hg2 + Hg3). The partition with the second best ASAP-score found support for three putative species corresponding to (i) all specimens from Sefidroud at IM, (ii) Hg1, and (iii) Hg2 + Hg3. In addition, if the K2P distance species threshold value (≥2%; Ward, 2009) is accepted, IM, and Hg1 + Hg2 + Hg3 would represent two separate species ( Table 2).
A haplotype network of the south Caspian samples is shown in Figure 4a. Similar to the ML tree, two main subclades separated by 11 fixed mutational steps are evident. The 14 samples from Sefidroud at IM corresponding to the smaller subclade on Figure 3 F I G U R E 3 ML phylogeny of Ponticola based on mitochondrial COI (108 individuals), with focus on the P. syrman group. New sequences produced in the context of this study are shown within the colored clades. The new limits set for P. iranicus and P. patimari on this tree are based on this phylogeny and the results derived from the species delimitation analyses (i.e., SP, ASAP, and 2% K2P new species threshold; the initial assignment of individuals and samples are presented in the text and Table 1 and Hg2 are deviated from Hg3 by six and two mutational steps, respectively ( Figure 4a).
Based on this phylogenetic, species delimitation and phylogeographic framework, and based on the type localities of P. iranicus and P. patimari, we revise the presently known distributional ranges of the two species (Figure 4b), which is a prerequisite for the following analyses. Ponticola iranicus is distributed in the upper Sefidroud subbasin including IM, while the rest of distribution (i.e., KA, SH, SI, KY, PO, TO, CH, NO, and KH) mitochondrially belongs to P. patimari: Hg1 is distributed in the eastern part, Hg2 is distributed in the central part, and Hg3 with a western-central distribution.
F I G U R E 4 (a) MJ haplotype network for the mitochondrial haplotypes observed in P. iranicus and P. patimari. The circle area is proportional to the frequency of each haplotype. Each color indicates a sampling site (as Figure 1), each soft line connecting neighboring haplotypes represents a single mutational step, and the lines with hatch marks/numbers indicate multiple mutational steps between sampled haplotypes. Small black circles correspond to missing/hypothetical haplotypes. (b) Revised distributional ranges for P. iranicus and P. patimari, and spatial distributions of the main P. patimari lineages (i.e., Hg1, Hg2, and Hg3). The left star and the right star on the map represent the type localities of P. iranicus (upper Sefidroud sub-basin, Tutkabon Stream, 36°50.756′N, 49°35.021′E) and P. patimari (Kheirud, 36°36'46.0"N 51°34'03.0"E), respectively.
TA B L E 2 Average K2P genetic distances (%, below the diagonal) between the freshwater lineages of the P. syrman group based on mitochondrial COI (see Figure 3). Abbreviation: SD, standard deviation (above the diagonal).

| Population differentiation and genetic structure within P. patimari
Excluding TO, NO, KH, and KA due to small sample sizes, the range of pairwise Fst values between the P. patimari samples was moderate to high and statistically significant (p < .001), ranging from 0.365 between two geographically close samples KY/PO, to 0.913 between two geographically distant samples SH/CH (Table 4). Since the P. patimari samples were distributed along an east-west axis, we expected geographic isolation to contribute to genetic affinities.
Concordantly, when all P. patimari localities were analyzed, an overall significant positive correlation was found between geographical and genetic distance matrices (r = .71, p < .001). Furthermore, a second Mantel test without taking into account the samples with small sample size (i.e., TO, NO, KH, and KA), indicated that the positive correlation was still high and significant (r = .86, p < .001).
AMOVA for classification of all P. patimari localities based on three geographic regions (i.e., central, western, and eastern) revealed that 53.23% of total variation was significantly correlated with differences among the geographical regions, whereas interand intra-sample differences explained 25.56% and 21.21% of the variation, respectively (Table 5). When samples were classified according to their sub-basins as groups (i.e., Talesh, Anzali, Sefidroud, and Mazandaran), percentage of variance (0.45%) was lower than among samples within groups (75.63%), and not significant. A second regional-based AMOVA for the populations with large sample sizes (i.e., SH, SI, KY, PO, and CH) also indicated that a regional classification of populations (among groups variation: 49.01%) for the TA B L E 3 Basic parameters of genetic diversity for mitochondrial COI in the studied samples/lineages of P. iranicus and P. patimari based on the revised distributional ranges in this study

| Historical demography
TO, NO, KH, and KA were excluded from the locality-level analyses.
Estimates of Tajima's D were negative for IM (P. iranicus) and PO, while positive for CH, SH, SI, and KY, however, none of these values were statistically significant (p > .5; for the pooled data set of P. patimari (Hg1 + Hg2 + Hg3) revealed a multimodal mismatch distribution ( Figure 6). In haplogroup-level analysis, the non-unimodal MMDs also rejected the model of sudden expansion for the P. patimari haplogroups. Accordingly, because TA B L E 5 AMOVA for P. patimari (revised distribution) based on geographic region and sub-basin classification of samples  q, u, trp, y, as 1 -as 3 ) and four longitudinal (x 1 , x 2 , la 1 -la 2 ) rows including the axillary series; row x 1 long, parallel to above AOC, reaching to trp; row x 2 posterior to rows x 1 and trp, above posterior third of opercle and above pore τ; row z upward from the PC dorsal end, ventrally starting close to, but not exceeding pore γ; row y behind pore τ and below row x 2 ; three short transverse rows in oculoscapular groove between pores ρ and θ: row q anteriormost, close to pore ρ, extending ventrally of the oculoscapular groove and sometimes with one or two papillae dorsal of pore ρ; posteriormost row trp close to pore θ, extending dorsally and passing upwards the level of row 3.5.2 | Ponticola iranicus from IM, near to the species type locality, and P. patimari Hg2-Hg3 Head lateral-line patterns in P. iranicus (Figure 7c,d) from IM, and in other P. patimari haplogroups (i.e., Hg2 and Hg3; Figure 8) are quite similar to that of P. patimari from the species type locality (i.e., KH: Hg1). AOC, POC, and PC canals present, with pores σ, λ, κ, ω, α, β, and ρ; θ and τ; and γ, δ, and ε, respectively. Suborbital infraorbital neuromast organs in seven transverse rows, four rows (1-4) before and three rows (5s, 6s and 7) above and two rows (5i, 6i) below row b and row a absent, row 7 consists of a few to several papillae descending from before anterior oculoscapular pore α, rows 5i and 6i separated, with row 5i well behind anterior end of row b and row 6i close to posterior end of row b; dorsal rows o separated in dorsal midline; oculoscapular row z ends near to pore γ; row x 1 ending anteriorly behind pore β.

| General morphology of sagittae otolith inner face
The otoliths of P. iranicus from IM and P. patimari from NO (samples near to the species type localities) have parallelogram shapes, with marked preventral and posterodorsal projections (Figure 9).

| Otolith outline differences
A total of 213 specimens belonging to six samples (i.e., IM, CH, KY, PO, SH, and SI) were analyzed for their left sagittal otolith outline.
According to the ANOVA-like permutation test, the lengths of the three major radii were significantly (p < .01) different between all the P. iranicus (i.e., IM) and P. patimari (i.e., CH, KY, PO, SH, and SI) samples (Table 8).
The mean otolith outline based on the wavelet coefficients   Boxplots of canonical score distances with respect to variation among samples based on the wavelet coefficient indicate shape differences between SH and the other samples, but no significant differences observed among all of the other samples ( Figure 10b). A cluster analysis dendogram based on CAP1 and CAP2 and using the Euclidean distance (Table S4)  showing a clear differences between SH and the others (Figure 10d).
The overall classification success with a leave-one-out crossvalidation estimation was 43.7%, the highest classification success was achieved for SH (71.7%), followed by KY (54.5%) and IM (47.4%; Table 9). The Mantel test confirmed that there was no correlation between the otolith outline variation and the geographical distances among samples (r = .37, p > .05). The levels of otolith outline resemblance between samples were somewhat dependent on genetic distance (r = .54, p = .019). and PC 2 scores showed a better separation of P. iranicus from both P. patimari Hg1 and P. patimari Hg2 + Hg3 (Figure 11a). PCA for significant morphometric variables extracted from variance-covariance matrix revealed that PC 1, PC 2, and PC 3 accounted for 28.59%, 16.33%, and 9.57% of the variation, respectively, which gave a cumulative variation of 54.5% for the first three PCs. In order of importance, the most significant variables loading on PC 1 and PC 2 were D1III, HL, E, PO, D1I, Cl, TL, D1II, A1I, Pl, ULl, D2I, Vl, SN/D1, D2h
For the meristic variables, the discriminant function test using Wilks' λ statistic was significant (λ = 0.41, χ 2 = 67.85, p < .0001), indicating a relatively high degree of inter lineages variance and that the means of the discriminant scores for the three lineages are different. DFA produced two discriminant functions for the meristic variables, DF 1 accounted for 81.6% and DF 2 accounted for 18.4% of among group variability. Biplot of DFA showed a separation between P. iranicus and both P. patimari Hg1 and P. patimari Hg2 + Hg3  showed morphological separation between the three lineages, especially a clear separation between P. iranicus and both P. patimari Hg1 and P. patimari Hg2 + Hg3 (Figure 11d).
Functions at group centroids analysis using the size adjusted variables revealed an average of 5.77 for P. iranicus, and −2.84 and −1.04 for P. patimari Hg1 and Ponticola patimari Hg2 + Hg3, respectively, suggesting a clear differentiation between P. iranicus and both P. patimari Hg1 and P. patimari Hg2 + Hg3. Classification success rates were estimated for the three lineages (Table S7). The proportion of individuals correctly classified into their original lineages was 100%, 100%, and 96% in P. iranicus, P. patimari Hg1, and P. patimari Hg2 + Hg3 (4% classified into P. patimari Hg1), respectively, indicating a high rate of correct classification of specimens. Furthermore, cluster analysis for the meristic and morphometric measurements resulted in a dendrogram with two clusters: P. iranicus in one, and the two P. patimari lineages in the other (Figure 11e).

| Species diversity and distributions
Our molecular results based on mitochondrial COI highlighted several main features, providing support for two or three freshwater species in the P. syrman group. All freshwater samples of the  P. syrman group belong to a monophyletic clade with two main subclades: (i) a smaller subclade comprising all individuals collected from the upper Sefidroud sub-basin, where P. iranicus is also described (Vasil'eva et al., 2015), and (ii) a larger subclade with three geographically structured haplogroups (i.e., Hg1, Hg2, and Hg3), which includes specimens from the rest of the distribution. Hg1 contains all the specimens collected from the eastern localities, where the type locality P. patimari is located (Eagderi et al., 2020); Hg2 and Hg3 are sister groups with central and western-central distributions, respectively. The two main subclades are separated by 12 mutational steps and a 2.32% K2P genetic distance in COI, while the three haplogroups within the larger subclade are separated from one another by 2-6 mutational steps and 1.19%-1.91% K2P genetic divergence, i.e., less than the conventional threshold proposed for COI in fish (i.e., 2%; Ward, 2009) used as an indicator of distinct species. Genetic divergences between the three haplogroups of the larger subclade and between Hg1 and Hg2 + Hg3 are significantly smaller than those between most species of Ponticola (Eagderi et al., 2020;Zarei, Esmaeili, Kovačić, et al., 2022a), except for the one between P. kessleri and P. eurycephalus (1.2%), and those between several other freshwater endemic species: P. platyrostris and P. constructor (1.2%), P. platyrostris and P. rhodioni (1.5%), TA B L E 1 0 MANOVA for differences among the three main lineages (P. iranicus, P. patimari Hg1, and P. patimari Hg2 + Hg3) based four tests. P. constructor and P. rhodioni (1.8%), and P. rizensis and P. turani

F I G U R E 11
(1.5%). The SP and ASAP delimitation analyses indicated two distinct species, corresponding to each major subclade. The partition with the second best ASAP-score however, found support for three hypothetical species, corresponding to: (i) the smaller subclade from upper Sefidroud, (ii) Hg1, and (iii) Hg2 + Hg3. Nevertheless, a consensus among molecular, distributional, and morphological data, i.e., integration by congruence (Padial et al., 2010), is often used as the major argument for establishing species delineation (Kekkonen & Hebert, 2014;Puillandre et al., 2021).
On page 29, Eagderi et al. (2020) state that "P. patimari is distinguished from P. iranicus by having six suborbital transverse rows (vs. seven rows)", however, in the extended description of P. patimari on page 27, they confirm the presence of row 7, incorrectly labeled as row 6. In addition, based on the original description, one may infer that P. patimari is also distinguished from P. iranicus by absence of anterior oculoscapular pores ω and preopercular pores γ and ε (vs. present). Based on the examination of types and newly collected materials in this study, it became clear now the cephalic lateral line system of P. patimari was incorrectly described; we redescribed it in detail based on corrected morphological character states: infraorbital neuromast organs in seven transverse rows, four (1-4) before and three (5s, 6s and 7) above and two (5i, 6i) below row b; AOC, POC and PC canals present, with pores σ, λ, κ, ω, α, β, ρ and θ, τ and γ, δ and ε, respectively. This is also the typical pattern of head lateral line system of Ponticola species, although many of these species still miss head lateral-line system description and P. syrman has three transverse rows below hyomandibular row b, third one marked 5i' as duplication of row 5i in Miller (2003). This head lateral line system pattern is also common among all Neogobiini + Ponticolini clade (Zarei, Esmaeili, Schliewen, et al., 2021b) with some exceptions: Mesogobius spp.
and Neogobius bathybius have infraorbital neuromast organs in 8-10 and 8 transverse rows, respectively. Intraspecific deviations in head lateral line pattern occur and are common among the Ponto-Caspian gobiids, which therefore should not be taken as diagnostic features to describe new species. These deviations concern the following character states: (i) addition or deletion of pores in the course of canals [e.g., additional pores in the course of AOC or loss of pore δ in some specimens of P. kessleri (Ahnelt & Duchkowitsch, 2001), or presence of an additional pore between pores κ and ω in some specimens of P. gorlap from Sefidroud (Zarei, 2022)]; (ii) paired interorbital pores [e.g., paired pores λ in some specimens of P. iranicus from the upper Sefidroud sub-basin (this study) or in some specimens of P. kessleri and Proterorhinus sp. (Ahnelt & Duchkowitsch, 2001)]; and (iii) presence of additional row before or below row b on one side of the head [e.g., in some specimens of P. ratan from Odessa (Pinchuk et al., 2003a), P. hircaniaensis (Zarei, Esmaeili, Kovačić, et al., 2022a), and P. patimari Hg3 from Shafaroud (this study)]. It was already found that between sister species of Ponticola genus there is no recognizable differences in the head lateral-line system (Kovačić & Engin, 2008;Zarei, Esmaeili, Schliewen, et al., 2021b), and even the recent descriptions of new cryptic species among the entire Gobius-lineage struggled to find differences to sister species or to similar congeneric species in the head lateral-line system (e.g., Kovačić & Šanda, 2016).
Furthermore, based on our data, ranges for D2 branched rays (15-16, usually 16 vs. 15-17, usually  it has been pointed out that for the sister species or similar congeneric species supported by molecular data, additional effort should be invested in a search for useful diagnostic characters that can distinguish the species (Kovačić et al., 2021;Kovačić & Šanda, 2016).
The morphological results only partly reinforced the indicated molecular phylogenetic subdivisions: a. The cephalic lateral line system pattern did not show interspecies differences supporting the results of mtDNA-based phylogenetic work. As discussed above, it is already known that among gobiid species from large genera of Gobius-lineage, like Gobius or Ponticola, some species share cephalic lateral line system so similar that it is hard or it is not possible to find the differences between species and express them as diagnostic characters (see species short descriptions and cephalic lateral line system pattern illustrations in Miller, 1986).
b. The otolith shape analysis did not show congruence with the results of phylogenetic pattern. Otolith analysis has made significant contributions to the understanding of systematics and evolution of teleosts (e.g., Lombarte et al., 2018;Reichenbacher et al., 2007;Teimori et al., 2019). and Neogobius) and at higher taxonomic levels (Zarei, 2022). On the other hand, otolith outline analysis using wavelet transform showed notable phenotypic variation among samples not reflecting major phylogenetic subdivisions and spatial distributions of the main lineages. This results can be explained by the fact that the morphological attributes apart from certain intrinsic (phylo-and population genetic) factors may be influenced by additional extrinsic, e.g., environmental, effects (Charmantier & Garant, 2005;Klingenberg, 2010;Landaeta et al., 2019).
c. The other qualitative characters, including coloration that has high intraspecific variability in these species (Figure 12) Accordingly, the present results revised distributional ranges of P. iranicus and P. patimari, narrowing known distributional range of P. iranicus, and expanding distributional range of P. patimari (Eagderi et al., 2020;Vasil'eva et al., 2015; Figures 1 and 4b).

| Evolutionary history and genetic structure
Multiple glacial-interglacial cycles affected the Caspian Sea during the Pleistocene, causing periods of transgressions and regressions and extreme salinity oscillations (Chepalyga, 1984;Dumont, 1998;Zenkevich, 1963 (Klenova et al., 1962;Krijgsman et al., 2019;Svitoch, 2012). During this phase of high salinity and reduced habitat availability, an arid climate prevailed in the south Caspian region (Abdullayev, 2000;Tagieva et al., 2013;Zastrozhnov et al., 2020), and rivers draining into the Caspian Sea may have become seasonal or dry for short periods (Sands et al., 2019). According to our molecular clock based estimates, brackish water ancestors of the contemporary P. iranicus + P. patimari diverged from P. syrman ca. 0.9 Mya: The ancestor of P. iranicus + P. patimari adapted to freshwater, while "marine" P. syrman adapted to increased salinity of the Caspian Sea in the Tyurkyanian period, so the ecological gap and isolation between them also increased. A similar effect has been outlined for evolution of several other aquatic species in the south Caspian basin (e.g., Sands et al., 2019), including the radiation of Ponticola hircaniaensis (Zarei, Esmaeili, Kovačić, et al., 2022a) and a possible bottleneck in Proterorhinus nasalis (Zarei, Esmaeili, Schliewen, & Abbasi, 2022b (Veith et al., 2003), and the freshwater crab Potamon ibericum (de Bieberstein, 1808) (Parvizi et al., 2018).
The Tyurkyanian regression was followed by the Bakunian sealevel high stand up to about 50 m a.s.l. and a temporal extent of ca.
378-480 kya (Krijgsman et al., 2019), which was the largest transgression in the Pleistocene history of the Caspian Sea. The Bakunian Sea became considerably desalinated, and a warmer and humid climate persisted during this period (Abramova, 1974;Filippova, 1997;Svitoch et al., 1998;Yakhemovich et al., 1986). It is likely that the Bakunian increase in climatically suitable habitats as well as in precipitation and freshwater flow had caused habitat connectivity during this period, leading to downstream expansion into the eastern and western sub-basins and to a further radiation within the freshwater clade; according to this scenario, (i) first, separation of the pre-iranicus and pre-patimari clades in the early Bakunian; the population that remained the headwaters of the Sefidroud sub-basin eventually evolved into P. iranicus, (ii) divergence of the pre-patimari Hg1 and pre-patimari Hg2 + Hg3 took place in the middle Bakunian, and (iii) separation of the pre-patimari Hg2 and pre-patimari Hg3 developed in the late Bakunian. Subsequently, the Hg2 and Hg3 groups underwent further expansion and radiation within the central and western sub-basins during the early Khazarian stage. Our phylogeographic assessments suggest a transitional zone between these two groups occurring in the region of the lower Sefidroud River (i.e., KY). Therefore, we propose that Hg2 and Hg3 mainly occur east and west of the Sefidroud River, respectively. Finally, during the early Khazarian phase, subsequent diversifications of P. iranicus and P. patimari in the colonized rivers took place, which produced a bush-like genealogy characterised by large numbers of private haplotypes and rare cases of shared haplotypes (e.g., H16 in P. patimari). This pattern suggests that contemporary connectivity and gene flow between catchments (compared with the Bakunian stage) is limited.

| Conservation units and management propositions
The south Caspian freshwater habitats have been and continue to be threatened by the main causes of biodiversity loss (Esmaeili et al., 2015;Mousavi-Sabet, 2021). In this context, a careful and accurate taxonomic assessment of their ichthyofauna is a prerequisite for sketching effective conservation measures. The taxonomic revision of the P. syrman group has several implications for conservation: tions], so that limited conservation resources can be utilized optimally (Ryder, 1986).
Reciprocal monophyly in mtDNA, high estimates of genetic divergence, lack of gene flow, distinct spatial distributions, and morphometric differences suggest long-term isolation between the three genetic clusters, P. iranicus, P. patimari Hg1, and P. patimari Hg2 + Hg3, which are sufficiently differentiated to be considered ESUs.
ESUs may themselves be sub-structured at hierarchical levels generally recognized as MUs (sensu Moritz, 1994), another logical units for conservation efforts (Moritz, 1994). The P. patimari Hg2 + Hg3 ESU is comprised of two geographically structured subgroups, Hg2 and Hg3, with a contact zone in the lower Sefidroud sub-basin. The molecular and morphometric data indicate that individuals from the Hg2 and Hg3 subgroups comprise just one biological group; however, this P. patimari Hg2 + Hg3 ESU is best viewed as comprising two MUs, Hg2 and Hg3, according to the original definition of Moritz (1994). Ponticola patimari has a high level of genetic diversity if treated as only one panmictic group, however, species structured into ESUs and MUs, necessitate separate monitoring and management of each unit (Frankham et al., 2002).
The three ESUs are each subject to different and distinct threats stemming from differences between the area occupied by each lineage and levels of habitat degradation within these areas; thus, they are natural candidates for differentially optimized management strategies. Ponticola iranicus is particularly vulnerable as it is composed of a single ESU with a narrow distributional range associated with low mtDNA variability which makes it more sensitive to extrinsic changes and therefore a greater risk of extinction (Frankham et al., 2002). The Sefidroud Dam (= Manjil Dam), constructed upstream of the Sefidroud River to store water for irrigation and produce hydroelectric power, is the main threat to P. iranicus. The Ponticola patimari Hg2 + Hg3 ESU features a broad geographical distribution from PO to KA. From a conservation perspective, this ESU possesses relatively high levels of genetic diversity indicating that the population does not require immediate conservation attention.
In addition, the distribution range of this ESU harbours many available habitats for the species to occupy. The P. patimari Hg1 ESU, despite having wider distributional range, shows levels of mtDNA variability similar to the P. iranicus, most likely biased by poor sampling at TO, NO, and KH.
In Iran, temperature will rise and precipitation will decrease under the climate change (Dastorani & Poormohammadi, 2016;Rahimi et al., 2018). These changes will affect freshwater habitats and will make them less habitable for freshwater species. Modeling the impacts of climate change on distribution of P. iranicus sensu Vasil'eva et al. (2015) predicts that the species will significantly lose its suitable habitats in the upper Sefidroud sub-basin over the next few decades which may lead to the range reduction or range shift in P. iranicus sensu stricto, while suitable habitats will increase in the eastern regions which may lead to the range expansion of the P. patimari Hg1 ESU (Yousefi et al., 2020).