Phylogeography of the Poecilimon ampliatus species group (Orthoptera: Tettigoniidae) in the context of the Pleistocene glacial cycles and the origin of the only thely

Parthenogenetic lineages are known to rapidly colonize large areas that become available after glacial periods as parthenogenetic reproduction is beneficial over mating when the favorable season is very short. The only obligatory parthenogenetic species of the largest bush-cricket subfamily Phaneropterinae is Poecilimon intermedius. It belongs to the Anatolio-Balkan lineage Poecilimon ampliatus species group and in contrast has a remarkably broad distribution from Central Europe to China, following the pattern of geographical parthenogenesis. In this study we provide a systematic revision of the P. ampliatus group based on mitochondrial (ND2) and nuclear (ITS) phylogeny. We estimate divergence times by applying secondary calibration on the ND2 tree to test for congruence between recent splits in the group and the Pleistocene climatic oscillations. We use ecological niche modelling to analyze the ecological requirements of the parthenogenetic P. intermedius and its sexually reproducing sister species P. ampliatus. By projecting on the conditions during the Last Glacial Maximum we outline the suitable areas for both species during the glacial cycles and discuss range shifts in response to climate change. Based on all results we hypothesize that the drought-tolerant P. intermedius originated during the recent glaciations in the southwestern part of its current range and rapidly radiated in a northeastern direction. Its sister species P. ampliatus, which is adapted to higher levels of precipitation, remained in the western Balkans, where populations retreated to higher altitudes during warming.

Large distribution areas are reported for various asexual organisms (Law and Crespi 2002;Stenberg et al. 2003;Cosendai et al. 2013). Тhe phenomenon that asexual lineages tend to occupy much bigger territories than their closely related sexual lineages was described as 'geographical parthenogenesis' by Vandel (1928). Explanations of that pattern cover different aspects, such as the balance between abiotic and biotic factors (Glesener and Tilman 1978), or alternative models of origin of asexual lineages (Lynch 1984;Vrijenhoek 1984). Because mating is a time-consuming process, origin of parthenogenesis is also linked to shortened periods of favorable conditions when rapid reproduction is required (Fernandez et al. 2010(Fernandez et al. , 2012. Recent reviews lead towards synthesis of the theories and highlight the importance of geographical parthenogenesis for better understanding of sexuality (Hörandl 2006;Vrijenhoek and Parker 2009;Tilquin and Kokko 2016).
It was argued that the advantage of parthenogenetic lineages in the colonization of new areas, such as formerly glaciated territories, is a consequence of hybridization rather than of parthenogenesis itself (Kearney 2005). However, no evidence of hybridization was found in P. intermedius (Lehmann et al. 2011). It has a diploid chromosome set (Warchalowska-Sliwa et al. 1996) and is a rare example of a thelytokous parthenogenetic bush-cricket, also following the pattern of geographical parthenogenesis (Lehmann et al. 2011). Poecilimon intermedius is considered most closely related to P. ampliatus (Ullrich et al. 2010) and a recent spontaneous origin of the thelytokous parthenogenesis of the former (less than 200 kya) was suggested (Lehmann et al. 2011). The remarkably broad range of the asexual P. intermedius is in contrast with the small restricted ranges of its closest relatives: P ampliatus and P. ebneri, which occur on some mountain ranges on the Balkan Peninsula (Heller and Lehmann 2004;Chobanov et al. 2020).
Our study aims to explore the origin of the only known obligatory parthenogenetic member of the richest bush-cricket subfamily Phaneropterinae (over 2,600 species; Mugleston et al. 2018) through reconstructing the molecular phylogeny of the P. ampliatus group. We test the hypothesis of a recent origin of the parthenogenetic lineage P. intermedius, triggered by the rapidly changing environmental conditions during the Pleistocene (e.g. Raymo 2005, 2007). We estimate divergence times and relate recent lineage splits to Pleistocene climatic events. Additionally, we use ecological niche modeling to outline possible dispersal routes and range shifts of the parthenogenetic P. intermedius and its closely related sexually reproducing ally P. ampliatus, driven by the glacial-interglacial cycles. As a result, we propose an evolutionary scenario for the origin of P. intermedius and redefine the monophyletic P. ampliatus species group in a broad sense.

Sampling, laboratory procedures and dataset preparation
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 dehy-drogenase 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 ITS 2 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 AATTAAGCTA-ATGGGTTCATACCC and TW-N 1284 AYAGCTTT-GAARGYTATTAGTTT (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 (Co-donCode, 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.

Phylogenetic analyses and divergence time estimations
Maximum likelihood (ML) analyses were run using RAx-ML ver. 8.2.12 (Stamatakis 2014) on the CIPRES Science Gateway (Miller et al. 2010 (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).

Ecological niche modelling
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).
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).

Characteristics of the datasets
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.

Phylogeny
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.

Estimation of divergence times
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. intermedi-us+P. ampliatus and P. pechevi was estimated at 0.43 Ma. P. intermedius and P. ampliatus diverge 0.37 Ma.

Ecological niche modelling
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.
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 MI-ROC-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 MI-ROC-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.
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 (BIO 5) 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 P. ampliatus species group
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) (Ü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 andChobanov 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).

Evolution of the P. ampliatus species complex
Т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(Hewitt , 2004Taberlet 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. in-termedius+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 al- titudes 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.

Origin of Poecilimon inter mediusrange shifts and loss of the male in response to cli matic oscil la tions
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(Fernandez et al. , 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).
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.

Systematic inferences: taxonomic reconsiderations and group composition
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:  (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.