The Balkan Peninsula is considered one of the most important Mediterranean refugia during the Quaternary glacial periods (Hewitt 2000). Multiple isolations and reconnections to Anatolia and Europe during the Neogene may underlie the huge biodiversity of this area with high levels of species richness and endemism. The region of the Balkan Peninsula is treated as a hotspot of biodiversity (Blondel and Aronson 1999; Myers et al. 2000; Mittermeier et al. 2003). Several land connections and submergences during the Miocene (23-5.33 Mya) and Pliocene (5.33-2.58 Mya) influenced the later development of this region (Steininger and Rögl 1984; Dermitzakis 1990; Popov et al. 2004; Husemann et al. 2014; Previšić et al. 2014; Poulakakis et al. 2015; Simaiakis et al. 2017; Španiel et al. 2017; Gömöry et al. 2020).

The Balkan Peninsula is at the forefront of the orthopteran diversity in the Palaearctic with over 40% of all European bush-crickets recorded from this region and new species being constantly described (Heller et al. 1998; Hochkirch et al. 2016). With the present study, we focus on one of the largest Palaearctic orthopteran genera, Poecilimon, comprising 145 species divided into 18 species groups (Cigliano et al. 2022). Members of the genus are distributed from the Apennines to Western Siberia and Central Tian-Shan (Bey-Bienko 1954) with the highest number of endemic species concentrated in the Aegean and Pontic areas. All species of Poecilimon are short-winged and flightless with complex acoustic communication. Cyclic glaciations during the Pleistocene influenced the diversity of the genus causing rapid radiation and diversification (La Greca 1999; Kaya et al. 2015; Borissov and Chobanov 2020; Borissov et al. 2020, 2021).

The taxonomy and phylogenetic relationships within Poecilimon are mainly based on morphological and bioacoustic traits (e.g., Heller et al. 2006, 2011; Chobanov and Heller 2010; Ingrisch and Pavićević 2010; Kaya et al. 2012, 2018; Boztepe et al. 2013; Sevgili et al. 2018; Chobanov et al. 2020). Many species groups of this genus have been studied in terms of molecular phylogeny and biogeography (Boztepe et al. 2013; Kaya et al. 2015; Kaya 2018; Borissov et al. 2020, 2021) while one of the largest groups – the Poecilimon ornatus group, has only recently been considered (Kociński 2020; Kociński et al. 2021). This species group contains bush-crickets distributed mostly in mountainous areas from the South-Eastern Alps to the Carpathians and Peloponnese and an isolated spot in Ukraine. The latest findings using cytochrome c oxidase subunit I (COI) barcodes showed the monophyly of the P. ornatus group (Kociński 2020). However, there is still an unclear relationship among the taxa associated with the Poecilimon affinis complex in the P. ornatus group (Chobanov and Heller 2010; Kociński 2020; Kociński et al. 2021). Currently, the P. affinis complex includes P. nonveilleri, P. pseudornatus and five subspecies of P. affinis (P. a. affinis, P. a. hajlensis Karaman, 1974, P. a. serbicus Karaman, 1974, P. a. komareki Cejchan, 1957, P. a. dinaricus Ingrisch & Pavićević, 2010). Recent studies suggested extending this complex with P. hoelzeli Harz, 1966 and P. ornatus (Schmidt, 1850) (Kociński 2020; Kociński et al. 2021).

‘Species complex’ refers to a group of sibling species with similar morphology or identical populations that are reproductively isolated (Mayr 1963; Sigovini et al. 2016) or cryptic species, where the boundaries between taxa are morphologically indeterminate. ‘Species complex’ has also been defined as consisting of closely related taxa that are still waiting for critical revision to clarify their taxonomic status (Sigovini et al. 2016). Cryptic species were defined as “two or more distinct species that are erroneously classified (and hidden) under one species name” (Bickford et al. 2007). In this sense, the P. ornatus group constitutes one or more species complexes that need to be resolved using interdisciplinary research.

Molecular data and species delimitation methods have become very important tools to detect and delimit new species (Luo et al. 2018; Mendes et al. 2021). DNA sequence analysis has revolutionized the way of recognizing species (Hajibabaei et al. 2007; Taylor and Harris 2012) and helped to reveal the existence of cryptic species in many taxa (Knowlton 1993; Bickford et al. 2007; Scheffers et al. 2012). The cytochrome c oxidase subunit I (COI) gene is a commonly used marker, easy to amplify due to the availability of conserved primers, with a strong phylogenetic signal, used in taxonomy (Folmer et al. 1994; Simon et al. 1994, 2006; Spicer 1995; Zhang and Hewitt 1997; Goto and Kimura 2001; Remigio and Hebert 2003; Kjer et al. 2014; Wang et al. 2017; Jafari et al. 2019; Karmazina et al. 2020). This marker is successfully used in Orthoptera and treated as a DNA barcode (Lehmann et al. 2017; Kaya and Çıplak 2018; Kundu et al. 2020; Liu and He 2021; Şirin et al. 2021; Warchałowska-Śliwa et al. 2021). NADH dehydrogenase subunit 2 (ND2) shows a higher proportion of variable and parsimony-informative sites (PI) and a lower heterogeneity of the substitution index than COI (Cheng et al. 2018), which was confirmed in Isophya – a closely related genus to Poecilimon (Chobanov et al. 2017), and in Hematopoecilimon (Borissov and Chobanov 2020). The control region (CR) is mainly used to study phylogenetic relationships in closely related taxa (Amaral et al. 2016; Li and Liang 2018), successfully tested in Poecilimon (Eweleit et al. 2015; Borissov and Chobanov 2020). The internal transcribed spacer 1 (ITS1) region represents a useful marker for the analysis of relationships in closely related species of Orthoptera and for recognition of new species because of higher evolutionary rates leading to greater variability in both, nucleotide sequence and length (Hillis and Dixon 1991; Gu et al. 2020). In this study, we perform molecular analyses of taxa in the P. ornatus group using a combined dataset (COI, ND2, CR, and ITS1).

Our study aims to reconstruct the phylogenetic relationships among taxa in the P. ornatus group and to elucidate the position of taxa related to the P. affinis complex. We test the hypothesis of a recent origin and divergence of the taxa in the P. affinis complex from the rest of the species in the P. ornatus group. The estimated divergence times were applied to test the correlation between the evolutionary history of this group and paleogeographic events in the Balkan Peninsula. Additionally, phylogeographical biogeographic tools were used to check if speciation was affected by vicariances, dispersal, and/or extinction events.

The final alignment of the COI sequence results in 607 bp with 129 parsimony-informative sites and 196 variable sites. The CR (including the 12S rDNA gene containing A+T-rich region) consists of 446 bp with 188 parsimony-informative and 272 variable sites. ND2 sequences include 695 bp, among them 168 are parsimony-informative and 245 variable sites. The final alignment of ITS1 sequences consists of 465 bp with 70 parsimony-informative and 130 variable sites. The combined matrix data of COI, ND2, CR, ITS1 consists of 2213 bp and involved six outgroup species. The genetic mean distance for CO1 and ND2 among taxa from the P. affinis complex is 0.02, whereas among the rest of the species from the P. ornatus group – 0.1. For CR, the genetic mean distance among taxa from the P. affinis complex is 0.05, among the rest of the species from the P. ornatus group is 0.2. The genetic mean distance for ITS1 is 0.04 for taxa from the P. affinis complex, and 0.09 for the rest of the taxa from the P. ornatus group. The genetic distances between species from the P. affinis complex and the P. ornatus group for each marker (COI, ND2, CR, ITS1) are available in Table 4.

The results of the substitution saturation test for COI, ND2, and CR alignments are summarized in Table 5. Calculated P-values were significant for all gene alignments and Iss (index of substitution saturation) values were lower than Iss.c (critical index of substitution saturation) in all cases. No saturation of the phylogenetic signal was observed for the COI, ND2, and CR datasets.

The substitution one-parameter model Jukes–Cantor (JC) with Gamma Distribution (G) and Invariable site (I) was the best fit for the COI, ND2, CR and ITS1 data matrix.

The BI and ML phylogenetic trees showed the same topology (Fig. 4) and confirmed the monophyly of the P. ornatus group (posterior probability support, PP = 1.0; bootstrap support, BP = 100), whereas the P. affinis complex was paraphyletic as suggested in Kociński (2020). The first clade consists of P. nobilis, P. artedentatus and P. obesus. The second clade includes P. gracilis, P. jablanicensis, and P. soulion. Poecilimon gracilioides and P. pindos occupy the branches between the second and third clade. The third clade consists of the taxa from the P. affinis complex: P. affinis affinis, P. a. dinaricus, P. poecilus, P. a. komareki, P. a. serbicus, P. nonveilleri, P. a. hajlensis, P. rumijae, P. pseudornatus; and two additional species: P. hoelzeli and P. ornatus. Poecilimon a. affinis is the most diverse taxon among the complex, which supports recent studies (Kociński 2020; Kociński et al. 2021). Poecilimon a. affinis, from Rilski Manastir and the Rila Mts., seems to be a sister taxon to the remaining representatives of the P. affinis complex. Poecilimon rumijae forms a separate branch among the third clade, as does P. poecilus, which is treated as a synonym of P. a. affinis according to the current systematics (Cigliano et al. 2022). Specimens of P. pseudornatus are grouped regardless of their location. Moreover, the phylogenetic relationship between taxa does not correlate with their place of occurrence (Fig. 4 – Locality).

Five species delineation tests revealed different taxonomic schemes that disagreed on some points with each other and with the current taxonomic classification. As a result of the ASAP analysis (Fig. 4 – ASAP), a barcoding gap of about 2–10% was estimated. The pairwise distance gap approach (Fig. 4 – ASAP) identified from 2 to 43 hypothetical species. We chose the fifth ASAP-score (6.50) which provides the best-fit scenario at the threshold distance of 2.68% (JC69) with 9 hypothetical species. The maximum-likelihood approach (Fig. 4 – GMYC) defined 34 species under a single threshold and 26 under multiple thresholds. The pairwise distance gap approach (Fig. 4 – ABGD) with the default settings (X = 0.5) suggested 9 groups with prior intraspecific divergence (P) reaching 0.007, while 36 groups were defined with P ≤ 0.001. For bPTP (Fig. 4 – bPTP ML), we conducted two analyses based on BI and ML approaches. BI showed 52 species, whereas ML identified 9 groups or species. Thus, only ML was used in this study. ASAP, ABGD, and bPTP grouped species from the P. affinis complex, P. hoelzeli and P. ornatus into one species, whereas GMYC recognized 17 species among the complex.

The time estimation analysis dated the last common ancestor (LCA) of the P. ornatus group at 1.62 Mya with the following main lineage splits dated between 1.33 and 0.42 Mya (Fig. 2) during the Calabrian and Chibanian stage of the Pleistocene. The divergence of the P. affinis complex from P. pindos was dated at ca. 0.71 Mya during the Pleistocene (95% –confidence interval) based on the molecular clock analysis and a priori calibration. The LCA of the P. affinis complex was dated at ca. 0.42-0.02 Mya in the Late Pleistocene.

The distribution pattern of the P. ornatus group results in six dispersal and five vicariance events (Fig. 2). The LCA of the group was positioned in the AB area and the group evolved by a vicariant event and subsequent dispersal within the Southern (A) and Central (B) areas where local lineage splits occurred. The Central region also represents the main speciation and dispersal centre of the Poecilimon ornatus group. From here, the Poecilimon affinis complex-ancestor evolved by dispersal in two main directions – North-West and (North-)East, where local dispersal and vicariant events contributed to the recent evolutionary history of the complex. Within the crown lineages, though poorly resolved, worth mentioning as stepping-stone - dispersal taxa are Poecilimon hoelzeli – distributed at the border of the Central with the (North-) Eastern lineage, and Poecilimon pseudornatus, having quite a wide distribution in the Central and North-Western regions. There was no correlation between genetic mean distance and geographic pattern in the P. ornatus group (Mantel Test, R = 0.0469; p = 0,193).

The present study represents the first comprehensive attempt to reconstruct the molecular phylogeny of the Poecilimon ornatus group. The molecular results support the monophyly of the P. ornatus group, as suggested in recent studies, based on ITS1, ITS2, 16S rRNA, tRNA-Val, 12S rRNA (Ullrich et al. 2010; part of the taxa), and the COI gene (Kociński 2020).

The Control region is the most variable marker, as confirmed in the previous studies on Poecilimon (Eweleit et al. 2015; Borissov and Chobanov 2020). It shows the highest genetic mean distance between taxa from the P. affinis complex and the remaining species from the P. ornatus group. The Control region is a useful phylogenetic marker with the potential of providing better resolution than COI (Vila and Björklund 2004; Cheng et al. 2018). The number of variable and PI sites in ND2 is about 20% higher than in COI which is similar to the results provided for Isophya (Chobanov et al. 2017). However, the internal transcribed spacer 1 (ITS1) region contains the lowest number of variable and PI sites.

Poecilimon nobilis, P. artedentatus, and P. obesus form the sister clade to the remaining species of the group. The latter lineage is consistent with the morphological similarity of these three species (Chobanov and Heller 2010). The present data do not confirm that P. gracilis is the sister species to the remaining taxa of the P. ornatus group, as suggested in previous studies based on morphology, bioacoustics (Chobanov and Heller 2010) and molecular data (Ullrich et al. 2010; Kociński 2020). Poecilimon gracilis is morphologically similar to P. jablanicensis and occurs parapatrically with the latter (Chobanov and Heller 2010) which is a prerequisite for close relationships as supported by our molecular results, where these species occupy the same subclade with P. soulion (Fig. 4). The sister clade to the latter includes the lineages of P. gracilioides, P. pindos, and the clade richest in taxa forming the P. affinis complex (Chobanov and Heller 2010; Kociński 2020; Kociński et al. 2021). Poecilimon hoelzeli and P. ornatus are placed among the taxa of the complex. Thus, the P. affinis complex is paraphyletic when these two species are not included. This finding is consistent with the previous studies (Kociński 2020; Kociński et al. 2021). Poecilimon pseudornatus occupies one subclade, regardless of where it occurs (North Macedonia (MK): Jablanica Mt.; Montenegro (MN): Durmitor, Treshnievik, Vusanje, Hajla; Serbia (SR): Kamena Gora) (Figs 1, 2), which corresponds to the low morphological variability of the species (Kociński et al. 2021). We can notice a distant genetic relationship between P. a. komareki and P. rumijae, which contradicts the current systematics where P. rumijae is treated as a synonym of P. a. komareki (Cigliano et al. 2022). Moreover, the results based on the geometric morphometric method of male pronotum and ovipositor confirmed that P. rumijae and P. a. komareki may be separate taxa (Kociński et al. 2021). This assumption is in line with the opinion of Ingrisch and Pavićević (2010), regarding P. rumijae as a species of the P. ornatus group, comparing it to P. nonveilleri. Nevertheless, as discussed by Kociński et al. (2021), P. nonveilleri does not seem to be closely related to P. rumijae, while the shape of the cercus and tegmen, length of the stridulatory row and number of stridulatory teeth in P. affinis komareki and P. rumijae show great similarity. In addition, the third clade (P. affinis complex) shows very low genetic structuring and low genetic variation, with poor resolution between groups of different taxonomic level. Specimens of P. a. affinis from different localities (Bulgaria (BG): Pirin Mts., Bratiya, Osogovo, Kirilova Polyana, Rila Mts., Rilski Manastir) form separate subclades (Figs 1, 4). Our results were confirmed by a geometric morphometric analysis of the male tegmen, cercus, pronotum, and ovipositor, where P. a. affinis was the most diffuse taxon among the group (Kociński et al. 2021). The above data suggest an infraspecific division of some local populations of Poecilimon a. affinis and contradict the assumption that the variability within this taxon depends mostly on the altitude of occurrence (Chobanov and Heller 2010). Despite the genetic variability in P. a. affinis from different localities, the Mantel test suggested no association between genetic and geographic distances in this group. Our results, based on three species delimitation methods (ASAP, ABGD, bPTP) (Fig. 4), suggest to divide the P. ornatus group into nine potential species, which contradicts the morphological, bioacoustics (Chobanov and Heller 2010; Ingrisch and Pavićević 2010; Kociński et al. 2021), and earlier molecular data (Kociński 2020). On the other hand, GMYC analysis reveals 26 hypothetical species among the group. The discrepancy in the results of species delimitation may indicate a greater conservatism of ASAP, ABGD, and bPTP over GMYC, which shows lower efficiency in data sets at the genus than at higher levels (Magoga et al. 2021). Though species delimitation has been defined as a method that sometimes causes confusion about almost every aspect of the definition of the ‘species’ level (Stanton et al. 2019), the problem with delineating species’ boundaries at the tree top must be related to the low-level independent genetic differentiation of the third clade in our tree. Based on the recent lineage splits (Fig. 2) and the large number of taxa occurring over a significant geographic area (most of the central and northern part of the Balkan Peninsula reaching the Eastern Alps and Carpathians), we assume a recent contemporary allopatric origin of the taxa within the Poecilimon affinis complex. The latter may still be in the genetic “gray” zone of speciation, forming clines of a multitude of phenotypes with poor genetic structure (de Queiroz 1998). In conclusion, our results confirmed the existence of the P. affinis complex, though they failed at separating species.

Poecilimon consists of groups of poorly morphologically distinguishable units/taxa that have been subjected to a rapid diversification following the set of the Miocene and especially during the Plio-Pleistocene climatic cycles (Borissov et al. 2020). According to our molecular clock (Fig. 2), most speciation processes in the P. ornatus group occurred between the middle Pleistocene (ca. 1.62 Mya) and the beginning of the Holocene (ca. 0.01 Mya). The dating of LCA of the P. ornatus group (1.62 Mya) coincides with a significant global climate cooling, which was also connected with the expansion of cold climate-adapted fauna in the North Atlantic (Lisiecki and Raymo 2005). Though most taxa of the group tend to occur in humid mountain areas with cool climates, the first clade of the group involves two species occurring in the lowland and middle-mountain belts in the Southern biogeographical region (in Peloponnesos) (P. nobilis and P. artedentatus) and one species with a narrower temperature tolerance (P. obesus) occurring in the lowlands of the Southern and southern part of the Central region (Chobanov and Heller 2010). Thus, the first lineage split in the group may have happened as a result of isolation due to climate deterioration in the Central or Southern region of distribution of the group (S and W Balkans) and subsequent adaptation of new lineage(s) with northern distribution to a cooler climate.

The following major lineage splits fall within the period called the Middle Pleistocene transition when climate cycles gradually changed from 41- to 100-Ka periods. This switch started ca. 1.25 Mya and after interruption continued after 0.9 Mya to be established ca. 0.7 Mya (Lisiecki and Raymo 2005; Clark et al. 2006). Within this irregular repetition of warmer, colder, wetter and dryer periods of variable temperature and humidity amplitude, multiple range shifts, accompanied by isolation and extinction events were driven. Thus, species like Poecilimon jablanicensis may have evolved from its ancestor, P. gracilis, from small populations subjected to the severe climate being isolated at mountain ridges by dense forest belt. The latter pattern may be applied to the origin of P. pindos, P. gracilioides and P. soulion, which possibly due to a wider ecological tolerance and/or eco-graphic factors have spread to a few or more mountain ranges.

The so-called Mid-Brunhes Transition ca. 430 ka ago marks a sharp increase in the temperature amplitude of the Pleistocene climate cycles (Barth et al. 2018). This time corresponds to a thermal minimum (l.c.), preceded by a minimum in the solar radiation in Europe (Boryczka and Stopa-Boryczka 2004) and concurs with the cold Marine Isotope Stage MIS 12 (478-424 ka ago) that was followed by Glacial Termination with a very large magnitude (Lisiecki and Raymo 2005). The time to LCA of the Poecilimon affinis complex (Fig. 2) corresponds well with the Mid-Brunhes Transition and interestingly – with the results for the two major lineage splits of the Poecilimon ampliatus complex (see Borissov et al. 2021). The larger temperature amplitudes with colder glacials and a larger decrease in humidity should be the main trigger for dispersal, isolation (vicariance), extinction, and ecological adaptation in the Poecilimon affinis complex, similarly to many other animals (Hewitt 1996, 2000; Taberlet et al. 1998; Wallis et al. 2016). As the multitude of geographic taxa within the Poecilimon affinis complex shows an overall low genetic differentiation of similar scale and a wider distribution than the ancestral lineages of the Poecilimon ornatus group, its evolution should have been ruled by fast spreading within comparatively short climatically favorable periods during the last two glacial periods. During this vast expansion accompanied by versatile morpho-acoustic diversification, distinct ecological forms evolved, including both mountain specialists (e.g., geographic forms of P. affinis s.str.), ecologically tolerant species (P. ornatus, P. pseudornatus), and early-seasonal Mediterranean species (P. a. komareki, P. ‘rumijae’ – synonym of P. a. komareki).

The ancestor(s) of the Poecilimon affinis complex splits off from the rest of the P. ornatus group in the Pleistocene (ca. 0.71 Mya). The results of the molecular clock confirmed the need to extend the complex with two species: P. ornatus and P. hoelzeli. The P. affinis complex diverged into two lineages ca. 0.42 Mya. The first lineage consists of P. hoelzeli, P. pseudornatus, P. a. komareki, P. poecilus, P. rumijae, P. a. serbicus, P. nonveilleri, P. a. hajlensis, which are partly consistent with their biogeographical regions (Central and North-Western). The second lineage includes species from the Eastern (P. ornatus, P. a. affinis), and North-Western regions (P. a. dinaricus).

The present study generated additional evidence for the relationships within the P. ornatus group. Our results indicate that COI, ND2, CR, and ITS1 markers can be successfully used for phylogenetic analyses, supporting the previous studies on the phylogeny of Poecilimon. The presented results confirmed the monophyly of the P. ornatus group and the existence of the P. affinis complex containing two additional species: P. hoelzeli and P. ornatus. Using phylogenetic and time estimation analyses, biogeographic reconstruction, and available paleoclimatic data, we reveal the origin and evolutionary patterns of the Poecilimon ornatus group and shed light on the climate-driven complex evolution of the Poecilimon affinis complex. These young taxa were formed by speciation modulated by dispersal, vicariance, and extinction events, and directed towards phenotypic and ecological diversification.