Glyptothorax (Teleostei: Sisoridae) from the Middle East: An Integrated Molecular and Morphological Insight into Its Taxonomic Diversity

: The Glyptothorax species from the Middle East are taxonomically revised based on extensive geographic range and taxon sampling, tree topologies from mitochondrial COI and Cyt b and nuclear RAG2 markers (2532 bps), molecular species delimitation and genetic distance analyses of DNA sequences against morphometric and morphological characters. A majority-rule consensus based on conceptually different molecular species delimitation analyses combined with the Bayesian and maximum likelihood tree topologies considered all the name-bearing Iranian endemic clades of Glyptothorax , except for G. pallens (i.e., G. alidaeii , G. galaxias , G. hosseinpanahii , G. shapuri and G. silviae ) as a single molecular entity. We also lent our years of experience to the morphology of Iranian Glyptothorax populations and tried to perceive consistent morphological differences, but without success. Therefore, based on this integrated molecular and morphological study, we treat G. alidaeii , G. galaxias , G. hosseinpanahii and G. shapuri as conspeciﬁc with G. silviae . Furthermore, our molecular and morphological results conﬁrmed the ﬁrst record of G. cous in Iranian waters. The species G. armeniacus , G. cous , G. daemon , G. kurdistanicus , G. pallens , G. silviae and G. steindachneri are considered as valid species.


Introduction
Siluriformes is one of the most species-rich actinopterygian orders, with about 4131 valid species in 39 families, including the family Sisoridae [1]. Members of the sisorid catfish genus Glyptothorax Blyth, 1860 (Siluriformes: Sisoridae) are typically inhabited in fastflowing hill streams or faster-flowing reaches of larger rivers from Asia (in the Tigris and Euphrates River drainages of Türkiye, Syria, Iraq and Iran), eastward to the Yangtze River drainage and southward to South-East Asia [2,3]. It is the most species-rich and widely distributed genus in the family, with new species being discovered regularly [1]. Until 2021, no data were published on the taxonomic revision of this genus in the Middle East. Recently, Freyhof et al. [4] and Mousavi-Sabet et al. [5] revised the Glyptothorax of the Euphrates, Tigris and the Persis River tributaries (i.e., the Persian Gulf basin) based on the morphological and mtDNA COI barcode region and described six new species.
Here, based on our available data since 2009, we present morphological and molecular (mitochondrial COI and Cyt b and nuclear RAG2) data to (i) provide the first multilocus phylogeny of the genus Glyptothorax in the Middle East, (ii) investigate the validity of newly described species from Iran and also species-level diversity within the genus in the area and (iii) update the distributional ranges for the species and lineages based on new samplings and DNA sequence data. To recognize the number of valid species within the genus Glyptothorax in the region, we identify distinct lineages based on tree topologies genus in the area and (iii) update the distributional ranges for the species and lineages based on new samplings and DNA sequence data. To recognize the number of valid species within the genus Glyptothorax in the region, we identify distinct lineages based on tree topologies (mitochondrial COI sequence data) as 'molecular species'. We then seek morphological differences to differentiate the 'molecular species' hypothesized by the tree topology and molecular species delimitations. This publication was postponed by more than a decade as we tried multiple times to find morphological differences between the Iranian molecular clades recognized in this genus. Results should provide a primary understanding of the genus taxonomic diversity and distribution in the region, which are relevant for future taxonomic efforts and formulating conservation strategies and management proposals.

Sampling and Morphological Measurements
New samples of Glyptothorax were collected from 21 Iranian localities in the tributaries of the Persian Gulf basin (Figure 1; Table 1). Specimens were collected using an electro-fishing device. Clove oil solution (1%) was used as an anesthetic.  Table 1. Un-numbered localities refer to those from Freyhof et al. [4] and Mousavi-Sabet et al. [5]. The map was originally designed in DIVA-GIS 7.5 (https://www.diva-gis.org/). Morphometric measurements were made with a dial calliper and recorded to 0.1 mm. All measurements were made point to point, never by projections. Subunits of the head are presented as proportions of head length (HL). Head length and measurements of body parts are given as proportions of standard length (SL). Methods for counts and measurements follow Ng and Dodson [6], except for caudal fin length which was measured as the length of the largest ray of the upper lobe.
For molecular analysis, the right pectoral fin of some specimens was clipped and fixed in 96% ethanol, while the voucher specimens were kept in 10% formaldehyde and subsequently stored in 70% ethanol. These specimens are deposited in the Zoological Museum of Shiraz University, Collection of Biology Department (ZM-CBSU). Additional  Table 1. Un-numbered localities refer to those from Freyhof et al. [4] and Mousavi-Sabet et al. [5]. The map was originally designed in DIVA-GIS 7.5 (https://www.diva-gis.org/) (accessed on 28 February 2022).
Morphometric measurements were made with a dial calliper and recorded to 0.1 mm. All measurements were made point to point, never by projections. Subunits of the head are presented as proportions of head length (HL). Head length and measurements of body parts are given as proportions of standard length (SL). Methods for counts and measurements follow Ng and Dodson [6], except for caudal fin length which was measured as the length of the largest ray of the upper lobe.
For molecular analysis, the right pectoral fin of some specimens was clipped and fixed in 96% ethanol, while the voucher specimens were kept in 10% formaldehyde and subsequently stored in 70% ethanol. These specimens are deposited in the Zoological Museum of Shiraz University, Collection of Biology Department (ZM-CBSU). Additional Glyptothorax specimens (G. daemon and G. armeniacus) were also included from Türkiye (Table 1).
Abbreviations used. SL, standard length; HL, lateral head length. Collection codes. BMNH, Natural History Museum, London; CMNFI, Canadian Museum of Nature; ZFMK, Zoologisches Forschungsmuseum Alexander Koenig; ZM-CBSU, Zoological Museum of Shiraz University, Collection of Biology Department, Shiraz. Table 1. List of novel and previously studied Glyptotorax spp. specimens from the Middle East used in this study. Sequences used for the combined analysis (COI + Cyt b + RAG2) are marked with Greek numbers (I-XXIII). Source: A, Mousavi-Sabet et al. [5]; B, Freyhof et al. [4]; C, this study (GenBank numbers are given in Table S3); † a possible hybrid population (see Freyhof et al. [4]; ‡ new locality numbers refer to the new sampling sites from Iran, as shown in Figure 1.  Iran: Kermanshah prov., stream Zemkan 3 km north of Zamkan-e Olya mtDNA genes, cytochrome c oxidase subunit 1 (COI) barcode region and Cytochrome b (Cyt b) and one nuclear locus, Recombination Activating 2 (RAG2). Primers, reaction profiles and cycling conditions were modified from Jiang et al. [8] and Zarei et al. [9] (Table S1). PCR amplifications were conducted in a total volume of 25 µL (including 12.5 µL of a ready 2X Taq PCR Master Mix (ParstousTM), 0.5 µL of each primer (10 pmol/µL), 6 µL of the target DNA and 5.5 µL DNase-free distilled water). Automated Sanger sequencing of purified PCR products was performed using ABI Prism R BigDyeTM Terminator v. 3.1 chemistry (Applied Biosystems, Foster City, CA, USA) on an Applied BiosystemsTM ABI PRISM 3730xl by Macrogen, South Korea.

Additional Sequence Data
In addition to the newly determined DNA sequences from Iran (25 COI, 24 Cyt b and 17 RAG2 sequences) and Turkey (4 COI sequences), we used 64 previously published and georeferenced COI barcode sequences from 10 Glyptothorax species from Iran and Turkey (i.e., G. pallens, G. kurdistanicus, G. daemon, G. daemon, G. cous, G. armeniacus, G. galaxias, G. silviae, G. shapuri, G. alidaeii and G. hosseinpanahii) obtained from GenBank and ZFMK and derived from the works of Freyhof et al. [4] and Mousavi-Sabet et al. [5] (Table 1). For the multilocus phylogenetic analysis, also additional COI, Cyt b and RAG2 sequences of 21 Glyptothorax species from different Asian localities and derived from the works of Jiang et al. [8] and Sullivan et al. [10] were obtained from GenBank (Table S2). Bagarius yarrelli (Sisoridae) from Chen et al. [11] and Peng et al. [12] was used as outgroup.

Phylogenetic Analyses
BioEdit 7.0.4 [13] was used to read and edit the DNA chromatograms. All novel DNA sequences including 29 COI, 24 Cyt b and 17 RAG2 sequences were deposited in GenBank (accession numbers and details of the specimens are given in Table S3: COI: OP585111-OP585139; Cyt b: OP589359-OP589381; RAG2: OP589382-OP589398) and the collection sites of the Middle Eastern Glyptothorax samples used in this study are illustrated in Figure 1. The DNA sequences were aligned using the ClustalW algorithm implemented in MEGA 7 [14]. After having trimmed out the sequences tails, the alignment of the novel fragments and those downloaded from GenBank led to trimmed aligned fragments of 651 bp (COI), 989 bp (Cyt b) and 892 bp (RAG2). The nucleotide substitution models best fitting the COI, Cyt b and RAG2 sequences were selected using the Akaike information criterion (AIC) [15] in JModelTest 2.1.3 [16]. The saturation test of Xia et al. [17] in DAMBE 7 [18] was used to test the nucleotide substitution saturation in the sequences. The fragments of both the mitochondrial and nuclear-selected markers were concatenated in a single partitioned dataset (2532 bp). Phylogenetic analyses of the partitioned concatenated dataset including the fragments of the amplified DNA markers were conducted using the software packages MrBayes v. 3.2.6 [19] and RaxML 7.2.5 (10,000 bootstrap replicates) [20] to infer the phylogenetic relationships through Bayesian Inference (BI) and Maximum Likelihood (ML) analyses. As support measures for the nodes, bootstrap values (BPs) [21] were calculated with 1000 replicates in the ML trees, whereas the node posterior probability (PP) values were reported in the BI trees. In the BI analyses, two independent Markov Chain Monte Carlo analyses were performed with 30 million generations (25% as burn-in). Convergence in the analysis was reached (effective sample size (ESS) >200 in all the analyses performed).

Molecular Species Delimitations
To delineate "molecular species", we used four molecular species delimitation methods on the COI dataset: (i) two distance-based methods, the Automatic Barcode Gap Discovery (ABGD) [22] and Assemble Species by Automatic Partitioning method (ASAP) [23]; (ii) a network-based method, the reversed Statistical Parsimony (SP) [24]; and (iii) a topologybased method, the multiple-rate Bayesian Poisson Tree Process (mPTP) [25]. We tested the COI dataset on the ABGD webserver (https://bioinfo.mnhn.fr/abi/public/abgd/ (accessed on 28 February 2022)) with a combination of ABGD settings within a parameter range of P min = 0.001, P max = 0.1 and a gap width of 1.1, all for a total of 10 steps. ASAP was run through its web interface (https://bioinfo.mnhn.fr/abi/public/asap/asapweb.html (accessed on 28 February 2022)) using the Kimura (K80) ts/tv (=0.2). TCS 1.21 [26] was used to calculate a Statistical Parsimony (SP) network, using a 95% connection probability threshold to delineate molecular species. The mPTP server (https://mptp.h-its.org (accessed on 28 February 2022)) was used with a Bayesian topology produced in MrBayes as an input tree and the analysis was run under default settings. Estimates of genetic divergence between the Middle Eastern species/lineages of Glyptothorax were calculated using the K2P model implemented in MEGA 7. PopART 1.7 [27] was used to investigate the phylogenetic depth and evolutionary relationships among haplotypes based on the median-joining (MJ) method.

Phyleogenetic Placements and Molecular Species Delimitations
The final COI alignment of Glyptothorax from the Middle East included 29 novel (25 sequences from 21 Iranian localities and 4 sequences from 3 localities in Turkey) and 64 previously published sequences derived from Freyhof et al. [4] and Mousavi-Sabet et al. [5], trimmed to 651 bps. The nucleotide substitution pattern showed that the COI sequences have not reached substitution saturation (Iss < Iss.cS and Iss < Iss.cA) and are, therefore, well applicable for phylogenetic analysis (Table S4). Based on a General Time-Reversible model of sequence evolution with a gamma-distributed rate variation among sites (GTR+G), this dataset was analyzed using the BI and ML methods (Figure 2), resulting in a phylogeny with 11 major clades corresponding to 10 nominal species: G. pallens, G. kurdistanicus, G. armeniacus, G. galaxias, G. silviae, G. shapuri, G. alidaeii, G. hosseinpanahii, G. cous and G. daemon. Glyptothorax daemon is separated into two mitochondrial lineages (hereafter, G. daemon 1 and G. daemon 2). According to Freyhof et al. [4], G. daemon 2 may represent a hybrid population between G. cous and G. daemon 1. Most of the nodes within the genealogy are supported with high PB and BP values. The deepest split in this genealogy is between G. pallens, followed by G. kurdistanicus and the remainder of the species. The next diverging clade contains G. cous and the two lineages of G. daemon. All the Iranian endemic members of the genus Glyptothorax, except for G. pallens (i.e., G. galaxias, G. silviae, G. shapuri, G. alidaeii and G. hosseinpanahii) form a monophyletic clade and are more derived than other species. Glyptothorax armeniacus, an endemic catfish to the headwater streams in the Euphrates drainage (Türkiye), is the sister group of this Iranian endemic clade.
The nucleotide substitution pattern also confirmed that the Cyt b and RAG2 sequences (24 and 17 novel sequences from Iran, respectively) have not reached substitution saturation (Iss < Iss.cS and Iss < Iss.cA) and are, therefore, well applicable for phylogenetic analyses (Table S4). The best-fit nucleotide substitution models for these loci were GTR+I+G and SYM+G, respectively. The combined COI+Cyt b (Figure 3a) and COI+Cyt b+RAG2 (Figure 3b) trees for the Iranian species with high support values (PP and BP) retrieved similar results with the topology from the COI tree. The deepest split is between G. pallens and the remainder of the species; the next diverging clade contains G. kurdistanicus and G. cous, followed by G. galaxias as the sister group of the remainder of the species. Glyptothorax hosseinpanahii and G. alidaeii formed a clade sister to another clade containing G. silviae and G. shapuri.
The molecular trees (Figures 2 and 3) show that the new specimens collected from Iran at the Kohmareh Sorkhi stream and the Dalaki River belong to G. shapuri, at the Ab-e Aala and Maroun rivers belonging to G. silviae, at the Fahlian River belonging to G. hosseinpanahii, at the Khersaan, Beshar, Armand and Kashkan rivers and Bibi Seyedan (Isfahan) belonging to G. alidaeii, at the Gamasiab and Chardaval rivers, Khak Daneh (Isfahan) and Kata (Kohgiluyeh and Boyer-Ahmad) belonging to G. galaxias, at the Sirvan River belonging to G. kurdistanicus and G. pallens and at the Seimare River (Ilam) belonging to G. cous, which represents the first record of G. cous in Iranian waters. The updated distributional ranges of these species based on molecular data are depicted in Figure 1 The nucleotide substitution pattern also confirmed that the Cyt b and RAG2 sequences (24 and 17 novel sequences from Iran, respectively) have not reached substitution saturation (Iss < Iss.cS and Iss < Iss.cA) and are, therefore, well applicable for phylogenetic analyses (Table S4). The best-fit nucleotide substitution models for these loci were GTR+I+G and SYM+G, respectively. The combined COI+Cyt b (Figure 3a) and COI+Cyt   The four conceptually different molecular species delimitation methods revealed 5-12 molecular species in a COI dataset comprising 93 sequences, belonging to 10 nominal species (Figure 2). SP was the most conservative method, as it delineated five molecular species. Out of 10 nominal species, 3 species (i.e., G. pallens, G. kurdistanicus and G. armeniacus) were delimited through the SP tool. It grouped G. daemon (including G. daemon 1 and G. daemon 2) and G. cous and all members of the Iranian endemic clade (i.e., G. galaxias, G. silviae, G. shapuri, G. alidaeii and G. hosseinpanahii) into single molecular entities. The ABGD and ASAP methods similarly delineated nine molecular species. Out of 10 nominal species, 4 species (i.e., G. kurdistanicus, G. cous, G. armeniacus and G. hosseinpanahii) were delimited through these tools. Both analyses presented potential species-level diversity in G. pallens (one in the Alvand stream and the other in the Siravan River, the Golain and Zemkan streams) and G. daemon (one as G. daemon 1 and the other as G. daemon 2). Both methods pooled all species of the Iranian endemic clade, except for G. hosseinpanahii (i.e., G. galaxias, G. silviae, G. shapuri, G. alidaeii), into a single molecular species. The mPTP method delineated 12 molecular species. Out of 10 nominal species, 8 species (i.e., G. kurdistanicus, G. cous, G. armeniacus, G. galaxias, G. silviae, G. shapuri, G. alidaeii and G. hosseinpanahii) were delimited through this tool. Similar to ABGD and ASAP, mPTP presented potential species-level diversity in G. pallens and G. daemon. In addition, species delimitation based on the K2P genetic distance new species threshold value (i.e., 2%: [28]) only provided support for six molecular species (Table 2, Figure 2). Out of 10 nominal species, 3 species (i.e., G. pallens, G. kurdistanicus and G. armeniacus) were delimited through the K2P method. It lumped G. daemon 2 and G. cous and all members of the Iranian endemic clade into single molecular entities. The median-joining (MJ) haplotype networks for COI and Cyt b in Figure 4 show the number of mutational steps and phylogeographic depth between the observed haplotypes of Glyptothorax from the Middle East. They further show marked divergences within two species: (i) G. pallens: divergence between the Alvand stream and the Siravan River, the Golain and Zemkan streams haplotypes (deviated by nine mutational steps = 1.5% K2P genetic divergence in COI), (ii) G. galaxias: divergence between the upper Karun drainage (G. galaxias) and the upper Karkheh drainage (G. cf. galaxias) haplotypes (deviated by six and thirteen mutational steps in their COI and Cyt b sequences, respectively).

The Name-Bearing Iranian Endemic Clades
See Figures 5 and 6 for general appearance and Table 3 for morphometric data of studied populations. The holotype general morphology of G. silviae is illustrated in Figure  7. As the morphological description of new Iranian species was mainly based on the thoracic adhesive apparatus (e.g., anteromedial striae), length of maxillary, inner mandibular and outer mandibular barbels and the caudal fin shape, the mentioned characters and several others were studied. Based on the data given in Table 3 and Figures 5 and 6, these characters show inter-and intra-population variation in closely related species (i.e., G. alidaeii, G. galaxias, G. hosseinpanahii, G. shapuri and G. silviae).

The Name-Bearing Iranian Endemic Clades
See Figures 5 and 6 for general appearance and Table 3 for morphometric data of studied populations. The holotype general morphology of G. silviae is illustrated in Figure 7. As the morphological description of new Iranian species was mainly based on the thoracic adhesive apparatus (e.g., anteromedial striae), length of maxillary, inner mandibular and outer mandibular barbels and the caudal fin shape, the mentioned characters and several others were studied. Based on the data given in Table 3 and Figures 5 and 6, these characters show inter-and intra-population variation in closely related species (i.e., G. alidaeii, G. galaxias, G. hosseinpanahii, G. shapuri and G. silviae).        Diagnosis: Glyptothorax silviae is a species with a large variation in morphological characteristics, such as length of barbels, shape of thoracic adhesive apparatus, caudal peduncle index, adipose-fin length and distance between the base of the last dorsal-fin ray and adipose-fin origin, which are used in the original descriptions of the four newly described species, G. alidaeii, G. galaxias, G. hosseinpanahii and G. shapuri. Based on the original description, G. silviae is distinguished from G. alidaeii and G. shapuri by having maxillary barbel as long as the head (90-110% HL vs. 76-89% HL); inner mandibular barbel 45-48% HL (vs. 19-33%); outer mandibular barbel 67-74% HL (vs. 42-49%). It is also distinguished from G. galaxias and G. hosseinpanahii by having long and numerous (vs. short or absent) anteromedial striae in thoracic adhesive apparatus. Our data do not confirm these Diagnosis: Glyptothorax silviae is a species with a large variation in morphological characteristics, such as length of barbels, shape of thoracic adhesive apparatus, caudal peduncle index, adipose-fin length and distance between the base of the last dorsal-fin ray and adipose-fin origin, which are used in the original descriptions of the four newly described species, G. alidaeii, G. galaxias, G. hosseinpanahii and G. shapuri. Based on the original description, G. silviae is distinguished from G. alidaeii and G. shapuri by having maxillary barbel as long as the head (90-110% HL vs. 76-89% HL); inner mandibular barbel 45-48% HL (vs. 19-33%); outer mandibular barbel 67-74% HL (vs. 42-49%). It is also distinguished from G. galaxias and G. hosseinpanahii by having long and numerous (vs. short or absent) anteromedial striae in thoracic adhesive apparatus. Our data do not confirm these diagnostic features, G. silviae vs. G. alidaeii and G. shapuri: (maxillary barbel 103-117% HL vs. 78-121% HL); (inner mandibular barbell 27-37% HL vs. 23-39% HL) and (outer mandibular barbell 54-65% HL vs. 42-67% HL). There is also a variety in the size and number of anteromedial striae in thoracic adhesive apparatus of G. silviae, G. galaxias and G. hosseinpanahii (see Figure 5a,d,e). Therefore, our data confirm the synonymy of G. alidaeii, G. galaxias, G. hosseinpanahii and G. shapuri with G. silviae.
Glyptothorax silviae is generally distinct in having long, or approximately as wide as long, thoracic adhesive apparatus (vs. short in G. kurdistanicus); head and flank without tubercles (vs. present in G. armeniacus and G. cous) and warts (vs. present in G. daemon); upper head, back and flank with few or many dark-brown spots and/or blotches (vs. in G. pallens) and medial pit of thoracic adhesive apparatus without striae (vs. present in G. steindachneri). 3a-Thoracic adhesive apparatus wider than long, as wide as long in juveniles (0.7-0.9-times longer than wide); adipose-fin short, its length 0.6-1.0-times larger than the distance between the base of last dorsal-fin ray and adipose-fin origin . . . . . . . . . . . . G. kurdistanicus. 3b-Thoracic adhesive apparatus longer than wide or as wide as long (1.0-1.6-times longer than wide); adipose-fin usually long, its length 1.1-3.6-times larger than the distance between the base of last dorsal-fin ray and adipose-fin origin . . . . . . . . . . . . . . . . G. silviae.

Glyptothorax cous (Linnaeus, 1766)
Silurus cous Linnaeus, 1766:504 ( Figure 8). Diagnosis. Glyptothorax cous is distinct in having large mouth width (37.4-40.1% HL), large dark-brown blotches and numerous cloudy spots on the body. Thoracic adhesive apparatus is almost as wide as long (1.0-1.1-times longer than wide), poorly delineated at its posterior margin, only partly situated on a shallow, horseshoe-shaped swelling, swelling absent in many adult individuals; head, back and flank usually have large, bony, striated and elongated tubercles (absent in some individuals).
General morphology. See Figure 8 for the general appearance of two collected specimens from the Seimareh River, Iran, and Figure 9 for the general morphology of the type specimen; morphometric data are given in Table 3. Body spindle-shaped, pre-ventral profile flat and pre-dorsal profile slightly arched; head depressed, with bluntly pointed snout tip and sub-terminal mouth and fleshy lips; upper jaw longer than lower jaw, eyes small; thorax and anterior portion of abdomen flattened ventrally with a conspicuous rhomboidal-shaped thoracic adhesive apparatus with a central pit, which is open caudally; barbels four pairs, the nasal reaches posterior margin of orbit and longer than inner mandibular, maxillary with broad bases, when straightened, extending beyond posterior end of pectoral fin base, outer mandibular reaches middle of pectoral fin base; dorsal fin is small and situated closer to the snout tip than to caudal fin, its posterior margin slightly curved, with two spines and 6 1 2 branched rays, spines smooth with no serrations; pectoral fin extending to vertical of origin of dorsal fin-base, with one spine and 8 1 2 branched rays, spine strong and broad with posterior serrations; pelvic fin-origin at level of behind posterior end of dorsal base and extending to the anus, with one spine and 5 1 2 branched rays; anterior margin of adipose dorsal gently sloping and its origin is over or slightly anterior to the origin of anal fin; anus with two spines and 7 1 2 -8 1 2 branched rays; caudal is quite deeply forked, lower lobe slightly longer; thoracic adhesive apparatus index (L/W) 1.0-1.1; caudal peduncle index (L/D) 1.9-2.0; adipose fin length/distance between dorsal and adipose fin 1.98-2.04; head blunt, spade-shaped, 31.0-31.4% SL. Distribution. Glyptothorax cous is known from the Tigris and Euphrates River drainages in Iraq and Turkey. It is also known from the Seimareh River (Karkheh tributary), Cham Namesht, Persian Gulf basin in Iran, as a new record (Figure 1).

Discussion
The Iranian Glyptothorax species are taxonomically revised based on extensive geographic range and taxon sampling, tree topologies from mitochondrial COI and Cyt b and nuclear RAG2 markers (2532 bps), molecular species delimitation and genetic distance Distribution. Glyptothorax cous is known from the Tigris and Euphrates River drainages in Iraq and Turkey. It is also known from the Seimareh River (Karkheh tributary), Cham Namesht, Persian Gulf basin in Iran, as a new record (Figure 1).

Discussion
The Iranian Glyptothorax species are taxonomically revised based on extensive geographic range and taxon sampling, tree topologies from mitochondrial COI and Cyt b and nuclear RAG2 markers (2532 bps), molecular species delimitation and genetic distance (K2P) analyses of DNA sequences against morphometric and morphological characters. All the Iranian endemic samples of Glyptothorax, following congruently well-supported BI and ML-based entities, formed six named molecular clades, i.e., G. pallens, G. galaxias, G. silviae, G. shapuri, G. alidaeii and G. hosseinpanahii. Only two groups were detected by the SP and K2P species delimitation approaches, i.e., G. pallens and the other Iranian endemic clades (i.e., G. galaxias, G. silviae, G. shapuri, G. alidaeii and G. hosseinpanahii) as the second group. However, a third group, G. hosseinpanahii, was detected by the ABGD and ASAP methods. Considering it as an independent species would render the rest of the second group a paraphyletic species. Therefore, the four molecular species delimitation methods combined with the BI and ML topologies (i.e., a majority rule consensus) consider all the Iranian endemic clades of Glyptothorax, except for G. pallens, as a single molecular entity. However, as there are not-and cannot be-any agreed universal thresholds for species delineation based on molecular data, we, as most other ichthyologists, usually adopt an iterative taxonomic approach [29]. This approach asks for the formulation of a testable species hypothesis and its testing by independent methods. If independent datasets (i.e., molecular data and morphology) agree in distinguishing two populations, these are recognized as valid species. We utilized our years of experience with the morphology of Iranian Glyptothorax populations and tried our best to perceive consistent morphological differences, but without success. We discuss this point in detail below.
Based on our available data to GF, Mousavi-Sabet et al. [5] described five new species of Glyptothorax from Iran. Glyptothorax alidaeii, G. galaxias, G. hosseinpanahii and G. shapuri are described based on low-genetic distances (<2%) and the overlapping morphological characters. However, due to the overlapping, we failed to distinguish these species based on the given morphometric characters. The morphological description of new species was mainly based on the thoracic adhesive apparatus (e.g., anteromedial striae), length of maxillary, inner mandibular and outer mandibular barbels and the caudal fin shape; however, these characters show inter-and intra-population variation in closely related species (i.e., G. alidaeii, G. galaxias, G. hosseinpanahii, G. shapuri and G. silviae). According to Mousavi-Sabet et al. [5], G. silviae is distinguished from G. hosseinpanahii by having the thoracic adhesive apparatus strongly elevated (vs. moderately elevated), without or with few short anteromedial striae (vs. many and long). However, as illustrated in Figure 5, the thoracic adhesive apparatus shows variation in the populations. According to them, G. alidaeii is distinguished from G. shapuri by having the thoracic adhesive apparatus moderately elevated (vs. strongly elevated), which is not a good taxonomic character state to distinguish these two species.
The mental disc was used as a significant character for taxonomic and systematic inferences within the genus Garra [34,35]. However, in the case of the Tashan cave barb Garra tashanensis, inhabiting a small cave in southwest Iran, there are two forms of fishes in terms of mental discs (sucking mouth disc), including disc-bearing and disc-less fishes. However, based on most species delimitation algorithms, both forms are members of a single taxonomic unit with genetic distances of 0.3-0.8% [36]. Hence, the observed mental disc variation was not inferred to be a taxonomic character or a consequence of feature displacement. Instead, it can be considered to be a case of character release to diversify among ecological niches in the limited subterranean habitat or a case of relaxed selection, which should be clarified in detail in follow-up ecological and population genetic studies [35].
The low genetic variations between G. alidaeii, G. galaxias, G. hosseinpanahii, G. shapuri and G. silviae might also be explained by ecological niche variation in the closely placed river systems, the geological events, such as the Late Pleistocene glacial and recent isolation of Iranian drainage basins [37].
Based on this integrated molecular and morphological study, we treat G. alidaeii, G. galaxias, G. hosseinpanahii and G. shapuri as the junior synonyms of G. silviae. The species, G. armeniacus, G. cous, G. daemon, G. kurdistanicus, G. pallens, G. silviae and G. steindachneri are considered as valid species.
Iranian inland water habitats have been and continue to be threatened by the main causes of biodiversity loss. The taxonomic revision of the Iranian endemic group of Glyptothorax here has several implications for conservation: (i) it confirms the presence of two valid species in this group as the main units of conservation (i.e., G. pallens and G. silviae), (ii) it documents a narrow distributional range and diversity for G. pallens, but a wider distributional range and mitochondrial variability for G. silviae and (iii) it allows for delimiting conservation units below species level, so that limited conservation resources can be utilized optimally [38]. The molecular and morphological data indicate that G. silviae comprises just one evolutionary significant unit; however, according to the definition of Moritz [39], this evolutionary unit is best viewed as comprising five management units corresponding to each of the identified subclades. Glyptothorax silviae has a high level of genetic diversity if treated as only one panmictic group; however, species structured into management units necessitate separate monitoring and management of each unit [40,41].
Supplementary Materials: The following supporting information can be downloaded at: https: //www.mdpi.com/article/10.3390/d14100884/s1, Tables S1: Primers used to amplify and their sequences; Table S2: List of different Glyptothorax species from outside the Middle East used in this study; Table S3: GenBank accession numbers (seq. ID/GenBank number) for the new Glyptothorax material examined from the Middle East (details of the specimens are given in Table 1); Table S4: Results of the Xia's nucleotide substitution saturation test for COI, Cyt b and RAG2 of the studied Glyptothorax specimens.