Toad-headed agamas of the genus Phrynocephalus Kaup, 1825, are distributed from south-eastern Europe and southwest Asia (including the Middle East and Arabian Peninsula) through Middle Asia to Central Asia (northern and central China and Mongolia). This taxonomically complicated genus currently contains up to 32 species (Uetz and Hošek 2016). The secret toad-headed agama, or Phrynocephalus mystaceus (Pallas, 1776), is one of the largest representatives of the genus, and is easily distinguished from all other congeners by a pair of large fringed cutaneous folds at the mouth angles. It is a specialized psammophilous species that inhabits sand dunes from Caspian region of the south-eastern part of European Russia in the west to the Ili valley in eastern Kazakhstan and western China in the east, and from Kazakhstan in the north through Middle Asia to northeastern Iran in the south (Bannikov et al. 1977; Zhao and Adler 1993; Anderson 1999; Ananjeva et al. 2004; Molavi et al. 2014; Fig. 1).

Phrynocephalus mystaceus was shown to have a high level of anatomical variability (Ananjeva “1986” 1987), which, together with its unique karyotype (Zeng et al. 1997), has led to its uncertain taxonomic classification at the generic level. Eichwald (1831) proposed the new generic name Megalochilus Eichwald, 1831 for Ph. mystaceus, which was synonymized with the genus Saccostoma by Fitzinger (1843). Ananjeva (“1986” 1987) restored the monotypic genus Megalochilus, but such taxonomic change was contradicted by Golubev and Sattorov (1992), as they argued that the differences proposed by Ananjeva were too slight to warrant a separate genus status. Molecular phylogenetic analyses based on mtDNA markers failed to resolve the phylogenetic position of Ph. mystaceus (Pang et al. 2003), which led Barabanov and Ananjeva (2007) to consider Megalochilus as a junior synonym of Phrynocephalus. However, a recent phylogeny based on the analysis of RAG1 nuDNA gene indicated Ph. mystaceus as a sister lineage with respect to all other examined Phrynocephalus species (Melville et al. 2009). Further study with better taxon sampling based on mtDNA data suggested that Ph. mystaceus is a member of the “core” Phrynocephalus clade and is associated with Ph. axillaris (Solovyeva et al. 2014). The most recent study proposed to consider Megalochilus as a subgenus of the genus Phrynocephalus (Solovyeva et al. 2014).

There was little consensus in the understanding of intraspecific taxonomy of Ph. mystaceus. Krassowky (1932) was the first to split Ph. mystaceus into two subspecies: European nominative subspecies Ph. m. mystaceus (Pallas, 1776) and Middle-Asian subspecies Ph. m. galli Krassowsky, 1932. This taxonomic classification was supported by subsequent studies of Soviet herpetologists (Shibanov 1941; Terentjev and Chernov 1949; Khonyakina 1961). However, morphometric studies by Vel’dre (1964a, 1964b) suggested that it is impossible to distinguish geographical races within Ph. mystaceus due to its high morphological variability among populations. Consequently, Ananjeva (1987 “1986”) suggested to upgrade the Middle-Asian subspecies Ph. m. galli to full species status and recognized a distinct subspeciesin Daghestan (Megalochilus mystaceus dagestanica in Ananjeva et al. 1987 “1986”). Semenov and Shenbrot (1990) analyzed morphological and chromatic differentiation of Ph. mystaceus from “Semirechye” (an area east of lake Balkhash in Eastern Kazakhstan), and suggested that this area is inhabited by a distinct subspecies, Ph. mystaceus aurantiacocaudatus Semenov et Shenbrot, 1990, which differs from the Middle Asian subspecies Ph. m. galli by its bright orange-red coloration of the ventral surface of the tail in young specimens (versus lemon-yellow coloration in other subspecies). However, Ph. mystaceus aurantiacocaudatus was synonymized with Ph. m. galli by Barabanov and Ananjeva (2007) without any discussion.

In summary, three subspecies of Phrynocephalus mystaceus are recognized in recent literature (see Barabanov and Ananjeva 2007):

1. Ph. m. mystaceus (Pallas, 1776), that inhabits eastern Ciscaucasia (eastern part of Chechen Republic, Daghestan, Kalmykia), Caspian region (southern part of Astrakhan Region, east of the Volga-Ural Sands; introduced to the Apsheron Peninsula, Azerbaijan) and northwestern Kazakhstan (Ananjeva et al. 2004). Terra typica restricta: Ryn-Peski (Ryn Sands), Ural Region, northwestern Kazakhstan (Barabanov and Ananjeva 2007). This form includes Megalochilus mystaceus dagestanica Ananjeva, “1986” 1987, described from Kumtorkala, Daghestan, Russia, as a junior synonym.

2. Ph. mystaceus galli Krassowsky, 1932, that inhabits Transcaspian Region and Middle Asia from Turkmenistan, Uzbekistan, Kazakhstan, to northeastern and eastern Iran and adjacent areas of Afghanistan (Anderson 1999; Ananjeva et al. 2004). Terra typica: Repetek station, Lebapsky District, Turkmenistan (Barabanov and Ananjeva 2007). Based on its distribution, this subspecies is supposed to inhabit north-eastern Iran (Anderson 1999).

3. Ph. mystaceus aurantiacocaudatus Semenov & Shenbrot, 1990, known from eastern Kazakhstan and western China (Ili River Valley in Xinjiang). Terra typica: 70 km north northwest of Ushtobe, Eastern Kazakhstan. Regarded as a junior synonym of Ph. mystaceus galli by Barabanov and Ananjeva (2007), however, without any justification.

It is notable that all previous works on geographic variations of Ph. mystaceus omitted populations from the southernmost edge of its range, Iran and Afghanistan, from the analyses. Morphological characterization and analysis of distribution of Ph. mystaceus in Iran was carried out by Anderson (1999) and Molavi et al. (2014). Anderson (1999) examined specimens from Iran and Uzbekistan, and proposed that Iranian populations demonstrate intermediate morphology between Ph. m. galli and Ph. mystaceus. Molavi et al. (2014), based on a study of seven specimens from Semnan Province, repeated earlier conclusions by Anderson (1999) and suggested that further investigation of both morphological and molecular characters are required to clarify the taxonomic status of Iranian Ph. mystaceus populations.

The recent analysis of phylogenetic relationships within the genus Phrynocephalus based on four mitochondrial genes revealed a remarkable divergence between Ph. mystaceus samples from Iran and Middle Asia (Solovyeva et al. 2014). Based on these results the Iranian population was tentatively indicated as a putative new subspecies Ph. mystaceus ssp. In the present study, we provide a detailed analysis of both morphological and genetic variation of Ph. mystaceus across its range and confirm deep differentiation between the population from Khorasan Province of Iran and other populations in the species range. The currently recognized subspecies of Ph. mystaceus are reviewed and a new subspecies from Khorasan Province is described, based on both molecular and morphological features.

Sampling. Historical collections of the Zoological Museum of Lomonosov Moscow State University (ZMMU) were examined, in total, 70 adult and subadult specimens of all currently recognized subspecies (Appendix 1). In addition, type specimens of Ph. mystaceus galli (lectotype, ZMMU R-6413) and Ph. mystaceus aurantiacocaudatus (holotype, ZMMU R-6412) were also examined. Sampling was carried out in the Khorasan Province of Iran in April of 2005, April of 2006, May and June of 2009, and May of 2010. Specimens from Iran were obtained through the collaboration with the Zoological Museum of International Center for Science, High Technology and Environmental Sciences (ICSTZM; Kerman, Iran; MOU no. 158/2010). Tissue samples from 31 Ph. mystaceus specimens were used in molecular analyses, and their geographic distribution is shown in Fig. 1. Details on museum IDs and localities of origin for each sample are summarized in Table 1.

Molecular analyses. Mitochondrial DNA COI gene (cytochrome oxidase c subunit I) fragment, 654 b. p. in length was analyzed. Muscle and skin tissues were disintegrated with Proteinase K and total genomic DNA was extracted using a standard phenol-chloroform extraction protocol followed by ethanol precipitation of DNA (Sambrook et al. 1989). PCR amplification was performed using MyCycler BioRad under conditions described by Ivanova et al. (2006). Standard pair of primers was used: VF1d (5'-TTCTCAACCAACCACAARGAYATYGG-3') and VR1d (5'-TAGACTTCTGGGTGGCCRAARAAYCA-3') or Rep-COI-F (5'-TNTTMTCAACNAACCACAAAGA-3') and Rep-COI-R (5'-ACTTCTGGRTGKCCAAARAATCA-3'). PCR reaction volume was 20 µl and it contained ca. 100 ng of template DNA, 0.3 pM/µl of each PCR primer, 1xTaq-buffer with 25 mM of MgCl₂ (Silex, Moscow Russia), 0.2 mM dNTPs, and 1 unit of Taq-polymerase (Silex, Moscow Russia; 5 units/µl). The results of the amplification were examined using electrophoresis in 1% agarose gel in presence of ethidium bromide. The length of the obtained fragments was 680 bp. We included two sequences of Ph. mystaceus from western China available from Genbank (NC022131 and KC578685; see Chen et al. 2014) in the analyses. Samples of Ph. melanurus (ZMMU R-12328, GenBank AN MF567976) and Trapelus sanguinolentus (ZMMU R-12709, GenBank AN KF691668) were used as outgroups.

Sequences were aligned using Seqman 5.06 and checked using BioEdit Sequence Alignment Editor 7.1.3.0 (Hall 1999). All sequences were deposited in GenBank (see Table 1 for all voucher information, with corresponding GenBank accession numbers).

Mean uncorrected p-distances and sequences characteristics were calculated using MEGA 6 (Tamura et al. 2011). Phylogenetic analyses were conducted using Treefinder (Jobb et al. 2011) and MrBayes 3.1.2 (Huelsenbeck and Ronquist 2001; Ronquist and Huelsenbeck 2003) software.

PartitionFinder v1.0.1 (Lanfear et al. 2012) was used to estimate the optimal evolutionary models for Bayesian inference analysis. The preferred model for COI alignment was HKY + G for two partitions (codon position 1 and 2 vs. codon position 3) as suggested by the Akaike information criterion (AIC). Bayesian phylogenetic analysis was performed using MrBayes v.3.1.2 (Ronquist and Huelsenbeck 2003) with two simultaneous runs, each with four chains, for 20 million generations, 2 million generations were cut as burn in. The convergence of the runs was checked to make sure that the effective sample sizes (ESS) were all above 200 by examining the likelihood plots using TRACER v.1.5 (Rambaut and Drummond 2007).

The Maximum Likelihood (ML) analysis was conducted using Treefinder (Jobb et al. 2011). Each dataset was divided into three partitions according to codon positions; for each partition the best fitting substitution model was selected using the AIC in Treefinder. For ML-analysis we used 1000 pseudoreplics (BS) and Expected Likelihood Weights (ELW).

Confidence in tree topology was tested by using non-parametric bootstrap analysis (Felsenstein 1985) with 1000 replicates and posterior probability (PP) for Bayesian inference (BA) in MrBayes 3.1.2 (Huelsenbeck and Ronquist 2001). Branches with bootstrap values of 70% or higher and posterior probabilities values over 0.95 were regarded as sufficiently resolved (Huelsenbeck and Hillis 1993).

Morphological analyses. Pholidosis was examined and morphometrics acquired for 79 individuals in four groups of Ph. mystaceus, including 20 specimens of nominative subspecies Ph. m. mystaceus, seven specimens from Khorasan Province of Iran, 32 specimens of Ph. m. aurantiacocaudatus from Eastern Kazakhstan, and 20 specimens of Ph. m. galli from Middle Asia (Appendix 1). In order to take into account sexual dimorphism, males (n = 26) and females (n = 44) were analyzed separately.

Morphological characteristics and the methods for their measurement are generally the same as in the study by Solovyeva et al. (2012). The following measurements and scalation counts were used: (1) snout-vent length (SVL); (2) tail length (TL); (3) SVL/TL ratio; (4) number of flat supralabials anterior to angular enlarged spine-like supralabial scales (SLbA); (5) total number of flat supralabials from tip of snout to insertion of cutaneous fold at mouth angle (SL); (6) relative length of the dark distal part of the tail to the total tail length (in ventral aspect, calculated as TL-black/TL ratio); (7) number of scales surrounding subnasal from below (SSbNb); (8) subnasal in contact with medial side of supranasal (vs. subnasal not in contact with medial side of supranasal) (SbN-SpN); (9) supranasal edges nostril dorsally along the full length of nostril (vs. supranasal edges nostril dorsally along only half of nostril length) (SpN); (10) height of supranasal is less than or equal to height of subnasal (vs. height of supranasal exceeds height of subnasal) (hSpN SbN); (11) number of scale rows that separate subnasal and labial scales (SbN-L); (12) longitudinal row of white scales in supraorbital area outlined by continuous black lines (or intermitted) (WS&BL); (13) number of small rows of scales between anterior (2d and 3d) inframandibulars and large rows of scales under infralabial scales – 1-2 or 2-3 (aIMd-IL); (14) number of scales that underlay enlarged spiny scales on edge of cutaneous fold at mouth angle (SuSSCF); (15) number of small granular scales between posteriormost supralabial and insertion of cutaneous fold at mouth angle (pSL-CF); (16) number of flat infralabials anterior to angular enlarged spine-like infralabial scales (ILbA); (17) total number of infralabials from tip of snout to insertion of cutaneous fold at mouth angle (IL); (18) number of subdigital lamellae under toe III (SLIII); (19) number of enlarged triangular scales on lateral fringes of toe III (FrIII); (20) number of subdigital lamellae under toe IV (SLIV); (21) number of enlarged triangular scales on lateral fringes of toe IV (FrIV). Characteristics 18-21 (SLIII, FrIII, SLIV, FrIV) were registered with no regard to the sex of the individual. Following standard measurements were additionally taken for holotype and paratypes: head height (HH); head length (HL, measured on ventral side from snout tip to gular fold); head width (HW, measured at broadest part of head excluding cutaneous folds); pileus width (PW). Measurements were taken using a digital caliper and rounded to the nearest 0.1 mm.

Box-and-whiskers-plots and values of descriptive statistics were calculated using R (R Core Team, 2013). The Mann-Whitney test of independent series was used to determine the differences between the pairs of subspecies (with confidence level of p ≤ 0.05). Principal Components Analysis (PCA) was performed using R (R Core Team, 2013) to visualize morphological variation between Khorasan specimens of Ph. mystaceus and specimens from other populations.

Sequence characteristics. The sequenced fragments from 31 Ph. mystaceus specimens were up to 654 b.p. in length, among which 577 sites were identified as conservative, 74 as variable and 59 as potentially parsimony-informative. Nucleotide frequencies were equal to: 30.2% (A), 27.5% (T/U), 27.9% (C), and 14.4% (G). The transition-transversion bias (R) was estimated to be 6.574 (all data given for in-group only).

Phylogenetic analysis. The results of phylogenetic analysis are presented in Fig. 2. BI and ML yielded trees that show essentially similar topologies. All analyses reveal the presence of two reciprocally monophyletic clades within Ph. mystaceus. The first clade consists of Irainan Ph. mystaceus ssp. from Khorasan Province (node support values are 1.0/86; hereafter given for BI PP/ ML BS; clade I on Fig. 2). The second clade includes all other Ph. mystaceus populations from Middle and Central Asia and Caspian Region (1.0/99; clade II on Fig. 2). Further phylogenetic structure within the second clade of non-Iranian Ph. mystaceus is poorly resolved. Populations from the eastern part of the range including Eastern Kazakhstan and Xinjiang (China) that correspond to the Ph. m. “aurantiacocaudatus” occupy basal position in the clade II, but are not monophyletic and fall into three poorly differentiated subclades: from the environs of Kapchagai (subclade A; 0.90/82), Ili River Valley (subclade B; from Zharkent to Xinjiang; 1.0/95), and the environs of Lake Balkhash (subclade C; 0.98/-) (see Fig. 2). Phylogenetic positions of two samples from Eastern Kazakhstan (ZMMU NAP-05510 and ZMMU R-12518-2) are not resolved. All other populations from Middle Asia (Kazakhstan and Uzbekistan – Fig. 2D) and Caspian Region (Astrakhan Province, Russia – Fig. 2E) form a significantly monophyletic clade (1.0/95), which is deeply nested within the basal differentiation of East Kazakhstan Ph. m. “aurantiacocaudatus” clades (see Fig. 2), rendering the latter taxon paraphyletic. The Middle Asian – Caspian clade (D + E) corresponds to the nominative subspecies Ph. mystaceus mystaceus and also includes populations previously classified as Ph. mystaceus galli (Aral Sea Region, Uzbekistan and Tajikistan – Fig. 2D). Populations of Ph. mystaceus mystaceus and Ph. mystaceus “galli” are mixed with each other without any clear structure (see Fig. 2 D1, D2, E1, E2).

Genetic distances. Uncorrected genetic p-distances between and within clades of Ph. mystaceus are shown in Table 2. The p-distances within the Middle Asian – Caspian clade of Ph. mystaceus, including comparisons between different lineages of Ph. m. “aurantiacocaudatus” and between Ph. m. “aurantiacocaudatus” and Ph. m. mystaceus are quite low (0.55–0.88% and 1.56–1.87%, respectively), which is less than intraspecific genetic distances for COI for some other species of Phrynocephalus (e.g. see Solovyeva et al. 2011 for Ph. helioscopus). However, p-distances between Ph. mystaceus ssp. from Khorasan Province and all other groups of Ph. mystaceus are very high (6.84–7.28%), they even exceed interspecific genetic distances for COI gene reported for certain species of Phrynocephalus (Solovyeva et al. 2014). This data clearly suggest a deep divergence between Ph. mystaceus populations from Khorasan Province and populations from the rest of the range of the species.

Morphology. Our study supports the results of previous researchers that indicated very high morphological variation in the absence of consistent morphological variation patterns that could delimit recognized subspecies in Middle Asian populations of Ph. mystaceus (Vel’dre 1964a, 1964b; Semenov and Shenbrot 1990). Most characteristics, including body size, were uninformative for distinguishing subspecies and local populations of Ph. mystaceus. Only four morphological characteristics showed consistent differences between Iranian and Middle Asian/Caspian populations of Ph. mystaceus, including SLIV, FrIII, SL, and TL-black/TL. Specifically, SLIV was lower in the population from Khorasan Razavi Province (N = 7) than in other subspecies of Ph. mystaceus (differences are significant; p = 0.000 for comparison with Ph. h. “aurantiacocaudatus”, N = 32; p = 0.000 for comparison with Ph. h. “galli”, N = 20; p = 0.087 for comparison with Ph. m. mystaceus sensu stricto, N = 20; for measurement ranges see Table 3). FrIII was also significantly lower in Khorasan population (N = 7) than in Ph. m. mystaceus sensu stricto (p = 0.000 N = 20). SL was also lower in Khorasan population (N = 7) than in other subspecies (p = 0.007 for comparison with Ph. m. “aurantiacocaudatus”, N = 32; p = 0.001 for comparison with Ph. m. mystaceus sensu stricto, N = 20; p = 0.050 for comparison with Ph. m. “galli”, N = 20). Finally, the dark distal part of the tail (TL-black/TL) was relatively longer in the Khorasan population (differences are significant; p = 0.023, for comparison with Ph. m. “aurantiacocaudatus”, N = 32; p = 0.000 for comparison with Ph. m. mystaceus sensu stricto, N = 20; p = 0.001 for comparison with Ph. m. “galli”, N = 20). Morphological comparison of four geographical population groups that correspond to the subspecies “mystaceus sensu stricto”, “galli”, “aurantiacocaudatus” and the “Khorasan population” for the diagnostic morphological characteristics mentioned above is given in Fig. 3. Other characteristics with p-values for pairwise comparisons <0.05 showed significant overlap of values between subspecies and cannot be reliably used in diagnostics; p-values for morphological characteristics for pairwise comparisons are summarized in Appendix 2. Standard measurements of Ph. mystaceus ssp. from Khorasan Province are presented in Table 4.

Comparison of Khorasan Ph. mystaceus ssp. population with other populations of Ph. mystaceus from Middle Asia and Caspian region (data for “mystaceus sensu stricto”, “galli”, “aurantiacocaudatus” combined together) also demonstrated significant differences for many traits with p<0.05 (Appendix 2), however, values for most of them were overlapping. The six of the characters with the minimal overlap were the following: SVL/TL, TL-black/TL, SL, ILbA, SLIII, SLIV (see Fig. 4 for details).

PCA showed differences between Khorasan population and other Ph. mystaceus populations, although these two groups are slightly overlapping with two Khorasan specimens falling into the Ph. m. mystaceus sensu lato area (Fig. 5). PCA failed to reveal any clear structuring within the Middle Asian / Caspian populations of Ph. mystaceus.