Does the Dzungarian racerunner (Eremias dzungarica Orlova, Poyarkov, Chirikova, Nazarov, Munkhbaatar,Munkhbayar & Terbish, 2017) occur in China? Species delimitation and identification with DNA barcoding and morphometric analyses

TheEremias multiocellata-przewalskiispecies complex is a viviparous group in the genusEremias, and a well-known representative of taxonomically complicated taxa. Within this complex, a new species -E. dzungarica(Orlova et al., 2017)- has been described recently from western Mongolia and eastern Kazakhstan, with an apparent distribution gap in northwestern China. In this study, we used an integrative taxonomic framework to address whetherE. dzungaricaindeed occurs in China. Thirty specimens previously classified asE. multiocellatawere collected in eastern Kazakhstan and the adjacent Altay region in China. The cytochromecoxidase I (COI) barcodes were sequenced and compiled with those from Orlova et al. (2017) and analyzed with the standard and diverse barcoding techniques. We detected an absence of a barcoding gap in this complex, which indicates potential cryptic species inEremiassp. 3 with high intraspecific diversity and multiple recently evolved species in Clade A.Both BIN and GMYC suggested an unrealistically large number of species (23 and 26, respectively), while ABGD,mPTP and BPP indicated a more conservative number of species (10, 12, and 15, respectively), largely concordant with the previously defined species-level lineages according to phylogenetic trees. Based on molecular phylogeny and morphological examination, all 30 individuals collected in this study were reliably identified asE. dzungarica- a distinct species - confirming the occurrence of this species in the Altay region, Xinjiang, China. Potentially owing to the larger sample size in this study, our morphological analyses revealed many inconsistencies with the original descriptions ofE.dzungarica, which were primarily associated with sexual dimorphism and a broader range of values for various traits.

The rapid development of DNA barcoding (Hebert et al.,2003) has facilitated the successful application of a standardized short mitochondrial gene fragment,COI, to most species level identifications (e.g., excluding plants), species discovery and global biodiversity assessment (DeSalle &Goldstein, 2019; Yang et al., 2020). DNA barcoding is particularly helpful for phylogenetic and taxonomic inference in species groups that have considerable morphological conservatism or ambiguity (e.g., Hofmann et al., 2019; Oba et al., 2015; Xu et al., 2020; Zhang et al., 2018).

Traditional taxonomy is mainly based on morphological characters, and hence can easily lead to misidentification as a result of phenotypic plasticity, cryptic species or different morphologies at different life history stages (Bickford et al.,2007; Lee, 2004; Rock et al., 2008). Taxonomy relying on DNA barcodes alone is also unrealistic, as mitochondrial genes have many inherent biases and limitations in species delimitation, e.g., those associated with maternal inheritance,reduced population size, inconsistent mutation rate, or evolutionary processes such as purifying selection (Blair &Bryson, 2017; Pino-Bodas et al., 2013; Rubinoff et al., 2006).With the advance of barcoding techniques, however, the evidence inferred from DNA barcoding can guide further targeted morphological analysis, and the use of multiple lines of evidence, such as nuclear loci, geographical and ecological data to make more robust inferences about the species boundaries under the rubric of integrative taxonomy (Damm et al., 2010; Dayrat, 2005; Miller, 2007; Padial et al., 2010; Will et al., 2005).

TheEremias multiocellata-przewalskiispecies complex comprises a natural group of viviparous species in the genusEremias(Guo et al., 2011; Orlova et al., 2017). The taxonomy of this species complex has been historically confusing due to the vast phenotypic variation within and among species, as well as the conservation of morphological characters in closely related species (Eremchenko et al., 1992; Eremchenko &Panfilov, 1999). So far, as many as 18 species/subspecies have been proposed across its wide distribution range covering Kyrgyzstan, eastern Kazakhstan, northern China,Mongolia, and southern Tuva Republic of Russia (Orlova et al., 2017; and references therein). Among these species is the newly delimited Dzungarian racerunner,Eremias dzungarica(Orlova et al., 2017). In addition to the apparent molecular and morphological deviations from congeners as described in Orlova et al. (2017), this species is characterized by a habitat preference for rocky hills and gravel ravines (“rocky form"coined in Orlova et al. (2017)) in western Mongolia at high elevations (2 400-2 600 m above sea level (a.s.l.)). However,it can also penetrate into low-altitude sandy dune areas in eastern Kazakhstan (400-1 000 m a.s.l.). As such, it remains unclear as to whetherE. dzungaricaoccurs in the vast territories of the northern Junggar Depression in Xinjiang,China, between western Mongolia and eastern Kazakhstan.To date, four occurrences of so-calledE. multiocellata(the multio-cellated racerunner) have been recorded from only two regions in the northern Junggar Depression: one in the Tacheng region and three in the Altay region reported in Zhao(1999) and Tao et al. (2018), respectively. As suggested by Orlova et al. (2017),E. dzungaricamay have been considered asE. multiocellatain China, hence it is possible that these reported populations are in factE. dzungarica, despite the lack of morphological and molecular data.

Orlova et al. (2017) for the first time utilized the DNA barcoding sequences (COI) to infer phylogenetic relationships and propose putative species in theE. multiocellataprzewalskiispecies complex based on mitochondrial lineages,incomplete morphological identification characters (e.g., no voucher specimens from certain lineages) and geographic distributions. However, they did not utilize more rigorous barcoding techniques to deeply explore the distribution of genetic distance and to test species boundaries in this species complex. To determine whetherE. dzungaricaoccurs in China, we sequenced the DNA barcodingCOIfragments and performed morphological measurement of the 30 purportedE.multiocellataindividuals collected from seven locations in eastern Kazakhstan and the adjacent western Altay region,Xinjiang, China (Figure 1A; Supplementary Table S1).Subsequently, sequences of the species complex from Orlova et al. (2017) were compiled and analyzed with diverse commonly used barcoding methods to explore the intra- and interspecific genetic distance patterns and reassess the species status of the species-level lineages proposed by Orlova et al. (2017). Finally, we explicitly identified the taxonomic status of the “multi-ocellated racerunners” collected from Kazakhstan and China with molecular and morphological data. Detailed methods are available in the Supplementary Materials and Methods.

Sequences from Orlova et al. (2017) had different lengths;most of the sequences (82.3%) were 651 bp, however some were 617-619 bp with missing data located at both ends of the sequences. To accommodate the majority of these sequences, a dataset of 651 sites was generated. A total of 81 haplotypes were determined, including four for the outgroups.235 sites were variable, and 190 were phylogeny-informative.

The general trends of Kimura 2-parameter genetic distance indicated an increment from lower to higher taxonomic levels(Supplementary Figure S1). However, the genetic distances at the species and genus level overlapped at a low frequency(Supplementary Figure S1), indicating no apparent barcoding gap between intra- and interspecific distances. Intraspecific distances ranged from 0 to 6.18% (mean±standard deviation(SD): 1.16%±0.91%, Supplementary Table S2); only Eremias sp. 3 exhibited variations over 3.0%, with a total frequency lower than 1.35% (Supplementary Figure S1). Interspecific divergences were highly variable (mean±SD: 9.09%±2.54%,Supplementary Table S2), ranging from extremely low between Eremias sp. 2 and E. cf. reticulata (0.77%) to remarkably high between Eremias sp. 3 and E. stummeri(13.19%); only species within Clade A (Figure 1B) exhibited variations ＜2.0%, with a total frequency lower than 0.62%(Supplementary Figure S1). Barcoding gap analysis indicated that the maximum intraspecific distance of each species was not always higher than the minimum distance to its nearest neighbor (Supplementary Figure S2). This evidence suggests the absence of a barcode gap. Three species (i.e., Eremias sp. 3, E. przewalskii and E. dzungarica) had lower distances to their nearest neighbor than their maximum intraspecific distances (Supplementary Figure S2 and Table S3). Of these,E. dzungarica exhibited moderately high maximum intraspecific distance (2.99%), while the interspecific distance to its nearest neighbor (E. cf. buechneri) was lower (1.87%;Supplementary Table S3).

Figure 1 Collection sites of Eremias multiocellata-przewalskii species complex samples and phylogenetic relationships and species delimitation

The phylogenetic trees reconstructed with three different methods (i.e., Bayesian inference (BI), Neighbor-Joining (NJ)and Maximum Likelihood (ML)) resulted in nearly consistent topologies (Figure 1B; Supplementary Figures S3, S4), which is also congruent with that in Orlova et al. (2017). Most clades(representative of species-level lineages) were recovered with high support from all analyses, except forEremiassp. 3 in which the Bayesian posterior probability was only moderate(BPP=0.84); the bootstrap support proportions (BSP) in NJ and ML were similarly moderate as well (74 and 59,respectively). The monophyly of theE. multiocellataprzewalskiispecies complex recovered here was remarkably lower (BPP=0.88) than in Orlova et al. (2017; BPP=0.97), but the NJ and ML trees recovered significantly high support for monophyly (BSP=98 and 95, respectively). More importantly,the 30 representative specimens (represented by the light green branches on the phylogenetic trees; Figure 1B and Supplementary Figures S3, S4) previously identified asE.multiocellatawere explicitly nested withinE. dzungaricawith strong support (BPP=1.0; BSP=98 and 96 in NJ and ML,respectively), indicating that these taxa could be allocated to the recently described new species.

The distance (BIN)- and tree (GMYC)-based methods both suggested an unrealistically large number of putative species.The application of BIN in our barcoding dataset identified 23 operational taxonomic units (OTUs), of which only six were taxonomically concordant with previously defined specieslevel lineages in the BI/NJ/ML gene trees (Figure 1B). The single-threshold GMYC model indicated a multiple species scenario with strong statistical support (P＜0.047). The GMYC analysis results - which identified as many as 26 entities -were largely incongruent with the BI/NJ/ML gene trees regarding the number of species-level lineages (Figure 1B).However, consistent with the results of genetic distance analysis based on predefined species, no apparent barcoding gap was detected in ABGD analysis based on pairwise comparisons of sequences across the dataset (Supplementary Figure S5A). ABGD has been suggested to be a powerful tool to partition the barcoding datasets into putative species, even when intra- and interspecific genetic distances overlap(Puillandre et al., 2012). Our results support this perspective,as ABGD conservatively delimited 10 putative species in the initial and recursive partitions (i.e., 10 initial and 10-12 recursive partitions with prior intraspecific divergences, which varied from 0.1% to 0.93%; Supplementary Figure S5B), half of which were taxonomically concordant (Figure 1B). BPP analyses suggested 14 or 15 species with relatively low(BPP=0.238) and high (BPP=0.735) support, respectively.Within the 15 genetic clusters, four (i.e.,Eremiassp. 2,E.cf.reticulata,E.cf.buechneriandE. yarkandensis) were delimited with relatively high support (BPP＞0.86 and ＜0.9) and others with even stronger support (BPP＞0.95), whileE.cf.buechneriandE. yarkandensismay form one putative species with extremely low support (BPP=0.13). These 15 delimited genetic clusters were largely consistent with the species-level lineages in the BI/NJ/ML gene trees, except thatEremiassp. 3 was delimited into three genetic clusters, consistent with the results in the ABGD analysis (Figure 1B). The mPTP model exhibited relatively conservative species delimitation for our dataset with highly uneven sampling. Twelve putative species were suggested by the mPTP, but only four of them were consistent with the species-level lineages in the BI/NJ/ML gene trees (Figure 1B). Both the ABGD and BPP indicated three consistent cryptic species, and the mPTP suggested two putative species inEremiassp. 3, whereas both the BIN and GMYC split this species into as many as five different genetic structures. One singleton cryptic species (Hap 58) from population 36 in eastern Inner Mongolia (Figure 1A;Supplementary Table S1) was recognized in all analyses,indicating that its taxonomic status deserves particular attention in future studies with more robust data. Moreover,the cryptic diversity inEremiassp. 3 was also supported with high intraspecific pairwise divergence (3.49%-6.18%) and large morphological variations that may have been overlooked by Orlova et al. (2017) according to our field observations in Inner Mongolia.

Despite the low interspecific pairwise genetic distances(＜2%) among the lineages in Clade A (including seven species-level lineages in the BI/NJ/ML gene trees; Figure 1B),Orlova et al. (2017) considered these as distinct species for the following reasons. First, these lineages in Clade A were apparently morphologically differentiated (exceptE.przewalskiiandE.cf.przewalskii). Second, these lineages covered separate geographic distributions. Similarly, our BPP analysis also suggests that all lineages in Clade A are putative species, partially supported by the BIN and GMYC analyses that congruently suggest the species status forE.cf.

buechneri,E. yarkandensisandE.cf.przewalskii(Figure 1B).Therefore, we propose that the divergence among the lineages in Clade A is indicative of multiple recent speciation events. However, Orlova et al. (2017) did not propose any interpretation for the controversial status ofE.cf.przewalskii,which morphologically resemblesE. przewalskii, although in our studyE.cf.przewalskiiwas not nested withinE.przewalskiiand did not form a sister taxon toE. przewalskiiin the phylogenetic trees (Figure 1B). There are two possible explanations for this discrepancy. Firstly,E.cf.przewalskiimay be a cryptic species that has evolved morphological traits similar toE. przewalskii. Alternatively,E.cf.przewalskiimay in fact beE. przewalskii, and the mitochondria of this population(site 29) had been replaced by an unknown species from theE. multiocellata-przewalskiispecies complex. With current data, however, we cannot rule out either of these hypotheses.As such, more thorough fieldwork and rigorous morphometric analyses and additional molecular data (e.g., nuclear loci) are needed for future taxonomic and evolutionary hypothesis testing for this species group.

The status ofE. dzungaricais the most contradictory in our species delimitation analyses. Although BPP and mPTP explicitly suggested its species status, BIN and GMYC split it into many different genetic clusters while ABGD merged it withE.cf.buechneriandE. yarkandensisas a single putative species. Given that the monophyly of most genetic clusters inE. dzungaricadelimited by BIN and GMYC was not strongly supported (BSP＜70; BPP＜0.9), the intraspecific differentiation inE. dzungarica(relatively high; maximum intraspecific divergence of 2.99%) may not have been high enough to form an independent evolutionary lineage. On the other hand, the relationships amongE. dzungarica,E.cf.buechneriandE.yarkandensiswere not resolved in the BI and ML trees(Supplementary Figures S3, S4). Thus, the ABGD results may be affected by the reciprocally nearest species among them(Supplementary Table S3). In addition, the BIN, GMYC and BPP analyses consistently suggestedE.cf.buechneriandE.yarkandensisas independent species (Figure 1B); these results, in turn, suggest thatE. dzungaricashould be allocated as a distinct species. Finally, the morphological concordance among the populations ofE. dzungaricalends further support for its distinct species status.

The range of values for most metric and meristic traits of the specimens from this study was largely consistent with that of theE. dzungaricaindividuals in Orlova et al. (2017). The only exception was the metric trait Dist.P.fm, with a larger range of values found in our study (3.52-6.34 mm) than that reported in Orlova et al. (2017) (1.60-2.50 mm; Supplementary Tables S4.1, S4.2). Moreover, statistical tests indicated that the specimens in this study were significantly different than those in Orlova et al. (2017) in four traits (Lab.total.R, Ventr.,Lam.subdig. and P.fm.L; Supplementary Table S5). However,the individuals in this study can further be morphologically identified asE. dzungaricabased on the following combination of characters: single frontonasal; two prefrontals; subocular shield not in contact with mouth margin, in touch with 6th-8th supralabials; 3-5 subocular shields; three pairs of nasals;subnasal not in contact with rostral shield, located above 1st to 3rd supralabials; two loreal shields, except one individual from Kazakhstan (Voucher No. KZL98; Supplementary Table S1) with single loreal at either side of head; 5-6 submaxillary shields at right or left side; first three pairs of submaxillary shields in contact with each other, no or minor split between last contacting pair of submaxillary shields; last submaxillary shield in contact with infralabials in certain individuals (26.7%);supraoculars separated completely from supraciliary shields by single row of granular scales in some individuals (60%), but partly in contact with supraciliary in other individuals due to deficiency of certain granular scales in the row; supraoculars in contact with frontal and frontoparietals without granular scales between them; 3-4 scales from distal femoral pore to knee; one or two explicitly enlarged shields in precloacal(preanal) area; background of dorsum and head dorsal surface grayish-brown; head dorsal surface with many (usually in males) or few (usually in females) random irregular black blotches; black blotches on ventral flanks forming two regular longitudinal rows, ventral sides near black blotch rows with sparse yellowish spots in some individuals (usually males)(Figure 2).

It should be noted that there are also some inconstancies between the morphological descriptions reported here and the original descriptions ofE. dzungaricain Orlova et al. (2017).One of the most important findings in this study is that the sexual dimorphism related to the dorsal coloration pattern reported in Orlova et al. (2017) may be unreliable. Males are distinguished by a bright green-yellowish coloration at the third row of ocelli, while not all females lack their bright color at the third row of ocelli, and many individuals even show the same well-developed ocelli with bright color as males (Figure 2A1,A2, A4). In fact, based on our field observations, we suspect that these dorsolateral ocelli coloration patterns may be related to age, for the following reasons. Firstly, both males and females displayed the bright ocelli before adulthood.Secondly, males may retain and develop the bright ocelli to attract females for mating, while females may gradually lose them since the bright ocelli may have no benefits for them.Thirdly, losing the bright color in adult females may help them avoid predators. Lastly, for the specimens examined in this study, all males had two apparent rows of whitish ocelli near the mid-dorsum (Figure 2B1-B3), while the whitish ocelli in the first row near the mid-dorsum was blurry in some females(Figure 2A1-A3). Taken together, these phenomena may indicate another unreliable sexually dimorphic trait related to the dorsal patterning inE. dzungarica.

Our morphometric analyses, however, showed apparent sexual dimorphism in body size, with males usually larger than females (Supplementary Tables S4.1, S4.2). Although one metric trait (Dist.P.fm) was not sexually dimorphic, it was closely related to the number of scales between two femoral pore series (scal.f.p) instead of body size. While the values of Dist.P.fm measured in Orlova et al. (2017) were substantially lower than those in this study, there is no explicit evidence to explain why such significant deviations occur. Moreover,Orlova et al. (2017) did not report sexual dimorphism in any of the meristic traits ofE. dzungarica, whereas we found a single meristic trait (Sq.c.cd) with significant sexual dimorphism. We also observed other inconsistencies in one individual (Voucher No. KZL98, Supplementary Table S1) from population 56 in Kazakhstan (Figure 1A; Supplementary Table S1) with a single loreal shield on both sides of the head and a broader range of values for many different traits (Supplementary Tables S4.1-S4.3). In general, these findings could be attributed to individual variation, like the fusion of two loreal shields into one, and/or to the larger sample size in this study than that in Orlova et al. (2017).

Our sampling sites in China are located in the western Altay region in Xinjiang, which is close to the known occurrences ofE. dzungaricain eastern Kazakhstan (Orlova et al., 2017).There is still a large degree of uncertainty about the existence of this species between the western Altay region and the known occurrences ofE. dzungaricain western Mongolia near the Altay Mountains (Figure 1A). For example, Tao et al.(2018) reported three occurrences ofE. multiocellatain the Altay region, one located in the west, close to our sampling localities, and two in the central and eastern parts. Although the precise taxonomic assignment of these populations is unknown, we suspect they could be allocated toE.dzungarica, given that this species may have a continuous distribution range from eastern Kazakhstan to western Mongolia. The habitats ofE. dzungaricaare described in Orlova et al. (2017) as rocky hills and ravines at elevations up to 2 400-2 600 m a.s.l. in Mongolia, and sandy dunes(400-600 m a.s.l.) and occasional rocky outcrops (1 000 m a.s.l.) at lower elevations in Kazakhstan. Consistent with these habitat descriptions, the populations from Kazakhstan sampled here were also associated with sandy dunes (～400 m a.s.l.). The habitats of individuals sampled from the western Altay region in China were more diverse, and included sandy dunes (sites 58, 59 and 61; 420-580 m a.s.l.; Figure 1A) at similar elevations, rocky outcrops (site 62; ～620 m a.s.l.), and rocky ravines (site 60; ～1 200 m a.s.l.; Figure 1A) at higher elevations. Considering thatE. dzungaricacan inhabit a wide range of altitudes, future fieldwork on this species should be conducted in the vast low-elevation territories from eastern Kazakhstan to western Mongolia, as well as the high-elevation territories in the Altay Mountain areas that span all three countries.

Figure 2 General view of wild specimens from field recordings in eastern Kazakhstan and western Altay region, China

SCIENTIFIC FIELD SURVEY PERMISSION INFORMATION

Permission for field surveys was granted by the Forestry Department and National Reserves of China and Committee for Forestry and Hunting of the Ministry of Environmental Protection of Kazakhstan.

SUPPLEMENTARY DATA

Supplementary data to this article can be found online.

COMPETING INTERESTS

The authors declare that they have no competing interests.

AUTHORS’ CONTRIBUTIONS

X.G.G and J.L.L. designed the study. J.L.L., T.N.D., M.A.C.,X.G., and D.J.L. collected specimens in the field. J.L.L., X.G.,and D.J.L. performed molecular experiments. J.L.L. and M.A.C. measured the specimens. J.L.L. performed data analyses. J.L.L., X.G.G and T.N.D. wrote and revised the manuscript. All authors read and approved the final version of the manuscript.

ACKNOWLEDGEMENTS

We thank three anonymous referees for insightful comments and Dr. Jacquelin DeFaveri for language editing.