Pleistocene phylogeography and cryptic diversity of a tiger beetle, Calomera littoralis, in North-Eastern Mediterranean and Pontic regions inferred from mitochondrial COI gene sequences

Sample collection and identification

In total, 169 imagines of Calomera littoralis were collected with entomological hand net on 43 sites on the Mediterranean coasts of the Balkan Peninsula, Crete and Turkey as well as on the northern and western coast of the Black and Azov Seas, in the years 2009–2012 (Fig. 1 and Table 1). At a site the material was fixed in 96% ethanol for DNA preservation. Taxonomic identification of the collected material followed Mandl (1981)

Molecular data analysis

First, all the obtained sequences were positively verified as Calomera DNA using GenBankBLASTn searches (Altschul et al., 1990). They were then edited and assembled with Clustal W algorithm (Chenna et al., 2003) using Bio Edit © 7.2.5. The resulting alignment was 697 bp long with no gaps, and composed of 169 COI sequences. The sequence data and trace files were uploaded to BOLD and subsequently also to GenBank (accession numbers KU905171–KU905339).

Pairwise Kimura 2-parameter (K2p) distances between sequences were estimated using Mega 6.2 (Tamura et al., 2013). Haplotypes were retrieved using Dna Sp v5 (Librado & Rozas, 2009). Phylogenetic relationships between the haplotypes were visualised with phylogenetic network computed using the neighbour-net algorithm and uncorrected p-distances in SplitsTree ver. 4.13.1 (Huson & Bryant, 2006).

To test for presence of distinct operational taxonomic units (OTUs) that may represent potential cryptic species/subspecies in the sequenced pool of individuals we used the Automatic Barcode Gap Discovery (ABGD) procedure (Puillandre et al., 2012). The default value of 0.001 was used as the minimum allowed intraspecific distance. The maximum allowed intraspecific distance was set to Pmax = 0.03 and 0.06, as both threshold values have been already used in literature to delimit insect species (Hebert et al., 2003; Hebert, Ratnasingham & DeWaard, 2003). We applied the K2p model sequence correction, which is a standard for barcode analyses (Hebert et al., 2003). We used primary partitions as a principal for group definition for they are usually stable over a wider range of prior values, minimise the number of false positive (over split species) and are usually close to the number of groups described by taxonomists (Puillandre et al., 2012).

To reveal the temporal framework for the divergence of the OTUs (potential cryptic species) defined within Calomera littoralis, the time calibrated phylogeny was reconstructed in BEAST, version 1.8.1 (Drummond et al., 2012). A COI sequence of Calomera lugens aphrodisia Baudi di Selve 1864 from GenBank (accession number KC963733) was used as an outgroup. This analysis was performed on a reduced dataset, containing only the most distant haplotypes from each OTU. Hasegawa–Kishino–Yano (HKY) model of evolution, selected as best-fitting to our dataset in MEGA 6.2, and coalescent model were set as tree priors. The strict clock with rate 0.0115, widely used for phylogenetic studies upon insects, was applied for the analyses (Brower, 1994). Five runs of 20 M iterations of Markov chain Monte Carlo (MCMC) sampled each 2000 iterations were performed. The runs were examined using Tracer v 1.6 and all sampled parameters achieve sufficient effective sample sizes (ESS > 200). Tree files were combined using Log-Combiner 1.8.1 (Drummond et al., 2012), with removal of the non-stationary 20% burn-in phase. The maximum clade credibility tree was generated using TreeAnnotator 1.8.1 (Drummond et al., 2012).

To provide insight into historical demography, i.e., the temporal changes of the effective population size of Calomera littoralis in the studied region, we performed Bayesian Skyline Plot (BSP) analysis (Drummond et al., 2005) in BEAST, version 1.8.1 (Drummond et al., 2012). Separate analysis was performed for each of the two phylogenetic lineages revealed in our study (see ‘Results’). The Northern Lineage was represented by 84 individuals from 22 localities, while the Southern Lineage was represented by 85 individuals from 32 localities. The HKY+I model of evolution was used as the best fitting model in case of the Eastern Lineage, while TN93+I was used in case of the Western Lineage. Two runs of MCMC, 20 M iterations long sampled each 2000 iterations, were performed. In both cases the runs were examined using Tracer v 1.6 (Drummond et al., 2012) and all sampled parameters achieved sufficient effective sample sizes (ESS > 200).

Two models of population expansion, demographic and spatial, were examined using mismatch distribution analysis (Slatkin & Hudson, 1991; Rogers & Harpending, 1992) and Tajima’s D neutrality test (Tajima, 1989). Analyses were performed for the COI groups, using Arlequin 3.5.1.3 (Excoffier & Lischer, 2010) with 1,000 replicates.

Morphometric data analysis

To test whether variation of morphometric traits reflects presence of two genetically divergent lineages (potential cryptic species), measurements of eight body parameters (Fig. 2) were taken from all the 69 males and 100 females used previously for the molecular analyses: 1, right mandible length (RML); 2, length of head (LH); 3, width of head (WH); 4, pronotum length (PL); 5, maximum pronotum width (MPW); 6, elytra length (EL); 7, maximum elytra width (MEW); and 8, total body length (TBL). The Principal Component Analysis (PCA) was performed separately for each sex (Fig. 3). To test for significance (p < 0.01) of morphological differences (separately for males and females) between the two divergent lineages one-way ANOSIM Pairwise Test was performed. All the above statistical analyses were done with PRIMER 6 software (Clarke & Gorley, 2006).

Molecular data

A total of 81 haplotypes were identified in the dataset composed of 169 individuals from 43 sites from the Mediterranean and the Pontic areas (Table 1). The phylogenetic network illustrating phylogenetic relationships among haplotypes (Fig. 4) uncovered presence of two distinct haplotype groups (phylogenetic lineages). The first group, from now on defined as southern lineage, includes 36 haplotypes present all over the studied range including the Balkan Peninsula and the Pontic area. The other group, from now on defined as northern lineage, is composed of 45 haplotypes present exclusively along the north-western coast of the Black Sea. The mean K2p genetic distance between both groups of haplotypes is relatively high (0.039, SD 0.007). Both variants of the ABGD analysis resulted in partitioning of the dataset into two OTUs, that may represent distinct operational taxonomic units—potential cryptic species or subspecies within Calomera littoralis in the studied area (Fig. 3A).

The Bayesian time-calibrated reconstruction of phylogeny shows that the two lineages split at ca. 2 Ma, i.e., in early Pleistocene (Fig. 5A). Results of the BSP analyses showing the temporal changes of the effective population size suggests that both lineages experienced rapid population growth that has started ca. 0.15 Ma (Fig. 5B). In both cases, a small decline in effective population size may be observed in most recent times (<0.05 Ma). Results of the mismatch analysis show that both lineages are currently in the stage of both demographic and spatial expansion (Fig. 5C). Interestingly, geographical distribution of both lineages shows that the spatial expansion of southern lineage was efficient enough to spread eastwards into the Black Sea and colonise effectively the north-western Black Sea coast. The northern lineage has spread only in the Pontic region.

Morphometric data

The results of PCA and ANOSIM revealed no differences in the analysed morphometric traits between the southern and the northern lineages, neither in males nor in females (Fig. 3B). In PCA (Fig. 3B), a very weak gradient (R = 0.03) could be seen in case of female body length. Females from the northern lineage clade were slightly larger than those from the southern one (body length; ANOSIM Pairwise Tests p = 0.03).

Cryptic diversity of Calomera littoralis

Known as very important hotspot of biodiversity, endemicity and cryptic diversity (e.g., Myers et al., 2000; Kryštufek & Reed, 2004; Blondel et al., 2010; Huemer & Timossi, 2014; Previšić et al., 2014; Caković et al., 2015), the southern Europe holds also the most diverse tiger beetle fauna in the entire Palearctic realm (Jaskuła, 2011). Presence of cryptic diversity was already pointed out for Cicindela hybrida in the Mediterranean (Cardoso, Serrano & Vogler, 2009) as well as for several species of tiger beetles occurring in other parts of the world (Vogler & Pearson, 1996; López-López, Hudson & Galián, 2012). Thus, existence of well-defined OTUs within Calomera littoralis is not surprising in the studied area. The level of divergence, 0.04 K2p distance, between the northern and the southern lineage is similar as those found between species of tiger beetles in other studies (e.g., Cardoso & Vogler, 2005; López-López, Abdul Aziz & Galián, 2015). Interestingly, we could not detect any conclusive morphological differences between the two lineages based on the multivariate analysis of eight morphometric traits. It must be mentioned that three subspecies of Calomera littoralis, described on the basis of morphology, were reported from the studied area: C. l. nemoralis from all the studied Balkan countries, Crete, Moldova, western Ukraine and western Turkey; C. l. conjunctaepustulata (Doktouroff, 1887) from the Azov Sea area; C. l. winkleri (Mandl, 1934) from Crete and the coastal zone of southern Turkey (Werner, 1991; Putchkov & Matalin, 2003; Avgin & Özdikmen, 2007). However, the morphological differences between the subspecies, such as body size, maculation of elytra and shape of aedeagus, are poorly defined and did not allow the identification of the studies material further than to the species level. Unfortunately, we had no opportunity to study the topotypical material—Provence, France, is locus typicus for C. l. nemoralis, Tibet for C. l. conjunctaepustulata, and Cyprus for C. l. winkleri. Thus, we cannot exclude a possibility that the two lineages we found in our material overlap with any of the above mentioned subspecies. However, only a further taxonomic revision combining more phenotypic traits, including e.g., cuticle ultrastructure, with several, mitochondrial and nuclear DNA data, could help to resolve this problem. Until such revision is done, we propose to use the tentative name “Calomera littoralis complex” for populations from the studied area.

Phylogeography of Calomera littoralis

Occurrence of C. littoralis in Europe is restricted predominantly to marine shorelines with sandy beaches and salt marshes as main habitats (e.g., Franzen, 2006; Jaskuła, 2011; Serrano, 2013). In the eastern Mediterranean it is distributed continuously all along the Adriatic and Aegean coasts, Turkish Straits and the Black Sea coastline (Cassola & Jaskuła, 2004; Jaskuła, Peśić & Pavicević, 2005; Franzen, 2006; Jaskuła, 2007a; Jaskuła, 2007b). However, pronounced genetic structure with two divergent operational taxonomic units (OTUs) implies prolonged spatial isolation in the evolutionary history of this species. The observed level of divergence indicates that this isolation initiated an allopatric speciation. Their present distribution i.e., sympatry in the Pontic region reveals secondary contact of the already divergent lineages in this area. The Bayesian time-calibrated reconstruction of phylogeny shows that split between these OTUs begun in early Pleistocene. This coincides with beginning of recurrent glaciations resulting in eustatic sea level changes and climate aridisation that ever since dominated the global climate and landscape/habitat distribution (Fagan, 2009). In the Mediterranean and in the Pontic region such global effects overlaid and strengthen the local effects of tectonic plate collision leading to Alpine orogeny, i.e., local land uplift and subsidence resulting in isostatic sea level changes, salinity fluctuations from freshwater to fully marine and habitat mosaicism (Stanley & Blanpied, 1980). For example, during that time the connections of Pontic basin to Mediterranean Sea was lost and regained for more than a dozen times (Kerey et al., 2004; Badertscher et al., 2011). A profound impact of these events on the evolution and, hence, distribution of local both aquatic (Audzijonyte, Daneliya & Vainola, 2006; Nahavandi et al., 2013) and terrestrial taxa (e.g., Böhme et al., 2007; Ferchaud et al., 2012). We can assume that in case of C. littoralis, a halophilic species bound to coastal habitats, sea level fluctuations would significantly affect its distribution. The 2 Ma divergence time for C. littoralis OTUs derived from our data coincides with one particular disconnection of the Mediterranean and Pontic basins. At that time, from ca. 2 to ca. 1.5 Ma, the Meothic Sea, one of several predecessors of the Black Sea, turned into the predominantly freshwater Pontos Sea/Lake (Grinevetsky et al., 2015). This surely broke the formerly continuous stretch of coastal habitats connecting the two basins and thus, could be an effective barrier leading to split of C. littoralis population into the allopatric southern and northern lineages. Their detailed history is impossible to unravel, yet results of BSP analyses reconstructing past changes in effective population size indicate that both lineages started their demographic expansions at ca. 0.15 Ma. This coincides with the terminal stage of MIS-6, i.e., Wartanian/Saalian glaciation, and beginning of MIS-5e, i.e., Eemian interglacial (Lisiecki & Raymo, 2005; Marks, 2011). The latter was characterized by warmer climate and sea level higher by 6–9 m in comparison to Holocene (Kopp et al., 2009; Dutton & Lambeck, 2012). In result, a wide connection between Mediterranean and the Pontic basin was re-established and the coastal habitats extended again, enabling exchange of faunas. Due to a deficiency of local studies, it is hard to compare our results to evolutionary history of any other terrestrial taxa in the area. However, a wealth of studies showing very similar spatiotemporal scenario in animal taxa comes from the coastal regions of the Gulf of Mexico and the adjacent Atlantic coast (summarised by Avise, 1992). During Pleistocene, Cuba was connected with a land bridge to the Florida Peninsula what lead to divergence of populations of several terrestrial and aquatic animals, including also a local tiger beetle species Cicindela dorsalis Say, 1817 (Vogler & DeSalle, 1993). Interestingly enough, however, according to our results both lineages are until now in the stage of demographic and spatial expansion, only the southern one has crossed the present Turkish straits. This asymmetry is hard to explain. Another interesting fact is that the isolation of Pontic basin from Mediterranean during the following Weichselian glaciation did not have probably any effect on the demography and phylogeography of the species. Based on the mitochondrial DNA marker only we cannot also conclude, whether the secondary contact of the divergent lineages effected in hybridization and or introgression. Answering this question requires employment of nuclear marker, what leaves a space for the future studies—much wider in terms of geographic coverage and molecular markers used.

Concluding, we have demonstrated that Pleistocene glaciations and associated sea level changes in the Mediterranean/Pontic region had a profound effect on the genetic diversity and distribution of widely distributed coastal insect species, generating some level of cryptic diversity. Our case study casts more light on the evolutionary relationships between populations of terrestrial animals inhabiting both the Mediterranean and Black Sea shorelines—a phenomenon that is still weakly explored in literature.