DOI: 10.5061/dryad.j9kd51cgp
title: Origin, evolution and systematics of the genus Poecilimon (Orthoptera: Tettigoniidae) – an outburst of diversification in the Aegean area

description: 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.
description: 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).
 <strong>Dataset preparation and phylogenetic
 analysis</strong> 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). <strong>Molecular clock calibration and
 divergence times</strong> 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
 <em>Poecilimon</em> found on the island of Crete is the
 endemic species <em>P</em>. <em>cretensis</em>
 Werner, which allows the assumption that the split of
 <em>P</em>. <em>cretensis</em> from its relatives
 in the <em>P</em>. <em>jonicus</em>-group was
 triggered by the isolation of Crete (see Borissov et al., 2020). Similar
 assumption is reasonable for the split of <em>P</em>.
 <em>deplanatus</em> Brunner von Wattenwyl, endemic to the
 island of Karpathos, from its relatives of the subgenus
 <em>Hamatopoecilimon</em> 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.