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 a connection with the Euxinian (Pontian) Basin (Meulenkamp 1977; Dermitzakis 1990; Steininger and Rögl 1984; Jolivet et al. 2006; Popov et al. 2004, 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. 2005a, 2005b; Douris 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).

Our study focuses on a group of sibling bush cricket species with a disjunct distribution between the Apennines, the western and southern Balkans, Crete, and southwestern Anatolia that form a poorly resolved monophyletic clade according to Ullrich et al. (2010 – 11 taxa studied). This group was first arranged under Poecilimon jonicus group based on considerable similarities in morphology and song structure by Heller (1984, 1988, 2004), while later, some of the taxa were transferred to a group of their own, P. inflatus group (Kaya et al. 2018). The external appearance of some of the studied taxa is presented in Fig. 1. The following 17 taxa (10 species), which Sevgili et al. (2018) defined as superspecies group, are considered in our study (ranges of species and the subspecies of P. jonicus are shown in Fig. 2):

  • Western and southern Balkans and Apennines, P. jonicus group sensu stricto: (1) P. jonicus superbus (Fischer, 1853), mainland Italy; (2) P. j. jonicus (Fieber, 1853), northwestern Greece, North Macedonia, and Albania, reaching as north as Montenegro; (3) P. j. lobulatus Willemse, 1982, restricted to a small area in western Greece; (4) P. j. tessellatus (Fischer, 1853), the Peloponnesus peninsula; (5) P. werneri Ramme, 1933, a small coastal area in Albania, southwestern edge of the Greek mainland, northwestern Peloponnesus; (6) P. laevissimus (Fischer, 1853), the Peloponnesus and Italy (Sicily), though the latter occurrence was suggested to be a result of human introduction (K.-G. Heller, pers. comm. in Hochkirch et al. 2016); and (7) P. erimanthos Willemse & Heller, 1992, endemic to the Erymanthos mountain and its surroundings in the Peloponnesus.

  • Crete, P. cretensis Werner, 1903—endemic to the island of Crete. The species was earlier associated with the P. jonicus group (Heller 1988), while a recent study points out morphological similarities with the P. inflatus species group, distributed in southwestern Anatolia (Kaya et al. 2018).

  • Anatolia, Poecilimon inflatus group (as defined by Kaya et al. 2018): (1) Poecilimon inflatus Brunner von Wattenwyl, 1891 with two subspecies—P. inflatus inflatus and P. inflatus lyciae Kaya & Çıplak, 2018—has a narrow range in southwestern Anatolia; (2) P. martinae Heller 2004 with two subspecies—P. m. martinae and P. m. tlos Heller, 2004—occurs in Antalya and Muğla provinces (SW Anatolia); (3) Poecilimon antalyaensis (Karabag, 1975) with three subspecies—P. a. antalyaensis (Karabag, 1975), P. a. myrae Chobanov & Heller, 2018, and P. a. anemurium Kaya, Chobanov & Çıplak, 2018—distributed along the western Mediterranean coast of Anatolia; (4) Poecilimon isopterus Kaya & Chobanov, 2018, with a small range in Antalya region; and (5) Poecilimon bilgeri Karabag, 1953, extreme southwestern Anatolia.

Fig. 1
figure 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)

Full size image
Fig. 2
figure 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

Full size image

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

Material and methods

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. 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 – CAACACTTATTTTGATTCTTTGG (forward) and TL-2-N-3014 – TCTAATGCATTAATCTGCCATCTTA (reverse) were used with modifications to match orthopterans (Simon et al. 1994). ND2 was amplified with TM-J210 AATTAAGCTAATGGGTTCATACCC (forward) and TW-N1284 AYAGCTTTGAARGYTATTAGTTT (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 (Kumar et al. 2018). 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.

Maximum likelihood (ML) analyses were conducted using RAxML (Stamatakis 2006) on the T-rex online platform (http://www.trex.uqam.ca/index.php?action=raxml&project=trex) (Boc et al. 2012). Support of nodes was acquired with 1000 bootstrap resampling. Bayesian inference (BI) analyses were accomplished in Mr. Bayes version 3.2.5 (Ronquist and Huelsenbeck 2003; Ronquist et al. 2005) with four simulations of Markov chains and 4 × 106 generations sampling each 100th tree. Stationary distribution of the MCMC parameters was checked with Tracer ver. 1.7.1 (Rambaut et al. 2018). The first 25% of trees were excluded as burnin. Results were visualized in FigTree 1.4.3 (http://tree.bio.ed.ac.uk/software/figtree/).

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 well-documented 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 × 106 generations sampling every 1 × 103 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 (Rambaut and 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, 2018). 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, 24-bit 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.

Results

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.

Fig. 3
figure 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

Full size image

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

Fig. 4
figure 4

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 (P, Poecilimon; L., Leptophyes; I., Isophya). a Combined maximum likelihood and Bayesian inference phylogenetic trees based on a 1659 bp alignment of the COI+ND2 mitochondrial dataset and b Bayesian inference phylogenetic tree based on a 1835 bp alignment of the 16S rRNA+12S rRNA mitochondrial dataset (Ullrich et al. 2010). Node values at branches show node support: above, ML bootstrap proportions; below, BI posterior probabilities; asterisk (*) indicates branches resolved by one of the analyses

Full size image

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

Fig. 5
figure 5

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

Full size image

Estimation of divergence times

The BEAST chronograms of the “unlinked” datasets of COI and ND2 calibrated by 0.0177 substitution/site/Ma (Fig. 6a) suggested that P. jonicus group sensu lato shares the last common ancestor 7.90/7.35 (COI/ND2) Ma ago. The TMRCAs for the internal clades in P. jonicus group sensu lato by the same chronograms (COI/ND2) are as follows: (i) 5.57/4.00 Ma for P. antalyaensis + P. isopterus, (ii) 7.33/6.66 Ma for Cretan + Balkan clades, (iii) 5.89/* Ma (COI) for P. inflatus + all Balkan species (COI chronogram only), (iv) 6.63/4.47 Ma for the three taxa in Peloponnesus, (v) */6.13 Ma for P. inflatus + (P. werneri + P. jonicus s.str.) clade (ND2 chronogram only), and (vi) 4.96/5.11 Ma for P. werneri + P. jonicus (except P. j. tessellatus). The chronograms of COI/ND2 estimated the TMRCA of Barbitistini as 10.89/11.19 and that of the members of Poecilimon as 9.55/9.63 Ma.

Fig. 6
figure 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

Full size image

The BEAST chronogram and RelTime timetree of the ND2+COI calibrated by the separation of Crete in Tortonian as 10.6 Ma produced significantly extreme TMRCAs contradicting the rest of our results and published analyses for Poecilimon and Isophya (compare Chobanov et al. 2017; for additional comments see Discussion); thus, we do not consider them and those are not shown here. The second chronogram of ND2 + COI, calibrated by separation of the Cretan lineage ca. 6–5.2 Ma ago, produced ancestral times largely congruent with the RelTime chronogram and rate calibrated analyses. The BEAST and RelTime (RTC) chronograms for ND2 + COI calibrated by the Cretan lineage split as 5.6 ± 0.4 Ma agreed well in node ages (combined results are shown in Fig. 6b). The TMRCAs (BEAST/RTC) are as follow: (i) Barbitistini 10.24/* Ma, (ii) Leptophyes + Poecilimon 9.12/*, (iii) Poecilimon clade 7.81/7.79 Ma, (iv) P. jonicus group sensu lato 6.05/6.13 Ma, (v) P. antalyaensis + P. isopterus 3.52/3.78 Ma, (vi) Cretan + Balkan clades 5.57/5.60 Ma (calibration node), (vii) Anatolian P. inflatus + mainland Balkan species 5.20/5.11 Ma, (viii) P. inflatus + (P. werneri + P. jonicus) clade 4.93/4.91 Ma, (ix) the clade including three species in Peloponnesus 4.44/4.47 Ma, and (x) P. werneri + P. jonicus 4.3/4.26 Ma, (xi) P. laevissimus + P. tessellatus ca. 1.73/1.71 Ma. Times to terminal nodes corresponding to intraspecific lineage splits were dated later than 0.9 Ma.

On the 16S rRNA + 12S rRNA tree BEAST calculated the TMRCAs of nodes as follow (Fig. 6c): (i) Barbitistini 9.34 Ma, (ii) Poecilimon clade 7.91 Ma, (iii) P. jonicus group sensu lato (calibration node) 6.24 Ma, (iv) P. jonicus group after the split of P. antalyaensis 5.81 Ma, (v) Anatolian P. inflatus + P. martinae + Cretan + Balkan species 5.81 Ma, (vi) P. cretensis + (P. inflatus + P. martinae) 5.27 Ma, (vii) the Balkan clade 5.30 Ma, (viii) the Southern Balkan clade 4.35 Ma, (ix) the northern Balkan clade 4.58 Ma, (x) the divergence of Italian P. jonicus superbus 3.87 Ma, (xi) P. inflatus + P. martinae 2.39 Ma, and (xii) P. laevissimus + P. tessellatus 2.06 Ma. Times to terminal nodes corresponding to intraspecific lineage splits were dated later than 0.97 Ma.

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.

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

Fig. 7
figure 7

Land-sea configuration in the Aegean area in different periods of the Neogene based on reconstructions by Popov et al. (2004) (with details taken from Creutzburg 1963 and Dermitzakis 1990) with historical biogeography hypothesis for the Messinian. a Late Serravallian (12–11 Ma), a possible time for the origin and early divergence of Barbitistini; b late Tortonian (8.5–7 Ma), transgression of the Mid-Aegean area in Tortonian corresponds to TMRCA of Poecilimon and Isophya; and c late Messinian (6–5.6 Ma), the ancestor of Poecilimon jonicus group emerged on the Anatolian side of Aegean land in Messinian; the regression of Mid Aegean area during the MSC provided corridors for dispersal to Crete and the Balkans

Full size image

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 was thought to have recently evolved (since Miocene; Çı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)].

Fig. 8
figure 8

Land-sea configuration in the Aegean area in different periods of the Neogene based on reconstructions by Popov et al. (2004) (with details taken from Creutzburg 1963 and Dermitzakis 1990) with historical biogeography hypotheses for the evolution of the P. jonicus group. a Mid-late Pliocene (ca. 3.5 Ma), transgression of the Aegean area at the end of Messinian/beginning of Pliocene led to vicariant speciations in Anatolia, Crete, and the Balkans and b late Pliocene–Mid Pleistocene (ca. 3.4–0.4 Ma), fluctuating levels of the Mediterranean Sea, small-scale geotectonic events and climatic cycles leading to frequent transformations in the land-sea configuration contributed to the recent history of Poecilimon jonicus group

Full size image

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, 2006; Simaiakis 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.

Fig. 9
figure 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

Full size image

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

  • Poecilimon antalyaensis anemurium Kaya, Chobanov & Çıplak, 2018

  • P. antalyaensis antalyaensis (Karabag, 1975)

  • P. antalyaensis myrae Chobanov & Heller, 2018

  • P. cretensis Werner, 1903

  • P. erimanthos Willemse & Heller, 1992

  • P. inflatus inflatus Brunner von Wattenwyl, 1891

  • P. inflatus lyciae Kaya & Çıplak, 2018

  • P. isopterus Kaya & Chobanov, 2018

  • P. jonicus jonicus (Fieber, 1853)

  • P. jonicus lobulatus Willemse, 1982

  • P. laevissimus (Fischer, 1853)

  • P. martinae martinae Heller, 2004

  • P. martinae tlos Heller, 2004

  • P. superbus (Fischer, 1853), stat. rev.

  • P. tessellatus (Fischer, 1853), stat. rev.

  • P. werneri Ramme, 1933

Conclusion

The dynamic geological history of the Aegean contributed to complex evolutionary patterns in terrestrial animals. The flightlessness in bush crickets was an additional trigger to isolation and thus speciation. Specializing in allopatry, the lineages did not require prompt changes in morphology, and sexual selection could have induced slow change. Thus, comparatively small to the eye morphological and behavioral differences corresponded to long-term speciation in isolation. Our molecular phylogenetic analyses provided a different viewpoint to the problem and suggested robust phylogeny of the Poecilimon jonicus group in a broad sense in parallel revealing a strong link between its evolution and the geomorphological history of the Aegean.