Molecular phylogeny of the Palaearctic butterfly genus Pseudophilotes (Lepidoptera: Lycaenidae) with focus on the Sardinian endemic P. barbagiae

Background: The Palaearctic butterfly genus Pseudophilotes Beuret, 1958 (Lycaenidae: Polyommatinae), that today occurs in North Africa and in Eurasia, includes ten described species with various distribution ranges, including endemics such as the Sardinian P. barbagiae. Phylogenetic relationships among these species are largely unresolved. In the present study, we analysed 158 specimens, representing seven species out of ten described in the genus, from widely distributed sites throughout the western Palaearctic region, using nuclear markers (28S rRNA, wingless, Internal Transcribed Spacer 2 and Elongation Factor 1α) and the barcoding region of the mitochondrial cytochrome c oxidase subunit 1 gene to investigate if the current taxonomic entities match the phylogenetic pattern. Further, we attempt to infer the geographic origin of the genus Pseudophilotes and estimate the timing of its radiations, including the split of the Sardinian endemic P. barbagiae. Results: Maximum Likelihood and Bayesian inference analyses confirmed the monophyly of the genus Pseudophilotes and clearly supported the closer affinity of P. barbagiae to the species assemblage of P. baton, P. vicrama and P. panoptes as opposed to P. abencerragus. The currently accepted species P. baton, P. vicrama and P. panoptes turned out to be weakly differentiated from each other, while P. bavius and P. fatma emerged as highly distinct and formed a well supported clade. The split between the lineage comprising bavius and fatma (sometimes treated as a distinct genus, Rubrapterus) with Salvia species as larval host plants, and the remaining Pseudophilotes that utilize Thymus and other Lamiaceae (but not Salvia), dates back to about 4.9 million years ago (Mya). Conclusions: Our results show that the last common ancestor of the genus probably lived in the Messinian period (5.33–7.25 Mya). At species level, they support the current taxonomy of the genus, although P. panoptes, P. baton and P. vicrama display complex patterns based on phylogeographic relationships inferred from mtDNA. The Sardinian endemic P. barbagiae turned out to be a young endemic, but clearly with European instead of North African origin and evolved through allopatric isolation on the island of Sardinia only about 0.74 Mya.


Background
Pseudophilotes Beuret, 1958 (Lycaenidae: Polyommatinae) is a Palaearctic genus of blue butterflies mainly distributed in Eurasia, but also in North Africa and the Levante (Fig. 1a). The phylogenetic position of the genus within the Glaucopsyche section (Polyommatinae) has been lately inferred by Ugelvig et al. [1]. They found that the genera Pseudophilotes and Scolitantides are clearly separated, although they had been lumped together by other studies [2,3]. This separation is explained by the affiliation with different larval host plants: Pseudophilotes with Lamiaceae, and Scolitantides with Crassulaceae. a b Fig. 1 a Sampling localities and approximate geographic distribution areas of the Pseudophilotes species considered in the present study. Sampling localities for the species included in this study are showed with differently coloured circles, diamonds, and triangles respectively. Geographic distribution areas are redrawn after [2,3,26,[74][75][76][77] and shaded in colours referring to the sample localities. The map was prepared using Quantum GIS 2.8.2 based on a map from Natural Earth (www.naturalearthdata.com). b Median-Joining Network of Pseudophilotes COI sequences. The circle size is proportional to haplotype frequency and the different colours refer to the different species used in the study (cf. a). Numbers of mutations between haplotypes are shown at the connections, except for single or double substitutions. Connections creating loops that are less probable according to frequency and geographical criteria are indicated in dashed lines. Haplotype codes match those in Additional file 1. A1 -A3: haplogroups including P. baton, P. vicrama and P. panoptes At the species level, the same study [1] showed that the Bavius Blue, Pseudophilotes bavius forms a sister clade to the rest of Pseudophilotes, which is also supported by its association to Salvia (Lamiaceae), whereas species in the remaining clades mainly feed on Thymus (Lamiaceae). Recent efforts to resolve the evolution of Pseudophilotes based on the COI barcoding region [4] indicated a closer relationship of the Sardinian Blue, P. barbagiae to the False Baton Blue, P. abencerragus than to the Baton Blue, P. baton. Further studies based on mtDNA [5] showed P. baton and the Eastern Baton Blue, P. vicrama to display complex patterns. Species delimitation in the genus remains thus a work-in-progress.
Species such as the Sardinian Blue, Pseudophilotes barbagiae [6], with highly restricted distributions and taxa representing unique lineages within a species flock can be suitable models to investigate the processes of differentiation and speciation and often provide the key for answers to broader-scale biogeographic questions [7]. In this respect, the study of islands where many endemic species or subspecies evolved has provided fundamental insights into the relationships between geographical patterns and evolutionary processes [8,9]. P. barbagiae has been described as distinct from the continental and Corsican P. baton by characteristics of the male genitalia and wing markings. Its divergence has been hypothetically linked to the geographic isolation of the Sardo-Corsican plate [10], but this assumption has never been tested on the basis of molecular data. Gerace [11] hypothesized that P. barbagiae probably diverged from P. baton as a result of the separation of Corsica and Sardinia due to important climatic changes, and later P. barbagiae became isolated in mountainous areas.
Due to its geographic position and geological history, the island of Sardinia is characterized by a remarkable richness of endemic species [12] and represents one of the most prominent biodiversity hotspots in the Mediterranean basin [13]. The question how its endemic fauna was formed and when it originated has been previously investigated using molecular analyses for a variety of vertebrate [14][15][16] and invertebrate animals [9,[17][18][19][20][21]. In all these cases, phylogenetic patterns could be related to several aspects of the geological history of the Mediterranean basin [11,12], but the specific evolutionary histories turned out to vary distinctly between taxa. About 5 Mya ago, during the Messinian period, the closure of the Strait of Gibraltar led to an almost complete desiccation of the Mediterranean Sea, and to a subsequent formation of salt deposits in all deep basins [22]. This event, called the Messinian Salinity Crisis (MSC), was often advocated as one of the main drivers for the emergence of a distinctive fauna and flora in the Mediterranean. The newly exposed land corridors permitted faunal exchange between southern Europe and the Maghreb region (Northern Africa) [23,24], as well as between Sardinia and northern Italy and southern France.
When the connection to the Atlantic ocean was re-established, the massive influx of water brought about a drastic change in marine flora and fauna, and at the same time terminated the connections between the Tyrrhenian islands and the Eurasian or African continental mainland, respectively, which had as a consequence the differentiation of many distinctly younger island species (e.g. butterflies, cave-beetles, mayflies and lizards). Land bridges between Sardinia and Corsica and the mainland continued to be created well into the Pleistocene due to glacial/interglacial sea-level oscillations [25]. During the Würm glaciation (115,000-11,700 years ago), Sardinia could have been in contact with the mainland via the island of Elba as the sea level was up to 120 m lower than today [11,12] and with the neighbouring Corsica.
In the present study, we analysed four nuclear genes (28S ribosomal [28S], wingless [wg], Internal Transcribed Spacer 2 [ITS2] and Elongation Factor 1α [EF1α]) and the barcoding region of the mitochondrial cytochrome c oxidase subunit 1 (COI) gene for seven out of ten Pseudophilotes species, from widely distributed sites throughout the Palaearctic region, including for the first time the Sardinian endemic P. barbagiae, as well as the African P. fatma. We aimed to (i) reconstruct the phylogeny of the genus Pseudophilotes and check if the obtained clades correspond to current taxonomic units, (ii) to trace the geographical origin of the genus since extant distributions of Pseudophilotes species cover both Eurasia and North Africa, and (iii) to shed light on the historical biogeographic processes that determined the present-day isolation of the Sardinian endemic P. barbagiae.

Taxon sampling and molecular procedures
We examined 96 Pseudophilotes individuals representing seven species out of ten described from sites across Europe and North Africa (Maghreb region) ( Fig. 1a and Additional file 1). We were unable to obtain sequences from P. sinaicus, P. jordanicus and P. jacuticus. We included sequences from the Eurasian species Scolitantides orion Pallas, 1771 in our analysis, since the Pseudophilotes species complex is frequently attributed to represent just a clade within a more inclusive genus Scolitantides Hübner, (1819) [2,26].
Specimens had either been dried or stored in 99% ethanol after collection. Subsequently, samples were stored at − 20°C until DNA extraction. The selected mitochondrial (mt) and nuclear (nc) genetic markers (COI, 28S, wg, ITS2, EF1α) were amplified with varying success. When necessary, new specific PCR primers were used (see Additional file 2), as the specimens had been collected between the years 2000 and 2014; older samples performed often worse in Polymerase Chain Reaction (PCR) than the more recent ones.
Extraction and amplification of mt and nc DNA was conducted in two different laboratories, the Canadian Centre for DNA Barcoding (CCDB) (University of Guelph, Guelph, Canada) and the Department of Botany and Biodiversity Research (University of Vienna, Austria), according to the following procedures. In the CCDB lab, total DNA was extracted from one leg of each individual on Biomek FX liquid handling robot using a semi-automated DNA extraction protocol [27] on glass fiber plates (PALL Acroprep 96 with 3 μm GF membrane over 0.2 μm Bioinert membrane). DNA was eluted in 35-40 μl of ddH 2 0 pre-warmed to 56°C and stored at − 80°C. PCR and sequencing for mtDNA were carried out following standard procedures for Lepidoptera [28,29]. A fragment of 658 bp of the mt COI gene corresponding to the barcoding region was obtained. In the Vienna lab, total DNA was extracted from legs of butterflies homogenized with ceramic beads using a Precellys 24 homogenizer set to 5000 min-1 for 2x20s. DNA extraction was performed according to the standard protocol supplied with the Analytik Jena innuPREP DNA Micro Kit (Analytik Jena AG). Four nuclear gene fragments were amplified using the Fermentas PCR system. PCR reactions were set up as follows: 2.5 μl of 10× (NH4)2SO4 PCR buffer, 2.5 μl MgCl2 (25 mM/l), 0.5 μl dNTPs (10 mM/μl), 0.5 μl of each primer (10 pg/μl), 0.2 μl BSA, 1 μl template DNA, 0.2 μl Taq polymerase and filled to 25 μl with PCR-grade H2O. PCR reactions were purified by digestion with shrimp alkaline phosphatase and exonuclease for 15 min at 37°C followed by 15 min at 80°C for enzyme deactivation. Sequencing reactions were set up with 1 μl ABI BigDye 3.1 (Applied Biosystems, Carlsbad, CA, USA), 1 μl primer, 1.5 μl sequencing buffer, 1 μl template DNA and filled to 10 μl with PCR grade H2O. Primer names, references, primer sequences as well as respective annealing temperature and time are shown in Additional file 2. Sequences were obtained by using an ABI capillary sequencer following the manufacturer's recommendations. All gene fragments were sequenced in both directions and they were edited and assembled using CodonCode Aligner 6.0.2. All sequences were submitted to GenBank (Accession Numbers in Additional file 1) and in BOLD system repository (dataset name: DS-PSEUDOPH; DOI: dx.doi.org/10.5883/DS-PSEUDOPH).

Data set compilation and tree reconstruction procedures
To reconstruct the phylogenetic relationships of Pseudophilotes, the 96 specimens sequenced in our labs were complemented by 62 publicly available specimens of Pseudophilotes and 53 outgroup species of the subfamily Polyommatinae and two species of the subfamily Lycaeninae, mainly selected from Vila et al. [30] (Additional file 1).
All mt and nc protein coding genes were aligned using CodonCode Aligner 6.0.2. The nc ribosomal 28S rRNA gene fragment and the nc ITS2 sequences were aligned with the software MAFFT v7.245, using the iterative refinement method G-INS-i [31] via the MAFFT online server (http://mafft.cbrc.jp/alignment/server/). As the nc ITS2 region consists of highly variable sections, its alignment remained partly ambiguous. We therefore used the software Aliscore v.2.0 [32] to identify ambiguously aligned or randomly similar sections (RSS) within the nc ITS2 alignment. Aliscore was applied with the following conditions: gaps were treated as ambiguous characters, a window size of six positions was selected and the maximum number of pairs was compared to check for replications. Taken together, the resulting alignments of all mt and nc target fragments consisted of 3929 bp.
Haplotype diversity and the number of variable sites of mt and nc target fragments were calculated using DnaSP 5.0 [33]. To represent the genealogical relationships among mtDNA haplotypes we calculated a Median Joining (MJ) network using NETWORK 5 [34].
Pseudophilotes relationships were assessed via Maximum Likelihood (ML), including all mt and nc target fragments (3929 bp), and using the tree reconstruction method IQ-TREE v.1.61 [35]. Prior to analysis, the data set was subdivided into eleven partitions, according to codon positions of the protein coding genes and the single ribosomal genes. To reduce over-parameterization and increase model fit, these partitions were further merged if possible and for each of the remaining partitions, the best fitting model was estimated within IQ-TREE, which resembles the partitioning algorithm of the software PartitionFinder v2.0 [36]. Node support was assessed by ultrafast bootstrap (UFBoot) analysis [37] with 1000 replications.

Multilocus coalescent species delimination
We tested the status of all Pseudophilotes species according to current taxonomy with the multilocus Bayesian Phylogenetics and Phylogeography software (BPP v4.0; [38]). Within BPP we employed the joint species delimitation and species-tree inference approach, for which the species tree inferred from the ML tree reconstruction was used as a guide tree. Prior to these analyses, we excluded all outgroups except Glaucopsyche, Philotes and Scolitantides, which resulted in a condensed data set, comprising 169 specimens and 3253 bp. BPP now implements inverse gamma priors for the parameters theta (θ) and tau (τ). This allows integrating out θ estimates analytically, which helps with moves between trees and delimitation models. We used an inverse gamma prior with mean 0.002 for θ and 0.001 for τ. The analyses were run for 200,000 generations with a burn-in of 10,000 generations and a sample frequency of five and were run three times to confirm consistency between runs.

Divergence times estimations
Divergence times of the species were further assessed using a species tree approach on the condensed data set in *BEAST 1.8.2. [39]. Each Pseudophilotes specimen in the dataset was assigned to one of the potential species delimited with a posterior probability of > 0.95 in BPP. The condensed dataset was further partitioned into the five gene fragments, which were subsequently loaded into PartitionFinder v2.0 to find the best fitting model scheme for each gene. In *BEAST, all substitution, clock and tree models were unlinked. We assigned Brower's rate of 0.0115 substitutions/site/million years [40] to the COI partition and estimated the rates of the other partitions. A lognormal relaxed clock with uncorrelated rates was assigned to each clock model. The Tree model was set to "Speciation: Yule Process". *BEAST analyses were run on the CIPRES Science Gateway v.3.3 [41], consisting of 150 million generations with sampling every 10,000 generations. The first 75 million generations were discarded as burn-in and convergence and mixing of the parameters were checked in Tracer 1.6 (http://tree.bio.ed.ac.uk/software/tracer/). Post-burnin samples were subsequently used to construct a maximum clade credibility tree in TreeAnnotator [42].

Tests of demographic equilibrium and population expansion in the mtDNA haplogroup
To test demographic equilibrium or events of recent expansion in the genus, different sets of mtDNA sequences were selected (Table 1) on a geographical basis and taking into account the results of the previous phylogenetic analysis. Demographic equilibrium was tested by calculating F S [43] and R 2 [44] parameters, which have both been shown to provide sensitive signals of historical population expansions [44]. Arlequin 3.0 [45] and DnaSP 5.0 [33] were used to compute F S (P < 0.02) and R 2 (P < 0.05) respectively, and to test their statistical significance by simulating random samples (10,000 replicates) under the null hypothesis of selective neutrality and constant population size using coalescent algorithms (both modified from Hudson [46]). Expected mismatched distributions and parameters of sudden expansion τ = 2 μ t (τ, sudden demographic expansion parameter; t, true time since population expansion; μ is the substitution rate per gene) were calculated using Arlequin 3.0 by a generalized least-squares approach [47], under models of pure demographic expansion and spatial expansion [48,49]. The probability of the data under the given model was assessed by the goodness-of-fit test implemented in Arlequin 3.0. Parameter confidence limits were calculated in Arlequin 3.0 through a parametric bootstrap (1000 simulated random samples).

Biogeographical analyses
We applied the R software package BioGeoBEARS v.0.21 [50,51] to infer the biogeographical history of Pseudophilotes and to test the support for its origin either in Eurasia or in North Africa. BioGeoBEARS uses the Lagrange DEC model [52,53], as well as likelihood interpretations of the DIVA model [54] and the BayArea model [55]. For each of these biogeographical models, it further implements a parameter describing founderevent speciation (+j). All model setups were compared in a statistical framework allowing the selection of the best-fitting model. The analyses were conducted with the BEAST maximum clade credibility tree (MCC) from the *BEAST analyses, excluding the outgroup species Glaucopsyche spp. This was necessary, as the relationships within this genus are still ambiguous [1] and therefore the unclear primary distribution of Glaucopsyche (Nearctic vs. Palaearctic) would have blurred the biogeographical reconstruction of Pseudophilotes. The following areas were designated (see Fig. 2c): Eurasia (A), Iberian Peninsula (B), North Africa (C), Sardinia (D), and North America (E). The maximum number of ancestral areas was set to three in both sets of analyses.

Phylogenetic reconstruction and coalescent species delimitation
The ML tree reconstruction (see Additional file 3) recovered a clade consisting of Philotes sonorensis and c a b Scolitantides orion as the closest relative to Pseudophilotes (Bootstrap support, BS = 98). Our results further support the monophyly of the genus Pseudophilotes (BS = 96) recognizing two well-supported main clades. One of them comprised P. bavius (distributed from the Balkan Peninsula and Romania eastwards to Turkey, southern Russia and Iran), and the Maghrebian endemic species P. fatma (BS = 100). In the other clade (BS = 100), P. abencerragus (distributed from Israel across the Maghreb into the Iberian Peninsula) was recovered as sister to all remaining species. This latter clade (BS = 96) was split into the Sardinian P. barbagiae (BS = 100) and a clade comprising P. vicrama, P. baton and P. panoptes (BS = 90). Although this clade shows three distinct groups, all species accepted in current taxonomy appeared polyphyletic with some members of P. panoptes scattered in the P. baton clade and members of P. baton appearing in P. vicrama. In contrast, the species delimitation analyses in BPP recovered all different taxonomical units of Pseudophilotes as species with highest posterior probabilities (pP = 1.0).

Divergence time estimations
The relationships inferred with ML were congruent with the results of the *BEAST analyses (Fig. 2a), except for the position of P. vicrama, which appeared as sister group to P. panoptes in the ML analyses (although with low support, Additional file 3) and as sister group to P. baton in the *BEAST analyses (well supported, Fig. 2a). The chronogram in Fig. 2a shows the median ages and 95% highest probability density (95% HPDs) computed from the post burn-in posterior topologies. We estimated the origin of the genus Pseudophilotes at 5.71 Mya (95% HPD: 3.24-9.13 Mya). The split between P. bavius and P. fatma was estimated at 1.65 Mya (95% HPD: 0.79-3.13 Mya). The split between P. abencerragus and the remaining species was estimated at 1.25 Mya (95% HPD: 0.59-2.24 Mya) and finally the emergence of P. barbagiae was estimated at 0.74 Mya (95% HPD: 0.32-1.32 Mya).

Geographic distribution of mtDNA haplotypes
The mtDNA haplotype network (Fig. 1b)  The thirteen sequences of P. barbagiae displayed two haplotypes that were connected to a central unresolved loop with P. abencerragus. Eight and nine mutational steps joined P. barbagiae with P. abencerragus from North Africa and Europe respectively, while thirteen mutational steps joined P. barbagiae with P. abencerragus from Israel. In contrast, 16 mutational steps joined P. barbagiae with the closest haplotype of the A1-A3 haplogroup (including P. panoptes, P. baton and P. vicrama), a result not reflected in the molecular analyses using both mt and nc markers, where P. barbagiae was recovered as sister to the clade including P. panoptes, P. baton and P. vicrama (Fig. 1a, Additional file 3).

Tests of demographic equilibrium
According to both the Fs and R 2 parameters, calculated for six mtDNA sequence sets (Table 1), the null hypothesis of constant population size could be rejected for two phylogeographic units (set in bold in Table 1): P. abencerragus (number of mtDNA sequences (N) = 23) and haplogroup A2 (N = 40) including P. vicrama sequences and P. baton from Corsica and southern Italy (Fig. 1b). The mismatched distribution of these sequence sets was also examined according to a sudden expansion model. Goodness of fit tests indicated significant deviations from expected distributions only for P. abencerragus that showed a multimodal and ragged shape, revealing demographic equilibrium or a stable population [53]. Conversely, for the haplogroup (A2) the parameter τ = 2 μ t could be used to estimate the time (t) elapsed from population expansion. According to Brower's mutation rate [40], the demographic expansion of the haplogroup A2 could be traced back to about 41 thousand years ago (kya) (0−207 kya).

Biogeographical analyses
The results of the ancestral area reconstructions in Bio-GeoBEARS are summarized in Tables 2 and 3. All models, which implemented founder-event speciation parameter (+j), were significantly better supported than their more simple counterparts. The AIC weights further showed the DEC + j model to achieve the highest relative probability (70%; see Table 2). However, the geographical origin of Pseudophilotes remained unclear, as the DEC + j model indicated several potential area combinations, none of these significant (Fig. 2b). Also the ancestral area for the P. abencerragus group (Pseudophilotes s. str.) remained unclear, with Eurasia+Iberia+Africa having the highest probability (60%; see Table 3). More details on the biogeographic modelling are given in the discussion.

Discussion
The taxonomical status of the genus Pseudophilotes and its species Mitochondrial and nuclear phylogenetic analyses (Figs. 1b, 2a, Additional file 3) confirmed the findings of Ugelvig et al. [1] that Pseudophilotes and Scolitantides are closely related, but clearly separated from each other. It therefore seems justified to maintain Pseudophilotes and Scolitantides as distinct genera (e.g. [56]), despite the fact that other studies have joined them [2,3]. Our analyses further corroborate the dichotomy between the clade of the species P. bavius and P. fatma, (equivalent to what has been proposed as the genus Rubrapterus Korshunov, 1987), and Pseudophilotes s. str. This evolutionary split between the Rubrapterus clade and the remaining Pseudophilotes follows larval food plant preferences. P. bavius and fatma feed on Salvia spp. [57,58], whereas Lamiaceae genera other than Salvia (mostly Thymus spp., but also Cleonia, Calamintha, and rarely Cuscuta as a parasite on thymes) are preferred by all other species [3,59]. Within P. bavius, several subspecies have been described in Europe (e.g. hungarica Diószeghy, 1913; macedonicus Schulte, 1958; casimiri Hemming, 1932; see for example Tshikolovets [3]), but our dataset did not reveal any clear spatial genetic differentiation below the species level. P. fatma has sometimes been treated as a subspecies of P. bavius [60,61], but given morphological [62] and genetic differences it is probably best regarded as a distinct species. Among the Thymus-feeders, the southern-mediterranean species P. abencerragus is clearly distinct from all other species. Within P. abencerragus, the sequences obtained from North African and Israeli specimens appear to be derived from those of Iberian individuals, reflecting previously proposed subspecies (North Africa: colonarium Turati, 1927; Israel: nabateus Graves, 1925 and amelia Hemming, 1927). However the support of these relationships is weak, prohibiting final statements on their taxonomic status. The Sardinian endemic P. barbagiae belongs to the baton-group but is clearly separated from the clade comprising P. baton, P. panoptes and P. vicrama. Although species delimitation analyses support the status of all of the latter (Fig. 2a), our tree reconstruction analyses (Additional file 3) and the mtDNA haplotype network (Fig. 1b) indicate a more complex pattern representing a challenging taxonomic case that requires further study to clarify species boundaries and their precise geographic distribution, as well as possible hybrid zones.

Biogeographical history of the genus Pseudophilotes
Our results are consistent with other studies on Polyommatinae [63] indicating that the most recent common ancestor of the genus Pseudophilotes probably lived in the Messinian period . This ancestor had already colonized Lamiaceae as larval hosts, a synapomorphic character differentiating Pseudophilotes from Scolitantides and Philotes; affiliation with Lamiaceae is a rather rare trait in the Polyommatinae worldwide [64,65].
During the Messinian, the closure of the Strait of Gibraltar led to an almost complete desiccation of the Mediterranean Sea, permitting massive faunal exchange between southern Europe and the Maghreb region [23,24]. Within Pseudophilotes s. str., the ancestor of P. abencerragus split into the extant P. abencerragus that occurs in southwestern Europe and from northwestern Africa eastwards through to southern Israel, and the ancestor of P.  baton, P. panoptes, P. vicrama and P. barbagiae. Although climate modelling analyses [66] showed differences in climate requirements that possibly suggest different evolutionary histories for the more western P. baton and its eastern vicariant P. vicrama, the two species are parapatric in Central Europe (Fig. 1a) and all P. baton analysed from populations in Corsica share their mtDNA haplotype (H 32 ) with P. vicrama. The large geographical distance of the Corsican populations from the extant contact zone in eastern central Europe prompts some degree of historical introgression between these two species. This hypothesis is supported both from the distribution of the H 32 haplotype ( Fig. 1a; Additional file 1) and the genetic traces of demographic expansion in P. vicrama (haplogroup A2: Fig. 1b [71] this plant species is not threatened. Our analyses of combined mt and nc markers show that P. barbagiae appears more closely related to the clade including P. baton, P. vicrama and P. panoptes than to P. abencerragus. On the other hand, the haplotype network (Fig. 1b) depicts a somewhat different picture. Here, P. barbagiae results more closely related to P. abencerragus than to P. baton. Such discrepancies between analyses based solely on the COI region and others based on combined data, including nuclear markers, are quite common, since the barcoding region by itself is not sufficient for drawing phylogenetic conclusions and mitochondrial inheritance often follows different pathways that do not necessarily reflect the evolutionary history of species. Our hypothesis is that during the Pleistocene, when Sardinia was connected with Corsica and the mainland due to glacial/interglacial sea-level oscillations [25], gene flow occurred among baton-group European populations and Sardinia. Afterwards, gene flow among these became interrupted and from this founder population P. barbagiae evolved into an endemic species about 0.74 Mya. Interestingly, the mtDNA of P. baton from Corsica is of the vicrama-type (H 32 ) shared with individuals of P. vicrama distributed from Finland to Iran (Fig. 1b;  Additional file 1). An expansion from Africa, as it has been hypothesized for two other Sardinian endemic butterflies (Maniola nurag: Kreuzinger et al. [72] and Papilio hospiton: Pittaway et al. [73]), does not seem a plausible scenario, despite the smaller geographical distance between Sardinia and the African continent.

Conclusions
This study confirmed the monophyly of the genus Pseudophilotes and shed further light on the validity of its component species, as described by current taxonomy. We found evidence for a deeper split between the Salvia feeders (P. bavius and P. fatma) and the Thymus feeders (Pseudophilotes s. str.). However, the precise species boundaries and possible signals of introgression and hybridization within the baton-group further merit closer scrutiny. Finally, P. barbagiae turned out to be a young endemic, but clearly with European instead of North African origin.