Material was collected during various field trips in the period 2014–2020. Collected specimens were conserved in absolute ethanol and kept at -20 C to preserve DNA. A total of 11 taxa from the P. ampliatus group sensu Ullrich et al. (2010) were sampled (Supplementary file 1: Locations and accession numbers).

Total DNA was isolated from the hind femur applying a salt-ethanol extraction protocol (e.g. Aljanabi and Martinez 1997). Mitochondrially encoded NADH dehydrogenase subunit 2 (ND2) was amplified through PCR. Sequences from the nuclear internal transcribed spacers (ITS1, 2), published by Ullrich et al. (2010), were accessed from GenBank. New sequences of ITS1 and ITS2 of P. pechevi were obtained in this study to completely represent the P. ampliatus group sensu Heller and Lehmann 2004). PCR was performed in a 25 µL volume using Thermo Scientific DreamTaq Hot Start Master Mix, following the instructions of the manufacturer. ND2 was amplified using the primers TM-J210 AATTAAGCTAATGGGTTCATACCC and TW-N 1284 AYAGCTTTGAARGYTATTAGTTT (Simon et al. 2006) with the thermal cycling following Chobanov et al. (2017). The ITS1+5.8S rRNa+ITS2 fragment was amplified with primers 18S–28S: TAGAGGAAGTAAAAGTCG and 28S–18S: GCTTAAATTCAGCGG (Weekers et al. 2001), applying the protocol by Ullrich et al. (2010) to ensure consistency of the results.

Visualization, trimming and assembly of sequences were performed with CodonCode Aligner v. 8.0.2 (CodonCode, Dedham, MA, USA). Alignments were prepared in Mega X (Kumar et al. 2018). The ITS dataset included mostly sequences from Ullrich et al. (2010). The 5.8S rDNA fragment was removed from the final dataset, resulting in an ITS1+ITS2 matrix. The protein-coding ND2 dataset was checked for stop codons and all datasets were tested for sequence saturation with DAMBE (Xia et al. 2003; Xia 2018). Best substitution models for all datasets defined using PartitionFinder ver. 2.1.1 (Lanfear et al. 2017), dividing the datasets of the protein-coding genes by codon positions, were implemented in the subsequent analyses.

Maximum likelihood (ML) analyses were run using RAxML ver. 8.2.12 (Stamatakis 2014) on the CIPRES Science Gateway (Miller et al. 2010). Node support was obtained through 1000 bootstrap resampling. Bayesian inference (BI) analyses were accomplished in Mr. Bayes v. 3.2.5 (Ronquist et al. 2012). For each dataset four simulations of Markov chains and 2 x 10⁶ generations were run, sampling 1 of every 100 trees. The MCMC parameters were examined in Tracer ver. 1.7.1 (Rambaut et al. 2018). The first 25% of trees were discarded as burn-in.

Choosing calibrations for molecular dating could be challenging, especially when fossil data is unavailable. A recent study estimated divergence dates in Poecilimon, using the isolation of Crete as a source of age constraint (Borissov et al. 2020). The split of the easternmost lineage of the only Cretan species P. cretensis Werner, 1903 was estimated at 0.8 Ma, possibly reflecting its allopatric evolution due to former disconnection of the easternmost part of Crete. Here we use this dating as a secondary calibration to date recent divergence times in the P. ampliatus species group. P. cretensis was included in the analyses based on ND2 and the age of the eastern lineage (Kotsounari) was constrained at 0.8 Ma (SD=0.2). Normal prior distribution was set to minimize error accumulation (Forest 2009; Hipsley and Müller 2014).

BEAST 2.6.3 (Bouckaert et al. 2019) was run on the CIPRES science gateway (Miller et al. 2010) with a Yule process and a total of 2 × 10⁸ generations, sampling every 1000 trees. ESS was monitored in Tracer v. 1.7.1 (Rambaut et al. 2018). Maximum clade credibility tree was constructed using TreeAnnotator (Drummond and Rambaut 2007) and visualized with FigTree (http://tree.bio.ed.ac.uk/software/figtree).

The two sister species Poecilimon ampliatus and P. intermedius diverged recently (Lehmann et al. 2010; Ullrich et al. 2010), although they strikingly differ in distributional area. As we hypothesize that the lineage split was triggered by a rapid climate transformation, we aim to model the environmental requirements of the two species (Soberón 2007; Soberón and Nakamura 2009; Junker et al. 2019) and compare their response to abiotic variables during the Last Glacial Maximum (LGM).

The area accessible via dispersal (M) for a particular taxon is an important component of niche modeling and needs to be carefully estimated, considering landscape changes over time, evolutionary history and relations with other species (Barve et al. 2011). The westernmost reported occurrence of P. intermedius is in Austria (Panrok 2010; Lechner and Zuna-Kratky 2017), while the easternmost record was recently reported from western China (Wu and Liu 2019). We assume that some areas with similar conditions beyond these points could have also been available for colonization. As a northern border of the study area, we chose the 60th parallel north since no literature record beyond this line has been found. The southern border of the study area was extended to include the Balkans and Anatolia covering any possible glacial refugia in these areas. Thus, for calibration of the model of P. intermedius most of Eurasia was included, reaching from the Atlantic coast to Lake Baikal (35ºN–60ºN; 5ºW–35ºE).

The sexually reproducing P. ampliatus is mainly distributed in the western Balkans reaching in an isolated locality Romania to the north (Heller and Lehmann 2004; Iorgu et al. 2008). The Apennine Peninsula was included since it could have been accessible to the species during the Late Pleistocene, when sea levels were significantly lower (Correggiari et al. 1996). The study area for P. ampliatus covered the Apennine and Balkan Peninsula, delineated by the Carpathians to the north and the Aegean Sea to the south (40ºN–47ºN; 9ºE–43ºE).

Occurrence data was collected from literature (Bey-Bienko 1954; Childebaev and Storozhenko 2004; Heller and Lehmann 2004; Kenyeres and Bauer 2008; Panrok 2010; Holuša et al. 2013; Sergeev et al. 2018; Szövényi et al. 2018; Wu and Liu 2019) and the Orthoptera Species File (OSF) (Cigliano et al. 2021). Two datasets were downloaded from the Global Biodiversity Information Facility (GBIF): https://doi.org/10.15468/dl.xmtv9b for Poecilimon intermedius and https://doi.org/10.15468/dl.6ngfu6 for Poecilimon ampliatus to compile the literature record with new localities. Some records found in literature were collected more than a century ago and their precise locality is uncertain or they have not been confirmed recently, although these records help modeling the realized niche of the species. Imprecise localities from literature were geographically referred to the closest and least modified steppe habitat based on similarities with observed habitats of the species using Google Earth (www.google.com/earth).

The occurrence datasets were thinned with the R package spThin (Aiello-Lammens et al. 2015) to ensure that the minimum distance between every two points is 5 km, thus avoiding spatial clustering and matching the resolution of the environmental variables (Sillero and Barbosa 2021).

Layers with the 19 bioclimatic variables available at Worldclim v. 1.4 (Hjimans et al. 2005; www.worldclim.org) were downloaded at 2.5 minutes resolution, together with the downscaled data from two alternative General Circulation Models (GCM) for the LGM: The Community Climate System Model (CCSM4; Gent et al. 2011) and Model for Interdisciplinary Research On Climate Earth System Model (MIROC-ESM; Watanabe et al. 2011). Four variables, namely mean diurnal range (BIO2) , isothermality (BIO3), precipitation of driest month (BIO14) and precipitation seasonality (BIO15), that were reported to show high variability between different GCMs (Varela et al. 2015), were removed to decrease uncertainty of the LGM projections.

The presence-only algorithm implemented in Maxent v. 3.4.3 (Phillips et al. 2017) was used for modeling the current conditions and projecting to the LGM. To select among different levels of model complexity we estimated AICc of 48 combinations of feature classes and regularization multiplier values with the package ENMeval (Muscarella et al. 2014). These included combinations of all feature classes used by Maxent and regularization multiplier values ranging from 0.5 to 4 with an increment of 0.5. The combination with the lowest AICc value (Warren and Seifert 2011) was used to adjust the Maxent settings. To account for the very large study area of P. intermedius, we increased the number of random background points to 20000 to better represent the environment. For P. ampliatus we kept the default Maxent setting of 10000 random background points. For Poecilimon intermedius the lowest AICc was obtained using linear, quadratic and hinge features with a regularization multiplier of 1. For P. ampliatus the model with the lowest AICc included hinge features only and a regularization multiplier of 2.5. Cross-validation with 10 replicates was run for both species. The average value for each cell across the 10 replicates was used to visualize the results. Models were evaluated using area under the receiver operating curve (AUC).

Figure 1. 

European representatives of the Poecilimon ampliatus species-group A,B: P. ampliatus (male, female), Montenegro, Durmitor National Park, 1500 m; C,D: P. pechevi (male, female), Bulgaria/North Macedonia border, Vlakhina Mt., 1900 m; E,F: P. ebneri (male, female), Bulgaria, Belassitsa Mt., 1850 m; G: P. klisuriensis (male), North Macedonia, Pelister Mt., Gjavato Pass, 1100 m; H: P. intermedius (female), Russia, Saratov, 300 m; I,J: P. marmaraensis (male, female), Bulgaria, E Stara Planina Mts, 950 m.

Figure 2. 

Map of the sampled taxa of the Poecilimon ampliatus species group.

The ITS1+ITS2 matrix (ITS hereafter) contained 27 sequences. The total length of the alignment was 708 bp, including 170 variable sites, 88 parsimony-informative sites and 12 gaps. The ND2 dataset included 29 DNA sequences and consisted of 936 bp from which 481 variable and 391 parsimony-informative sites. No signs of significant saturation or NUMTs were found.

Nuclear data. ML analyzes provided poor bootstrap support, though generally agreed with the BI topology (results not shown). The ITS BI-tree (Fig. 3A) reaffirmed the monophyly of the P. ampliatus group sensu lato as suggested by Ullrich et al. (2010). Poecilimon amissus+P. marmaraensis form a strongly supported clade, basal for P. ampliatus group s.l. The tree did not resolve the relationships between P. doga, P. excisus+P. luschani species group, and a clade consisting of the following species: P. ampliatus, P. armeniacus, P. davisi, P. ebneri, P. haydari, P. intermedius, P. klisuriensis, P. pechevi. From the latter clade P. armeniacus+P. haydari forms the basal branch. The following topology was highly supported for the rest: P. davisi+(P. klisuriensis+(P. pechevi+(P. intermedius+P. ampliatus)+P. ebneri)). Although the tree did not resolve some of the nodes closest to the tip, populations of P. ebneri, previously described as P. klisuriensis (Willemse 1982), are clearly separated from the other related taxa (Fig. 3A). BI analyses based on ITS suggested that P. intermedius is most closely related to P. ampliatus.

Figure 3. 

Bayesian inference phylogenetic trees of the Poecilimon ampliatus species group. A: ITS; B: ND2. Posterior probabilities are represented as numbers next to nodes. Colored branches correspond to the Poecilimon ampliatus species complex (yellow), the rest of P. ampliatus group (blue) and outgroups (black), respectively.

Mitochondrial data. The ND2 BI-tree showed the best node support (Fig. 3B) confirming the P. ampliatus group sensu Ullrich et al. (2010) and our ITS phylogeny, with P. amissus+P. marmaraensis representing its basal clade. P. davisi was arranged within the P. ampliatus group. The ND2 phylogeny fully supported the following monophyletic clade at the tip of the tree: (((P. ebneri+P. klisuriensis)+(P. intermedius+P. ampliatus+P. pechevi))+(P. orbelicus+P. haydari)). P. karakushi Ünal, 2003 (P. syriacus group according to Heller et al. 2008; later transferred by Ünal 2010 to the P. minutus group but see the position of P. doga in Fig. 3A as well as the trees by Ullrich et al. 2010), formed a sister clade to the P. ampliatus lineage. P. glandifer arranged outside the P. ampliatus species group distantly grouping with P. brunneri (Frivaldzsky, 1868).

According to our divergence times estimation, all splits within P. ampliatus s.l. happened during the Pliocene and Pleistocene (Fig. 4) (highest posterior density intervals are given in Supplementary material 2: List of 95 % HPD intervals). The P. ampliatus s.l. group shared a common ancestor 4.09 Ma when P. amissus+P. marmaraensis diverged. Divergence of P. karakushi was dated much earlier, 5.62 Ma. The time of most recent common ancestor (TMRCA) of P. luschani group (represented by P. orbelicus) and P. armeniacus group (represented by P. haydari) was estimated at 2.22 Ma. TMRCA of the P. ampliatus complex (P. ampliatus, P. ebneri, P. intermedius, P. klisuriensis, P. pechevi) is estimated as 2.32 Ma. These taxa shared most recent common ancestor with P. luschani+P. armeniacus 2.96 Ma. The divergence between P. klisuriensis and P. ebneri was dated at 0.45 Ma. The split between P. intermedius+P. ampliatus and P. pechevi was estimated at 0.43 Ma. P. intermedius and P. ampliatus diverge 0.37 Ma.

Figure 4. 

Mitochondrial chronogram based on the ND2 phylogenetic reconstruction of the Poecilimon ampliatus species group. Blue line indicates climate fluctuations over time (by Robert A. Rohde, based on data from Petit et al. 1999; Lisiecki and Raymo 2005). Blue arrow indicates the switch from 41 to 100 kyr glacials. Red arrow indicates the Mid-Brunhes Transition. Black numbers before the nodes correspond to time estimation in millions of years. Numbers in white circles at the branches refer to Supplementary material 2, where they correspond to the 95% HPD intervals.

After the thinning procedure a total of 58 localities for P. intermedius (Fig. 5A) and 32 for P. ampliatus (Fig. 6A) were retained. The AICc and AUC values for each combination of model parameters are given in Supplementary material 3: Model parameters. The combinations with the lowest AICc were used to set the parameters of Maxent and for replication.

Figure 5. 

Ecological niche model of P. intermedius in geographic space. A: occurrence points; B: current conditions; C: LGM projection on CCSM4; D: LGM projection on MIROC-ESM. Scale bar represents 1000 km.

Figure 6. 

Ecological niche model of Poecilimon ampliatus in geographic space. A: occurrence points; B: current conditions; C: LGM projection on CCSM4; D: LGM projection on MIROC-ESM. Scale bar represents 250 km.

The variables with the highest contribution to the model of P. intermedius were precipitation of warmest quarter (BIO18) and precipitation of driest quarter (BIO17). The highest gain from a single variable was obtained from BIO17 and maximum temperature of warmest month (BIO5) (see jackknife test Fig. 7C). The model highlighted suitable conditions for P. intermedius in Central Europe and Asia, generally beyond the 45^th parallel north (Fig. 5B). The projections on the LGM conditions showed a significant difference between CCSM4 and MIROC-ESM (Fig. 5C, D). The projection on the CCSM4 selected three sharply delineated regions with suitable conditions for the species during the LGM: 1) a small area between the Balkan and the Apennine Peninsula, including parts of the Adriatic Sea that were subaerial at the time; 2) an area northwest of the Black Sea; 3) an area north of the Caspian Sea. The projection on the MIROC-ESM model (Fig. 5D) outlined broader suitable areas, covering almost the whole Apennine Peninsula, a significant part of the southeastern Balkans and parts of northwestern Anatolia and the northern Caucasus but with low probability of presence.

Figure 7. 

Importance of environmental variables. A: response of Poecilimon intermedius to precipitation of driest quarter (BIO17); B: Response of P. ampliatus to BIO17; C: Jackknife test of regularized training gain of the P. intermedius model; D: Jackknife test of regularized training gain of the P. ampliatus model.

The most contributing variable in the P. ampliatus model was precipitation of driest quarter (BIO17) with 65 % relative contribution to the model. BIO17 also obtained the highest gain when used in isolation (Fig. 7D). Among the others, maximum temperature of warmest month (BIO5) showed a relative contribution of 16% to the model. The most suitable conditions for the species were found in the mountain ranges of the Western Balkans (Fig. 6B). The two LGM simulations showed significant disagreement when compared. The projection on CCSM4 showed a significant increase in the suitable area during the LGM compared to the current conditions including a large area on the Apennine Peninsula (Fig. 6C). MIROC-ESM highlighted much smaller areas with suitable conditions during the LGM, generally resulting in a much lower probability of presence (Fig. 6D).

The phylogenetic reconstructions presented in this study support the P. ampliatus species group in a broad sense as defined by Ullrich et al. (2010). This clade, hereafter referred to as P. ampliatus species group, includes most representatives of the following four morphological species groups, currently listed in Cigliano et al. (2021): P. ampliatus (sensu Heller and Lehmann 2004 and partly sensu Ünal 2010), P. armeniacus, P. davisi, and P. luschani (sensu Ünal 2010 and Boztepe et al. 2013). In addition, P. doga, included within the P. minutus group (Ünal 2010), clearly groups here. As at least two members of the latter show either morphological (P. minutus Karabag, 1975) or both morphological and genetic (P. karakushi Ünal, 2003) affinities with the P. syriacus group (Heller et al. 2008; Ullrich et al. 2010; Chobanov et al. 2020; this paper), the P. minutus group is polyphyletic. The remaining taxon within the latter, P. solus Ünal, 2010, has unclear relationships. The species groups P. ampliatus and P. syriacus are here supported by the deep split between the former and P. karakushi, and thus their early divergence (Fig. 3A, B and Fig. 4).

The ITS phylogeny suggests that the P. armeniacus species group after Ünal (2010) is polyphyletic and P. excisus is related to the P. luschani group sensu Boztepe et al. (2013). The latter is considered here as the P. orbelicus species complex. The monophyly of the P. orbelicus species complex is strongly supported (Fig. 3A). At least a few taxa closely related to P. armeniacus (P. guichardi, P. harveyi, P. haydari) fit well within the P. ampliatus species group (see Ullrich et al. 2010 and Chobanov et al. 2020) and thus, though polyphyletic, the composition of P. armeniacus species group fits the P. ampliatus group s.l. with a few taxa requiring additional study.

As a result, the Poecilimon ampliatus groupings sensu Heller and Lehmann (2004) and Ünal (2010) represent polyphyletic lineages with P. marmaraensis+P. amissus being the basal clade of the group while P. ataturki and P. glandifer are not directly related to that group (Fig. 3A and comments in Chobanov et al. 2020). In addition, strong support was obtained for a clade at the tip of the tree that comprises five taxa, namely P. ampliatus, P. ebneri, P. intermedius, P. klisuriensis and P. pechevi. This cluster will be referred to as the P. ampliatus species complex hereafter. Genetic distinction of P. klisuriensis from P. ebneri was supported by both the nuclear and mitochondrial phylogenies with the ITS-phylogeny suggesting its basal position within the P. ampliatus species complex. Thus, P. klisuriensis may be considered a distinct species as suggested by the stable cercus shape and ecological preferences of all sampled populations (Willemse 1982 and own observations). Though it occurs sympatrically with P. ebneri, the populations of both taxa were never found together with P. klisuriensis occurring in the low-mountain forest belt appearing early in the season and P. ebneri occupying mountain summits over the tree-line later in the season.

The nuclear phylogeny strongly supports that P. davisi is a sister taxon to the P. ampliatus species complex, which was not supported by the mitochondrial phylogeny (compare Fig. 3A and B).

Тhe Pleistocene climatic oscillations caused well-studied dramatic range shifts and shaped genetic diversity and speciation of numerous lineages on a world scale (Hewitt 1996, 2004; Taberlet et al. 1998; Wallis et al. 2016). This “orbitally forced species’ range dynamics” is reported to select against specialization and for dispersal ability (Dynesius and Jansson 2000). According to our mitochondrial chronogram (Fig. 4) members of the P. ampliatus species complex share a common ancestor 2.3 Ma (Middle Gelasian), and all subsequent splits fall into the Pleistocene. The closest relatives of the P. ampliatus complex occur in Anatolia (P. davisi, P. armeniacus, P. haydari and most members of the P. orbelicus complex) with a single species found on the Balkan Peninsula (P. orbelicus). Broad terrestrial connections between northwestern Anatolia and the Balkans existed during Pliocene and Pleistocene times (Elmas 2003), providing corridors for multiple dispersal events (e.g., Chobanov et al. 2017). Divergence of P. orbelicus, a member of the P. orbelicus complex, distributed in southern Bulgaria and northern Greece, was estimated at ca. 1 Ma (Kaya et al. 2015). The higher diversity and the larger distributional area of the P. ampliatus complex suggests that its ancestor reached the Balkans earlier, probably during the late Pliocene.

Two splits at the tip of the tree:1) between P. ebneri and P. klisuriensis and 2) between P. pechevi and P. intermedius+P. ampliatus, were estimated at 0.40–0.45 Ma (also matching the 95 % HPD intervals). These splits coincide with the Mid-Brunhes Transition, an increase in the amplitude of climatic cycles ca. 0.43 Ma, that resulted in significantly warmer interglacials (Jansen et al. 1986; Candy et al. 2010; Barth et al. 2018). This rapid transition forced populations, adapted to colder and more humid conditions, to retreat towards higher altitudes, thus isolating mountainous populations (e.g. Berger et al. 2010). Our ecological niche modelling outlined the major role of precipitation in limiting the distribution of P. ampliatus (Fig. 7B. D). Comparison between the current distribution and LGM projections supports a shift towards higher altitudes during interglacial (Fig. 6B) and expansion during glacial periods (Fig. 6C), which we believe is a common model for all mountainous species in the group including P. ebneri and P. pechevi. On the other hand, the habitat of P. klisuriensis, similar to that of P. marmaraensis in Europe—lush forest meadows in the mid-altitude belt, is also defined in terms of high humidity, though the ecological niche of the species is not only spatially but also temporarily restricted as those meadows dry out with the beginning of the summer. Therefore, the populations of the species must have retracted towards the remaining humid habitats during dry glacial periods.

Molecular dating showed that P. intermedius, P. ampliatus and P. pechevi shared a common ancestor 0.43 Ma. The split between P. ampliatus and P. intermedius was estimated at 0.37 Ma (Fig. 4). A detailed study of the parthenogenesis in P. intermedius ruled out hybridization or bacteria-induced parthenogenesis (Lehmann et al. 2011). Sexually reproducing species, that have sister parthenogenetic lineages, tend to have broader niches and distributional ranges than other relatives (van der Kooi et al. 2017). Poecilimon ampliatus has the largest distributional area among the members of the P. ampliatus species complex and its range extends northwards in comparison to all representatives of the group, this serving as additional evidence that the latter might be the closest sexual relative of P. intermedius.

Though obligately thelytokous lineages are not rare in insects, these usually have a recent origin and are not expected to have long evolutionary history (Bell 1982; Normark 2003). Pleistocene climatic shifts are reported to favor very rapid speciation, guided by sexual selection, forming reproductive barriers, without accumulating significant DNA substitutions (Knowles 2000). Our ITS phylogeny did not solve relationships between most lineages in the P. ampliatus complex, which indicates recent simultaneous splits in the group. The latter is supported by weak morphological differentiation and similar song patterns between the members of the P. ampliatus complex (Willemse 1982; Heller and Lehmann 2004; Chobanov et al. 2020).

Decay of sexual traits and mating behavior is expected and has been reported in parthenogenetic females (Carson et al. 1982; Kraaijeveld et al. 2009). Such decay in the auditory mating-response was documented in P. intermedius (Lehmann et al. 2007). This species shows significant vestigialization of auditory structures, compared to other representatives of Poecilimon (Strauß et al. 2014). Apparently, reduction of auditory function could be fast, as a result of rapid specialization, as is the case with P. jablanicensis (see Chobanov and Heller 2010). In contrast, a stick insect from New Zealand, Clitarchus hookeri (White, 1846), has both sexual and asexual populations with males being rare or absent in marginal, previously glaciated areas, where parthenogenetic females yet maintain sexual traits (Nakano et al 2019). One possible explanation for this phenomenon is the recent origin of C. hookeri, less than 40 kya (Morgan-Richards et al. 2019). Thus, in order to lose basic sexual traits, P. intermedius should have existed way longer. However, sexual trait decay in parthenogenetic lineages do not depend on divergence age only and could be rather fast (see Schwander et al. 2013).

Parthenogenetic lineages are known to rapidly colonize large areas that become available after glacial periods (Stenberg et al. 2003; Kearney 2005). One example is the parthenogenetic bush-cricket Saga pedo (Pallas, 1771), which resembles P. intermedius in distribution range, extremely larger than the range of its closest relatives (Kolics et al. 2012), though the former is pentaploid (Dutrillaux et al. 2009) and thus possibly of a hybrid origin. Asexually reproducing organisms have important advantages in colonizing areas with relatively harsh conditions compared to their sexual relatives. For instance, parthenogenetic reproduction is beneficial over mating when the favorable season is very short (Fernandez et al. 2010, 2012), as asexual reproduction saves time and resources that are usually spent for mate searching and courtship performance. Besides, severe climatic conditions could cause a massive decrease in population size, thus restricting mating success and leading to inbreeding in sexual populations, while asexual populations are less affected (Haag and Ebert 2004).

The main characteristic of our ecological niche model of P. intermedius is that the projection on the LGM conditions significantly reduces suitable area (Fig. 5C, D). Projections on the two general circulation models, CCSM4 and MIROC-ESM, did not show coherent results (Figs 5, 6). Different GCMs show low agreement in some continental areas (Varela et. al. 2015). For P. intermedius MIROC-ESM (Fig. 5D) predicts several small disconnected suitable areas on the Apennine Peninsula, and the south-eastern Balkans during the LGM. The Eastern Balkans retained dry conditions until mid-Holocene (Wright et al. 2003), which could have favored the radiation of P. intermedius (see Fig. 7A). Yet, the wide distribution of Artemisia-chenopod semi-deserts even at high altitudes in the Eastern Balkans (Wright et al. 2003), suggests very dry conditions inappropriate for P. intermedius, as suggested by the CCSM4 model (Fig. 5C), excluding occurrence of the species from all the Balkans. In addition, access to this southern area could have been restricted by geomorphological and biotic factors. Two 874 bp DNA sequences of ITS1+5.8S+ITS2 from P. intermedius from two remote areas (Saratov, Russia and Czech Republic) (Supplementary material 1: Localities and accession numbers)—over 2000 km apart, were identical, indicating recent expansion from a small population that had undergone a single loss of sexual reproduction. One of the GCMs – MIROC-ESM predicted an extremely limited refugial zone for P. ampliatus in the Dinarides (Fig. 6D), which is in strong disagreement with the result from CCSM4 (Fig. 6C) and does not correspond to our phylogeographic implications and the current distribution of the species (e.g. isolated populations in Romania and Serbia). Therefore, we further opt for the CCSM4 models.

Precipitation variables were the most informative for the models of both species. For P. ampliatus precipitation of driest quarter (BIO17) was the most important variable. The projection of the P. ampliatus ecological niche model on the CCSM4 demonstrated that during glaciations suitable habitats for the species extended on a large area including lowland and mountain areas over the Dinarides, while populations retreated to higher latitudes in response to warming (compare Fig. 6B and C). This allows us to speculate that P. ampliatus had a significantly wider distribution during the Pleistocene and the populations in Romania and eastern Serbia (Heller and Lehmann 2004) are remnants from this former range. At the same time, favorable habitats of P. intermedius have significantly shrunk to restricted areas north of the Adriatic, Black and Caspian Sea (Fig. 5C). Hence, during the LGM, suitable areas for P. intermedius and P. ampliatus were in close proximity in the Western Balkans. Comparison between the response to precipitation (BIO17) of P. ampliatus and P. intermedius demonstrates that the latter is significantly more tolerant to drought (Fig. 7A, B) and even demanding a certain low amount of precipitation during winter (Fig. 7A). A sharp detachment of the ecological niches of both species based on the levels of precipitation in the driest quarter of the year below and over ca. 200 mm may be exemplified based on their current modelled distribution (Fig. 8).

Figure 8. 

Models of current distribution of Poecilimon intermedius and P. ampliatus mapped in relation to precipitation of driest quarter (BIO17). A: color-coded map of BIO17 (Hijmans et al. 2005); B: territories with suitable conditions for P. intermedius; C: territories with suitable conditions for P. ampliatus; D: overlay of the mapped models of P. intermedius and P. ampliatus in the context of BIO 17. Scale bar represents 1000 km.

All the above speculations call for a comparatively recent climate-driven origin of P. intermedius out of a high humidity-dependent ancestor. The latter might have been an isolated population of P. ampliatus or a common ancestor of both taxa, subjected to isolation during climate and habitat deterioration and population decline. Based on dating the main lineage splits within the P. ampliatus species complex, shrinking of the eco-niche of P. ampliatus during interglacial periods, and poor genetic distinction between distant populations of P. intermedius, we suggest that the current vast population of P. intermedius evolved in a single evolutionary event ca. 0.4 Ma, shortly after the Mid-Brunhes Transition. The warmer interglacials that followed were suitable for P. intermedius to colonize the vast territories to the east. The distance on land between the Dinarides and China is ca. 5000 km which gives an average colonization speed of ca. 0.01 km per year. Colonization was rapid during interglacials followed by retreat to refugia during glacials. Thus, the expansion of the species must have had a pulsate stepping stone character with successive cycles forward and back. Even though considering that the expansion periods were shorter than retraction, a significantly larger speed (>>10 meters/year) combined with dispersal from the refugial stepping stone-populations allowed this vast ‘migration’. According to our ecological niche models P. intermedius suffered a significant population crisis during the Last Glacial Maximum and is suspected to be in a current expansion following the expansion of open habitats, increase of temperature and decrease of humidity.

Poecilimon klisuriensis was initially described as a species, differing from P. ebneri in the shape of the male cercus tip, length of the subgenital plate and the shape of the lower valve of the ovipositor (Willemse 1982). It was later synonymized and regarded as a form of P. ebneri (Heller and Lehmann 2004). However, our phylogenetic analyses (supporting result by Ullrich et al. 2010) show considerable genetic distinction comparable with the genetic distances between well-outlined species within this group (Fig. 3A, B). Poecilimon klisuriensis populations seem to express two stable independent characters – ecological requirements and morphology of cerci and female ovipositor lower valve. Populations showing these morpho-characters were always found at lower altitudes and develop earlier in the season, though they fit within the geographical range of P. ebneri. Yet, specimens representing morphotypes of P. ebneri or klisuriensis do not occur together. As a result, we propose re-establishing the species status of P. klisuriensis.

The present study also supports the species status of P. pechevi that was recently questioned based on its similarity with P. ebneri (Lemonnier-Darcemont and Darcemont 2020).

A general systematic implication of the presented phylogenies is the monophyly of the group as basically proposed by Ullrich et al. 2010. Herewith we confirm the monophyletic origin of the Poecilimon ampliatus species group sensu lato including the following 27 taxa (24 species) ordered into complexes by an approximate phylogenetic order:

Poecilimon amissus species complex

Poecilimon amissus Brunner von Wattenwyl, 1878

Poecilimon marmaraensis marmaraensis Naskrecki, 1991

Poecilimon marmaraensis nalbanti Ünal, 2005

Poecilimon orbelicus species complex (= P. luschani species group sensu Boztepe, Kaya & Çiplak 2013)

Poecilimon egrigozi Ünal, 2005

Poecilimon helleri Boztepe, Kaya & Çiplak, 2013

Poecilimon ledereri Ramme, 1933

Poecilimon luschani birandi Karabag, 1950

Poecilimon luschani chobanovi Boztepe, Kaya & Çiplak, 2013

Poecilimon luschani luschani Ramme, 1933

Poecilimon orbelicus Pančić, 1883

Poecilimon tuncayi Karabag, 1953

Poecilimon armeniacus species complex (partly sensu Ünal, 2010; still in need of revision)

Poecilimon armeniacus (Uvarov, 1921)

Poecilimon eskishehirensis Ünal, 2003

Poecilimon ferwillemsei Ünal, 2010

Poecilimon harveyi Karabag, 1964

Poecilimon haydari Ramme, 1951

Poecilimon guichardi Karabag, 1964

Poecilimon inopinatus Ünal, 2010

Poecilimon karabagi (Ramme, 1942)

Poecilimon ampliatus species complex

Poecilimon ampliatus Brunner von Wattenwyl, 1878;

Poecilimon ebneri Ramme, 1933

Poecilimon intermedius (Fieber, 1853);

Poecilimon klisuriensis Willemse, 1982, stat. rev.;

Poecilimon pechevi Andreeva, 1978

Species with distant relationships (ungrouped):

Poecilimon davisi Karabag, 1953

Poecilimon doga Ünal, 2004;

Poecilimon excisus Karabag, 1950.