Evolution of Poecilimon jonicus group (Orthoptera: Tettigoniidae): a history linked to the Aegean Neogene paleogeography

The Aegean archipelago is among the largest on Earth with astonishing biodiversity within Europe. Its territory underwent a massive geotectonic transformation in Neogene that resulted in multitude of changes in land-sea configuration and disintegrated the formerly united Aegean land to a complicated mainland-archipelago system. Therefore, it represents an excellent laboratory for studying evolution of terrestrial fauna. In the present study, we use a model group of flightless bush crickets with annual reproduction cycle—Poecilimon jonicus species group—to trace correlation of lineage diversification with the known paleogeographic events in the Aegean area. The group belongs to the hyperdiverse genus Poecilimon and has a disjunct distribution along the Hellenic arc from southwestern Anatolia through Crete to the western Balkans and the Apennines. To test our hypothesis, we inferred phylogenetic relationships of the P. jonicus group sensu lato using a nuclear fragment covering two spacers of the ribosomal cistron (ITS1 + ITS2). To study intra-group phylogeny, we compared mitochondrial phylogenies based on two matrices—(1) a concatenated ND2 and COI dataset of 1656 bp and (2) a 16S rRNA + 12S rRNA dataset of 1835 bp. As a second step, we estimated divergence times applying Bayesian approach with BEAST and a relative rate framework with RelTime on the mitochondrial matrices. We compare trees calibrated based on evolutionary rates and tectonic events and discuss radiation scenarios in concordance with known paleogeographic events in the Aegean area. Our results revealed robust phylogeny of the Poecilimon jonicus group and confirmed a strong link between its evolution and the Aegean paleogeography. The phylogenetic relationships of the group supported reconsideration of its systematics.


Introduction
Archipelagoes provide rich and comparatively simple experimental sets for studying the evolution of terrestrial organisms due to their small size, isolated biotas, and possibilities for a more precise dating of vicariant events in comparison with mainland systems (Shaw and Gillespie 2016). In addition, isolated islands maintain flightless fauna due to a combination of factors selecting against unwanted spread and excessive energy and water loss. Vogler and Timmermans (2012) state that "the lack of dispersal may not only result in deeper genetic subdivision, but also in a speed-up of molecular rate of change"; hence, micropterous island populations may be subjected to accelerated diversification. A combination of a diverse group with fast reproduction and low mobility occurring on a diverse archipelago-mainland system, therefore, is an excellent object for phylogeographic studies.
In the Mediterranean, the Aegean represents one of the largest archipelagos on Earth with more than 7500 islands and islets (Triantis and Mylonas 2009), making it especially interesting for phylogeographic studies. The post-Oligocene collision of the African and Arabian tectonic plates with Eurasia led to the disintegration of the formerly united Aegean landmass, consequently forming the current land-sea configuration (Creutzburg 1963;Steininger and Rögl 1984;Dermitzakis 1990). In late Miocene, tectonic events formed basins and grabens, which were filled by the Tortonian transgression of the Mediterranean Sea, resulting in the first formation of the Aegean Sea that established Electronic supplementary material The online version of this article (https://doi.org/10.1007/s13127-020-00466-9) contains supplementary material, which is available to authorized users. a connection with the Euxinian (Pontian) Basin (Meulenkamp 1977;Dermitzakis 1990;Steininger and Rögl 1984;Jolivet et al. 2006;Popov et al. 2004Popov et al. , 2010. This marine barrier was largely interrupted in the Messinian in two phases-first one (ca. 5.8 Ma ago) is characterized with a smaller-scale sea drop, while in the second (5.6 Ma ago), the basin was drastically affected (Popov et al. 2004) by the desiccation of the Mediterranean resulting in an over 1500 m drop of the sea level (Dermitzakis 1990;Jolivet et al. 2006;Loget et al. 2006). Re-flooding of the Aegean by the beginning of the Pliocene strongly affected the early evolution of terrestrial biota in this area, while later (Plio-Pleistocene) vicariant and dispersal events were ruled by smaller-scale land-sea configuration changes due to climatic or paleogeographic processes, mostly within the areas separated by the three main Aegean basins and the so-called mid-Aegean trench (reviewed by Poulakakis et al. 2014).
Interdisciplinary studies aiming phylogeographic assumptions accumulated a significant amount of data for various animal groups in the Aegean area (Gantenbein and Keightley 2004;Çıplak 2004;Kasapidis et al. 2005;Poulakakis et al. 2005aPoulakakis et al. , 2005bDouris et al. 2007;Simaiakis and Mylonas 2008;Çıplak et al. 2010;Papadopoulou et al. 2010). The long history of the Aegean archipelago and availability of paleogeographic, palynological, and paleontological data for this area (e.g., reviews by Dermitzakis 1990; Anastasakis et al. 2006) provide reliable calibration points for divergence time estimations and allow calibration with mutation rates and hypothesizing about major speciation events (Lymberakis et al. 2007;Papadopoulou et al. 2010). As a result, the number of phylogeographic studies in the area recently increased significantly (for a review, see Poulakakis et al. 2014).
Recent studies affirmed Orthoptera as a marker group for phylogeographic reconstructions and testing speciation processes (Çıplak 2004;Çıplak et al. 2010;Korkmaz et al. 2014;Allegrucci et al. 2017;Chobanov et al. 2017;Kaya and Çıplak 2017). The vast diversity of tribus Barbitistini (Ensifera, Tettigoniidae) in the Aegean area (with ca. 10 to 20% of all Orthoptera in the region) indicates the latter as the possible center of origin of Barbitistini. Poecilimon Fischer, 1853 is the largest genus within this tribus with remarkable diversity in the Eastern Mediterranean (La Greca 1999;Çıplak 2004). This genus represents a group of micropterous bush crickets with low mobility and complicated acoustic communication system that has obviously been subjected to a rapid diversification following the set of the Miocene and especially during the Plio-Pleistocene climatic cycles (compare Çıplak 2004;Kaya et al. 2015;Chobanov et al. 2017). The latter has led to a vast variety of closely related lineages, forming groups of poorly distinguishable morphologically, yet behaviorally and ecologically specialized taxa, most of them currently grouped according to morpho-acoustic and/or genetic similarity (e.g., Heller 1984;Ullrich et al. 2010;Chobanov and Heller 2010;Heller et al. 2011).
The main objectives of the present study are to, first, reconstruct phylogenetic relationships of P. jonicus group sensu lato and consider its systematics in the light of new knowledge and, second, estimate divergence times for the group members and test if radiation steps in their evolutionary history correlate with the known paleogeographic events in the Aegean area. To test these hypotheses, we used sequence data from one nuclear DNA dataset covering two spacers of the ribosomal cistron (ITS1 and ITS2-internal transcribed spacer 1 and 2, after removing the 5.8S rDNA in-between) and two mitochondrial datasets-ND2 + COI (nicotinamide adenine dinucleotide dehydrogenase subunit 2 and cytochrome C oxidase subunit I) and 16S rRNA + 12S rRNA (16S rRNA and 12S rRNA fragments after removing the tRNA-Val in-between).

Sampling and available data
Sampling for the study was carried by the authors in Greece, Albania, Turkey, and the Republic of North Macedonia between 2014 and 2018. Altogether 12 taxa were sampled, all representing species of the P. jonicus and P. inflatus groups. Fig. 1 General appearance of representatives of the Poecilimon jonicus group. a, b P. jonicus lobulatus (a male, b female); c P. werneri (female); d, e P. erimanthos (d male, e female); f P. inflatus lyciae (male); g, h P. tessellatus (g male, h female); i P. laevissimus (male); j, k P. cretensis (j male, k female); l, m P. isopterus (l male, m female); n, o P. antalyaensis anemuruim (n male, o female) Besides, six outgroup taxa were included representing four distinct lineages of Poecilimon and two related genera of Barbitistini-Isophya and Leptophyes. In addition to newly obtained sequences for this study, all available published sequences from the species groups P. jonicus and P. inflatus, as well as outgroup taxa, best supplementing our samples, were accessed from GenBank (https://www.ncbi.nlm.nih.gov/ genbank/). In order to widen the taxonomic scope of the study, we also used the mitochondrial dataset (16S rRNA-tRNA-Val-12S rRNA fragment) by Ullrich et al. (2010) available from GenBank. Origin of samples, gene fragments studied, and GenBank accession numbers for the downloaded and newly obtained unique haplotypes are shown in Supplementary material A.

DNA extraction and marker amplification
Total DNA was isolated from hind leg muscle tissue, applying standard protocol for ethanol precipitation. Diluted genomic DNA was used as a template for amplification of two mitochondrial markers-a fragment from the 3' end of the ND2 and COI. For the COI fragment, primers C1-J-2183 -CAACACTT ATTTTGATTCTTTGG (forward) and TL-2-N-3014 -TCTAATGCATTAATCTGCCATCTTA (reverse) were used Fig. 2 Distribution of the members of the Poecilimon jonicus group sensu lato-ranges of all species currently considered within the groups P. jonicus and P. inflatus including subspecies of P. jonicus shown with different colors as marked in the plates with modifications to match orthopterans (Simon et al. 1994). ND2 was amplified with TM-J210 AATTAAGCTAATGG GTTCATACCC (forward) and TW-N1284 AYAGCTTT GAARGYTATTAGTTT (reverse) (Simon et al. 2006). ITS1 and ITS2 together with the 5.8S rDNA gene were amplified with primers TAGAGGAAGTAAAAGTCG (forward) and GCTTAAATTCAGCGG (reverse) resulting in an ITS1-5.8S-ITS2 fragment (Weekers et al. 2001). The PCR was performed using Supreme NZYTaq II 2x Colourless Master Mix (NZYTech, Lda. -Genes and Enzymes) following the manufacturer's instructions. Temperature cycling was held following Chobanov et al. (2017), with adaptations for hot-start PCR and slight adjustments. For the COI fragment, an initial step at 95°C was held for 5 min., followed by 35 cycles, including denaturation (94°C for 15 s), annealing (51°C for 20 s), and elongation (70°C for 90 s). Final elongation step was performed at 72°C for 15 min. ND2 was amplified with initial step at 95°C for 5 min, followed by 35 cycles of denaturation (95°C for 50 s), annealing (51°C for 40 s), and elongation (72°C for 80 s), with a final elongation step at 72°C for 15 min. For the ITS fragment, the protocol by Ullrich et al. (2010) was applied. Presence of amplification products with certain length was confirmed electrophoretically after 1 h at 80 V in 2% agarose gel containing GelRed™ gel stain (Biotium, Inc., Hayward, CA, USA). Standard sequencing reaction was performed from both 5' and 3' ends with the corresponding primers by Macrogen Europe (Macrogen, Inc., Amsterdam, Netherlands).

Sequence preparation and phylogenetic analyses
Obtained sequences were trimmed, assembled, and visually checked using CodonCode Aligner version 8.0.2 (CodonCode, Dedham, MA, USA). Sequence alignments were performed in MEGA X . Unique haplotypes were determined using DnaSP ver. 5.10 (Librado and Rozas 2009). For the protein-coding sequences, alignments were checked for stop codons to avoid numts. Due to differences in the length of sequences and/or position of primers, we were not able to use all available from GenBank COI sequences, while ND2 sequences of Poecilimon are mostly missing. Therefore, to reflect full taxonomic composition as used in the ITS phylogeny, we prepared additional mitochondrial matrix from available 16S rRNA-tRNA-Val-12S rRNA sequences (Ullrich et al. 2010).
Sequences of the 16S rRNA-tRNA-Val-12S rRNA fragment by Ullrich et al. (2010) available from GenBank were edited by removing the tRNA, and a final matrix of the 16S rRNA and 12S rRNA was prepared. The 5.8S rDNA fragment of the ITS1-5.8S-ITS2 alignments was removed, resulting in a ITS1 + ITS2 matrix, further in the text referred to as ITS. Own ND2 and COI sequences complemented with corresponding available published DNA fragments were concatenated in a ND2 + COI matrix. As a result, three matrices-ITS, COI + ND2, and 16S rRNA + 12S rRNA-were obtained. The matrices of protein-coding genes were checked for saturation, and alignments were divided by codon position (1st, 2nd, 3rd) using DAMBE 7.0.39 (Xia et al. 2003;Xia 2018). Substitution models for all codon positions of each molecular marker were estimated with jModelTest v. 2.1.7 (Darriba et al. 2012) and PartitionFinder ver. 2.1.1 (Lanfear et al. 2016). Best models were selected under the corrected Akaike information criterion (AICc). Estimated parameters and respective substitution models were implemented in subsequent phylogenetic analyses.

Estimation of divergence times
Prior to time estimation analyses in BEAST, we performed a maximum likelihood clock test in MEGA X applied to the mitochondrial matrices. The null hypothesis of an equal evolutionary rate was rejected at a 5% significance level (P < 0.05). Thus, for time estimation analyses, we applied uncorrelated lognormal relaxed clock (Drummond et al. 2006). Priors on tree topology were set as obtained from ML and BI trees produced from respective datasets, and taxa for which monophyly was strongly supported (above 90%) were marked as monophyletic.
To estimate divergence times for the ND2 + COI phylogeny, we followed two calibration strategies. First, we calibrated the BEAST analysis using a substitution rate of 0.0177 site/ Ma for COI, estimated for Aegean tenebrionid beetles and calculated by the age of the mid-Aegean trench (Papadopoulou et al. 2010), which has already been used for bush crickets (Chobanov et al. 2017). Two separate partitions-COI and ND2-were analyzed "unlinked," and the COI substitution rate was set as a prior, while the substitution rate for ND2 was estimated as "linked." Second, in order to reduce possible bias caused by substitution rate variations between lineages, we applied neutral geotectonic calibration on an internal branch of our phylogenies. The welldocumented dating of the separations of Crete from the mainland (1) in Tortonian (10.6 Ma ago; Dermitzakis and Papanikolaou 1981;Dermitzakis 1990) and (2) at the end of Messinian (5.6-5.2;Meulenkamp 1985;Dermitzakis 1990;Anastasakis et al. 2006;Loget et al. 2006) were used as two alternative hypothetical calibration points as they were thought to provide efficient calibration assuming that lineages of terrestrial animals were fully isolated on the island (Beerli et al. 1996;Lymberakis et al. 2007). Therefore, we conducted BEAST analyses for the ND2 + COI dataset setting the divergence of P. cretensis to 5.6 ± 0.4 and at 10.6 ± 0.4 Ma.
BEAST analyses were run using a Yule speciation process and MCMC chain for 100 × 10 6 generations sampling every 1 × 10 3 generations in BEAST 2.5.2 (Bouckaert et al. 2019). BEAST was run on the CIPRES Science Gateway webserver (Miller et al. 2010). Convergence to stationary and the effective sample size of model parameters were checked via Tracer ver. 1.7.1 (Rambaut et al. 2018). Maximum clade credibility trees were established with TreeAnnotator Drummond 2002-2019) to summarize the BEAST output, discarding the initial 10% of the trees as burn-in.
Alternatively, the RelTime method, implemented in Mega X, was applied on the ND2 + COI dataset to infer a timetree (Tamura et al. 2012. Tree topology was set according to our BI and ML results. One constraint at the branching off of P. cretensis was set within the range between 5.2 and 6.0 Ma, and alternatively between 10.2 and 11.0 Ma. The Tamura-Nei (Tamura and Nei 1993) substitution model with 5 gamma evolutionary categories was applied.
Molecular clock calibration for the 16S rRNA + 12S rRNA phylogeny using substitution rates would be dubious at this stage. On the other hand, dating the split of P. cretensis would have been irrelevant due to unresolved branching (see Results). Therefore, we calibrated the tree according to the results from the paleo-geotectonic dating of the ND2 + COI tree for a highly supported in both phylogenies branch, basal for the studied group.

Bioacoustics
Male acoustic signals and partly female responses are described in detail for the studied taxa (e.g., Heller 1984;Kaya et al. 2018). To exemplify specific male calling songs, we used our own recordings, performed using Pettersson D500 external microphone connected to ZOOM H2 (96 kHz, 24-bit sampling frequency) or Tascam DR-680MKII (192 kHz, 24bit sampling frequency) digital recorders and visualized with Audacity 2.1.2. (https://www.audacityteam.org/). Own data were compared with recordings available from OSF online (Cigliano et al. 2020) and Fauna d'Italia (Massa et al. 2012). Details are given in Fig. 9.

Data characteristics and genetic distances
The ITS dataset combined from own and published sequences consisted of 665 bp with gaps (630 bp without gaps), including 214 variable and 122 parsimony informative sites, and involved 16 ingroup and 6 outroup taxa (respectively 28 ingroup and 11 outroup haplotypes). The ITS dataset from Ullrich et al. (2010), involving all P. jonicus taxa sensu lato and same outgroup taxa (except for the available taxa from Leptophyes and Isophya), included 12 ingroup and 6 outgroup taxa (16 ingroup and 10 outgroup haplotypes). The latter consisted of 717 bp with gaps (679 bp without gaps), including 212 variable and 112 parsimony informative sites. The aligned concatenated mt-dataset of ND2 + COI contained 1656 bp and involved 25 unique haplotypes from 11 ingroup and six outgroup taxa. Included COI fragment had a length of 705 bp with 268 variable and 237 parsimony informative sites, and the ND2 fragment had 951 bp with 548 variable and 443 parsimony informative sites. The 16S rRNA + 12S rRNA dataset had a final length of 1835 bp, involved 25 unique haplotypes from 11 ingroup and six outgroup taxa, and included 72 indels, 666 variable sites, and 506 parsimony informative sites. No numt signs were detected in alignments of protein-coding sequences. Saturation tests applied did not show signs of significant saturation. Best substitution models for each matrix and for codon positions of the protein-coding genes are given in Supplementary material B.

Phylogenetic inferences
The BI trees obtained from the combined (Fig. 3a) and published ( Fig. 3b; Ullrich et al. 2010) ITS matrices well supported the monophyly of the group including the Balkan, Cretan, and Anatolian lineages, except for Poecilimon bilgeri that grouped with P. pergamicus Brunner von Wattenwyl, 1891, positioned closer to the basal branch of Poecilimon. The support for the monophyly of the group rose when lowering the sample size (Fig. 3b), though one haplotype described as P. j. jonicus 031 grouped with P. brunneri (Frivaldzsky, 1868) and P. zwicki Ramme, 1939. The ML analysis provided the same topology; yet, with poor support for some branches, therefore it is not shown here. The analyses failed to resolve the ingroup relationships, grouping only sequences from the same taxa, including subspecies.
The ML and BI analyses of the ND2 + COI matrix supported the group monophyly while also providing good resolution and very high support for most branches not only at a group level but also within groups (Fig. 4a). P. bilgeri, P. martinae, and P. jonicus superbus were missing from the analyses due to a lack of or very short sequences available. The ML and BI trees of the concatenated matrix placed Isophya + Leptophyes as a sister clade to Poecilimon. Within Poecilimon, the outgroup taxa formed a sister group to P. jonicus sensu lato, as follows: P. hamatus Brunner von Wattenwyl, 1878 + (P. pergamicus + (P. zwicki + P. brunneri)).
The ML tree based on the 16S rRNA + 12S rRNA dataset did not show good resolution and support (tree not shown). BI tree generally agrees with the published maximum parsimony tree (Ullrich et al. 2010) (Fig. 4b). P. bilgeri groups with P. pergamicus within the branch [(P. pergamicus + P. bilgeri) + (P. zwicki + P. brunneri)], sister to P. hamatus + P. jonicus group s.l., while the monophyly of P. jonicus group is fully supported. Anatolian P. antalyaensis (its sister taxon, P. isopterus, was not presented in the dataset) branches out first. Support for the following grouping was poor, except for the well-outlined groups of P. inflatus + P. martinae, P. erimanthos + (P. jonicus tessellatus + P. laevissimus), and P. werneri + (P. j. superbus + P. j. jonicus). The latter groups reflect the terminal grouping in the ND2 + COI tree.
Phylo-spatial dimension of the P. jonicus group based on combined data from the phylogenetic trees is summarized in Fig. 5 (mostly based on Fig. 4a; for relative branch lengths and branch support compare with Fig. 4a and b). The Anatolian taxa P. antalyaensis and P. isopterus form a clade that branches off first within P. jonicus group. The Cretan lineage (P. cretensis) branches off next (Fig. 4a) or belong to the poorly resolved crown clade (Fig. 4b). Based on the ND2 + COI phylogeny, the crown clade is formed by the Balkan lineages together with the Anatolian P. inflatus. Within the crown clade, two Balkan lineages may be outlined as follows. The Northern Balkan clade consists of P. werneri + (P. j. superbus + (P. j. jonicus + P. j. lobulatus)). The Southern Balkan clade includes the taxa occurring on the Fig. 3 Phylogenetic relationships of the members of the Poecilimon jonicus group sensu lato with four outgroups from the genus Poecilimon and two outgroups from additional genera of Barbitistini (P, Poecilimon; L., Leptophyes; I., Isophya) based on Bayesian inference analysis of ITS1+ITS2 sequences. a 665 bp alignment-tree of sequences obtained in this study and from GenBank (Ullrich et al. 2010) and b 717 bp alignment-tree of sequences available from Ullrich et al. (2010). Node support is shown at (below) resolved branches Peloponnesus Peninsula, namely, P. erimanthos + (P. jonicus tessellatus + P. laevissimus). The third group in the crown clade consists of the Anatolian taxa P. inflatus and P. martinae, whose relationships with either of the Balkan clades remain uncertain.

Discussion
Results of the present phylogenetic and time estimation analyses can be evaluated from different aspects such as the evolution of Barbitistini, Poecilimon, and P. jonicus and P. inflatus groups. In the analyses, the number of taxa from Barbitistini and Poecilimon is limited; yet, support to monophyly of the genus is the first implication to be mentioned (Ullrich et al. 2010;Chobanov et al. 2017). In addition, present data offer significant clarification about the evolutionary history of the lineages concerned and some implications concerning the paleogeographic history of the Aegean area.
Five out of seven time estimation analyses produced largely agreeing ages for the resulting chronograms. The BEAST chronograms of COI and ND2 datasets calibrated by the substitution rate of the COI (Papadopoulou et al. 2010), the RelTime chronogram of COI+ND2, and the BEAST chronograms of COI+ ND2 and 16S rRNA+12S rRNA, calibrated by the separation of Crete in post-Messinian period, resulted in very similar TMRCAs for almost all nodes (Fig. 6a, b, c). Small, though, for basal nodes, significant differences were observed between analyses calibrated for substitution rate versus those calibrated for paleogeographic events. Yet, these differences largely overlap if consider 95% HPD intervals (Supplementary material C). This is not unexpected as the rate of evolution may vary depending on lineage, time, and geography (Brower 1994;Shapiro et al. 2006;Percy et al. 2004;Papadopoulou et al. 2010;Kaya and Çıplak 2016). Received TMRCAs allow us to make robust statements about the age of the ancestor, the place of origin, and the historical factors which triggered the evolution of Barbitistini. Phylo-spatial dimension of the P. jonicus group based on the mtDNA phylogeny (data combined from COI+ND2 and 16S rRNA+12S rRNA datasets). Blue, "Anatolian" lineage; green, "Cretan" lineage; orange, "Balkan" lineage including Poecilimon inflatus Members of Poecilimon-the most speciose genus of the tribus-shared last common ancestor ca. 9.6-7.8 Ma, a period well corresponding to the first transgression of the Aegean area in Tortonian. Similarly, the last common ancestor of Isophya, the second largest lineage in the tribus, was dated at 8.3-8.8 Ma (Chobanov et al. 2017), again corresponding to the same geotectonic event. The split of Poecilimon from Isophya + Leptophyes was dated ca. 14-12 Ma (Chobanov et al. unpublished data), a time prior to Tortonian that well corresponds to current time estimations. Thus, we conclude Fig. 6 Timetrees for the P. jonicus group based on different phylogenies and time calibrations with node labels referring to the time estimation in millions of years. a BEAST maximum clade credibility tree inferred from COI and ND2 partitions calibrated for the substitution rate of COI; b BEAST chronogram and RelTime timetree inferred from the COI+ND2 phylogeny calibrated for the last separation of Crete (values above nodes, BEAST chronogram estimations; values below nodes, RelTime timetree estimations); c BEAST chronogram inferred from the 16S rRNA+12S rRNA phylogeny calibrated for the TMRCA (interval) of P. jonicus group (node 12) as suggested from the COI+ND2 chronogram. Different topologies of the trees are marked with "*" for the COI partitions and "**"for the ND2 partition on (a) and different topology of the BEAST and RelTime trees on (b). Numbers in white circles correspond to data in Supplementary material C. Light blue vertical band marks the onset of the northern hemisphere glaciation and dark blue lines mark major climatic switches in the Pleistocene. Uniform gray areas below the trees show united Aegean landmass (left) and reunification of large areas in Messinian (right), while striped areas show disintegrated state of the Aegean that Barbitistini radiated from a common ancestor in Serravallian/Langhian and the marine transgression of the Aegean area in Tortonian (Fig. 7b) is the main event that triggered early lineage divergence. At present, members of the tribus are mainly ranged in derivatives of the old Aegean Plate. The latter inferences support Çıplak (2004) in his assumption that the tribus radiated from an Aegean lineage in Langhian when the Aegean land represented united subcontinent in the Tethys Sea (Fig. 7a). The transgression of the Aegean as an evolution driving factor has been observed in several other lineages (Douris et al. 2007;Çıplak et al. 2010;Papadopoulou et al. 2010;Parmakelis et al. 2013). On the other hand, the contemporary ranges of both Poecilimon and Isophya largely overlap; yet, Isophya is not represented in the southern and southwestern part of the Aegean (most of mainland Greece, Crete and the majority of the Aegean islands), while Poecilimon occurs with a number of taxa all over the Aegean area. Thus, as a working hypothesis to be tested, we assume that Isophya represents the main northern and Poecilimon represents the main southern lineage of Barbitistini (see also Chobanov et al. 2017).
Tribus Barbitistini was claimed to be an impressive example of diversification within a limited spatial and temporal scale (Kaya et al. 2015;Grzywacz et al. 2018). The group  Çıplak 2004;Chobanov et al. 2017) and is currently known with ca. 300 micropterous species concentrated in the Aegean and Pontic area (Cigliano et al. 2020). Chobanov et al. (2017) argued that genus Isophya evolved after the middle Miocene climatic transition from an "ephemeral" ancestor that developed fast within a short favorable season with appropriate humidity and temperatures. The loss of flight is beneficent to organisms that need to reduce energy costs in an ecologically harsh environment and is typical for insects of the dryer, cooler, and windier climates in comparison with tropical and temperate zones (e.g., Roff 1990;Grzywacz et al. 2018). Flightlessness on the other hand can promote diversification not only by contributing to isolation but also through acceleration of the mutation rates (Ikeda et al. 2012). Thus, a combination of geological (disintegration of the Aegean area) and climatic factors (the Middle Miocene Climatic Transition) may have initiated and accelerated the early evolution of Barbitistini.

Evolution of Poecilimon jonicus group
According to mitochondrial chronograms (Fig. 6a-c), P. jonicus group sensu lato shares a common ancestor 8-6 Ma ago, a time well corresponding to the Messinian period of the Miocene. As the Anatolian clade (P. antalyaensis + P. isopterus) constitutes the basal branch and the Cretan (P. cretensis) the subsequent one, Anatolian mainland seems to be the place of origin of the group ancestors. Hence, it can be assumed that desiccation in the Mediterranean during the Messinian salinity crisis (MSC) (5.96-5.33 Ma), especially in the Aegean area (e.g., Steininger and Rögl 1984;Dermitzakis 1990), provided distribution corridor for dispersal of an ancestral stock over the present distribution range of the group. However, the poorly resolved relationships between the Cretan, Balkan, and P. inflatus + P. martinae lineage suggest an early radiation within Anatolia or along the Aegean landmass connecting the Balkans and Anatolia during the MSC (Fig. 7c). The following post-Messinian transgression of the Mediterranean and Aegean area ca. 5.3-5.2 Ma (Jolivet et al. 2006) acted as a vicariant event resulting in four geographically isolated stocks-ancestral P. inflatus + P. martinae in Anatolia and P. cretensis in Crete, the Northern, and Southern Balkan group in the southern Balkan mainland and on the Peloponnesos, respectively (Fig. 8a). In that case, the expected phylogeny of the group will be [(P. antalyaensis + P. isopterus) + (P. cretensis + (P. inflatus + P. martinae) + Northern Balkan clade + Southern Balkan clade)].
During Messinian, the level of the Mediterranean dropped after its isolation from the Atlantic. As a result of the MSC, the Aegean landmass was largely united again, and terrestrial migrations between the Balkan mainland and Crete were possible and documented (e.g., Dermitzakis 1990;Lymberakis et al. 2007). With the end of the Messinian 5.33 Ma ago and the Zanclean refill of the Mediterranean basin the Balkans, Anatolia and Crete were disconnected again ca. 5.2 Ma ago (Meulenkamp 1985;Dermitzakis 1990). We believe these two periods are of major importance for dispersal, early speciation, and subsequent vicariance leading to major lineage splits within the Poecilimon jonicus group (Figs. 7c, 8a). Within this time frame (5.8/5.6 to 5.2 Ma), the Cretan lineage was established and shortly after that the Balkan lineage splits into a Northern and a Southern group. Roughly at the same time, the P. inflatus branch (present P. inflatus + P. martinae) was established.
In Pliocene, the land-sea configuration further complicated with the development of many islands and peninsulas, as well as vast inland freshwater areas (marshes and lakes) (Popov et al. 2004) (Fig. 8a, b). The latter may have contributed to the basal lineage splits within the Northern Balkan group (P. werneri, P. jonicus superbus), Southern Balkan group (P. erimanthos), and in the Anatolian group (P. antalyaensis and P. isopterus) due to vicariant events. Specifically, the lack of P. jonicus superbus in most of Puglia (southeastern corner of Italy) (Fig. 2), though this area represented an island at the time this taxon originated, may possibly be explained as a result of recent climatic factors.
Latest speciation events are dated to the Pleistocene, starting with the onset of the northern hemisphere glaciation (2.5-2.6 Ma ago), when climate cycles with alternating cooler drier and warmer more humid periods dominated the area. Those climate events have been suggested as triggers for recent lineage splits in other Barbitistini (Kaya et al. 2015;Chobanov et al. 2017). Pleistocene climate was dominated by two major climatic switches-first ca. 1.4 Ma, when 41-Ka cold-warm cycles become constant, and second, 0.8 Ma, when 100-Ka cycles established (Lisiecki and Raymo 2007). Those periods of cooling and warming were connected not only with climate changes but also with sea level drops (glacials) and rises (interglacials). Thus, thermophilous animals dependent on modest humidity retreated to refuges and were isolated on islands during stadials (cold dry periods) and regained territories during interstadials (warm humid periods). Interesting correlation may be observed between the intraspecific divergence of the Cretan lineage (P. cretensis) and the second Pleistocene climatic switch (Fig. 6b). Patterns of the evolutionary history of Cretan lineages based on early vicariance are documented in various groups (Parmakelis et al. 2005(Parmakelis et al. , 2006Simaiakis and Mylonas 2008). Populations on the island were completely isolated after the end of the Messinian. Yet, the current appearance of Crete as a united landmass was established only recently, while many smaller island configurations existed during the Pliocene (Creutzburg 1963;Dermitzakis 1990;Douris et al. 1998;Welter-Schultes 2000) (Fig. 8a, b). Relatively high genetic distances between sampled populations of P. cretensis could be a reflection of vicariant events during the Pleistocene as a result of isolation of the easternmost Cretan territory by lakes/brackish waters around that time.

Implications for the systematics of Poecilimon jonicus group
Clear implications about the systematics of P. jonicus group sensu lato may be outlined from the present study. Both nuclear and mitochondrial DNA sequences support a clear monophyletic outline of the group including the present P. jonicus and P. inflatus species groups and excluding P. bilgeri. Although the group is monophyletic, its basal internal phylogeny is not doubtlessly resolved. Though both mitochondrial trees agree for the ancestral position of the Anatolian P. antalyaensis (+ P. isopterus), the Cretan, Northern Balkan, Southern Balkan, and the Anatolian P. inflatus + P. martinae clades have unresolved relationships. The phylogenetic position of P. inflatus is especially intriguing, as it appears as an internal branch of the Balkan lineages with very high support in the ND2 + COI phylogeny, at the same time occurring at the opposite side of the Aegean Sea. Yet, based on the long history of the discussed lineages and contradictions in the trees, independent evolution of each genome is possible (Bernardo et al. 2019). The instability at the base of the P. jonicus group tree may express a real situation. If a contemporary split of the Balkan, Cretan, and Anatolian lineages by the end of Messinian is assumed, then a basal polytomy or basal instability will be expected.
All taxa currently recognized in the group seem to represent distinct phylogenetic units, and the position of certain lineages allows for accepting a higher number of species. For instance, P. tessellatus groups with P. laevissimus and is clearly distinct from P. jonicus, with which it was recently united; yet, P. jonicus jonicus and P. j. lobulatus seem representing distinct units with a subspecies rank. Poecilimon jonicus superbus, occurring in the Apennines, shows large genetic distances from P. j. jonicus and P. jonicus lobulatus, exceeding those between P. inflatus and P. martinae and comparable with those between the well morpho-acoustically distinguished P. antalyaensis and P. isopterus. The divergence between haplotypes of P. cretensis is prominent, and its genetic subdivision extends those of P. jonicus and P. antalyaensis (for P. a. antalyaensis and P. a. myrae). On the other hand, genetic distances between subspecies of P. inflatus and P. antalyaensis (subspecies of P. martinae not included in the study) (see Kaya et al. 2018) support their status.
Actual taxonomic changes arising from our results are the following. Poecilimon jonicus tessellatus was initially described as a distinct species, and later, it was synonymized with P. jonicus and obtained subspecies status (Heller 1988). Although similarities with P. jonicus are obvious, the molecular phylogeny of the group points out that P. j. tessellatus is part of the Southern Balkan group and has not common evolutionary history with P. jonicus being relative to P. laevissimus. Divergence between P. j. jonicus and P. j. tessellatus is also observed in the distinct song patterns (compare Fig. 9). While our data suggests that P. j. tessellatus characterizes with a stable complex song pattern with an isolated short main syllable and two complex after parts ( Fig. 9a) (compare the stridulatory movement by Heller 1988 -Abb. 36 H), possibly rarely involving only the main part (Heller 1988 -Abb. 36 G), own and published data for other subspecies of P. jonicus suggest mostly a simple song pattern with a long main syllable (compare the stridulatory movement by Heller 1988 -Abb. 36 E, F). Yet, single after parts are present in some recordings of P. jonicus (Fig.  9e, g). On the other hand, though the song of P. j. superbus (also described as a distinct species) is almost identical to that of P. j. jonicus, we revealed large genetic distances between them, suggesting 3.1 to 4.7 Myr of isolation (Figs. 6c,8a). Considering the fact that closely related allopatric bush crickets usually retain a conservative song pattern (Heller 2006) and that the genetic differentiation of P. j. superbus is comparable with that of well-distinguished sympatric taxa within the group, we may conclude for its species status.
Concluding from the above discussion, we propose the following systematic and nomenclatural changes. Poecilimon bilgeri, currently related to P. inflatus (Kaya et al. 2018), does not belong to the Poecilimon jonicus group sensu lato and is relative to Poecilimon pergamicus (this study and own unpublished data). Poecilimon jonicus group represents a monophyletic lineage, whose basal clades diversified within a short time frame and thus (at the present state of knowledge) cannot be divided to less than four or five taxonomic complexes; therefore, the existence of Poecilimon inflatus group (Kaya et al. 2018) is not justified. Hence, the current assemblage of the Poecilimon jonicus group (stat. rev.) is as follows (in alphabetical order): Fig. 9 Examples of the male calling song presented in a time frame of 500 ms with oscillograms (a, c, e, f, g, h) and spectrograms (b, d) of taxa currently considered subspecies of Poecilimon jonicus. a P. j. tessellatus, Greece, Peloponnesos, Kallithea, 22.4°C, daily recording, own data; b same; c P. j. lobulatus, mainland Greece, Amphilochia, 21°C, daily recording, own data; d same; e P. j. jonicus, Albania, Kolonje district, Qafa e Qarrit pass, 1.5 km W of Pepellash, 1200 m a.s.l, lab recording, 30°C, ZOOM-H4 digital recorder, sample rate 96 kHz, G. Puskas, source: OSF online; f P. j. jonicus, Albania, S of Vlorë, N of Qeparo vill., ca. 25°C, daily recording, own data; g P. j. superbus, Italy, Liguria, Orsomarso, Valle Fiume Argentino, 25°C, source: Massa et al. 2012; and h same, different part of the recording Poecilimon antalyaensis anemurium Kaya, Chobanov & Funding This study was supported by the National Science Fund (MES) of Bulgariaproject DN11/14-18.12.2017 to Dragan Chobanov and by the Program for Supporting Young Scientists of the Bulgarian Academy of Sciences -Project 17-86/27.08.2017 to Simeon Borissov. The material was collected with the support by The Teodore J. Cohn Research Fund (OSF), Deutsche Gesselschaft für Orthopterologie (grants to Simeon Borissov), and a joint project 'Phylogenetic relationships in the group Barbitistini (Orthoptera, Tettigoniidae, Phaneropterinae) based on genetic studies' of D. Chobanov and B. Çıplak. The Scientific and Technological Research Council of Turkey (TÜBITAK) (postdoctoral research grant to Dragan Chobanov) supported the acquisition of sequences from Anatolian specimens.
Data availability The datasets generated and/or analyzed during the current study will be freely available after July 1, 2020. Edited and aligned sequences from this study are referred to GenBank (https://www.ncbi. nlm.nih.gov/genbank/) accession numbers given in Supplementary material A. Phylogenetic and timetrees will be available from DRYAD (https://datadryad.org/stash). All published and unpublished data are also available from the corresponding author on reasonable request.
Open Access This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if changes were made. The images or other third party material in this article are included in the article's Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article's Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this licence, visit http://creativecommons.org/licenses/by/4.0/.