Phylogeography and morphological evolution of Pseudechiniscus (Heterotardigrada: Echiniscidae)

Tardigrades constitute a micrometazoan phylum usually considered as taxonomically challenging and therefore difficult for biogeographic analyses. The genus Pseudechiniscus, the second most speciose member of the family Echiniscidae, is commonly regarded as a particularly difficult taxon for studying due to its rarity and homogenous sculpturing of the dorsal plates. Recently, wide geographic ranges for some representatives of this genus and a new hypothesis on the subgeneric classification have been suggested. In order to test these hypotheses, we sequenced 65 Pseudechiniscus populations extracted from samples collected in 19 countries distributed on 5 continents, representing the Neotropical, Afrotropical, Holarctic, and Oriental realms. The deep subdivision of the genus into the cosmopolitan suillus-facettalis clade and the mostly tropical-Gondwanan novaezeelandiae clade is demonstrated. Meridioniscus subgen. nov. is erected to accommodate the species belonging to the novaezeelandiae lineage characterised by dactyloid cephalic papillae that are typical for the great majority of echiniscids (in contrast to pseudohemispherical papillae in the suillus-facettalis clade, corresponding to the subgenus Pseudechiniscus). Moreover, the evolution of morphological traits (striae between dorsal pillars, projections on the pseudosegmental plate IV’, ventral sculpturing pattern) crucial in the Pseudechiniscus taxonomy is reconstructed. Furthermore, broad distributions are emphasised as characteristic of some taxa. Finally, the Malay Archipelago and Indochina are argued to be the place of origin and extensive radiation of Pseudechiniscus.

Tardigrades represent a group of miniaturised panarthropods 1 , which is recognised particularly for their abilities to enter cryptobiosis when facing difficult or even extreme environmental conditions 2 . Their relationships with Arthropoda and Onychophora are a subject of long-standing debate 3 , although the sister position with respect to these two lineages is the most popular hypothesis 4 . Internal tardigrade relationships are enigmatic likewise, despite the fact that the effort to decipher their phylogeny has been on an increase in the last decade [5][6][7][8] . The insufficient understanding of evolutionary relationships between tardigrade species (Darwinian shortfall) is not the only gap in our knowledge 9 , as it results from the scarcity of described tardigrade species 10 (Linnean shortfall). The other parallel problem is that tardigrade distribution on the globe is vastly unknown 11 (Wallacean shortfall), and many old records are actual misidentifications caused by poor understanding of how to dissect intraspecific and interspecific variability [12][13][14] . Currently, tardigradologists aim at unravelling biodiversity patterns within the phylum with the application of modern molecular tools, i.a. [15][16][17] .
One of the most taxonomically absorbing conundrums of recent years concerns the genus Pseudechiniscus 18 , the echiniscid genus extremely scarce in morphologically informative traits 19 after the erection of Acanthechiniscus to accommodate the distinct evolutionary history and morphology of some species earlier attributed to Pseudechiniscus 20 . The recent disclosure of the molecular diversity within Pseudechiniscus 21 , with an overview of morphological criteria [21][22][23][24] and an integrative re-description of the nominal species Pseudechiniscus suillus 23 constituted the turning point in the classification of this taxon. Perhaps the two most significant findings of these studies were the discovery of the crucial taxonomic role of ventral sculpturing patterns [22][23][24] , independently shown also for Hypechiniscus 25 , and the distinguishing of two phylogenetic lineages concordant with the shape of the cephalic papillae (secondary clavae) 21 , which is a unique state in the entire family, as other echiniscid genera are homogeneous in this criterion 19 . The suillus-facettalis lineage exhibits pseudohemispherical (illusively dome-shaped) papillae (e.g. P. insolitus or P. jubatus, Fig. 1A-B), whereas members of the novaezeelandiae line exhibit dactyloid papillae (e.g. P. novaezeelandiae or P. juanitae, Fig. 1C-D). Cesari et al. 21 implied that these lines correspond with potential subgenera, simultaneously noting that the genetic evidence was insufficient at

Results
Phylogeny and taxonomy. Overall, both Bayesian and Maximum Likelihood approaches gave wellresolved and largely congruent phylogenies (Fig. 2). The monophyly of Pseudechiniscus received a maximal support, and the genus was clearly divided in two large, deeply divergent clades, coherent with the suillus-facettalis ( Fig. 1A-B) and novaezeelandiae morphotypes (Fig. 1C-D); see identical division in the COI phylogeny (Supplementary Material 2). Notably, branches in the latter are visibly longer than in the suillus-facettalis line. Moreover, the only inconsistent topology occurred in the novaezeelandiae line, in which ML recovered the following arrangement of species: ((((P. cf. angelusalas (((P. sp. 6 ((P. cf. saltensis (P. sp. 7 + P. sp. 8)))), but with weak support. On the basis of ITS-1 tree, the number of species estimated by bPTP varied between 23 and 43, with i.e. 1.00 for BI and 100 for ML, are indicated by asterisks. Support values for tip nodes (within intra-specific phylogenetic structure) are not shown for simplicity. Topologies of all consensus trees were identical, with the exception of a single clade, in which ML gave incongruent results (marked with hashtags). The scale refers to the BEAST consensus tree. Values in round brackets and in the superscript at each species name refer to the support they received in the bPTP delimitation method (maximum values, i.e. 1.00, are indicated by asterisks www.nature.com/scientificreports/ the average of 30. Both ML and the simple heuristic search indicated the existence of 26 species (Fig. 2), which was congruent with classical taxonomic identification. Remarks: Tumanov 22 indicated dubious character of P. clavatus. Not only the shape of both primary and secondary clavae is atypical, but also the pseudosegmental plate IV' is unidentifiable in the original description based on juvenile specimens 33 . Therefore, we ascertain that P. clavatus is a nomen dubium and this name should not be used in modern literature. We are in agreement with Tumanov regarding the status of P. jiroveci sp. dub. and P. marinae. Contrarily to Grobys et al. 23 , we are of the opinion that P. megacephalus does not belong to a different genus. The alleged mushroom-shaped cephalic papillae with a peduncle-like bases are unknown in any other echiniscid species 19 , making this trait most likely artefactual rather than autapomorphic 34 . Consequently, P. megacephalus is designated as a nomen dubium. Pseudechiniscus pulcher should be included within Antechiniscus 22 . P. transsylvanicus has long cirri C and the body ca. 350 μm long 35 , that is much larger than the typical length for Pseudechiniscus, which can reach up to 250 μm, but usually is below 200 μm. Moreover, cirriform trunk appendages are absent in Pseudechiniscus s.s. 22 , thus P. transsylvanicus almost certainly belongs in a different echiniscid genus, but the currently available data preclude its transfer. The same is for P. bispinosus characterised by two rigid, long spines C 36 . Pseudechiniscus shilinensis has an inadequate description preventing its taxonomic identification 37 , and is designated here as a nomen dubium. Finally, the following species are unassignable to subgenera as their descriptions do not specify the cephalic papillae shape: P. dicrani and P. marinae 38,39 , thus they are tentatively retained in the subgenus Pseudechiniscus until type or fresh material is available.
Diagnosis: As for the nominal subgenus, except the cephalic papillae which are dactyloid, and attached to the body cuticle only at their bases. Striae present, usually large and well-developed.
Etymology: From Latin meridionalis = southern + ending -iscus derived from Echiniscus, the first established echiniscid genus, literally meaning "an echiniscid from the South". The name refers to the geographic origin of the known species belonging to the new subgenus, that most likely has its roots in the Southern Hemisphere and tropical/subtropical regions of the world.
Remarks: Pseudechiniscus angelusalas, P. dastychi and P. indistinctus were incorrectly assigned to the suillusfacettalis line by Roszkowska et al. 24,41 . The original description of P. yunnanensis does not specify the shape of cephalic papillae 42 , but the development of evident striae clearly indicates its affinity within Meridioniscus subgen. nov.
The estimation of the number of Pseudechiniscus species to be described by the year 2050 was performed using the exponential model, by applying the best-fit curve to data on cumulative species richness. Curve was described by the formula y = 0.5791e 0.2491x , where x signified the subsequent decades since the description of P. (P.) suillus in 1853 until now. The diagram showing the predicted increment in the number of Pseudechiniscus species is shown in Fig. 3.
Morphological evolution and the evolution of reproductive mode. Firstly, we assessed the phylogenetic signal using 1000 most credible BEAST trees, taking four key taxonomic traits and reproductive mode into consideration: (1) the shape of cephalic papillae, (2) the type of ventral ornamentation pattern, (3) the development of projections on the posterior margin of the pseudosegmental plate IV' , (4) the development of striae and (5) the presence of males in a population. Of these, p-value for Pagel's λ was > > 0.05 for traits 3. and 5., thus they were discarded from further analyses. Subsequently, BayesTraits implemented in RASP was used to reconstruct ancestral states along phylogenetic lineages of Pseudechiniscus. The analysis showed that the ancestor of Pseudechiniscus had elongated (dactyloid) cephalic papillae, therefore the pseudohemispherical papillae of the nominal subgenus should be considered autapomorphic (

Discussion
Systematic and morphological studies with the application of phylogenetic and statistical methods on meiofaunal phyla, represented e.g. by tardigrades 43 , are notoriously laborious and challenging. The finding of abundant populations of Pseudechiniscus allowed for the confirmation of hypotheses 21,22 that its two evolutionary lineages correspond with some morphological traits and can be concurrently elevated to the rank of subgenera. Among the examined traits, the shape of cephalic papillae was recovered as the most conservative character, followed by the development of striae, and the most labile traits, namely the type of ventral ornamentation (Fig. 4) 47,48 , and Echiniscoididae have strongly modified papillae that are comparatively large with respect to the size of the head and are hemispherical in shape [49][50][51] . Not excluding the homologous character of the shape of the cephalic papillae in some Echiniscidae and Echiniscoididae, we are inclined to hypothesise that hemispherical and pseudohemispherical papillae are derived states, present only in some lineages of the clade Echiniscoidea. As noted by Tumanov 22 , the presence of striae in the dorsal sculpturing correlates with phylogeny. Striae are a derived state present in Meridioniscus subgen. nov. and only rarely occur in the subgenus Pseudechiniscus (Fig. 4C). However, some of the examined species, for which multipopulation data are available, exhibit intraspecific variability in the development of striae (e.g. P. (P.) sp. 9 and P. (M.) cf. angelusalas, see Supplementary Material 2). Their homoplasious convergent nature is supported by the independent development in several not directly related lineages within Echiniscidae, e.g. in Cornechiniscus 44 , Stellariscus 52 , or some Echiniscus species, such as E. rackae 53 . Importantly, the ancestor of the genus Pseudechiniscus did not exhibit striae, thus its dorsal sculpturing probably resembled that in Hypechiniscus 25 . Such a type of sculpturing, composed solely of endocuticular pillars, is plesiomorphic for Echiniscidae in general. Consequently, we recognise the development of pores or a sponge layer as advanced characters that evolved convergently in the Bryodelphax and Echiniscus lineages 44 . www.nature.com/scientificreports/ Ventral ornamentation pattern only recently was found to be pivotal in the taxonomy of Pseudechiniscus 22-24,41 , albeit a reticulated design was reported in many species long before these analyses 54 . Although taxonomically informative, this trait is not easily identifiable and requires an optical equipment of high quality and properly fixed specimens to comprehensively describe the belts of pillars and connective points. Moreover, inter-population variability was observed in some widely distributed species, such as Pseudechiniscus (M.) cf. angelusalas, P. (M.) sp. 4 or P. (P.) cf. ehrenbergi, evidenced particularly in the intensity of sculpturing in the intersegmental areas between the subsequent pairs of limbs. Moreover, our populations of P. (P.) suillus exhibited reduced ventral sculpturing, whereas the redescription provides an intricate pattern in this species 23 . The absence of ventral ornamentation is a derived character in the genus Pseudechiniscus, so far observed only in P. (M.) quadrilobatus (replaced by a system of epicuticular thickenings known in Cornechiniscus 44 ). This species has strongly sclerotised armour, and the examination of the type material of P. alberti, also having unusually for the genus sclerotised plates 55 , would be desirable to adjudicate whether the reduction of ventral sculpturing is correlated with the strengthening of the dorsal armour.
Projections, in the form of teeth or lobes, occur on the posterior margin of pseudosegmental plate IV' in both Pseudechiniscus subgenera. The variability in the development of the projections varies between species, as some  (Fig. 4), indicating spurless claws as autapomorphic. Noteworthy, the Javanese P. (P.) bidenticulatus also has spurless claws 59 , suggesting its affinity to this clade and a potential transfer to Meridioniscus. Internal claws with primary spurs are plesiomorphic to Echiniscidae, and the loss of spurs evolved several times in different genera, e.g. in Cornechiniscus 44 . Irrespectively of this, Pseudechiniscus exhibits a tendency towards miniaturisation of spurs and their further merging with the claw branch, especially in the Asian lineages 41 .
Our analyses show that Pseudechiniscus biogeography is complex and requires a detailed interpretation. Our initial assumption was in accordance with the 'everything is everywhere, but environment selects' hypothesis 60,61 , i.e. that the majority of species may have wide geographic distributions, but are limited by climate, and if exceeding a single biogeographic region 62 , then they should occur in realms with similar habitats (e.g. Oriental -Afrotropical -Neotropical). Premises fur such zero hypothesis are numerous: (1) we have already demonstrated a pantropical distribution for an echiniscid species 14 , which we expect to be the only pattern of a widespread geographic range in the tropical region (in other words, e.g. a combination Nearctic-Afrotropical should be considered as rather non-parsimonious and unlikely); (2) a wide distribution was shown for Pseudechiniscus (P.) ehrenbergi, which inhabits both the Western (Italy) and the Eastern (Mongolia) Palaearctic 21,24 ; (3) Meridioniscus subgen. nov. was preliminarily identified in our samples from Chile, Colombia, French Guyana, Kenya, Tanzania and India (not sequenced due to the scarcity of individuals; see also an unidentified species from Argentina in Tumanov 22 ), but never in a much richer material from the Holarctic (see below for the case of P. (M.) indistinctus); contrarily, individuals of the subgenus Pseudechiniscus were found worldwide, and particularly frequently in Europe, including Iceland, Germany, Sweden, Poland, Croatia, or Italy (not sequenced for identical reasons as above). These observations implied a certain degree of regionalisation at least for one of the evolutionary lineages of the genus Pseudechiniscus; (4) much better studied micrometazoans, such as mites, are known to be geographically regionalised, with a low proportion of pantropical and cosmopolitan species inhabiting a given region (up to only 10-15% of the entire fauna 63,64 ); similar values were reported for collembolans (20%) 65 , but see the contrasting evidence for gastrotrich 66  have wide geographic ranges (21%). Thus, it seems that regional endemism prevails over pantropical or cosmopolitan ranges of Pseudechiniscus spp. In other words, only two of the analysed species, P. (P.) cf. ehrenbergi, P. (M.) cf. angelusalas, exhibit genetically verified geographic ranges that are in agreement with the 'everything is everywhere, but environment selects' hypothesis (in the case of species with preferences for tropical habitats, pantropical distributions are predicted by the hypothesis).
The analyses indicated the West Palaearctic-Oriental (Indomalayan) origin of Pseudechiniscus, which is not surprising given the geographically biased dataset. It is probable that sequencing Antarctic (post-Gondwanan) 16 and Australian representatives of Pseudechiniscus may affect the inference about the ancestral geographic range. Although Pseudechiniscus (M.) dastychi and P. (M.) titianae are autochthonous for Antarctica and there are some barcodes available for these species 20,24 , they are non-homologous with our DNA sequences, thus we were not able to include them in the reconstructions. However, it is Australia, where the crucial phyletic lineages of the genus are most likely to be found. This is so, because ancient and highly diversified lineages are known within many animal groups in Australia [68][69][70][71] . Unluckily, the Australian tardigrade fauna is almost unknown, although there are at least few Pseudechiniscus species representing both subgenera present on this continent 72 .
The West Palaearctic and Oriental regions were recovered as ancestral for the subgenus Pseudechiniscus and for Meridioniscus subgen. nov., respectively (Figs. 6-7). Moreover, the first subgenus exhibits a cosmopolitan distribution, whereas the second is primarily tropical/subtropical-Gondwanan (although we do not preclude its presence in e.g. the Mediterranean Palaearctic). A puzzling exception to this pattern is the genetically verified presence of P. (M.) indistinctus in the southern part of Norway (western Scandinavian Peninsula) 24 , but the tropical origin of Meridioniscus subgen. nov. does not exclude dispersal to the Holarctic. Placing this species in the generic phylogeny will allow for the clarification of its affinities. Comparatively longer branches in Meridioniscus subgen. nov. may be caused by two factors: (1) undersampling of lineages, or (2) a relatively slower rate of speciation and extinction, accompanied by morphological stasis, suggesting bradytely (arrested evolution) in this evolutionary line 73 . This hypothesis finds support in the main retained plesiomorphy of the subgenus -dactyloid cephalic papillae (Fig. 4A).
We stress that ancestral area reconstructions using BioGeoBEARS (Fig. 7) indicated the Oriental region to be the centre of dispersal of Meridioniscus subgen. nov. lineages (four points of dispersal), which stays in agreement with the wide geographic distributions of P. (M.) quadrilobatus (potential dispersal to Neotropics from tropical Asia) and P. (M.) cf. angelusalas (potential dispersal to the Afrotropic from tropical Asia). Among the sequenced 25 species, only one representative of the subgenus Pseudechiniscus exhibited a similarly wide tropical distribution (P. (P.) cf. ehrenbergi). This means that only ca. 12% of species in our dataset show signs of pantropical ranges, and we anticipate this fraction to be similar in many other cosmopolitan tardigrade genera (based also on the echiniscid-rich material from South Africa, data in preparation). The analyses predicted that P. (M.) cf. angelusalas and P. (P.) cf. ehrenbergi are likely to be pantropical (which resembles the analyses conducted on Echiniscus lineatus 14 ), with the majority of suitable areas present in the subtropical and tropical zones, and highly suitable areas in Western Europe for P. (P.) cf. ehrenbergi. This species most likely inhabits also the Palaearctic, as indicated by COI data, placing our sequences with those published previously 21,24 in one clade (Supplementary Material 2). We also stress the evident importance of mountain ranges throughout the modelled geographic www.nature.com/scientificreports/ ranges (Fig. 8) for the presence of both taxa, which is in line with the known pattern of tardigrade diversity and endemism increasing to some extent in relation to the elevation 74 (also observed in many other animal groups 75 ).
Pseudechiniscus is the second most speciose echiniscid genus after Echiniscus 23 , and its taxonomy was for long considered disorganised and distrustful. For example, a considerable fraction of literature reports of Pseudechiniscus constitute the records of P. (P.) suillus and P. (P.) facettalis 11 , two allegedly cosmopolitan species. Paradoxically, it is likely that they both have restricted geographic ranges, as their loci typici are at high elevations in the mountains of West Palaearctic 23,32 and in Greenland 76 , i.e. places that favour typical stenotherms rather than eurytopic species. Our model (Fig. 3) predicts that by the mid of twentieth century, the number of described Pseudechiniscus species will be likely at least doubled, approaching the currently known number of Echiniscus species. In fact, it would not be surprising if it is actually Pseudechiniscus that will have the status of the most speciose echiniscid, but testing this supposition requires strict taxonomic regime, that is redescriptions of "old" taxa, e.g. P. (P.) facettalis or P. (M.) conifer, joined with descriptions of new species always associated with DNA barcodes. This will help to avoid synonymies and will also allow for supplementation of generic phylogenetic tree, as well as for more detailed biogeographic analyses and testing of the hypotheses put forward in this study. Beside of the abovementioned Australia, we foresee that sampling in the Oriental region may bring further improvement to our understanding of the natural history of Pseudechiniscus. Recently, Borneo and Indochina, i.e. areas that arose at the border between Laurasia and Gondwana, were argued as crucial for Oriental biodiversity 77 . This goes in tandem with what we discovered for Pseudechiniscus, thus tardigradologists interested in taxonomy and phylogeny of the genus are likely to benefit from focusing on the Indopacific area.

Material and methods
Sampling, microscopy and imaging. Pseudechiniscus specimens were extracted from 65 samples collected throughout the world (Supplementary Table 1). Each population was divided into at least two groups destined for DNA sequencing and light microscopy analyses. If a population was sufficiently abundant, specimens were preserved also for scanning electron microscopy. Permanent slides were made from some individuals using Hoyer's medium and later examined under an Olympus BX53 light microscope with phase contrast (PCM), associated with an Olympus DP74 digital camera. Supplementary figures illustrating dorsal/ventral sculpturing of females and males (provided always in this order) were assembled in Corel Photo-Paint X6, ver. 16.4.1.1281. For structures that could not be satisfactorily focused in a single light microscope photograph, a stack of 2-6 images was taken with an equidistance of ca. 0.2 mm and assembled manually into a single deep-focus image.
Cephalic papillae (secondary clavae) morphotypes. Although all species in the genus Pseudechiniscus have elongated cephalic papillae, the subgenera Pseudechiniscus and Meridioniscus subgen. nov. differ in the orientation and attachment of the papillae to the cuticle of the body, which results in their very much different appearance under LCM (Fig. 9) 22 . Specifically, papillae in Meridioniscus, as in the great majority of Echiniscidae, are attached to the body only by their base and they appear finger-like under both SEM and LCM (Fig. 9A-D). However, in the subgenus Pseudechiniscus, the papillae are attached to the head also along their side (Figs. 9E-H). This peculiar longitudinal attachment is clearly visible in SEM (Fig. 9F,H), but most often makes these papillae appear in LCM as hemispherical (dome-shaped) 21 (Fig. 9E,G). Thus, under LCM, cephalic papillae in the subgenus Pseudechiniscus resemble those in Cornechiniscus (Fig. 9I-J) and Mopsechiniscus. However, the papillae in the two latter genera are flattened and widened from the apex, thus they are indeed hemispherical (dome-shaped) and are visible as such in both SEM and LCM. Thus, in order to avoid confusion when describing the shape of cephalic papillae in the subgenus Pseudechiniscus, and to differentiate them from the truly hemispherical (dome-shaped) papillae in Cornechiniscus (Fig. 9I-J) and Mopsechiniscus 45 , we propose to term them as pseudohemispherical.
Comparative material. We borrowed and examined the following type specimens of Pseudechiniscus species: P. bartkei, P. gullii, P. insolitus, P. jubatus, P. nataliae, P. quadrilobatus, P. ramazzottii, P. santomensis, P. spinerectus. They were used to determine their subgeneric affinity and to formulate taxonomic notes (Supplementary Material 2). www.nature.com/scientificreports/ Estimating species richness. In order to predict species richness in the genus, described Pseudechiniscus species were categorised into cumulative species richness increasing since the description of P. suillus in 1853 32 . Then, a logarithmic best-fit curve was assigned to data to infer a probable progress in the taxonomic descriptions of new taxa. We decided to use this simple method in order to obtain conservative estimates of potential species richness 78 . Genotyping. Prior to DNA extraction, all specimens were examined in vivo under PCM. Individual DNA extractions were made from animals and/or exuviae following recent protocols 79,80 . Three nuclear DNA markers were sequenced: 18S rRNA, 28S rRNA and ITS-1. All fragments were amplified using the primers listed in Supplementary Table 3. We attempted to sequence the mitochondrial COI, but obtained good quality chromatograms only for ca. one third of the analysed populations (see Supplementary Table 3 for all used primer combinations). Nevertheless, GenBank accession numbers for obtained COI sequences are presented alongside the nuclear markers in Supplementary  92 . In MrBayes, random starting trees were used and the analysis was run for ten million generations, sampling the Markov chain every 1000 generations. An average standard deviation of split frequencies of < 0.01 was used as a guide to ensure the two independent analyses had converged. The program Tracer v.1.6 93 was then used to ensure Markov chains had reached stationarity and to determine the correct "burn-in" for the analysis, which was the first 10% of generations. The effective sample size values were greater than 200, and a consensus tree was obtained after summarising the resulting topologies and discarding the "burn-in". The original matrix was also analysed using BEAST in order to obtain a set of Bayesian phylogenetic trees needed for biogeographic analyses. Four clock and tree prior combinations were chosen and ran in parallel: (a) random local clock 94 with the coalescent tree prior, (b) random local clock with the speciation: Yule process as the tree prior, (c) strict clock 95  www.nature.com/scientificreports/ the use of Maxent algorithm, ver. 3.4.1 103 , with the kuenm package 104 , in R 105 . To eliminate a potential bias of clustered occurrences, the datasets were filtered so that there was only one record per a 5 arc-min cell for each species. Thus, only 15 occurrence records for P. (M.) cf. angelusalas and 12 records for P. (P.) cf. ehrenbergi were used in modelling (see Supplementary Material 5 for the full list of records). We used the bioclimatic variables available in WorldClim version 2.1 106 [http:// www. world clim. com/ versi on2], with a 5-arc-minute resolution, as environmental variables for Maxent modelling. Out of 19 available variables, we excluded those that combined temperature and precipitation (bio8, bio9, bio18 and bio19), because they displayed artificial discontinuities between adjacent grid cells in some areas, which could introduce artefacts to modelling 107 . From the remaining 15 bioclimatic variables, we removed the highly correlated ones (Spearman rank correlation value <|0.85|) and selected those with the highest importance based on jack-knife procedure (regularised training gain) for each species (Supplementary Tables 6-7). Finally, we used six variables in the ENM: bio2 (mean diurnal range), bio4 (temperature seasonality), bio10 (mean temperature of the warmest quarter), bio12 (annual precipitation), bio14 (precipitation of the driest month) and bio15 (precipitation seasonality). Based on these variables and randomly selected 60% of occurrence records, we created 255 candidate models for each species combining 17 values of regularisation multiplier (0.1-1.0 at intervals 0.1, 2-6, 8 and 10) and 15 combinations of four feature classes (linear = l, quadratic = q, product = p, and hinge = h). Then, all candidate models were evaluated based on the partial ROC approach 108 and predictive power of model based on omission rates 109 , using the test data subset (40% of species records not included in training data). Statistically significant models with omission rates ≤ 10% were selected as the best models. The 10 best models for P. (M.) cf. angelusalas and 4 for the P. (P.) cf. ehrenbergi were selected according to these criteria (Supplementary Table 8). For each parameter setting selected as the best, we created 10 bootstrap replicates of models with complemental log-log (cloglog) output format 103 . Models were calibrated and projected using the whole world as a training area and all runs were set with 500 iterations and 10,000 background points. Final maps of potential distribution of both species were created in QGIS 110 by averaging all best models. The map showing similarities in the distribution of both species is provided in Supplementary Material 9.

Data availability
Data generated or analysed during this study are included in the published article and its supplementary information files. All unique sequences are deposited in GenBank (i.e. if several specimens from one population shared an identical haplotype, only one sequence was uploaded). www.nature.com/scientificreports/