Specimens were collected at ten sampling locations by ultraviolet (UV) light detection at night and rock-rolling during the day. Newly collected material, transferred to 75–96% ethanol, is deposited in the American Museum of Natural History (AMNH), New York, U.S.A., and the Zoological Museum of Ferdowsi University of Mashhad (ZMFUM), Iran. Tissue samples are deposited in the Ambrose Monell Cryocollection (AMCC) at the AMNH. Type material of H. saulcyi (not examined) is deposited in the Museum National d’Histoire Naturelle (MNHN), Paris, France, and the Zoologisches Museum, Universität Hamburg (ZMUH), Germany.

Morphological characters of adult specimens, adopted from Stahnke (1970) and Lamoral (1979), were recorded with a > 0.02 mm Olympus micrometer eyepiece (OSM-4) attached to an optical Olympus SZ40 stereomicroscope. Morphological terminology follows Sissom (1990). Abbreviations for morphometric and meristic characters and ratios used in the statistical analyses and taxonomic descriptions are as follows:

Measurements.PeL: pedipalp length; FL: pedipalp femur length; FW: pedipalp femur width; PaL: pedipalp patella length; PaW: pedipalp patella width; ChL: pedipalp chela length; ML: pedipalp chela manus length; MW: pedipalp chela manus width; MD: pedipalp chela manus depth; MFL: pedipalp chela movable finger length; CAW: carapace anterior width; CPW: carapace posterior width; CL: carapace length; CX: distance between anterior margin of carapace and anterior edge of median ocelli; CY: distance between anterior edge of median ocelli and posterior margin of carapace; MsL: mesosoma length; TIIIL: mesosomal tergite III length; TIIIW: mesosomal tergite III width; MtIL, MtIIL, MtIIIL, MtIVL, MtVL: length of metasomal segments; MtIIW, MtIIIW, MtIVW, MtVW: width of metasomal segments I–V, respectively; MtID, MtIID, MtIIID, MtIVD, MtVD: depth of metasomal segments I–V, respectively; TL: telson length; TW: telson width; TD: telson depth; MtTL: metasoma and telson length; BL: total length; PL: pecten length; DBP: distance between pectines; PDO: distance between prolateral margins of ocelli; RDO: distance between retrolateral margins of ocelli.

Meristic characters: MFDR: count of median denticle subrows on pedipalp chela movable finger; PDR: sinistral and dextral pectinal tooth counts (male, female).

Morphometric ratios.FL/W: pedipalp femur length : width; PaL/W: pedipalp patella length : width; ChL/ML: pedipalp chela length : chela manus length; CAW/PW: carapace anterior width : posterior width; CL/AW: carapace length : anterior width; CL/PW: carapace length : posterior width; CX/CY: distance between anterior margin of carapace and median ocelli : distance between median ocelli and posterior margin of carapace; TIIIL/W: mesosomal tergite III length : width; MtIL/W, MtIIL/W, MtIIIL/W, MtIVL/W, MtVL/W: metasomal segments I–V length : width, respectively; MtIL/D, MtIIL/D, MtIIIL/D, MtIVL/D, MtVL/D: metasomal segments I–V length : depth, respectively; TL/W: telson length : width; TL/D: telson length : depth: TW/D: telson vesicle width : depth.

Distribution maps were created using DIVA-GIS 7.5 by overlaying point locality records of sampling locations on spatial layers depicting political boundaries and topography (elevation) at 2.5 arc-minutes altitude (Hijmans et al. 2004; Fig. 1).

The ingroup comprised 12 samples from nine locations in Iran and one in Iraq, rooted on the outgroup, Hottentotta schach (Birula, 1905), from Iran.

DNA was extracted from muscle tissue using the Favorgene (Taipei, Taiwan) DNA Extraction kit or the Qiagen (Hilden, Germany) DNeasy Blood and Tissue Kit. The nuclear and mitochondrial markers were amplified using standard primers (Prendini et al. 2005; Loria & Prendini 2020; Prendini & Loria 2020; Prendini et al. 2021).

A variety of optimization protocols were employed in the polymerase chain reaction (PCR) (Prendini et al. 2005, 2021; Azghadi et al. 2014; Barahoei et al. 2022; Loria & Prendini 2020; Prendini & Loria 2020). Each PCR reaction contained 12.5 μL Ampliqon (Odense, Denmark) ready master mix, 7.5 μL deionized H₂O, 1 μL (= 1 pmol) of each primer and 3 μL of DNA template. Sequencing was conducted using ABI Big Dye terminator chemistry on an ABI Prism 3700 (Applied Biosystems, Foster City, CA, USA). Sequences were edited and assembled using SEQUENCHER v. 4.5.6 (GeneCodes Corporation, Ann Arbor, MI, U.S.A.).

Fifty-two new sequences were generated for the study. Sequence data are deposited in GenBank (http://www.ncbi.nlm.nih.gov) with the following accession numbers (Table 1): 28S: PP133558–PP133570; 12S: PP133532–PP133544; 16S: PP133545–PP133557; COI: PP133849–PP133859 and PP135557–PP135558.

Sequences were aligned using MUSCLE (Edgar 2004). The aligned, concatenated dataset was 2009 nucleotide base-pairs (bp) in length comprising 513 bp (28S), 343 bp (12S), 496 bp (16S), and 657 bp (COI). In the 28S alignment, 510 (99.4%) positions were conserved, 3 (1.36%) variable, and 2 (0.38%) parsimony informative. In the 12S alignment, 257 (74.9%) positions were conserved, 85 (24.78%) variable, and 55 (14.5%) parsimony informative. In the 16S alignment, 386 (77.82%) positions were conserved, 104 (20.9%) variable, and 63 (12.7%) parsimony informative. In the COI alignment, 535 (81.4%) positions were conserved, 122 (18.56%) variable, and 96 (12.78%) parsimony informative. Kimura 2-parameter pairwise genetic distances were calculated for each marker using MEGA7 (Kumar et al. 2016), and ExcaliBAR (Aliabadian et al. 2014).

Bayesian Inference (BI) and Maximum Likelihood (ML) methods were applied to the concatenated data. BI was performed using the Markov Chain Monte Carlo method in MrBayes v. 3.2.2 (Huelsenbeck and Ronquist 2001). The dataset was partitioned into mitochondrial and nuclear subsets. The number of generations was set to 50 × 10⁶, and a tree sampled every 1000^th generation. The average standard deviation of split frequencies of the four simultaneous and independent runs performed was applied to resolve the stationary point of likelihoods (Ronquist & Huelsenbeck 2003). Bayesian tree and posterior probabilities were calculated by majority rule consensus after burning off all pre-asymptotic topologies.

ML analysis was performed in RAxML v. 1.3 (Silvestro & Michalak 2012) with 1000 rapid bootstraps. The best-fitting models of nucleotide substitution were estimated for each marker using the Akaike Information Criterion (Akaike 1973) in jModeltest v.2.1.10 (Posada 2008), and a tree constructed under GTR+I+G.

Nodes with posterior probabilities greater than 95% (Huelsenbeck and Ronquist 2001) and bootstrap support values greater than 70% (Hillis and Bull 1993) were considered strongly supported.

Following Barahoei et al. (2022), 41 linear distances (measurements, in mm) and three meristic characters were recorded on 66 adult specimens representing 10 populations deposited in ZMFUM (Table 2; Fig. 1A). Twenty-one ratios were calculated from the measurement data. The normality of variables was checked using the Shapiro-Wilk test and graphing a Q-Q plot, confirming that all met the normality assumption. A Welch’s two sample t-test was performed to test for significant differences between the two divergent clades of H. saulcyi sensu lato.

The morphological variation of the populations studied was assessed by calculating the two main aspects of body form, i.e., size and shape. The size of each specimen (hereafter called the overall size) was computed as the square root of the sum of all squared variables (Sundberg 1996; Navarro et al. 2004). Shape variables were computed as the log shape ratio of the morphometric variables (Eldredge 1972; Navarro et al. 2004). These shape data were size corrected and, because size may be influenced by environment and developmental conditions (Lira et al. 2020), are important for comparing the body shape of populations in multivariate space.

Pairwise correlation between the shape variables was estimated to ensure character independence. MIIID was removed from the multivariate statistical analyses due to its high correlation (r² > 0.95) with MIID and MtIVD.

In order to evaluate whether the populations are morphologically different from one another and whether there is significant sexual dimorphism in either size or shape, type II two-way ANOVA and MANOVA were performed on overall size and shape data, respectively, using the factors of location and sex. The same analyses were carried out on the factors of species and sex, to examine the size difference between species. When ANOVA revealed a significant difference in size, the posthoc Tukey’s Honestly Significant Difference (HSD) test was performed for pairwise comparison among populations. A box plot was also performed on the overall size to visualize the size differences in males and females of the populations studied.

As sexual dimorphism in shape was significant and the sex ratio of the specimens studied for each population was unequal, the effect of sex was corrected by extending the Burnaby (1966) procedure, i.e., sending the shape variation in an orthogonal space against the sexually dimorphic shape. The phenotypic differences between populations were analyzed via a linear discriminant analysis (LDA) on the sex-corrected shape variables using population factor and via a minimum spanning tree on the Mahalanobis distances between populations. In order to mitigate against the effects of small sample size, populations represented by fewer than five specimens were omitted from the LDA and projected onto the discriminant axes a posteriori (sensu Claude 2008).

All statistical analyses of the morphometric data were conducted using R (R Development Core Team 2023). The MASS package (Venables & Ripley 2002) was used for LDA; the ape package (Paradis and Schliep 2019) for the minimum spanning tree, the car package (Fox and Weisberg 2019) for ANOVA II; the Agricolae package (De Mendiburu & Yaseen 2020) for the Tukey HSD test; the ggplot2 package (Wickham 2016) for box plots; and the Picante package (Kembel et al. 2010) for pairwise correlation among characters.

The tree topologies obtained by analysis of the concatenated aligned data with BI and ML were entirely congruent (Fig. 2). Both recovered two well supported clades (Fig. 2). Clade A, comprising samples from Aleshtar, Borojerd and Iraq, represents the typical form of H. saulcyi and Clade B, comprising samples from Ilam, Izeh, Masjedsoleiman and Andimeshk represents a distinct species, described below as H. hatamtiorumsp. nov. Separate analyses of the aligned 28S, 12S, 16S and COI data with BI recovered similar tree topologies to those obtained by analysis of the concatenated dataset.

The average K2P genetic distance between the two clades, corresponding to H. saulcyi (Clade A) and H. hatamtiorumsp. nov. (Clade B), ranged from 10.0–17.8% for 12S, 11.2–14.6% for 16S, and 9.8–11.7% for COI (Table 3). As expected, K2P distances were much lower (0.2–0.6%) for the nuclear 28S marker than for the mitochondrial markers. Values for intrapopulational genetic distances ranged from 1.1–7.8% for 12S, 0.6–5%, for 16S and 1.6–6.8% for COI. The maximum intrapopulation genetic distances for all mitochondrial markers were observed in H. hatamtiorumsp. nov.

The type II two-way ANOVA on overall size for species × sex and population × sex, followed by the posthoc Tukey’s HSD test suggested a significant size difference among populations (F = 6.17; d.f. (effect) = 9; d.f. (residual) = 48; p < 10^–5), but was non-significant (F = 0.002; d.f. (effect) = 1; d.f. (residual) = 62; p > 0.96) between the two species. This suggests that size may be under strong selection within species which could be due to habitat loss and transformation (Miyashita et al. 1998; Hillaert et al. 2018, Penell et al. 2018; Salomão et al. 2018; Lira et al. 2020). Sexual dimorphism in size was significant either considering species (F = 4.02; d.f. (effect) = 1; d.f. (residual) = 62; p = 0.049) or population (F = 5.79; d.f. (effect) = 1; d.f. (residual) = 48; p = 0.02) as the main factor, suggesting that size should be under slight sexual selection in H. saulcyi sensu lato. The non-significant interaction between species and sex (F = 2.13; p > 0.14) and between population and sex (F = 1.73; p > 0.12) implies that the evolution of sexual size dimorphism is similar in H. saulcyi and H. hatamtiorumsp. nov. and in their respective populations.

The results of Tukey’s HSD test on the overall size of the different populations confirmed that H. saulcyi (Clade A) is generally larger at Borojerd and smaller at Sahneh than at other localities whereas H. hatamtiorumsp. nov. is generally larger at Masjedsoleiman, and Poldokhtar than at Ilam and, especially, Alborz. The size variation among populations is shown in Fig. 3.

The type II MANOVA on shape data revealed significant differences between populations (F = 2.55; p < 10^–10) and between sexes (F = 14.95; p < 10^–9) but no significant difference for the interaction thereof (F = 0.92; p > 0.66). This suggests that shape may be under strong selection both across and within species, but that sexual dimorphism in shape is similar for both species and their populations. The t-test indicated significant differences among H. saulcyi and H. hatamtiorumsp. nov. in 17 characters (Table 4).

In the linear discriminant analysis (LDA) on population factor using sex-corrected shape variables (Fig. 4), the first two components (LD1 and LD2) explained 36.5% and 22% of the variance, respectively. The populations comprising Clades A and B formed distinct phenetic clusters, and LD1 supported the divergence between H. saulcyi (Clade A) and H. hatamtiorumsp. nov. (Clade B). The population of Poldokhtar, geographically close to the populations of H. hatamtiorumsp. nov., and for which no molecular data were available, clustered more closely to the populations of H. saulcyi. A linear discriminant analysis on species factor was therefore performed, taking only the populations with molecular data support into account and projecting all populations on the discriminant axis, a posteriori. The results (not shown) suggest that the population from Poldokhtar is conspecific with the populations of H. hatamtiorumsp. nov. from Alborz, Ilam and Masjedsoleiman, and not with H. saulcyi. This finding was also confirmed by the minimum spanning tree on the Mahalanobis distance among populations (Fig. 4).