New insights into the genetic diversity of the Balkan bush-crickets of the Poecilimon ornatus group (Orthoptera: Tettigoniidae)

the Abstract The Balkan Peninsula is treated as a hotspot of biodiversity with over 40% of European bush-crickets occurring there. Poecilimon Fischer, 1853 is one of the largest Palaearctic orthopteran genera containing several species groups. One of them is the Poecilimon ornatus group (Schmidt, 1850) with 13 species and 5 subspecies. Among the group, the Poecilimon affinis complex is designated as consisting of P. pseudornatus Ingrisch & Pavićević, 2010, P. nonveilleri Ingrisch & Pavićević, 2010, and five subspecies of P. affinis (Frivaldszky, 1868). The aim of this study is to reconstruct the phylogenetic relationships among taxa of the P. ornatus group and to elucidate the position of taxa related to the P. affinis complex. Molecular phylogeny supported the monophyly of the P. ornatus group and showed that their ancestor probably originated in the southern Balkans. The underlying processes are thought to be six dispersals and five vicariance events linked to geological events and climate changes in the Pleistocene. The species delimitation analysis showed mostly nine hypothetical species among the group.


Introduction
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;Mitter-meier 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; 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 shortwinged 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 et al. , 2021. The taxonomy and phylogenetic relationships within Poecilimon are mainly based on morphological and bioacoustic traits (e.g., Heller et al. 2006Heller et al. , 2011Chobanov and Heller 2010;Ingrisch and Pavićević 2010;Kaya et al. 2012Kaya et al. , 2018Boztepe 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. , 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 andP. 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. 1994Simon et al. , 2006Spicer 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 . 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;. 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.

Taxon sampling
A total of 74 specimens from 34 populations representing 19 formerly recognized taxa of the Poecilimon ornatus group were used in this study (Table 1). Six outgroup species were selected representing three other species

Molecular laboratory procedure
DNA was extracted from hind leg-muscle tissue using the NucleoSpin tissue kit (Macherey-Nagel, Germany) according to the manufacturer's protocol. Genomic DNA was used for the amplification of three mitochondrial markers (COI, ND2, CR) and one nuclear marker (ITS1). The Polymerase chain reaction (PCR) primer pairs used in this study are included in Table 2. The amplification was performed in 25 µl reaction volume containing 12.5 µl 2x Phanta Max Master Mix (Vazyme, China), 10 mM dNTP mixture, 10 µM forward and reverse primers, 1-3 µl genomic DNA, and sterile deionized water. The PCR protocols used for amplification of COI, ND2, CR, and ITS1 are included in Table 3. All PCR products were purified using Exo-BAP Mix (EURx, Poland, following the standard protocol). The sequencing reaction was carried out in 10 µl reactions containing: 1.5 µl of sequencing buffer, 1.0 µl of BrilliantDye TM v3.1 Terminator Cycle Sequencing Kit (NimaGen, The Netherlands), 1.0 µl of primer (forward or reverse), 3.0 µl of the purified DNA and 3.5 µl of sterile water. The sequencing protocol was as follows: the initial melting step of 3 min at 94°C followed by 25 cycles of 10 s at 96°C, 5 s at 55°C and a final step of 90 s at 60°C. The obtained sequences were deposited in GenBank (www.ncbi.nlm.nih.gov/genbank) under the accession numbers provided in Table 1. Additionally, 85 DNA sequences were acquired from Gen-Bank. The nucleotide sequences were edited and aligned in CodonCode Aligner 9.0 (CodonCode Corporation; https://www.codoncode.com/aligner) with default parameters. All sequences were checked for stop-codons in MEGA 11 (Tamura et al. 2021), verified using BLAST of NCBI (http://blast.ncbi.nlm.nih.gov/Blast.cgi). Genetic distances were calculated using MEGA 11 (Tamura et al. 2021). The saturation of the nucleotide substitution was checked for CR, ND2, and two separate partitions of COI (with codon positions 1 + 2 and codon position 3) (Xia et al. 2003) through the substitution saturation test in DAMBE (Xia 2013). The partition homogeneity test (Farris et al. 1995) was conducted in PAUP (Swofford 2002) with 1000 replicates to determine whether all regions (COI, ND2, CR, ITS1) could be combined in a unique data matrix.

Phylogenetic analyses
To infer evolutionary relationships, two methods were used -Bayesian inference (BI) and maximum likelihood (ML). The substitution model of evolution was estimated in MrModeltest software (Nylander 2004) using the Akaike Information Criterion (AIC). MrBayes (Ronquist et al. 2012) was used to obtain the Bayesian tree (BI). Posterior probabilities were based on two independent Markov chain Monte Carlo (MCMC) runs, each composed of four chains (three heated chains and one cold chain). BI was performed for 6,000,000 generations, with a sampling of trees every 100 generations. The convergence of the analyses was validated by monitoring the likelihood values using Tracer (Rambaut et al. 2018). Maximum likelihood (ML) estimates of the phylogeny were conducted using IQ-TREE (Nguyen et al. 2015

Sequence-based species delimitation test
To detect independently evolved lineages, three different DNA sequence-based species delimitation approaches were chosen. The first approach was the general mixed Yule-coalescent (GMYC) model. It uses the maximum likelihood approach based on the prediction that independent evolution leads to the appearance of distinct genetic clusters (Fujisawa and Barraclough 2013). This approach was successfully used for detecting cryptic lineages (e.g., Pons et al. 2006;Jörger et al. 2012;Chobanov et al. 2017). The next approaches were the Automatic Barcode Gap Discovery (ABGD) and Assemble Species by Automatic Partitioning (ASAP). These methods use pairwise distances to group sequences into potential species based on detecting gaps in the variation between supposed in-tra-and interspecies groups (barcode thresholds) (Puillandre et al. 2012(Puillandre et al. , 2021. The last method was the Poisson Tree Processes (bPTP), which is mainly intended for delimiting species in single-locus molecular phylogenies (Zhang et al. 2013).

Estimation of divergence time and biogeographic analysis
To date the most recent common ancestor, the Bayesian approach with an MCMC integration was used in BEAST (Drummond et al. 2012) based on COI sequences. In order to follow the phylogenetic tree-topology, we have constrained monophyly for the well-supported clades of the P. ornatus group, while monophyly was not set for the branches within the P. affinis complex due to poor resolution. The analysis was run for 10,000,000 generations with sampling every 1,000 generations and a 10% burn-in. For time estimation analyses, an uncorrelated lognormal relaxed clock was applied (Drummond et al. 2006). The convergence to stationary distribution and the effective sample size of model parameters were checked using Tracer. The maximum clade credibility trees were built with TreeAnnotator (Drummond et al. 2012). In a recent study, divergence dates in Poecilimon were estimated based on the minimum time of isolation of Poecilimon cretensis, endemic to the island of Crete . As a result, an intraspecific lineage split between the easternmost and the other lineages of P. cretensis was estimated at 0.8 Ma, possibly reflecting former vicariant events as a result of the former disconnection of the easternmost part of Crete. The latter dating is here used as a secondary calibration to date recent divergence times in the P. ornatus species group.
Poecilimon 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) (see also Borissov et al. 2021). In order to infer the biogeographic history of the Poecilimon ornatus group, we first selected areas defined as centers of endemism. As most taxa concerned are regional endemics (occurring in a mountain range or a geographic outline of a few mountain ranges and/ or valleys) and only one species (Poecilimon jablanicensis Chobanov & Heller, 2010) is strictly a local endemic, the regions selected cover the geographical extent of a few sympatric taxa. Thus, wider distributed species may occur in more than one region. As a result, four biogeographical regions (Fig. 2, 3 Lunt et al. 1996 ND2 TM-J210 (F) TW-N1284 (R)

A A T T A A G C T A A T G G G T T C A T A C C C A Y A G C T T T G A A R G Y T A T T A G T T T Simon et al. 2006
CR SR-J14610 (F) T1-N18 (R)

Results
The  ic 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 hoelzelidistributed 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).

Discussion
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;. 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).  (Figs 1, 2), which corresponds to the low morphological variability of the species . 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 . 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 . 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 spe-cies 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 ). 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(Hewitt , 2000Taberlet 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).

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

Acknowledgements
We thank the Biology Students' Research Society (BSRS; Skopje, Republic of North Macedonia) and its 2017 Chair Marija Trencheva for the accommodation and logistic support, and Slobodan Ivković for the help in the field, during our collecting trips in North Macedonia.
This work was partly supported by a joint research project between the Bulgarian Academy of Sciences and the Polish Academy of Sciences (project Convergent evolution of polyphyletic bush-crickets (Orthoptera: Phaneropterinae): micropterism and speciation). DC was supported by Grant DN11/14-18.12.2017 from the National Science Fund (MES) of Bulgaria.