Research Article
Research Article
Phylogeography of the Poecilimon ampliatus species group (Orthoptera: Tettigoniidae) in the context of the Pleistocene glacial cycles and the origin of the only thelytokous parthenogenetic phaneropterine bush-cricket
expand article infoSimeon Borislavov Borissov, Georgi Hristov Hristov, Dragan Petrov Chobanov
‡ Institute of Biodiversity and Ecosystem Research, Bulgarian Academy of Sciences, Sofia, Bulgaria
Open Access


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.

Key words

ecological niche modelling, evolution, parthenogenesis, phylogeny, Poecilimon intermedius, systematics

1. Introduction

Poecilimon Fischer, 1853 is the most diverse genus of the Eurasian bush-cricket tribe Barbitistini with more than 140 currently described species (Cigliano et al. 2021). Representatives are flightless, phytophagous and show complex acoustic behavior. The exceptional diversity of the genus is explained by recent radiation and rapid diversification, driven by the Pleistocene cyclic glaciations (La Greca 1999; Kaya et al. 2015). Taxonomic studies define species groups within the genus, based on morpho-acoustic characters (e.g. Heller et al. 2006, 2011; Chobanov and Heller 2010; Kaya et al. 2012; Boztepe et al. 2013; Kaya et al. 2018). Molecular phylogeny supported monophyly of many species groups (Ullrich et al. 2010; Boztepe et al. 2013; Kaya 2018; Borissov et al. 2020), yet some systematic problems remain unsolved (Ullrich et al. 2010; Chobanov et al. 2020).

Diversity of Poecilimon is concentrated in the southern Balkans, Anatolia and the Caucasus with very few taxa spreading outside this range. Nevertheless, the parthenogenetic species P. intermedius (Fieber, 1853) has a remarkably broad range stretching from Central Europe to western China, reaching well over the 55th parallel north in Middle Asia (Bey-Bienko 1954; Kenyeres and Bauer 2008; Wu and Liu 2019). Heller and Lehmann (2004) provided evidence that no males of the latter have been found and defined the P. ampliatus species group including P. intermedius, P. amissus Brunner von Wattenwyl, 1878, P. ampliatus Brunner von Wattenwyl, 1878, P. ebneri Ramme, 1933, and P. marmaraensis Naskrecki, 1991. Additionally, two sibling species – P. glandifer Karabag, 1950, and P. ataturki Ünal, 1999, were included in the P. ampliatus group (Ünal 2010). In addition to the morpho-acoustic grouping, a molecular phylogeny published by Ullrich et al. (2010) suggested that the group in its original composition is polyphyletic unless the following taxa are included – P. davisi Karabag, 1953, P. doga Ünal, 2004, P. excisus Karabag, 1950, P. haydari Ramme, 1951, P. armeniacus (Uvarov, 1921), P. ledereri Ramme, 1933, P. lushani Ramme, 1933, P. orbelicus Pančić, 1883, P. tuncayi Karabag, 1953, P. birandi Karabag, 1950. The last five, together with P. egrigozi Ünal, 2005, and P. helleri Boztepe, Kaya and Çiplak, 2013, were later included in the P. luschani species group (Boztepe et al. 2013). Recently, Chobanov et al. (2020) provided details on the song pattern and morphology of some taxa of Poecilimon discussing the acoustic and morphological belonging of P. pechevi Andreeva, 1978, P. armeniacus, P. harveyi Karabag, 1964, P. guichardi Karabag, 1964, P. haydari, P. doga, and P. davisi to the P. ampliatus species group, while excluding P. glandifer and P. ataturki based on the amplitude-temporal characteristics of the male calling song and the peculiar structure of the stridulatory file.

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, 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. Lisiecki and 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.

2. Methods

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

2.2. Phylogenetic analyses and divergence time estimations

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 106 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 × 108 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 (

2.3. 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): for Poecilimon intermedius and 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 (

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

3. Results

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

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

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

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

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

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 45th 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).

4. Discussion

4.1. 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): 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).

4.2. 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, 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.

4.3. Origin of Poecilimon intermedius — range shifts and loss of the male in response to climatic oscillations

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.

4.4. 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:

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.

5. Acknowledgements

This study was supported by the National Science Fund (MES) of Bulgaria – project DN11/14–18.12.2017 to Dragan Chobanov (sampling, DNA, isolation, amplification and sequencing), and by the Bulgarian Ministry of Education and Science under the National Research Program “Young scientists and postdoctoral students” approved by DCM Nº577 / 17.08.2018 – a grant to Simeon Borissov (data processing and analyses). We thank the reviewers Klaus-Gerhard Heller and Gerlind Lehmann, as well as the editor Marianna Simões, for their valuable corrections, comments and contribution in improving the manuscript.

6. References

  • Aiello-Lammens ME, Boria RA, Radosavljevic A, Vilela B, Anderson RP (2015) spThin: an R package for spatial thinning of species occurrence records for use in ecological niche models. Ecography 38: 541–545.
  • Aljanabi S, Martinez I (1997) Universal and rapid salt-extraction of high quality genomic DNA for PCR- based techniques. Nucleic Acids Research 25: 4692–4693.
  • Barve N, Barve V, Jiménez-Valverde A, Lira-Noriega A, Maher SP, Peterson AT, Soberón J, Villalobos F (2011) The crucial role of the accessible area in ecological niche modeling and species distribution modeling. Ecological Modelling 222: 1810–1819.
  • Berger D, Chobanov DP, Mayer F (2010) Interglacial refugia and range shifts of the alpine grasshopper Stenobothrus cotticus (Orthoptera: Acrididae: Gomphocerinae). Organisms Diversity & Evolution 10: 123–133.
  • Bey-Bienko GY (1954) Phaneropterinae. Fauna of the USSR. Orthoptera 2 (2), 385 pp.
  • Borissov SB, Bobeva A, Çıplak B, Chobanov D (2020) Evolution of Poecilimon jonicus group (Orthoptera: Tettigoniidae): a history linked to the Aegean Neogene paleogeography. Organisms Diversity & Evolution 20: 803–819.
  • Bouckaert R, Vaughan TG, Barido-Sottani J, Duchêne S, Fourment M, Gavryushkina A, Heled J, Jones G, Kühnert D, De Maio N, Matschiner M, Mendes FK, Müller NF, Ogilvie HA, du Plessis L, Popinga A, Rambaut A, Rasmussen D, Siveroni I, Suchard MA, Wu C-H, Xie D, Zhang C, Stadler T, Drummond AJ (2019) BEAST 2.5: An advanced software platform for Bayesian evolutionary analysis. PLOS Computational Biology 15: e1006650.
  • Boztepe Z, Kaya S, Çiplak B (2013) Integrated systematics of the Poecilimon luschani species group (Orthoptera, Tettigoniidae): radiation as a chain of populations in a small heterogeneous area: Poecilimon luschani species group. Zoological Journal of the Linnean Society 169: 43–69.
  • Candy I, Coope GR, Lee JR, Parfitt SA, Preece RC, Rose J, Schreve DC (2010) Pronounced warmth during early Middle Pleistocene interglacials: Investigating the Mid-Brunhes Event in the British terrestrial sequence. Earth-Science Reviews 103: 183–196.
  • Childebaev MK, Storozhenko SY (2004) An annotated list of the long-horned orthopterans (Orthoptera, Ensifera) of Kazakhstan. Tethys Entomological Research 9: 213–228.
  • Chobanov DP, Heller K-G (2010) Revision of the Poecilimon ornatus group (Orthoptera: Phaneropteridae) with particular reference to the taxa in Bulgaria and Macedonia. European Journal of Entomology 107: 647–672.
  • Chobanov DP, Kaya S, Grzywacz B, Warchałowska-Śliwa E, Çıplak B (2017) The Anatolio-Balkan phylogeographic fault: a snapshot from the genus Isophya (Orthoptera, Tettigoniidae). Zoologica Scripta 46: 165–179.
  • Chobanov DP, Sevgili H, Heller K-G (2020) Bioacoustics of poorly known Poecilimon taxa (Insecta: Orthoptera: Tettigoniidae) with redescriptions of P. pechevi and P. stschelkanovzevi. Zootaxa 4890: 535–553.
  • Correggiari M, Trincardi A, Roveriand F (1996) Late Pleistocene and Holocene evolution of the north Adriatic Sea. Il Quaternario 9(2): 697–704.
  • Cosendai A-C, Wagner J, Ladinig U, Rosche C, Hörandl E (2013) Geographical parthenogenesis and population genetic structure in the alpine species Ranunculus kuepferi (Ranunculaceae). Heredity 110: 560–569.
  • Dutrillaux AM, Lemonnier-Darcemont M, Darcemont C, Krpac V, Fouchet P, Dutrillaux B (2009) Origin of the complex karyotype of the polyploid parthenogenetic grasshopper Saga pedo (Orthoptera: Tettigoniidae). European Journal of Entomology 106: 477–483.
  • Dynesius M, Jansson R (2000) Evolutionary consequences of changes in species’ geographical distributions driven by Milankovitch climate oscillations. Proceedings of the National Academy of Sciences 97: 9115–9120.
  • Elmas A (2003) Late Cenozoic tectonics and stratigraphy of northwestern Anatolia: the effects of the North Anatolian Fault to the region. International Journal of Earth Sciences 92: 380–396.
  • Fernández R, Almodóvar A, Novo M, Simancas B, Díaz Cosín DJ (2012) Adding complexity to the complex: New insights into the phylogeny, diversification and origin of parthenogenesis in the Aporrectodea caliginosa species complex (Oligochaeta, Lumbricidae). Molecular Phylogenetics and Evolution 64: 368–379.
  • Fernández R, Novo M, Gutiérrez M, Almodóvar A, Díaz Cosín DJ (2010) Life cycle and reproductive traits of the earthworm Aporrectodea trapezoides (Dugès, 1828) in laboratory cultures. Pedobiologia 53: 295–299.
  • Gent PR, Danabasoglu G, Donner LJ, Holland MM, Hunke EC, Jayne SR, Lawrence DM, Neale RB, Rasch PJ, Vertenstein M, Worley PH, Yang Z-L, Zhang M (2011) The Community Climate System Model Version 4. Journal of Climate 24: 4973–4991.
  • Glesener RR, Tilman D (1978) Sexuality and the components of environmental uncertainty: clues from geographic parthenogenesis in terrestrial animals. The American Naturalist 112: 659–673.
  • Haag CR, Ebert D (2004) A new hypothesis to explain geographic parthenogenesis. Annales Zoologici Fennici 539–544
  • Heller K-G, Korsunovskaya OS, Sevgili H, Zhantiev RD (2006) Bioacoustics and systematics of the Poecilimon heroicus-group (Orthoptera: Phaneropteridae: Barbitistinae). European Journal of Entomology 103(4): 853–865.
  • Heller K-G, Lehmann A (2004) Taxonomic revision of the European species of the Poecilimon ampliatus-group (Orthoptera Phaneropteridae). Memorie della Societa Entomologica Italiana 82(2): 403–422.
  • Heller K-G, Reinhold K (2014) Specimen data to the phylogenetic study of the tribe Barbitistini (Orthoptera: Tettigonioidea: Phaneropteridae) in Ullrich et al. 2010. Articulata 29(1): 79–92.
  • Heller K-G, Sevgili H, Reinhold K (2008) A re-assessment of the Poecilimon syriacus group (Orthoptera Tettigonioidea, Phaneropteridae) based on bioacoustics, morphology and molecular data. Insect Systematics & Evolution 39: 361–379.
  • Heller K-G, Willemse L, Odé B, Volleth M, Feist R, Reinhold K (2011) Bioacoustics and systematics of the Poecilimon hamatus group (Tettigonioidea: Phaneropteridae: Poecilimon: Hamatopoecilimon n. subg.). Journal of Orthoptera Research 20: 81–95.
  • Hewitt G (1996) Some genetic consequences of ice ages, and their role in divergence and speciation. Biological Journal of the Linnean Society 58: 247–276.
  • Hewitt G (2004) Genetic consequences of climatic oscillations in the Quaternary. Willis KJ, Bennett KD, Walker D (Eds). Philosophical Transactions of the Royal Society of London. Series B: Biological Sciences 359: 183–195.
  • Hijmans RJ, Cameron SE, Parra JL, Jones PG, Jarvis A (2005) Very high resolution interpolated climate surfaces for global land areas. International Journal of Climatology 25: 1965–1978.
  • Holuša J, Kočárek P, Vlk R, Marhoul P (2013) Annotated checklist of the grasshoppers and crickets (Orthoptera) of the Czech Republic. Zootaxa 3616: 437–460.
  • Iorgu I, Pisica E, Pai P, Lupu G, Iuan C (2008) Checklist of Romanian Orthoptera (Insecta) and their distribution by eco-regions. Travaux du Muséum National d’Histoire Naturelle “Grigore Antipa 51: 119–135.
  • Junker RR, Lechleitner MH, Kuppler J, Ohler L-M (2019) Interconnectedness of the Grinnellian and Eltonian Niche in Regional and Local Plant-Pollinator Communities. Frontiers in Plant Science 10: 1371.
  • Kaya S (2018) Poecilimon zonatus Bolívar (Orthoptera, Tettigoniidae) revisited: genetic data revealed two new species and one new subspecies. Trakya University Journal of Natural Sciences 19(1): 85–99.
  • Kaya S, Boztepe Z, Çiplak B (2015) Phylogeography of the Poecilimon luschani species group (Orthoptera, Tettigoniidae): a radiation strictly correlated with climatic transitions in the Pleistocene. Zoological Journal of the Linnean Society 173: 1–21.
  • Kaya S, Chobanov D, Heller K-G, Yahyaoğlu Ö (2018) Review of Poecilimon species with inflated pronotum: description of four new taxa within an acoustically diverse group. Zootaxa 4462: 451–482.
  • Kaya S, Çiplak B, Chobanov D, Heller K-G (2012) Poecilimon bosphoricus group (Orthoptera, Phaneropterinae): iteration of morpho-taxonomy by song characteristics. Zootaxa 3225: 1–71.
  • Kenyeres Z, Bauer N (2008) Habitat selection and daily activity of Poecilimon intermedius (Fieber, 1853) (Orthoptera: Phaneropteridae)–autecological studies in a typical habitat of the species (Hungary). Russian Entomological Journal, 17(3): 247–257.
  • Kolics B, Ács Z, Chobanov DP, Orci KM, Qiang LS, Kovács B, Kondorosy E, Decsi K, Taller J, Specziár A, Orbán L, Müller T (2012) Re-Visiting phylogenetic and taxonomic relationships in the genus Saga (Insecta: Orthoptera). PLoS ONE 7(8): e42229.
  • Kumar S, Stecher G, Li M, Knyaz C, Tamura K (2018) MEGA X: Molecular Evolutionary Genetics Analysis across computing platforms. Molecular Biology and Evolution 35: 1547–1549.
  • La Greca M (1999) Il contributo degli Ortotteri (Insecta) alla conoscenza della biogeografia dell’Anatolia: la componente gondwaniana. Biogeographia—The Journal of Integrative Biogeography 20(1).
  • Lanfear R, Frandsen PB, Wright AM, Senfeld T, Calcott B (2017) PartitionFinder 2: New methods for selecting partitioned models of evolution for molecular and morphological phylogenetic analyses. Molecular Biology and Evolution 34(3): 772–773.
  • Lechner L, Zuna-Kratky T (2017) Checkliste der Heuschrecken Österreichs. Denisia, 39, Kataloge des oberösterreichischen Landesmuseums Neue Serie 184: 181–192
  • Lehmann GUC, Siozios S, Bourtzis K, Reinhold K, Lehmann AW (2011) Thelytokous parthenogenesis and the heterogeneous decay of mating behaviours in a bushcricket (Orthopterida). Journal of Zoological Systematics and Evolutionary Research 49: 102–109.
  • Lehmann GUC, Strauß J, Lakes-Harlan R (2007) Listening when there is no sexual signalling? Maintenance of hearing in the asexual bushcricket Poecilimon intermedius. Journal of Comparative Physiology A 193: 537–545.
  • Lemonnier-Darcemont M., Darcemont C (2020) Presence of Poecilimon pechevi Andreeva, 1978 (Orthoptera, Phaneropterinae) on Orvilos Mountain, Greece. Articulata 35: 87–91.
  • Lisiecki LE, Raymo ME (2005) A Pliocene-Pleistocene stack of 57 globally distributed benthic δ 18 O records: Pliocene-Pleistocene benthic stack. Paleoceanography 20(1).
  • Lynch M (1984) Destabilizing hybridization, general-purpose genotypes and geographic parthenogenesis. The Quarterly Review of Biology 59: 257–290.
  • Miller MA, Pfeiffer W, Schwartz T (2010) Creating the CIPRES science gateway for inference of large phylogenetic trees. In: 2010 Gateway Computing Environments Workshop (GCE). IEEE, New Orleans, LA, USA, 1–8.
  • Morgan-Richards M, Langton-Myers SS, Trewick SA (2019) Loss and gain of sexual reproduction in the same stick insect. Molecular Ecology 28: 3929–3941.
  • Mugleston JD, Song H, Whiting MF (2013) A century of paraphyly: A molecular phylogeny of katydids (Orthoptera: Tettigoniidae) supports multiple origins of leaf-like wings. Molecular Phylogenetics and Evolution 69: 1120–1134.
  • Muscarella R, Galante PJ, Soley-Guardia M, Boria RA, Kass JM, Uriarte M, Anderson RP (2014) ENMeval: An R package for conducting spatially independent evaluations and estimating optimal model complexity for Maxent ecological niche models. Methods in Ecology and Evolution 5: 1198–1205.
  • Nakano M, Morgan-Richards M, Godfrey AJR, Clavijo McCormick A (2019) Parthenogenetic females of the stick insect Clitarchus hookeri maintain sexual traits. Insects 10: 202.
  • Panrok A (2010) Erstnachweis der Mittleren Buntschrecke Poecilimon intermedius (Fieber, 1853) (Orthoptera, Ensifera) in Österreich. Beiträge zur Entomofaunistik 11: 89–93.
  • Petit JR, Jouzel J, Raynaud D, Barkov NI, Barnola J-M, Basile I, Bender M, Chappellaz J, Davis M, Delaygue G, Delmotte M, Kotlyakov VM, Legrand M, Lipenkov VY, Lorius C, PÉpin L, Ritz C, Saltzman E, Stievenard M (1999) Climate and atmospheric history of the past 420,000 years from the Vostok ice core, Antarctica. Nature 399: 429–436.
  • Phillips SJ, Anderson RP, Dudík M, Schapire RE, Blair ME (2017) Opening the black box: an open-source release of Maxent. Ecography 40: 887–893.
  • Rambaut A, Drummond AJ, Xie D, Baele G, Suchard MA (2018) Posterior summarization in Bayesian phylogenetics using Tracer 1.7. Susko E (Ed.). Systematic Biology 67: 901–904.
  • Ronquist F, Teslenko M, van der Mark P, Ayres DL, Darling A, Höhna S, Larget B, Liu L, Suchard MA, Huelsenbeck JP (2012) MrBayes 3.2: Efficient bayesian phylogenetic inference and model choice across a large model space. Systematic Biology 61: 539–542.
  • Saupe EE, Barve V, Myers CE, Soberón J, Barve N, Hensz CM, Peterson AT, Owens HL, Lira-Noriega A (2012) Variation in niche and distribution model performance: The need for a priori assessment of key causal factors. Ecological Modelling 237–238: 11–22.
  • Schwander T, Crespi BJ, Gries R, Gries G (2013) Neutral and selection-driven decay of sexual traits in asexual stick insects. Proceedings of the Royal Society B: Biological Sciences 280: 20130823.
  • Sergeev MG (2018) An annotated check-list of Orthoptera of Tuva and adjacent regions. Part 1. Suborder Ensifera. Far Eastern entomologist 372: 1–24.
  • Simon C, Buckley TR, Frati F, Stewart JB, Beckenbach AT (2006) Incorporating molecular evolution into phylogenetic analysis, and a new compilation of conserved polymerase chain reaction primers for animal mitochondrial DNA. Annual Review of Ecology, Evolution, and Systematics 37: 545–579.
  • Soberón J, Nakamura M (2009) Niches and distributional areas: Concepts, methods, and assumptions. Proceedings of the National Academy of Sciences 106: 19644–19650.
  • Soberón J, Peterson AT (2005) Interpretation of models of fundamental ecological niches and species’ distributional areas. Biodiversity Informatics 2: 1–10.
  • Stenberg P, Lundmark M, Knutelski S, Saura A (2003) Evolution of clonality and polyploidy in a weevil system. Molecular Biology and Evolution 20: 1626–1632.
  • Strauß J, Lehmann AW, Lehmann GUC (2014) Sensory evolution of hearing in tettigoniids with differing communication systems. Journal of Evolutionary Biology 27: 200–213.
  • Szövényi G, Skejo J, Rebrina F, Tvrtković N, Puskás G (2018) First data on the Orthoptera diversity of Poštak Mountain and its surroundings (Croatia). Annales de la Société entomologique de France (N.S. ) 54: 559–569.
  • Tilquin A, Kokko H (2016) What does the geography of parthenogenesis teach us about sex? Philosophical Transactions of the Royal Society B: Biological Sciences 371: 20150538.
  • Ullrich B, Reinhold K, Niehuis O, Misof B (2010) Secondary structure and phylogenetic analysis of the internal transcribed spacers 1 and 2 of bush crickets (Orthoptera: Tettigoniidae: Barbitistini). Journal of Zoological Systematics and Evolutionary Research 48(3): 219–228.
  • Ünal M (2010) Phaneropterinae (Orthoptera: Tettigoniidae) from Turkey and the Middle East II. Transactions of the American Entomological Society 136: 125–183.
  • van der Kooi CJ, Matthey-Doret C, Schwander T (2017) Evolution and comparative ecology of parthenogenesis in haplodiploid arthropods. Evolution Letters 1: 304–316.
  • Vandel APM (1928) La parthénogenèse géographique: contribution á l’étude biologique et cytologique de la parthénogenèse naturelle. Bulletin biologique de la France et de la Belgique 62: 164–281.
  • Vrijenhoek RC (1984) Ecological differentiation among clones: The frozen niche variation model. In: Wöhrmann K, Loeschcke V (Eds) , Population Biology and Evolution. Proceedings in Life Sciences. Springer Berlin Heidelberg, 217–231.
  • Vrijenhoek RC, Parker ED (2009) Geographical parthenogenesis: general purpose genotypes and frozen niche variation. In: Schön I, Martens K, Dijk P (Eds) , Lost Sex. Springer Netherlands, Dordrecht, 99–131.
  • Warchalowska-Šliwa E, Bugrov AG, Maryanska-Nadacowska A (1996) Karyotypes and C-banding patterns of some species of Phaneropterinae (Orthoptera, Tettigonioidea). Folia biologica (Kraków) 44(1): 5–10.
  • Warren DL, Seifert SN (2011) Ecological niche modeling in Maxent: the importance of model complexity and the performance of model selection criteria. Ecological Applications 21: 335–342.
  • Watanabe S, Hajima T, Sudo K, Nagashima T, Takemura T, Okajima H, Nozawa T, Kawase H, Abe M, Yokohata T, Ise T, Sato H, Kato E, Takata K, Emori S, Kawamiya M (2011) MIROC-ESM: model description and basic results of CMIP5-20c3m experiments. Geoscientific Model Development Discussions 4: 1063–1128.
  • Weekers PHH, De Jonckheere JF, Dumont HJ (2001) Phylogenetic relationships inferred from ribosomal ITS sequences and biogeographic patterns in representatives of the genus Calopteryx (Insecta: Odonata) of the West Mediterranean and adjacent West European Zone. Molecular Phylogenetics and Evolution 20: 89–99.
  • Willemse F (1982) A survey of the Greek species of Poecilimon Fischer (Orthoptera, Ensifera, Phaneropterinae). Tijdschrift voor Entomologie 125(6): 155–203.
  • Wright HE, Ammann B, Stefanova I, Atanassova J, Margalitadze N, Wick L, Blyakharchuk T (2003) Late-glacial and Early-Holocene dry climates from the Balkan Peninsula to Southern Siberia. In: Tonkov S (Ed.) Aspects of palynology and palaeoecology. Pensoft Publishers, Sofia-Moscow, 127–136.
  • Xia X (2018) DAMBE7: New and improved tools for data analysis in molecular biology and evolution. Kumar S (Ed.). Molecular Biology and Evolution 35: 1550–1552.

Supplementary materials

Supplementary material 1 

File 1

Borissov SB, Hristov GH, Chobanov DP (2021)

Data type: .xlsx

Explanation note: Table S1. Localities and accession numbers of sequences used in the analyses.

This dataset is made available under the Open Database License ( The Open Database License (ODbL) is a license agreement intended to allow users to freely share, modify, and use this Dataset while maintaining this same freedom for others, provided that the original source and author(s) are credited.
Download file (13.67 kb)
Supplementary material 2 

File 2

Borissov SB, Hristov GH, Chobanov DP (2021)

Data type: .pdf

Explanation note: Table S2. List of 95 % HPD intervals of the ND2 chronogram. Node numbers correspond to the numbers on Fig. 4.

This dataset is made available under the Open Database License ( The Open Database License (ODbL) is a license agreement intended to allow users to freely share, modify, and use this Dataset while maintaining this same freedom for others, provided that the original source and author(s) are credited.
Download file (398.54 kb)
Supplementary material 3 

File 3

Borissov SB, Hristov GH, Chobanov DP (2021)

Data type: .xlsx

Explanation note: Table S3. Model parameters estimated with ENMeval (Muscarella et al. 2014).

This dataset is made available under the Open Database License ( The Open Database License (ODbL) is a license agreement intended to allow users to freely share, modify, and use this Dataset while maintaining this same freedom for others, provided that the original source and author(s) are credited.
Download file (16.47 kb)
login to comment