Tempo and mode of species diversification in Dolichopoda cave crickets (Orthoptera, Rhaphidophoridae)

Molecular Phylogenetics and Evolution 60 (2011) 108–121

Contents lists available at ScienceDirect

Molecular Phylogenetics and Evolution journal homepage: www.elsevier.com/locate/ympev

Tempo and mode of species diversification in Dolichopoda cave crickets (Orthoptera, Rhaphidophoridae) Giuliana Allegrucci a,⇑, Emiliano Trucchi a,b, Valerio Sbordoni a a b

Dept. of Biology, University of Rome Tor Vergata, Via della ricerca scientifica snc, 00173 Rome, Italy Centre for Ecological and Evolutionary Synthesis (CEES), Dept. of Biology, University of Oslo, P.O. Box 1066 Blindern, N-0316 Oslo, Norway

a r t i c l e

i n f o

Article history: Received 20 October 2010 Revised 28 March 2011 Accepted 7 April 2011 Available online 14 April 2011 Keywords: Lineage diversification Molecular rates Biogeography ABC analysis Cave crickets Dolichopoda

a b s t r a c t This study focuses on the phylogenetic relationships among ninety percent of known Dolichopoda species (44 out of 49); primarily a Mediterranean genus, distributed from eastern Pyrenees to Caucasus. A total of 2490 base pairs were sequenced corresponding to partial sequences of one nuclear (28SrRNA) and three mitochondrial genes (12S, 16S and COI). A relaxed molecular clock, inferred from Bayesian analysis was applied to estimate the divergence times between the lineages using well dated palaeoevents of the study areas. Molecular substitution rates per lineage per million years were also obtained for each analyzed gene. Based on the nearly complete species phylogeny, temporal patterns of diversification were analyzed using Lineage-Through-Time plots and diversification statistics. Alternative hypotheses about the colonization of present range by Dolichopoda species were tested by means of Approximate Bayesian Computation analysis. Results from this analysis carried out on the 90% of known Dolichopoda species confirmed the previous ones based on subgroups of species, suggesting the ABC analysis as a remarkable tool in biogeographic studies. Based on these results, the distribution of Dolichopoda species appears to have been shaped by the palaeogeographic and climatic events that occurred from Late Miocene up to the Plio-Pleistocene. Both vicariance and dispersal events appear to have influenced Dolichopoda species distributions, with many processes occurring in ancestral epigean populations before the invasion of the subterranean environment. Ó 2011 Elsevier Inc. All rights reserved.

1. Introduction Dolichopoda cave crickets (Orthoptera, Rhaphidophoridae) represent an excellent material to explore historical biogeographic patterns. A deep taxonomic research on this genus has resulted in about 50 relatively well defined species. Species delimitation has properly been confirmed by molecular systematics (Sbordoni et al., 2004; Allegrucci et al., 2005, 2009). The distributions of species are also well documented and encompass a wide geographic range across the Mediterranean region. Distributional patterns of the extant Dolichopoda species suggest a possible scenario for historical biogeography and the complex historical process of the Mediterranean region offers the possibility to refer to well-dated palaeoevents to calibrate the molecular clock of our robust phylogeny. During the last decade several studies have used molecular data to explore the biogeographic affinities in the Mediterranean region. Major palaeogeographical events such as the formation of Mid-Aegean trench, the Messinian salinity crisis and the Milankovitch climate oscillations since the late Tertiary have been advocated to ⇑ Corresponding author. Fax: +39 06 72595965. E-mail address: [email protected] (G. Allegrucci). 1055-7903/$ - see front matter Ó 2011 Elsevier Inc. All rights reserved. doi:10.1016/j.ympev.2011.04.002

explain the diverse evolutionary histories of Mediterranean taxa (Dermitzakis, 1990; Dermitzakis and Papanikolaou, 1981; Steininger and Rögl, 1984). Most of the studies carried out on the biogeography in the Mediterranean region concluded that a principally vicariant pattern of differentiation related to the geotectonics of the area has shaped the present-day distribution of the taxa studied (Weisrock et al., 2001; Gantenbein and Largiader, 2002; Parmakelis et al., 2005, 2006a,b; Kasapidis et al., 2005; Poulakakis et al., 2003, 2005). In this paper we analyze samples of ninety percent of known Dolichopoda species, collected across the entire range of the genus. This nearly complete species sampling allows a more accurate assessment of historical patterns and provides well-supported clades for assessing temporal trends in biogeography and diversification. Dolichopoda, a Mediterranean genus distributed from the eastern Pyrenees to Turkish Armenia and Caucasus mountains, belongs to the family Rhaphidophoridae (Orthoptera, Gryllacridoidea) which includes a large number of cave-adapted genera and species with worldwide distribution. Dolichopoda cave crickets species inhabit both natural and artificial caves. The latter are especially colonized in the northern part of the range where favorable climatic conditions occur and where individuals are frequently observed outside in moist or mesic woods. The dry hot

G. Allegrucci et al. / Molecular Phylogenetics and Evolution 60 (2011) 108–121

Mediterranean ‘‘macchia’’ which is mainly widespread in the south part of the range, represent a strong constraint for dispersal. Actually, the highest Dolichopoda diversity is found in insular and peninsular Greece where most of known species occur and where the Mediterranean ‘‘macchia’’ predominates. Results from previous papers based on molecular markers and carried out on the western and eastern species separately (Allegrucci et al., 2005, 2009) revealed a rather well resolved phylogeny at species level and suggested an eastern origin of Dolichopoda. In this paper, the historical biogeography of Dolichopoda genus was examined using the comprehensive species phylogeny. We used both published (Allegrucci et al., 2005, 2009) and original data. In particular, three more species from the most eastern part of the range (Anatolia and Caucasus mountains) have been added. Phylogeny was reconstructed on the basis of three mitochondrial genes (12S rRNA, 16S rRNA and COI), as in previous studies. Moreover, one nuclear gene (28S rRNA) was also considered. Patterns of dispersal were derived by divergence time estimates calculated assuming a relaxed molecular clock and based on Bayesian analysis (Drummond et al., 2006). Well-dated palaeoevents of the Mediterranean region have been used to calibrate the molecular clock and to obtain substitution rates for each of the four genes examined. We focused on three ‘‘evolutionary steps’’ of Dolichopoda history; (a) the colonization of Caucasus and Anatolia regions, (b) the colonization of insular and peninsular Greece, (c) the colonization of western Mediterranean. Testing of different alternative colonization scenario, for each step, was carried out by means of Approximate Bayesian Computation analysis (ABC; Beaumont et al., 2002; Tallmon et al., 2004). ABC method is based on the comparison of summary statistics computed on the observed data and on simulated datasets derived from an evolutionary model. The integration of ABC with software packages implemented for Coalescent-based molecular data simulation (Laval and Excoffier, 2004; Anderson et al., 2005) produces an extremely flexible framework that allows testing-complex evolutionary hypotheses. ABC is a relatively recent methodology that is more commonly used in population genetics and phylogeography (Csilléry et al., 2010). In this paper we assayed the use of this powerful analytical tool for testing models of deeper phylogenetic divergence (i.e. diversification at genus level). Considering species phylogenies, the population tree and the gene tree are expected to coincide, as a consequence of the long time that would lead to fixation of different lineages. However, lineage sorting can occur, with increasing probability as population sizes become larger or phylogenetic branches become shorter, thus representing a discrepancy between species and gene trees (Maddison, 1997). Incomplete lineage sorting can be a problem when analyzing recently diverged species. However, even in the case of an old but rapid burst of speciation, many short and deep internal branches could occur, thus increasing the probability that genetic lineages coalesce with lineages from distant species (Rosenberg and Tao, 2008). In these cases, phylogenetic methods that consider both nucleotide substitutions and the processes of population genetics could produce a more reliable species tree (Maddison and Knowles, 2006). Moreover, given the growing availability of multilocus and genomic data, there is a strong interest in methods based on a wide variety of approaches that can account for multispecies coalescent in estimating species tree (Degnan and Rosenberg, 2009; Liu et al., 2009). Model-choice methods like ABC that rely on a summary-statistic approach are considered to be easier to implement and more flexible than those based on likelihood analysis (Knowles, 2009), when testing alternative hypotheses. However, which summary statistics to be used to sufficiently distinguish between competing hypotheses and inform about parameters of interest remains problematic (Bertorelle et al., 2010). In our analyses, this aspect has been taken into consideration and discussed.

109

The nearly complete species sampling provides an appropriate framework to make reliable assessments of temporal patterns of diversification which are analyzed using Lineages-Through-Time (LTT) plots (Nee et al., 1992; Nee and May, 1994). Phylogenetic trees, especially those including all the living species in the taxonomic group of interest, provide an indirect record of the speciation events that have led to present-day species (Barraclough and Nee, 2001). The use of dated phylogenies allows visualizing fluctuations in the rate of diversification by plotting the number of lineages through time (Nee et al., 1992; Nee and May, 1994). Specific statistics, such as the gamma statistics (Pybus and Harvey, 2000) or birth–death likelihood test (Rabosky, 2006a) allow analyzing if these plots deviate from a constant rate of diversification. 2. Materials and methods 2.1. Taxon sampling Most of the species included in this study were previously analyzed in Allegrucci et al. (2005, 2009) (Table 1, Fig. 1). In this paper we present new data for species from Turkey and Caucasus Mountains region. In particular, we considered Dolichopoda noctivaga Di Russo and Rampini, 2007, from south-west Caucasus Mountains, Dolichopoda hyrcana Bey-Bienko, 1969, from eastern Caucasus and Dolichopoda lycia Galvagni, 2006, from south-western Anatolia (Table 1). This dataset contains 74 specimens including representative taxa of all four subgenera for a total of 44 species out of the 49 known species. Three additional taxa belonging to three different genera within the same family were used as outgroups (Troglophilus cavicola, Hadenoecus cumberlandicus and Euhadenoecus insolitus). Out of the five missing species, four are from the western part of the range and one is from Anatolia. However, they are from geographic areas widely represented from the species here considered. 2.2. Laboratory procedures DNA extraction, amplification and sequencing of mitochondrial genes (12S, 16S and COI) followed the same protocols as in Allegrucci et al., 2005, 2009. In this study the large subunit of the nuclear ribosomal DNA (28S rDNA) was also included. It was partially amplified and sequenced for all the 44 species of Dolichopoda here considered. In particular, a fragment of 580 base pairs, belonging to domains 3–5, was obtained using primers from Friedrich and Tautz (1997). Moreover, 12S gene was sequenced for the western species, previously analyzed using only 16S and COI genes (Allegrucci et al., 2005). Double-stranded amplifications were performed with a PerkinElmer-Cetus thermal cycler (PE Applied Biosystems, Foster City, CA, USA) in 25 ll reaction volume containing genomic DNA (10– 100 ng), 1.5 mM MgCl2, 2.5 mM of each dNTP, 0.5 lM primer, 1U of Amplitaq (Applied Biosystems, Foster City, CA, USA) and the buffer supplied by the manufacturer. 28S rDNA PCR was performed as follows: after a 2 min denaturation step at 95 °C, each cycle of PCR consisted of denaturation for 30 s at 95 °C, annealing for 30 s at 50 °C, and extension for 30 s at 72 °C. After 35 cycles, a 2 min extension step at 72 °C followed. Optimal cycling parameters for mitochondrial genes are the same as in Allegrucci et al. (2005). PCR products were purified using the ExoSAP digestion (GE Healthcare Europe, Munich, Germany), directly sequenced in both directions using the BigDye terminator ready-reaction kit, and resolved on ABI 3100 Genetic Analyzer (PE Applied Biosystems, Foster City, CA, USA), following the manufacturer’s protocols.

110

G. Allegrucci et al. / Molecular Phylogenetics and Evolution 60 (2011) 108–121

Table 1 Dolichopoda species and outgroup taxa included in this study. Localities

Code

Genebank Accession No.

Covoli di Veneto Cave, Veneto, Northern Italy Bat Cave, Carter Cave State Park, Carter Co., KY, USA

TRO BAT

12S: EF216946 16S: AY793564 COI: AY793624 28S: EF217003 12S: EF216948 16S: AY793562 COI: AY793592 28S: EF217004

Indian Grave Point Cave, De Kalb Co., TN, USA

IND

12S: EF216947 16S: AY793563 COI: AY793591 28S: EF217005

Sirach Cave, Eastern Pyrenees Brando Cave, Corsica Island

SIR BRA

Sisco Cave, Corsica Island

SIS

Valletto Cave, Corsica Island Sabara Cave, Corsica Island

VLT SAB

12S: JF826039 12S: JF826047 JF826069 12S: JF826048 JF826070 12S: JF826050 12S: JF826049

Corno Cave, Piemonte Pugnetto Cave, Piemonte

CON PUG

12S: JF826040 16S: AY793568 COI: AY793604/AY793605 28S:JF826062 12S: JF826041 16S: AY793569 COI: AY793601/AY793602/AY793603 28S:JF826063

Eastern-North and Central Italy D. laetitiae Poscola Cave, Veneto Diavolo cave, Tuscany Piane Cave, Umbria

PSC DIA GDP

12S: JF826054 16S: AY793581 COI: AY793611/AY793613 28S:JF826076 12S: JF826052 16S: AY793580 COI: AY793614/AY793615 28S: JF826074 12S: JF826053 16S: AY793582 COI: AY793610/AY793612 28S: JF826075

Western-Central and Southern Italy D. baccettii Punta degli Stretti Cave, Tuscany D. aegilion Campese Mine, Giglio Island, Tuscany D. schiavazzii Pipistrelli Cave, Tuscany Marciana Cave, Elba Island, Tuscany Fichino Cave, Tuscany D. muceddai Limbara Mount, Sardinia D. geniculata Valmarino Cave, Latium Fontanelle Cave, Campania Ischia cellars, Ischia Island, Campania D. capreensis San Michele Cave, Capri Island, Campania D. palpata Tremusa cave, Calabria

PST CAM ORS MRC FIC SAR VAL FON ISC CPR TRE

12S: 12S: 12S: 12S: 12S: 12S: 12S: 12S: 12S: 12S: 12S:

Velika Cave, Blato, Miljet Kod Solina Cave, Govedari, Miljet

VEL SOL

12S: EF216944 16S: EF216967 COI: EF217019 12S: EF216944 16S: EF216967 COI: EF217018 28S: EF216982

D. thasosensis D. graeca

Pozarska Mala Pestera, Loutrakiou, Pella Waterfall Cave, Edessa, Pella Apano Skala Cave, Naoussa, Imathia Saranda Outdate Cave, Naoussa, Imathia Aghlia Paraskevi Cave, Tembi Valley, Larissa Small caves, Kato Olimpos, Leptokaria, Kalipefki Drakotripa Cave, Panayia, Thasos Island, Kavala Perama Cave, Ioannina, Epiro

POZ EDE IZB NAU TEM OLI DRA PER

12S: 12S: 12S: 12S: 12S: 12S: 12S: 12S:

EF216939 16S: EF216969 COI: AY793637 EF216939 16S: EF216969 COI: AY793637 28S: EF217001 EF216943 16S: EF216973 COI: EF217031 EF216943 16S: EF216973 COI: EF2217032/3/4 28S: EF216990 EU887846 16S: EU887861 COI: EU887894/5/6 28S: EU887876 EU887845 16S: EU887860 COI: EU887891/2/3 28S: EU887875 EF216926 16S: EF216956 COI: EF217020 /1 28S: EF216983 EF216923 16S:EF216953 COI:EF217013 28S: EF216979

Ionian Islands D. steriotisi D. gasparoi D. ithakii D. patrizii D. pavesi

Antropograva Cave, Klimatia, Kerkira, Corfu Chirospilia Cave, Evghiros, Levkada Marmarospilia cave, Vathi, Ithaki Island Small cave, Petalas Drogarati Cave, Sami, Kefalonia Island

ANT CHI ITA PET SPI

12S: 12S: 12S: 12S: 12S:

EF216925 16S: EF216955 COI: EF217016/7 28S: EF216981 EF216920 16S: EF216950 COI: EF217008/9 28S: EF216976 EF216919 16S: EF216949 COI: EF217006/7 28S: EF216975 EU887847 16S: EU887862 COI: EU887898/9/ 900 28S: EU887877 EF216921 16S: EF216951 COI: EF217010/1 28S: EF216977

ORO AGH AND

12S: EF216922 16S: EF216952 COI: EF217012 28S: EF216978 12S: EF216924 16S: EF216954 COI: EF217014/5 28S: EF216980 12S: EU887848 16S: EU887863 COI: EU887878 28S: EU887878

KOR HER GLK PAN JOA TRI PKI KSA

12S: 12S: 12S: 12S: 12S: 12S: 12S: 12S:

Outgroup Troglophilus cavicola Hadenoecus cumberlandicus Euhadenoecus insolitus Ingroups Genus Dolichopoda Western-South France D. linderi D. bormansi

D. cyrnensis Western-North Italy D. ligustica

Croatia D. araneiformis Northern Greece D. remyi D. hussoni D. annae

Central-Western Greece D. giachinoi Megalospilio Cave, Monastirakion, Aitolo-Akarnania D. kiriakii Kiriaki Cave, Korifè, Aghlia Kiriaki, Parga D. lustriae Aghios Andreas Cave, Valtou M., Halkiopuli, Etolia Central-Eastern Greece D. sp.‘‘Parnaso’’ D. vandeli D. D. D. D.

insignis petrochilosi cassagnaui makrykapa

Korykion Andron Cave, Mount Parnitha, Attica Hermes Cave, Orkomenos, Dhionisos, Beotia Cave over Kopais Lake, Orkomenos, Beotia Panos Cave, Marathon, Athene, Attica Aghlia Joannis, Nea Pendeli, Athene, Attica Aghlia Triada Cave, Karistos, Eubea Island Paralia Kilidau Cave, Lamari, Eubea Island Paralia Pot Cave, Kao Seta, Aghlia Triada, Eubea Island

JF826046 JF826045 JF826044 JF826043 JF826042 JF826051 JF826055 JF826057 JF826058 JF826059 JF826060

16S: AY793567 COI: AY793598/AY793599 28S: JF826061 16S: AY793578 COI: AY793631/AY793632/AY793627 28S: 16S: AY793579 COI: AY793625/AY793626/AY793628 28S: 16S:AY793577 COI: AY793620/AY793621 28S:JF826072 16S: AY793576 COI: AY793618/AY793619 28S:JF826071

16S: 16S: 16S: 16S: 16S: 16S: 16S: 16S: 16S: 16S: 16S:

AY793571 AY793570 AY793573 AY793572 AY793574 AY793575 AY793583 AY793584 AY793585 AY793587 AY793588

COI: AY793639/AY793640 28S:JF826068 COI: AY793600 28S: JF826067 COI: AY793633 28S: JF826066 COI: AY793635 28S:JF826065 COI: AY793634/AY793636 28S: JF826064 COI: AY793629/AY793630 28S: JF826073 COI: AY793616/AY793617 28S: JF826077 COI:AY793594 28S: JF826079 COI: AY793595 28S: JF826080 COI: AY793606/AY793607 28S: JF826081 COI: AY793608/AY793609 28S: JF826082

EU887851 16S: EU887866 COI: EU887905/6 28S: EU887880 EF216932 16S: EF216962 COI: EF217039/40 28S: EF216992 EF216932 16S: EF216962 COI: EF217038 EF216938 16S: EF216968 COI: EF217054 28S: EF217000 EF216937 16S: EF216967 COI: EF217053/1/2 28S: EF216999 EF216931 16S: EF216961 COI: EF217035/6/7 28S: EF216991 EF216941 16S: EF216971 COI: EF217942 EF216941 16S: EF216971 COI: EF217041 28S: EF216993

G. Allegrucci et al. / Molecular Phylogenetics and Evolution 60 (2011) 108–121

111

Table 1 (continued)

Eastern Greek Islands D. naxia D. giulianae D. calidnae

Northern Peloponnesus D. matsakisi

Localities

Code

Genebank Accession No.

Zeus Cave, Filotas, Naxos Island, Cyklades Small Cave, Stauros, Naxos Island, Cyklades Moni Spilianis Cave, Pithagorion, Samos Island Seven Virgins Cave, Kalimnos, Kalimnos Island, Dodecanese Skalia cave, Skalia, Kalimnos Island, Dodecanese

ZEU STA SPS EPT

12S: 12S: 12S: 12S:

SKA

12S: EF216934 16S: EF216964 COI: EF217048 28S: EF216996

Analipsi Cave, Pititsa, Achaia Ton Limnon Cave, Kastri, Kalavrita, Achaia

ANA KASA KASB KEF

12S: EF216927 16S: EF216957 COI: EF217022/3 28S: EF216984 12S: EF216927 16S: EF216957 COI: EF217024 28S: EF216985

Kataphingi Cave, Selitsa, Messenia River cave of Glyfada, Dirou, Aeropolis, Laconia Aghia Varvara Cave, Parori, Sparta, Lakonia (Taygeto Mt) Taygeto Shelter,Parori, Sparta, Lakonia

KAT DIR VAR

12S: EF216940 16S: EF216970 COI: EF217045/6 28S: EF216994 12S: EF216940 16S: EF216970 COI: EF217043/4 12S: EU887850 16S: EU887865 COI: EU887904

TAY

12S: EU887849 16S: EU887864 COI: EU887902/3 28S: EU887879

Aghlia Paraskevi Cave, Skotinon, Iraklio Atzigano cave, Adrianos, Lassithi Dikteion Antron Cave,Psychro, Iraklio

PAR NIK DHI

12S: EF216942 16S: EF216972 COI: EF217027/8 28S: xxxxxxxx 12S: EF216942 16S: EF216972 COI: EF217030 28S: EF216989 12S: EF216930 16S: EF216960 COI: 217029 28S: EF216988

Karain Cave, Antalya, Turkey Insuyu cave, Burdur, Turkey Gedelma cave, Antalya, Turkey Kafkasor cave, Artvin, Turkey

12S: 12S: 12S: 12S:

Duzkoy cave, Artvin, Turkey Coruh Valley, Erzurum, Turkey

KAR INS TUR KAFb DUZ COR

Vorontzovskaya, Caucasus, Russia Golova Otapa Cave, Caucasus, Russia Talysh cave, Lenkoran, Trancaucasia eastern

VOR GOL LEN

12S: EF216945 16S: AY793566 COI: AY793622 28S: EF217002 12S: xxxxxxxx 16S: AY793565 COI: AY793623 12S: EU887856 16S: EU887871 COI: EU887856 28S: EU887885

D. sp. ‘‘Limnon’’

Ton Limnon Cave, Kastri, Kalavrita, Achaia

D. dalensi

Kefalovrisi Cave, Argos, Argolide

Southern Peloponnesus D. unicolor D. sp. ‘‘Taygeto’’

Crete Island D. paraskevi D. sp. ‘‘Crete’’ Turkey D. sbordonii D. lyciae D. noctivaga

Caucasus region D. euxina D. hyrcana

2.3. Sequence alignment and analysis Each gene fragment (12S, 16S, COI and 28S) was considered separately for the alignment. Sequences of 16S, 12S and 28S rDNA were aligned using CLUSTAL_X 1.81 (Thompson et al., 1997) with opening gap = 10 and extending gap = 0.10. The alignment, also checked by eye, did not require further improvement by considering a secondary structure model. Cytochrome oxidase I nucleotide sequences were assembled, aligned, and translated with CodonCode Aligner 2.0.2. The four genes were treated as separate data partitions for phylogenetic analyses. Congruence among partitions was assessed by the partition-homogeneity test [also named the incongruence length difference (ILD) test; Farris et al., 1994, 1995], as implemented in PAUP. All data sets were subsequently combined for a total evidence analysis. All phylogenetic analyses were performed using Bayesian inference as implemented by the software MrBayes 3.1b4 (Huelsenbeck and Ronquist, 2001). Mrmodeltest (Nylander, 2004) was used to perform hierarchical likelihood ratio test and calculate approximate AIC values of the nucleotide substitution models for each gene fragment. Indels were coded as binary characters and included as separate binary data partition in the analysis. The evolutionary model for this partition was a single F81-like model (Felsenstein, 1981) as implemented in MrBayes. For all Bayesian analyses, at least two simultaneous searches were conducted comprising four Markov chains started from a randomly chosen tree and run for 5000,000 generations, with sampling every 100 generations. The first 5000 trees were discarded as burn-in and posterior probabilities (PP) were calculated from

EU887853 16S: EU887868 COI: EU887909/10/11 28S: EU887882 EU887852 16S: EU887867 COI: EU887908 28S: EU887881 EF216935 16S: EF216965 COI: EF217049 28S: EF216997 EF216933 16S: EF216963 COI: EF217047 28S: EF216995

12S: EF216928 16S: EF216958 COI: EF217025 28S: EF216986 12S: EF216929 16S: EF216959 COI: EF217026 28S: EF216987

EF216936 16S: EF216966 COI: EF217050 28S: EF216998 EU887854 16S: EU887869 COI: EU887917 28S: EU887883 EU887855 16S: EU887870 COI: EU887918 28S: EU887884 EU887859 16S: EU887874 COI: EU887922 28S: EU887888

12S: EU887858 16S: EU887873 COI: EU887921 28S: EU887887 12S: EU887857 16S: EU887872 COI: EU887920 28S: EU887886

post burn-in trees. The following descriptors were assumed to indicate convergence on a common phylogenetic topology by separate Bayesian searches: similarity in log-likelihood scored at stationarity, similarity in consensus tree topologies and PP values for supported nodes, and a final average standard deviation of split frequencies (ASDSF) for simultaneous searches approaching zero.

2.4. Dating of the cladogenetic events A likelihood ratio test (LRT) was used to test the molecular clock hypothesis, according to Huelsenbeck and Crandall (1997). The LRT results suggested a heterogeneous rate variation for each gene and therefore a clocklike evolution of sequences could not be assumed. Consequently, dates of divergence were inferred using a relaxed molecular clock, following the uncorrelated relaxed lognormal clock as implemented in BEAST (v.1.5.4, Drummond et al., 2006; Drummond and Rambaut, 2007). This analysis was limited to species level and carried out on the concatenated matrix, partitioned by genes. A Yule or ‘‘pure birth’’ prior process was also applied to model speciation. The time to the most recent common ancestor (MRCA) between each clade was estimated under the models highlighted in MrModeltest (Nylander, 2004) for each gene. We did two independent runs with BEAST, each one for 10 million steps. One million steps were discarded as burn-in. Convergence to stationarity and effective sample size (ESS) of model parameters were assessed using TRACER 1.5 (Rambaut and Drummond, 2007). Samples from both independent runs were pooled after removing a 10% burn-in using Log Combiner 1.5.4 (Drummond and Rambaut, 2007).

112

G. Allegrucci et al. / Molecular Phylogenetics and Evolution 60 (2011) 108–121

Fig. 1. Range of Dolichopoda genus. Sampling sites are also indicated. Acronyms refer to the assayed populations as detailed in Table 1.

Well-dated events of geological vicariance were used to calibrate the molecular clock. In particular, we used two dated palaeogeographic events: (i) the separation of Giglio Island (an island in the Tuscan Archipelago, Thyrrenian Sea, Italy, inhabited by Dolichopoda aegilion) from the continent, occurred during a marine transgression that took place around the Calabria-Ionian transition, 0.7– 1.0 million years ago (Ma, Lipparini, 1976; Mazzanti, 1983); (ii) the isolation of Crete Island, inhabited by Dolichopoda paraskevi (Aegean Sea, Greece) from neighboring mainland, happened at the end of Messinian salinity crisis (5.33 Ma, Krijgsman et al., 1999). After this period, Crete became permanently isolated (Dermitzakis, 1990). These calibration points were implemented in BEAST using lognormal prior distribution ranging from 0.75 to 1.2 Ma (isolation of Giglio) and from 5 to 5.3 Ma, (isolation of Crete). Estimates of substitution rates were summarized in TRACER 1.5, multiplying the ‘‘mu’’ parameters of each partition by the overall rate to obtain the absolute rate for each partition. 2.5. Biogeographic analyses Based on the phylogeny, ancestral distributions were reconstructed using different analyses. First of all we used dispersal– vicariance analysis in the program DIVA (vers. 1.1, Ronquist, 1996). This method reconstructs ancestral distributions on a phylogeny from current distributions by optimizing dispersal and vicariance events using a cost matrix. We performed DIVA analysis based on the rooted Bayesian tree generated in BEAST and limited to the species level. The distribution of each species was classified in ten different geographic areas: (1) Pyrenees mountains region,

(2) Coastal Italy, (3) Inland Italy, (4) Sardinia and Corsica Islands, (5) Eastern Greece, (6) western Greece, (7) Peloponnesus, (8) Aegean islands, (9) South-western Turkey and (10) Caucasus Mountains region. The selection of these areas was based on the geographic distribution of considered species. The maximum number of ancestral areas was set to ten (the total number of areas in the analysis) and the analysis was carried out with and without outgroup. Taking advantage of the Approximate Bayesian Computation approach (ABC; Beaumont et al., 2002; Tallmon et al., 2004), we then focused on three ‘‘evolutionary steps’’ of Dolichopoda history implementing the ABCtoolbox pipeline (Wegmann et al., 2010) with BayeSSC, as the molecular simulation routine (Laval and Excoffier, 2004; Anderson et al., 2005). With the application of this methodology, the following competing hypotheses on uncertain phases of Dolichopoda diversification and spread were compared and the relative posterior support estimated: (1) The colonization of Anatolia and Caucaso (Clade H and I). We tested three alternative hypotheses concerning the relative ancestry of these species: (i) both clades originated from a single ancestral population that diverged at the base of the Dolichopoda phylogeny (AC-H1); (ii) Caucasian species diverged first (AC-H2); (iii) Anatolian species diverged first (AC-H3). (2) The colonization of peninsular and insular Greece. We tested three alternative hypotheses, previously analyzed in Allegrucci et al. (2009) by means of deep coalescence analysis. The comparison of our results with those of previous analyses allowed us to test the performance of ABC. In particular,

G. Allegrucci et al. / Molecular Phylogenetics and Evolution 60 (2011) 108–121

we considered: (i) a single colonization event from Anatolia, through the present-day Bosporus, and the colonization of the Aegean islands during the Messinian marine regression (GR-H1); (ii) the separation of the Hellenic and Aegean species from the Anatolian species, after the opening of the MidAegean Trench. These two groups colonized their present Greek range following two separate routes: some species expanded towards north, ascending the Hellenic Peninsula, while others expanded towards south and east, crossing the Aegean during the Messinian salinity crisis (GR-H2); (iii) the origin of Greek species from two independent Messinian dispersals from Anatolia (GR-H3). (3) The colonization of western Mediterranean. We tested the following hypotheses: (i) the diversification of inland and coastal groups from a single colonization event from the North of Italy (WM-H1); (ii) the diversification of inland and coastal groups from two different colonization events from the North of Italy (WM-H2); (iii) the diversification of inland and coastal groups from two different colonization events from the North and South-East of Italy, respectively (WM-H3). In order to find out which of these competing hypotheses would have the highest posterior probability support, we modeled three sets of three alternative hypotheses (3 hypotheses  3 evolutionary phases in Dolichopoda). As regards each evolutionary phase, a model-choice procedure based on the Approximate Bayesian Computation statistical framework was implemented as follows. 100,000 molecular datasets were simulated under each of the nine hypotheses using the ABCsampler software in the ABCtoolbox pipeline coupled with BayeSSC: each dataset consists of two blocks corresponding to (i) the concatenation of the three mitochondrial markers (COI, 12S and 16S) and (ii) the nuclear marker (28S). Each block was produced independently by a single ABCsampler iteration during which it sampled model parameter values from specified prior distributions and passed these values to BayeSSC, parsing two separate input files. These two files were pre-compiled with fixed values for locus-specific parameters, as regards length in base pairs (mt locus: 1883 bp; nuclear locus: 580 bp), nucleotide substitution pattern (mt locus: Ti/Tv = 5; Gamma distribution shape parameter = 0.5 and number of rate classes = 4; nuclear locus: Ti/ Tv = 5; no Gamma distribution was considered), and molecular substitution rate per generation for the whole fragment (mt locus: 0.000025; nuclear locus: 0.0000006.). Substitution rates were derived from the chronogram calculated in BEAST. As regard the mitochondrial genes (12S, 16S and COI) a weighted mean among the obtained values was considered. One individual was sampled in each simulated species and included in the statistical analyses. Prior distributions for model parameters were indicated in one model-specific estimation file common to both mitochondrial and nuclear block. It constrained the simulations according to specific values concerning the effective population size of each species and the timing of the speciation events (species tree node heights). Effective population size of each species was set constant through time and with values drawn at random from a uniform distribution (U: 100–50,000). Speciation events were modeled as historical events, considering the calibrated Bayesian phylogenetic tree (Fig. 3A). The effects of different prior distributions on node heights were studied by considering: (i) a truncated normal distribution (N-dist) based on Credible Region detected by Bayesian analysis and (ii) a uniform distribution (U-dist) with an arbitrarily very narrow range around mean node heights derived by the same Bayesian chronogram. These parameters were regulated on each other in order to avoid inversions in tree topology. Alternative topologies for the species trees were drawn accordingly to the alternative biogeo-

113

graphic hypotheses (Fig. 4). During each iteration of the simulation process, ABCsampler drew a set of parameter values from the prior distribution and parsed these values in the two locus-specific files specified above completing the input information for the simulation software. We focused on different portions of the phylogenetic tree whose changes could be considered as essential in order to distinguish among competing hypotheses; in our tests, we considered two different pruned trees: one in the case of Anatolian/Caucasian and Greek species, and the other one in the case of west Mediterranean species. Detailed description of the historical events modeled in each hypothesis is given in Supplementary Information. In order to summarize and compare the simulated data, we decided to employ the pairwise genetic distance as summary statistics. It can be considered as a raw index of genetic differentiation and it is strongly related with the tree topology. Considering each of the three evolutionary phases described above, different groups of species pairs, whose genetic distances may be differentially affected under competing models (i.e. species trees), were preliminary screened and tested as regards their effectiveness as summary statistics. Then, a selected group of species pairs (reported in the Supplementary Information), whose changes in the genetic distances were consistent with the differences in the source models, were selected and included in the summary statistic vector. Three different vectors were built, one for each of the three evolutionary phases considered in Dolichopoda diversification. Since the simulation software used here allows for calculation of the summary statistics we were interested in, its output was stored apart after each iteration. When a run for each competing model was completed, this output was passed to the R-Statistical Package in order to estimate their relative posterior support. Considering the observed dataset, the pairwise genetic differences (pdistance, as calculated on the simulated datasets) between the selected species pairs were calculated by Arlequin 3.1 (Excoffier et al., 2005). The vector of the summary statistics was used to estimate the Euclidean error, according to the formula ffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi ffi s p p 2 p p  j i j i þ þ    with: pi,j,. . . = ith, jth, . . . , pairwise differa a i

j

ence in observed dataset; pi;j;... ¼ ith; jth; . . ., pairwise difference in simulated dataset; ri;j;;... ¼ standard deviation of ith, jth, . . . , pairwise difference in simulated dataset. Posterior probability (PPABC) of each competing hypothesis was estimated through a weighted multinomial logistic regression procedure (Fagundes et al., 2007), using the calmod function as implemented in R-Statistical Package. Local regression was performed with threshold (e) values ranging from 0.0005 to 0.01 in order to check for consistency of the results.

2.6. Diversification analysis The tempo of diversification was assessed using the comprehensive Bayesian phylogenies of the genus Dolichopoda. Chronogram obtained by BEAST analysis was used to construct semilogarithmic LTT plots in GENIE ver. 3.0 (Pybus and Rambaut, 2002). To examine confidence in the estimated datings, confidence intervals (95%) were estimated using 2000 trees from the pool of converged Bayesian trees in BEAST analysis. Trends resulting were plotted together in the observed LTT plot. Two statistical analyses were used to test for significant departures from the constant speciation rate model (CR), the c-statistic (Pybus and Harvey, 2000) and the birth–death likelihood test (BDL, Rabosky, 2006a). To take into account incomplete sampling simulations were carried out using different numbers of missing taxa with the Monte Carlo Constant rates test (MCCR, Pybus and Harvey, 2000) in LASER (Rabosky, 2006b). The gamma statistic of Pybus and Harvey

114

G. Allegrucci et al. / Molecular Phylogenetics and Evolution 60 (2011) 108–121

(2000) relates to the distributions of internode distances through time and under the pure birth model follows a standard normal distribution with a mean of 0. A c value less than 1.645 rejects the pure birth model under a one-tailed test and provides support (a = 0.05) for the alternative hypothesis of slowdown. BDL analysis (Rabosky, 2006a) was carried out using LASER library in R (Rabosky, 2006b) to test the null hypothesis of rate-constancy against rate-variable models using AIC (Akaike, 1973). Significance of the change in AIC scores was determined by creating a distribution of AIC scores by simulating 5000 trees and using yuleSim in LASER (Rabosky, 2006b). The same number of taxa and the same speciation rate as those determined under the pure birth model were considered.

3. Results and discussion 3.1. Phylogenetic analysis MrModeltest (Nylander, 2004) indicated GTR + I + U (Lanave et al., 1984; Gu et al., 1995) and K80 + I (Kimura, 1980) as the best models of DNA substitution for the mitochondrial genes and the nuclear one, respectively. The phylogeny in Fig. 2 supports the major phylogenetic relationships previously demonstrated (Allegrucci et al., 2005, 2009; Martinsen et al., 2008) and it can be considered the best estimate of relationships in the genus Dolichopoda as it includes data from different independent genes and almost all known species. Some unresolved relationships within the western (Allegrucci et al., 2005) and eastern (Allegrucci et al., 2009) species are here well defined. Each lineage was well supported (PP P 0.90, Fig. 2) and corresponded to distinct geographical areas in the Mediterranean region: clades A and B are from western Mediterranean (Italy and France); clades C-G are from eastern Mediterranean (Greece, Anatolia and Caucasus region). In particular, clade C is from north-eastern Greece, clade D from north-western Greece and north Peloponnesus, clade E from western Greece, clade F from south-eastern Greece and south Peloponnesus, clade G from Aegean Islands. The new assayed species from Anatolia and Caucasus regions, considered in the present paper, form two clearly distinct clades (H-I) and are sister to the rest of Dolichopoda species although their relative relationships are not resolved (Fig. 2). Trees topology of each studied gene (Fig. 2) is rather similar to the topology obtained by combining all data sets. However, there is a conflict between COI and 12S/16S trees over the position of Dolichopoda thasosensis. As hypothesized in Allegrucci et al. (2009), the colonization of Greece by Dolichopoda species followed two different routes, probably originating in the easternmost part of the current range. D. thasosensis, endemic of Thasos island and belonging to the lineage currently occurring in western Mediterranean, northern Greece, Ionian islands, and north-eastern Peloponnesus, could share several characters both with the western Mediterranean species and with the Ionian ones, as revealed by the individual COI and 12S/16S gene trees. However, D. thasosensis is sister to the geographically closer Ionian species (PP = 1.00, Fig. 2) when combined dataset is considered. As far as the other taxa relationships are concerned we could observe the same topology, with some unresolved deep nodes. In particular, relations among A–D clades were always resolved both in mitochondrial ribosomal genes (12S and 16S rRNA) and COI, while relationships among E– H clades were resolved only in COI gene. Nuclear ribosomal gene (28S rRNA) produced a topology distinguishing only some clades although most of the taxa relationships were unresolved. This is probably due to the low level of polymorphism revealed in this gene. However, its presence in the combined dataset increased the resolution of the deep nodes in phylogeny.

3.2. Divergence times The data showed a non-clock like behavior, with a coefficient of variation of 0.417 (95% highest posterior density, HPD: 0.289– 0.545). Molecular substitution rate estimates for each of the four studied genes are reported in Table 2. COI gene showed the highest estimate, being the mean substitution rate equal to 1.6% per site/ lineage/million years. 12S and 16S rRNA showed lower estimates with values of 1.1% and 0.7%, respectively. Finally, a mean substitution rate of 0.06% was estimated for 28S rRNA. These estimates are close to those calculated in other groups of insects in the recent literature. In particular, Papadopoulou et al. (2010) found a mean substitution rate for COI equal to 1.68% for Aegean Coleoptera tenebrionids, Shapiro et al. (2006) equal to 1.5–2% for Orthoptera Tettigonidae, Pons and Vogler (2005) equal to 1.67% for Coleoptera Cicindelidae and Kiyoshi and Sota (2006) equal to 1.55% for Odonata Gomphidae. In contrast the mitochondrial rRNA genes evolve more slowly, being the rates equal to 0.54% for Aegean tenebrionids (Papadopoulou et al., 2010), to 0.23% for Coleoptera Crysomelidae (Gomez-Zurita et al., 2000) and to 0.38% for Coleoptera Cicindelidae (Pons and Vogler, 2005). Regarding the 28S nuclear gene, the obtained mean rate for Dolichopoda (0.06%) completely overlap the estimates from Papadopoulou et al. (2010) for Aegean tenebrionids. Dating estimates indicated that radiation of the ingroup started in the late Miocene and that the Messinian salinity crisis in the Mediterranean region has been the main palaeogeographic event responsible for the initial diversification of the genus Dolichopoda (Fig. 3a). Estimation of molecular rates is in principle a remarkable tool for investigating evolutionary history; however the reliability of such estimates strongly depends on the accuracy by which they are estimated and on the appropriateness of the calibration method. The choice of calibration events is therefore crucial to the accuracy of molecular dating. In the present case we chose calibration events that only partially astride the threshold (more recent than ca. 1 Ma, Ho and Larson, 2007; Ho et al., 2005; Gratton et al., 2008) sensible to the time-dependency of molecular rates. Then, the predicted speed-up of molecular rates should not affect our calibration. Following the complex geological history of the Mediterranean region, we could refer to two major geological events to calibrate our molecular clocks. The choice of these calibration points produced molecular substitution rates in agreement with those estimated in other groups of insects and remarkably close to those estimated for Orthoptera by Shapiro et al. (2006). Dating estimates obtained in the present paper are in excellent agreement with the previous ones (Allegrucci et al., 2005, 2009) obtained by using substitution rates already known for insects (Brower, 1994). Moreover, as far as the western Italian species are concerned, these dating estimates completely overlap those obtained by using different molecular markers (Allegrucci et al., 1992, 2005; Sbordoni et al., 1985, 1987). These observations reinforce our interpretation on the likely ancient scenario that shaped the diversification pattern in Dolichopoda genus. 3.3. Biogeographic analyses The best ancestral distribution reconstructions, using dispersal vicariance analysis (DIVA), were obtained using outgroup and no ‘‘maxareas’’ restrictions. The root node distributions are large and include almost all of the area occupied by the terminal species. However, nodes immediately following the root node clearly identify the distributions and the ancient splits (Fig. 5b). The DIVA analysis suggested that eight phylogenetically independent dispersal events have occurred in the direction of eastern to western Medi-

G. Allegrucci et al. / Molecular Phylogenetics and Evolution 60 (2011) 108–121

115

Fig. 2. Bayesian analysis for each gene separately and for the combined dataset. Scale refers to number of substitution per site. Posterior probability (PP) values P0.80 are reported on each node.

116

G. Allegrucci et al. / Molecular Phylogenetics and Evolution 60 (2011) 108–121

Table 2 Substitution rates per site/lineage/million years in Dolichopoda cave crickets. Gene

Rate ± SD

95% HPD Lower  Upper

COI 12S 16S 28S

1.62 ± 1.14 1.12 ± 7.65 6.63 ± 4.45 6.44 ± 4.06

1.32  2.02 7.73  1.42 4.83  8.53 3.74  9.34

terranean. The first dispersal event would have occurred toward Anatolia and eastern Greece from Caucasus region. The other dispersal events concerned the colonization of Peloponnesus, western Greece and western Mediterranean region. In particular, dispersion toward Peloponnesus occurred two times, independently, from both eastern Greece and North-western Greece. Dispersion toward mainland Italy and French Pyrenees occurred from eastern Greece. Dispersion to Corsica and Sardinia probably occurred more recently from mainland Italy to the two Islands (Fig. 5b). ABC analysis, carried out considering the three ‘‘evolutionary steps’’, mainly confirmed results from DIVA analysis and (Fig. 5a) suggested that: (1) Caucasian species (clade H) would be the most ancient ones (e = 0.005, U-dist: PPABC = 0.97; N-dist: PPABC = 0.99 for ACH2). (2) Greek species would originate from two independent Messinian colonization events from Anatolia (e = 0.005, U-dist and N-dist: PPABC close to 1 for the GR-H3), confirming results in Allegrucci et al. (2009). (3) The diversification of inland and coastal groups in western species was not completely disentangled (e = 0.005, U-dist: PPABC = 0.85 for WM-H2; N-dist: PPABC = 0.55 for WM-H1; PPABC = 0.45 for WM-H2), when the wider normal distribution range for tree node heights was applied. In that case, the time-scale of hypotheses WM-H1 and WM-H2 mainly overlapped and became hardly distinguishable. Employing different threshold values on the retained simulations for the local regression step did not substantially change our results. In this first application of ABC at phylogenetic level, we decided to employ the pairwise genetic distance, which is a basic and simple index of genetic differentiation, strongly related with the tree topology. The selection of a sufficient vector of summary statistics can likely be considered as the most challenging aspect of ABC analysis. Many efforts have been made in order to find a method that could maximize the information retained in summary statistics and maintaining, at the same time, a low level of dimensionality (e.g. Wegmann et al., 2009); however, advantages of these efforts are not fully demonstrated. As regards the application at phylogenetic level, a sufficient summary statistic has not yet been developed, though some applications of the topological distance index in an ABC framework seem promising and should be evaluated (Battagliero et al., 2010). Based on our results, it seems that ABC analysis can be considered as a useful tool to deal with phylogenetic as well as phylogeographic/population genetics issues. Its main advantage lies in the fact that it allows comparison of complex historical models in a more intuitive way. By this perspective, it can be considered as an actual improvement for the broad scientific community of biogeographers and phylogenetists. However, a long work of refinement by statistically-skilled researchers is needed to further implement and test different aspects of ABC analysis. Based on these results we attempted to reconstruct possible ancient scenarios explaining the present-day distribution of Dolichopoda species. ABC and DIVA analyses suggested a Caucasian origin

of this genus (Fig. 5a and b) with Caucasian species (north-eastern Anatolia) and Anatolian ones (south-western Anatolia) not forming a monophyletic group (Fig. 2). Divergences between these species appear to be consistent with the rise of Anatolian Plateau that uplifted 5–10 Ma due to acceleration of northward movement of the Arabian Plate (Quennell, 1984; Steininger and Rögl, 1984). Populations from the north margin of Anatolia dispersed westwards colonizing the Europe, while populations from the south margin of Anatolia were restricted to the Lycian Taurus. Such a geographic pattern has been also observed in Castanea sativa (Villani et al., 1999) that is part of mesophilous woods to which the ancestor of Dolichopoda was strictly linked. Moreover, our results appear to be consistent with those obtained for other animal taxa ecologically associated to mesophilous woods, presenting more or less disjunct distributions in the ponto-caucasian area. For example, Weisrock et al. (2001) suggested that the two species of the ‘‘true’’ salamander from south-western and north-eastern Anatolia have had separate historical biogeographic origins, following the complex geological history of Anatolia. As far as the colonization of Greece by Dolichopoda species is concerned, results from DIVA and ABC analyses confirmed outcomes in Allegrucci et al. (2009), indicating two colonization episodes and two different routes. The southern lineage likely arose from a trans-Aegean colonization during the Messinian salinity crisis (5.96–5.33 Myr ago). The northern lineage could be the result of dispersal from the north through the Balkan Peninsula. The opening of the Mid-Aegean trench could have promoted an initial diversification within the uprising Anatolian plateau, while the Messinian marine regression offered the conditions for a rapid dispersal through the whole Aegean-Hellenic region. As far as the colonization of western Mediterranean is concerned, results from ABC and DIVA analyses do not completely agree. In particular, ABC analysis did not clearly discriminate between WM-H1 and WM-H2 hypotheses, but based on previous studies (Sbordoni et al., 1985; Allegrucci et al., 1987, 1992, 1997) and on DIVA results, we support the hypothesis WM-H1 as biogeographically more reliable. Indeed, we argue that this result could be influenced by the degree of genetic differentiation of the Italian inland species, Dolichopoda geniculata and Dolichopoda laetitiae. Based on previous studies carried out on independent molecular data (Sbordoni et al., 1985; Allegrucci et al., 1987, 1992), we know that populations belonging to these two species are highly genetically structured with interspecific genetic distance values sometimes overlapping the interspecific ones. This result has also been confirmed by present data that do not allow distinguishing populations belonging to D. geniculata from those belonging to D. laetitiae. Indeed, while the cluster including both species is highly supported (PP = 1, Fig. 2), clades including populations from these two species are unresolved. We think that ABC analysis and especially the summary statistic vector could have been driven from this situation and it could be not adequate to model fine scale diversification processes as in this particular case. On this basis, we think that the western Mediterranean lineage (clades A and B in Fig. 2) probably arose from a single colonization event from the North of Italy during the Middle/Late Pliocene (Figs. 3a and 5). An ancestral group spread in the entire inland of Italian Peninsula and from this stock, groups of marginal populations in the Italian coasts could have become isolated at different times, due to the marked climatic and vegetation changes which have occurred since the Pliocene. This result completely agrees with previous studies based on allozyme variation in populations of Dolichopoda belonging to the different Italian species (Allegrucci et al., 1987, 1992, 2005; Sbordoni et al., 1985, 1987). The most western lineage currently occurring in the Pyrenees, in France, probably arose simultaneously to the peninsular Italian ancestral stock, during the Early Pliocene (Figs. 3a and 5).

G. Allegrucci et al. / Molecular Phylogenetics and Evolution 60 (2011) 108–121

117

A

B

Fig. 3. (A) Divergence times among Dolichopoda species inferred by Bayesian analysis using relaxed molecular clocks. Bars at the nodes represent the 95% highest posterior density (HPD) credibility interval. Black circles indicate nodes considered as calibration points. (B) Observed lineage through time plot (black line). Confidence intervals (95%), estimated using 2000 trees from the pool of converged Bayesian trees in BEAST analysis, are also indicated (gray lines).

118

G. Allegrucci et al. / Molecular Phylogenetics and Evolution 60 (2011) 108–121

A

B

C

Fig. 4. Test of alternative biogeographic hypotheses by ABC approach. Three alternative models of species diversification (pruned species trees) were proposed for each of the three evolutionary steps on which we focused. (A) Colonization of Caucasus and Turkey regions (AC-H1-3). (B) Colonization of Greece region (GR-H1-3). (C) Colonization of Western Mediterranean region (WM-H1-3). Node heights correspond to the mean value estimated by the Bayesian chronogram and node labels correspond to those reported in Supplementary Information.

3.4. Mode of diversification The lineage accumulation plot shows sharp increases in the number of lineages from the Late Miocene to Late Pliocene (Fig. 3b). Credibility intervals associated with dates of diversification are broad; nevertheless the date still fits within the general trend for major palaeogeographic changes and climatic shifts. The c-statistics showed a negative value (c = 4.03), rejecting the model of constant diversification rate (p < 0.0001) and suggested a slower net diversification rate in the late history of the genus Dolichopoda. The MCCR test of c-statistics using different number of missing taxa confirmed this result rejecting the null hypothesis of constant rate (p = 0.0001–0.0003). The BDL test (Table 3) showed a positive value of DAICRC (19.39) indicating that rate-variable models (in particular yule2rate) have the best fit to the data. This result was also confirmed by the distribution of simulations that rejected the null hypothesis of a constant diversification model at a = 0.05 showing p = 0 (Rabosky, 2006a). These patterns of rate slowdown are also apparent in several other taxa, such as reptiles (Rawlings et al., 2008; Dolman and Hu-

gall, 2008), birds (Zink and Slowinski, 1995; Weir, 2006) and beetles (Mckenna and Farrell, 2006) and have often been interpreted in terms of niche saturation following adaptive radiation. In the present case, the observed diversification pattern may reflect climate – driven environmental changes and/or palaeogeographic events. After a long-lasting period of relatively poor cladogenetic evolution, the early sudden diversification pattern of Dolichopoda (Fig 3b) could have been affected by palaeogeographic and palaeoclimatic changes occurred during the Messinian salinity crisis. As stated before, our estimates of divergence time suggest that this earlier cladogenesis may be associated with the uplift of the Anatolian Plateau (5–10, Quennell, 1984; Steininger and Rögl, 1984). The divergence of the major clades falls within the Messinian salinity crisis (Hsü, 1972; Hsü et al., 1977), 5.96–5.33 Ma (Krijgsman et al., 1999). The Messinian was a time of high rainfall and high sediment yield rates. This period, named the Zeit Wet Phase, stands in marked contrast to the arid conditions of the preceding Tortonian Stage (Griffin, 1999). The desiccation of the Mediterranean basin during the Messinian salinity crisis, probably allowed broad land connections between formerly separated Aegean landmasses.

G. Allegrucci et al. / Molecular Phylogenetics and Evolution 60 (2011) 108–121

119

A

B

Fig. 5. (A) Biogeographic history of Dolichopoda genus as inferred by ABC analysis. Supported hypotheses are reported for each evolutionary step (see text for details). From East to West: colonization of Anatolian and Caucasian region (AC-H3); colonization of Greece through two migration routes (GR-H3); colonization and diversification in the western Mediterranean (WM-H1). (B) Results from DIVA analysis. Dispersal events are indicated.

Table 3 Results of fitting diversification models to Dolichopoda species using likelihood from multi-rate variant of the Yule model. Model

Pure birth

Birth–death

DDL

DDX

yule2rate

Parameters

r1 = 0.352

r1 = 0.352 a=0

r1 = 0.977 k = 49.59

r1 = 1.389 x = 0.441

r1 = 0.584 r2 = 0.100 st = 1.351

Ln(L) AIC DAIC

39.217 76.435 19.390

39.217 74.435 21.390

48.858 93.717 2.109

42.589 81.178 14.648

50.912 95.826 0

DDL and DDX, logistic and exponential density-dependent speciation model; yule2rate, multi-rate variant of the Yule model; r1 and r2, net diversification rate (speciation event per million year); a, extinction fraction; st, time of rate shift (Mya); k, parameter in the logistic density-dependent model; x, parameter in the density-dependent exponential model; Ln(L), log-likelihood; AIC, Akaike information criterion; DAIC, change in AIC between model and the best-fit model.

Radiation of the genus Dolichopoda in this area may have been promoted by the combination of low sea level and relatively humid climate, occurring during the Messinian and favoring regional dispersal. The rise of sea level occurred at the end of the Messinian salinity crisis seems to have promoted the speciation of the main lineages. Finally, our results would indicate that diversification rate reached a plateau in the last million years. We do not emphasize such an outcome as we believe it may suffer from the following analytical bias. As the chronogram in Fig. 3a was constructed assuming a Yule process, which does not correctly model recent diversification processes, we did not include data at population

level in our analysis. In fact, independent molecular data obtained on the western species (Sbordoni et al., 1985; Allegrucci et al., 1992, 1997) suggested a considerable degree of genetic structuring among populations indicating that especially in the youngest western species diversification is still going on. 4. Conclusions We wish to emphasize two main issues of this study. First, this study is based on a robust taxon sample of cave crickets representative of the whole range of Dolichopoda, attaining 90% of the known species. Second, the study organisms have been the object of several researches developed across many years, involving several aspects of the life cycle (Carchini et al., 1983; Sbordoni et al., 1987), ecology and genetic studies at the population level (Sbordoni et al., 1985, 2000; Allegrucci et al., 1982, 1987, 1992, 1997; Cesaroni et al., 1997). These features add value and robustness to the scenarios analyzed and proposed in the present paper. On the basis of results so far discussed we could attempt to track the tempo and mode of species diversification and colonization in Dolichopoda cave crickets. Both phylogenetic reconstruction (Figs. 2 and 3a) and biogeographic analysis (DIVA and ABC; Fig. 5) suggested an eastern origin of Dolichopoda species, indicating a strongly supported phylogeographic pattern and suggesting that the temporal sequence of cladogenetic events is in good agreement with the geographic distribution of the studied taxa. Climatic events during the Plio/Pleistocene may be responsible of the speciation within each phylogeographic unit, principally driven by

120

G. Allegrucci et al. / Molecular Phylogenetics and Evolution 60 (2011) 108–121

vicariance events. However, also adaptation to cave life seems to have played an important role in this process. The sylvicolous ancestors of Dolichopoda might have used caves as refugia during the unfavorable climatic conditions, beginning their adaptation to subterranean habitat. Therefore, the current distribution of Dolichopoda can be explained by a combination of both vicariance and dispersal events, with many processes occurring in ancestral epigean populations before the invasion of the subterranean environment. Our application of ABC analysis at the interspecific level had also a testing purpose. This method has been mainly used at the intraspecific level or more recently in comparing population histories of co-distributed species (e.g. Fagundes et al., 2007; Carnaval et al., 2009). In our work, the use of different statistics together with ABC allowed a comparative test on the efficiency of this new methodology at higher taxonomic level in biogeographic issues. Results coming from different analyses (e.g. DIVA, Deep coalescences) provided in this paper as well as in Allegrucci et al. (2009) support the ABC approach. Therefore, the ABC approach could be used in those situations where taking into account incomplete lineage sorting maybe an important issue. Moreover, this can be achieved avoiding the time-consuming, and often not affordable, estimate of the likelihood (Beaumont et al., 2002). In the present case, ABC was used as a model-choice methodology and helped to describe the evolutionary pattern of the genus Dolichopoda. Our results encourage biogeographers to further test and implement ABC analysis in studying complex histories at the species and genus level. Further implementations need to be mainly directed towards the development of sufficient summary statistics. Although it still requires to be more widely tested, this approach can be considered as a promising analytical tool. Acknowledgments We are grateful to Donatella Cesaroni, Anna Fabiani and Gabriele Gentile who provided useful criticisms on a previous version of the paper and to Mauro Rampini for providing samples from Caucasus region. We also like to thank Paolo Gratton, Saverio Vicario and two anonymous reviewers for their help and constructive criticsms about ABC analysis. This research was financially supported by PRIN (the Italian Program for Relevant National Researches) grant to Valerio Sbordoni. Appendix A. Supplementary material Supplementary data associated with this article can be found, in the online version, at doi:10.1016/j.ympev.2011.04.002. References Akaike, H., 1973. Information theory and an extension of the maximum likelihood principle. In: Petrov, B.N., Csaki, F. (Eds.), Proceedings of 2nd International Symposium on Information Theory Budapest. Akademia Kiado, pp. 267–281. Allegrucci, G., Caccone, A., Cesaroni, D., Cobolli Sbordoni, M., De Matthaeis, E., Sbordoni, V., 1982. Natural and experimental interspecific hybridization between populations of Dolichopoda cave crickets. Experientia 38, 96–98. Allegrucci, G., Cesaroni, D., Sbordoni, V., 1987. Adaptation and speciation of Dolichopoda cave crickets (Orth. Rhaph.): geographic variation of morphometric indices and allozyme frequencies. Biol. J. Linn. Soc. 31, 151–160. Allegrucci, G., Caccone, A., Cesaroni, D., Sbordoni, V., 1992. Evolutionary divergence in Dolichopoda cave crickets: a comparison of single copy DNA hybridization data with allozymes and morphometric distances. J. Evol. Biol. 5, 121–148. Allegrucci, G., Minasi, M.G., Sbordoni, V., 1997. Patterns of gene flow and genetic structure in cave-dwelling crickets of the Tuscan endemic, Dolichopoda schiavazzii (Orthoptera, Rhaphidophoridae). Heredity 78, 665–673. Allegrucci, G., Todisco, V., Sbordoni, V., 2005. Molecular phylogeography of Dolichopoda cave crickets (Orthoptera, Rhaphidophoridae): a scenario suggested by mitochondrial DNA. Mol. Phylogenet. Evol. 37, 153–164. Allegrucci, G., Rampini, M., Gratton, P., Todisco, V., Sbordoni, V., 2009. Testing phylogenetic hypotheses for reconstructing the evolutionary history of

Dolichopoda cave crickets in the eastern Mediterranean. J. Biogeogr. 36, 1785– 1797. Anderson, C.N.K., Ramakrishnan, U., Chan, Y.L., et al., 2005. Serial SimCoal: a population genetic model for data from multiple populations and points in time. Bioinformatics 21, 1733–1734. Barraclough, T.G., Nee, S., 2001. Phylogenetics and speciation. Trends Ecol. Evol. 16, 391–399. Battagliero, S., Puglia, G., Vicario, S., Rubino, F., Scioscia, G., Leo, P., 2010. An Efficient Algorithm for Approximating Geodesic Distances in Tree Space IEEE/ACM Trans. Comp. Biol. Bioinformatics. IEEE Computer Society, p. 99. Beaumont, M.A., Zhang, W., Balding, D.J., 2002. Approximate bayesian computation in population genetics. Genetics 162, 2025–2035. Bertorelle, G., Benazzo, A., Mona, S., 2010. ABC as a flexible framework to estimate demography over space and time: some cons, many pros. Mol. Ecol. 19, 2609– 2625. Brower, A.V.Z., 1994. Rapid morphological radiation and convergence among races of the butterfly Heliconius hereto inferred from patterns of mitochondrial DNA evolution. Proc. Natl. Acad. Sci. USA 91, 6491–6495. Carchini, G., Rampini, M., Severini, C., Sbordoni, V., 1983. Population size estimates of four species of Dolichopoda in natural and artificial caves of Central Italy (Orth. Raph.). Mém. Biospéol. 10, 341–347. Carnaval, A., Hickerson, M.J., Haddad, C.F.B., Rodrigues, M.T., Moritz, C., 2009. Stability predicts genetic diversity in the Brazilian Atlantic Forest Hotspot. Science 323, 785–789. Cesaroni, D., Matarazzo, P., Allegrucci, G., Sbordoni, V., 1997. Comparing patterns of geographic variation in cave crickets by combining geostatistic methods and Mantel tests. J. Biogeogr. 24, 419–431. Csilléry, K., Blum, M.G., Gaggiotti, O.E., François, O., 2010. Approximate Bayesian Computation (ABC) in practice. TREE 25, 410–418. Degnan, J.H., Rosenberg, N.A., 2009. Gene tree discordance, phylogenetic inference and the multispecies coalescent. Trends Ecol. Evol. 24, 332–340. Dermitzakis, D.M., 1990. Paleogeography, geodynamic processes and event stratigraphy during the late Cenozoic of the Aegean area. In: International Symposium on Biogeographical Aspects of Insularity, Roma 1987. Acc. Naz. Lincei 85, 263–288. Dermitzakis, M.D., Papanikolaou, D.J., 1981. Paleogeography and geodynamics of the Aegean region during the Neogene. In: Seventh International Congress on Mediterranean Neogene, Athens. Ann. Geol. Pays Hellen 30, 245–289. Dolman, G., Hugall, F.H., 2008. Combined mitochondrial and nuclear data enhance resolution of a rapid radiation of Australian rainbow skinks (Scincidae: Carlia). Mol. Phylogenet. Evol. 49, 782–794. Drummond, A.J., Rambaut, A., 2007. BEAST v1.4.7. Bayesian Evolutionary Analysis Sampling Trees (computer program). . Drummond, A.J., Ho, S.H.W., Phillips, M.J., Rambaut, A., 2006. Relaxed phylogenetics and dating with confidence. PLoS Biol. 4, 699–710. Excoffier, L., Laval, G., Schneider, S., 2005. Arlequin ver. 3.0. An integrated software package for population genetics data analysis. Evol. Bioinform. 1, 47–50. Fagundes, N.J.R., Ray, N., Beaumont, M., Neuenschwander, S., Salzano, F.M., Bonatto, S.L., Excoffier, L., 2007. Statistical evaluation of alternative models of human evolution. Proc. Natl. Acad. Sci. USA 104, 17614–17619. Farris, J.S., Källersjö, M., Kluge, A.G., Bult, C., 1994. Testing significance of incongruence. Cladistics 10, 315–319. Farris, J.S., Källersjö, M., Kluge, A.G., Bult, C., 1995. Constructing a significance test for incongruence. Syst. Biol. 44, 570–572. Felsenstein, J., 1981. Evolutionary tree from DNA sequences: a maximum likelihood approach. J. Mol. Evol. 17, 368–376. Friedrich, M., Tautz, D., 1997. An episodic change of rDNA nucleotide substitution rate has occurred during the emergence of the insect order Diptera. Mol. Biol. Evol. 14, 644–653. Gantenbein, B., Largiader, C.R., 2002. Mesobuthus gibbosus (Scorpiones: Buthidae) on the island of Rhodes – hybridization between Ulysses’ stowaways and native scorpions? Mol. Ecol. 11, 925–938. Gomez-Zurita, J., Juan, C., Petitpierre, E., 2000. The evolutionary history of the genus Timarcha (Coleoptera, Chrysomelidae) inferred from mitochondrial COII gene and partial 16S rDNA sequences. Mol. Phylogenet. Evol. 14, 304–317. Gratton, P., Konopinski, M.K., Sbordoni, V., 2008. Pleistocene evolutionary history of the Clouded Apollo (Parnassius mnemosyne): genetic signatures of climate cycles and a ‘time-dependent’ mitochondrial substitution rate. Mol. Ecol. 17, 4248–4262. Griffin, D.L., 1999. The late Miocene climate of northeastern Africa: unravelling the signals in the sedimentary succession. J. Geol. Soc. 156, 817–826. Gu, X., Fu, Y.X., Li, W.H., 1995. Maximum likelihood estimation of heterogeneity of substitution rate among nucleotide sites. Mol. Biol. Evol. 12, 546–557. Ho, S.Y.W., Larson, G., 2007. Molecular clocks: when times are a-changin’. Trends Genet. 22, 79–83. Ho, S.Y.W., Phillips, M.J., Cooper, A., Drummond, A.J., 2005. Time dependency of molecular rate estimates and systematic overestimation of recent divergence times. Mol. Biol. Evol. 22, 1561–1568. Hsü, K.J., 1972. Origin of saline giants: a critical review after the discovery of the Mediterranean evaporates. Earth-Sci. Rev. 8, 371–396. Hsü, K.J., Montadert, L., Bernoulli, D., Cita, M.B., Erickson, A., Garrison, R.E., Kidd, R.B., Mèlierés, F., Müller, C., Wright, R., 1977. History of the Messinian salinity crisis. Nature 267, 399–403. Huelsenbeck, J.P., Crandall, K.A., 1997. Phylogeny estimation and hypothesis testing using maximum likelihood. Ann. Rev. Ecol. Syst. 28, 437–466.

G. Allegrucci et al. / Molecular Phylogenetics and Evolution 60 (2011) 108–121 Huelsenbeck, J.P., Ronquist, F., 2001. MrBayes: Bayesian inference of phylogenetic trees. Bioinformatics 17, 754–755. Kasapidis, P., Magoulas, A., Mylonas, M., Zouros, E., 2005. The phylogeography of the gecko Cyrtopodion kotschyi (Reptilia: Gekkonidae) in the Aegean archipelago. Mol. Phylogenet. Evol. 35, 612–623. Kimura, M., 1980. A simple method for estimating evolutionary rates of base substitutions through comparative studies of nucleotide sequences. J. Mol. Evol. 16, 111–120. Kiyoshi, T., Sota, T., 2006. Differentiation of the dragonfly genus Davidius (Odonata: Gomphidae) in Japan inferred from mitochondrial and nuclear gene genealogies. Zool. Sci. 23, 1–8. Knowles, L.L., 2009. Statistical phylogeography. Ann. Rev. Ecol. Evol. Syst. 40, 593– 612. Krijgsman, W., Hilgen, F.J., Ray, I., Sierro, F.J., Wilson, D.S., 1999. Chronology, causes and progression of the Messinian salinity crisis. Nature 400, 652–655. Lanave, C., Preparata, C., Saccone, C., Serio, G., 1984. A new method for calculating evolutionary substitution rates. J. Mol. Evol. 20, 86–93. Laval, G., Excoffier, L., 2004. Simcoal v. 2.0: A program to simulate genomic diversity over large recombining regions in a subdivided population with a complex history. Bioinformatics 20, 2485–2487. Lipparini, T., 1976. Per la storia del popolamento delle isole dell’Arcipelago Toscano (contributo geo-paleontologico) Lav. Soc. Ital. Biogeogr. 5, 13–25. Liu, L., Yu, L., Kubatko, L., Pearl, D.K., Edwards, S.V., 2009. Coalescent methods for estimating phylogenetic trees. Mol. Phylogenet. Evol. 53, 320–328. Maddison, W.P., 1997. Gene trees in species trees. Syst. Biol. 46, 523–536. Maddison, W.P., Knowles, L.L., 2006. Inferring phylogeny despite incomplete lineage sorting. Syst. Biol. 55, 21–30. Martinsen, L., Venanzetti, F., Johnsen, A., Bachmann, L., 2008. Phylogeography and mitochondrial DNA divergence in Dolichopoda cave crickets (Orthoptera, Rhaphidophoridae). Hereditas 146, 33–45. Mazzanti, R., 1983. Il punto sul Quaternario della fascia costiera dell’Arcipelago Toscano. Bollet. Soc. Geo. Ital. 52, 419–556. McKenna, D.D., Farrell, B.D., 2006. Tropical forests are both evolutionary cradles and museums of leaf beetle diversity. Proc. Natl. Acad. Sci. USA 103, 10947–10951. Nee, S., May, R.M., 1994. The reconstructed evolutionary process. Philos. Trans. Roy. Soc. London Ser. B 344, 305–311. Nee, S., Harvey, P.H., Mooers, A.Ø., 1992. Tempo and mode of evolution revealed from molecular phylogenies. Proc. Natl. Acad. Sci. USA 89, 8322–8326. Nylander, J.A.A., 2004. MrModeltest v.2. Program Distributed by the Author. Evolutionary Biology Centre, Uppsala University. Papadopoulou, A., Anastasiou, I., Vogler, A.P., 2010. Revisiting the insect mitochondrial molecular clock: the Mid-Aegean trench calibration. Mol. Biol. Evol. 27, 1659–1672. Parmakelis, A., Pfenninger, M., Spanos, L., Papagiannakis, G., Louis, C., Mylonas, M., 2005. Inference of a radiation in Mastus (Gastropoda, Pulmonata, Enidae) on the island of Crete. Evolution 59, 991–1005. Parmakelis, A., Stathi, I., Chatzaki, M., Simaiakis, S., Spanos, L., Louis, C., Mylonas, M., 2006a. Evolution of Mesobuthus gibbosus (Brullé, 1832) (Scorpiones: Buthidae) in the northeastern Mediterranean region. Mol. Ecol. 15, 2883–2894. Parmakelis, A., Stathi, I., Spanos, L., Louis, C., Mylonas, M., 2006b. Phylogeography of Iurus dufoureius (Brullé, 1832) (Scorpiones, Iuridae). J. Biogeogr. 33, 251–260. Pons, J., Vogler, A.P., 2005. Complex pattern of coalescence and fast evolution of a mitochondrial rRNA pseudogene in a recent radiation of tiger beetles. Mol. Biol. Evol. 22, 991–1000. Poulakakis, N., Lymberakis, P., Antoniou, A., Chalkia, D., Zouros, E., Mylonas, M., Valakos, E., 2003. Molecular phylogeny and biogeography of the wall-lizard Podarcis erhardii (Squamata: Lacertidae). Mol. Phylogenet. Evol. 28, 38–46. Poulakakis, N., Lymberakis, P., Tsigenopoulos, C.S., Magoulas, A., Mylonas, M., 2005. Phylogenetic relationships and evolutionary history of snake-eyed skink Ablepharus kitaibelii (Sauria: Scincidae). Mol. Phylogenet. Evol. 34, 245–256. Pybus, O.G., Harvey, P.H., 2000. Testing macro-evolutionary models using incomplete molecular phylogenies. Proc. Roy. Soc. London Ser. B 267, 2267– 2272. Pybus, O.G., Rambaut, A., 2002. GENIE: estimating demographic history from molecular phylogenies. Bioinformatics 18, 1404–1405.

121

Quennell, A.M., 1984. The western Arabia rift system. In: Dixon, J.E., Robertson, A.H.F. (Eds.), The Geological Evolution of the Eastern Mediterranean. Geol. Soc. Spec. Publ. No. 17. Blackwell Science, Oxford, pp. 775–778. Rabosky, D.L., 2006a. Likelihood methods for detecting temporal shifts in diversification rates. Evolution 60, 1152–1164. Rabosky, D.L., 2006b. LASER: a maximum likelihood toolkit for detecting temporal shifts in diversification rates. Evol. Bioinform. Online 2, 257–260. Rambaut, A., Drummond, A.J., 2007. TRACER v.1.5 (computer program). . Rawlings, L.H., Rabosky, D.L., Donnellan, S.C., Hutchinson, M.N., 2008. Python phylogenetics: inference from morphology and mitochondrial DNA. Biol. J. Linn. Soc. 93, 603–619. Ronquist, F., 1996. DIVA version 1.1. Computer program and manual available by anonymous FTP from Uppsala University (ftp.uu.se or ftp.systbot.uu.se). Rosenberg, N.A., Tao, R., 2008. Discordance of species trees with their most likely gene trees: the case of five taxa. Syst. Biol. 57, 131–140. Sbordoni, V., Allegrucci, G., Cesaroni, D., Cobolli Sbordoni, M., De Matthaeis, E., 1985. Genetic structure of populations and species of Dolichopoda cave crickets: evidence of peripatric divergence. In: Sbordoni, V. (Ed.), Genetics and ecology in contact zones of populations. Boll. Zool. 52, 95–114. Sbordoni, V., Allegrucci, G., Caccone, A., Carchini, G., Cesaroni, D., 1987. Microevolutionary studies in Dolichopodinae cave crickets. In: Baccetti, B. (Ed.), Evolutionary Biology of Orthopteroid Insects. Horwood Ltd. Publ., Chichester, UK, pp. 514–540. Sbordoni, V., Allegrucci, G., Cesaroni, D., 2000. Population genetic structure, speciation and evolutionary rates in cave dwelling organisms. In: Wilkens, H., Culver, D., Humphrey, B. (Eds.), Ecosystems of the World. Subterranean Ecosystems. Elsevier, Amsterdam, pp. 450–483 (Chapter 24). Sbordoni, V., Allegrucci, G., Todisco, V., 2004. Il genere Dolichopoda in Sardegna: filogenesi molecolare e ipotesi sull’evoluzione del popolamento. Studi Trent. Sci. Nat., Acta Biol 81, 103–111. Shapiro, L.H., Strazanac, J.S., Roderick, G.K., 2006. Molecular phylogeny of Banza (Orthoptera: Tettigoniidae), the endemic katydids of the Hawaiian Archipelago. Mol. Phylogenet. Evol. 41, 53–63. Steininger, F.F., Rögl, F., 1984. Paleogeography and palinspatic reconstruction of the Neogene of the Mediterranean and Paratethys. In: Dixon, J.E., Robertson, A.H.F. (Eds.), The Geological Evolution of the Eastern Mediterranean. Geol. Soc. Spec. Publ. No. 17. Blackwell Science, Oxford, pp. 659–668. Tallmon, D.A., Luikart, G., Beaumont, M.A., 2004. Comparative evaluation of a new effective population size estimator based on approximate Bayesian computation. Genetics 167, 977–988. Thompson, J.D., Gibson, T.J., Plewniak, F., Jeanmougin, F., Higgins, D.G., 1997. The CLUSTAL_X windows interface: flexible strategies for multiple sequence alignment aided by quality analysis tools. Nucleic Acids Res. 24, 4876–4882. Villani, F., Sansotta, A., Cherubini, M., Cesaroni, D., Sbordoni, V., 1999. Genetic structure of natural populations of Castanea sativa in Turkey: evidence of a hybrid zone. J. Evol. Biol. 12, 233–244. Wegmann, D., Leuenberger, C., Excoffier, L., 2009. Efficient approximate bayesian computation coupled with markov chain Monte Carlo without likelihood. Genetics 182, 1207–1218. Wegmann, D., Leuenberger, C., Neuenschwander, S., Excoffier, L., 2010. ABCtoolbox: a versatile toolkit for approximate Bayesian computations. BMC Bioinform. 11, 116. Weir, J.T., 2006. Divergent timing and patterns of species accumulation in lowland and highland neotropical birds. Evolution 60, 842–855. Weisrock, D.W., Macey, J.R., Ugurtas, I.H., Larson, A., Papenfuss, T.J., 2001. Molecular phylogenetics and historical biogeography among salamandrids of the ‘‘true’’ salamander clade: rapid branching of numerous highly divergent lineages in Mertensiella luschani associated with the rise of Anatolia. Mol. Phylogenet. Evol. 18, 434–448. Zink, R.M., Slowinski, J.B., 1995. Evidence from molecular systematic for decreased avian diversification in the Pleistocene epoch. Proc. Natl. Acad. Sci. USA 92, 5832–5835.