Our study focuses on the origin, dispersal patterns, evolutionary strategies and systematics of Poecilimon, the largest bush-cricket genus in the Palearctic with over 150 taxa described. We employ phylogenetic and divergence time estimation analyses based on multilocus sequence data (ND2+COI+12S+16S+ITS+28S), perform ancestral area reconstruction, and track the evolution of behavioural (evolution of sound communication) and morpho-physiological traits (body size and shape, and spermatophore size) in this genus. Based on our results, we propose a revised systematics of Poecilimon, including description of a new species, P. nivalis sp. nov., and hypothesise three stages in the evolution of Poecilimon. (1) In the early evolution of the genus in Tortonian, when open dry habitats appeared in the Eastern Mediterranean, diversification rates were low and speciation was possibly induced by vicariance and habitat fragmentation; physiology and morphology during this period retained their ancestral states but the evolution of main lineages may have been accompanied by behavioural specialisations. (2) Climate cooling and aridification during the Messinian induced dispersals and adaptation to new habitats, followed by physiological and behavioural adaptations; major clades formed or started diversifying. (3) Starting at the end of Messinian and continuing through the Plio- and Pleistocene, a few dispersal events from Anatolia to the Balkans took place and climatic oscillations were followed by allo- and parapatric divergence of habitat specialists, while ecological adaptations enhanced song diversity and led to morpho-physiological changes.

Material used in this study was sampled in the period 2014–2020 in the southern Balkans and Anatolia (including Albania, Bulgaria, Greece, Republic of North Macedonia, and Turkey) as shown in Table S1. Tissue samples (usually a hind leg) for DNA extraction were preserved in absolute ethanol. Total DNA was isolated from hind leg muscle, applying a standard protocol for ethanol precipitation (Aljanabi & Martinez, 1997) or using the Qiagen DNeasy® Blood & Tissue kit following the manufacturer’s instructions. Two fragments from mitochondrial and one from nuclear DNA were amplified as follows. The mitochondrially encoded NADH dehydrogenase subunit 2 (ND2) was amplified with standard primers, modified for orthopterans – forward TM-J210 (AATTAAGCTAATGGGTTCATACCC) and reverse TW-N1284 (AYAGCTTTGAARGYTATTAGTTT) (Simon et al., 2006). Partial cytochrome oxidase subunit I (COI) was amplified with the primer pair – forward C1-J-1751 (GGAGCTCCAGATATAGCTTTTCC) and reverse TL2-N-3014 (TCTAATGCATTAATCTGCCATCTTA) (Simon et al., 1994, 2006). The internal transcribed spacers 1 and 2 (ITS1 and ITS2) were partially amplified as a single fragment together with the 5.8S rDNA gene in-between using the primer pair forward TAGAGGAAGTAAAAGTCG and reverse GCTTAAATTCAGCGG (Weekers et al., 2001). Polymerase chain reactions were carried out in 25 µl volume using Thermo Scientific DreamTaq Hot Start Master Mix following the manufacturer’s protocol. Temperature cycling followed Chobanov et al. (2017) for the mitochondrial markers and Weekers et al. (2001) for the ITS1+5.8S+ITS2 fragment. Sequencing was carried out by Macrogen Europe B.V. (Amsterdam, The Netherlands). Chromatograms were processed and assembled with CodonCode Aligner version 8.0.2 (CodonCode, Dedham, MA, USA).

Dataset preparation and phylogenetic analysis

A total of 64 sequences of ND2, COI and ITS1+5.8S+ITS2 were obtained in the present study. In addition to the newly obtained data, DNA sequences of the latter fragments were available from own previously published works (Chobanov et al., 2017; Borissov & Chobanov, 2020; Borissov et al., 2020, 2021), while additional sequences for the corresponding taxa from the nuclear 28S rRNA gene and a mitochondrial fragment of 12S rRNA+tRNA-Val+16S RNA (e.g., Ullrich et al., 2010) were downloaded from GenBank. In order to cover more taxa, the entire nuclear dataset (ITS1+5.8S+ITS2) from Ullrich et al. (2010) was accessed from GenBank, compiled with other published and new data, and processed separately as described above. An overview of all DNA sequences used in this study is given in Table S1 with corresponding GenBank accession numbers. Unique haplotypes were selected using DAMBE ver. 7.0.39 (Xia, 2018). Multiple sequence alignments were performed in MEGA X (Kumar et al., 2018). All protein-coding sequences were checked for numt possibility (Kaya & Çıplak, 2018). The 5.8S rRNA gene in-between the ITSl and ITS2 (hereafter referred to as ITS) and the tRNA-Val sequence from the mitochondrial rRNA (hereafter 12S+16S) were removed from the final datasets through multiple alignment.

The phylogenetic analyses described below were performed on the ITS dataset and on a single concatenated matrix (ND2+COI+12S+16S+ITS+28S). Bayesian inference (BI) phylogenetic analyses were conducted using MrBayes v. 3.2.5. (Ronquist et al., 2012). Best substitution models were selected using PartitionFinder ver. 2.1.1 (Lanfear et al., 2016) under the corrected Akaike criterion. Partitions referred to codon positions for protein-coding fragments and a single partition for each of the other datasets (see Table S2). Four simulations of Markov chains with 4x106 generations were run sampling every 100 trees. Chain parameters were examined with Tracer v. 1.7.1 (Rambaut et al., 2018) to confirm stationarity and convergence. The first 25% of sampled trees were excluded as burn-in. Maximum likelihood analyses (ML) were performed with IQ-TREE (Nguyen et al., 2015) through the W-IQ-TREE online interface (Trifonopoulos et al., 2016). Substitution models were inferred by the incorporated ModelFinder (Kalyaanamoorthy et al., 2017), using partitions (Chernomor et al., 2016). Ultrafast bootstrap with 1000 replicates was used to obtain branch support (Hoang et al., 2017) together with two single branch tests – SH-aLRT with 1000 replicates (Guindon et al., 2010) and approximate Bayes (Anisimova et al., 2011).

Molecular clock calibration and divergence times

Biogeographic calibrations of the molecular clock are gaining popularity and provide alternatives to fossil calibrations, especially for organismal groups and geographical regions with insufficient fossil record (De Baets et al., 2016). The active geological past of the Eastern Mediterranean strongly affected lineage divergence and several geological events have been used to calibrate molecular phylogenies (Lymberakis et al., 2007; Papadopoulou et al., 2010; Trichas et al., 2020). One of the most widely accepted events is the reopening of the Gibraltar strait 5.33 Ma that caused a massive inflow of Atlantic waters (Garcia-Castellanos et al., 2009), thus rapidly submerging large parts of the Mediterranean. Consequently, in the Eastern Mediterranean, Crete lost its connection to the mainland and was never reconnected (Meulenkamp, 1985; Dermitzakis, 1990). The only representative of Poecilimon found on the island of Crete is the endemic species P. cretensis Werner, which allows the assumption that the split of P. cretensis from its relatives in the P. jonicus-group was triggered by the isolation of Crete (see Borissov et al., 2020). Similar assumption is reasonable for the split of P. deplanatus Brunner von Wattenwyl, endemic to the island of Karpathos, from its relatives of the subgenus Hamatopoecilimon Heller (Borissov & Chobanov, 2020). Karpathos was connected to Anatolia through Rhodes until the late Pleistocene (3–3.5 Ma) and was not reconnected to a major landmass ever since (Daams & Van Der Weerd, 1980). We used the latter isolation events as calibration points to set prior information in BEAUti for the BEAST v. 2.6.6 analysis (Bouckaert et al., 2019). Normal distribution was used for both calibration nodes with mean of (1) 5.6 Ma for Crete (SD=0.2) and (2) 3.3 Ma for Karpathos (SD=0.1). An alternative scenario of divergence as a result of the earlier documented insulation of Crete ca. 10 Ma (Dermitzakis and Papanikolaou, 1981; Dermitzakis, 1990) was ruled out since it produces significantly old divergence times that contradict independently estimated ages (Chobanov et al., 2017; Kaya et al., 2015; Çıplak et al., in press). Divergence times were estimated from the concatenated dataset partitioned as follows – ND2 (positions 1–3), COI (positions 1–3), ITS1+2, 28S, 16S+12S (Table S2). A lognormal relaxed clock was applied (Drummond et al., 2006). The BI tree produced in this study was used as a prior. BEAST was run under Yule speciation process using 1x108 generations and sampling every 1000 trees. Effective sample size and stationarity were assessed with Tracer v. 1.7.1 (Rambaut et al., 2018). Results were summarised with TreeAnnotator (Rambaut & Drummond, 2002–2019) discarding 10% of the trees as burn-in.