A molecular phylogeny of Polyommatus s. str. and Plebicula based on mitochondrial COI and nuclear ITS2 sequences (Lepidoptera: Lycaenidae)

The phylogenetic relationships of the subgenera Polyommatus and Plebicula, within the Palaearctic butterfly genus Polyommatus, were inferred from a combined analysis of the nuclear marker ITS2 and the barcoding section of the mitochondrial gene COI. Eight major clades were recovered within Polyommatus s. l., which correspond closely to subgenera based on traditional systematics and are of late Pliocene to early Pleistocene origin. Extraordinary chromosomal evolution occurred independently in three of these clades. The disputed position of several species formerly placed in the subgenus Plebicula is clarified. A group of Cen- tral Asian species (Bryna) was recovered as a monophyletic clade within Polyommatus s. str. The Kurdistanian endemic P. buzul- mavi appears as a sister species to P. icarus. P. celina replaces P. icarus in NW Africa and the Canary Islands, and split from the last common ancestor with P. icarus back in the early Pleistocene.


INTRODUCTION
Polyommatus Latreille, 1804 is a genus of blue butterflies, which is distributed throughout the Palaearctic region.Its type-species, Polyommatus icarus (Rottemburg, 1775), is the most widespread and common representative and was recently even introduced into Ontario in Canada (Hall, 2007).The delineation of Polyommatus has been contentious.Among the nomenclaturally available genus names listed in Eliot (1973) within his "Polyommatus section" several are currently often regarded as synonyms or subgenera of Polyommatus (e.g.Hesselbarth et al., 1995;Bálint & Johnson, 1997; see Table 1).These are Cyaniris Dalman, 1816, Bryna Evans, 1912, Meleageria de Sagarra, 1925, Agrodiaetus Hübner, 1822, Lysandra Hemming, 1933and Plebicula Higgins, 1969.The same applies to three further (sub-)genera, which were described in 1977: Neolysandra Koçak, Sublysandra Koçak and Paragrodiaetus Rose & Schurian.Kretania Beuret, 1959 was included within Polyommatus by Bálint & Johnson (1997), whereas most other authors (e.g.Hesselbarth et al., 1995) treat Kretania as a subgenus of Plebejus Kluk, 1802.Bálint & Johnson (1997) also include the genera Glabroculus Lvovsky, 1993 (= Elviria Zhdanko, 1994) and Rimisia Zhdanko, 1994within Polyommatus. Gorbunov (2001) synonymized Polyommatus with Plebejus, but this is not accepted by most current authors.It should be noted that no clear synapomorphic morphological features have been found that could be used to delineate the genus Polyommatus, although Bálint & Johnson (1997) tried to give a diagnosis.Also no phylogenetic analyses have been attempted based on morphological characters, and all currently available systematic treatments are mostly based on the intuition of the authors.Indicative of the state of the taxonomy within Polyommatini is the frustration expressed by Eliot (1973): "I have to admit complete failure in my efforts to find a satisfactory basis for subdividing this very large tribe into a few major natural groups.I have therefore fallen back on naming no less than thirty 'sections', many of them of no more than subsection or even generic worth." Recent molecular genetic studies (Wiemers, 2003;Kandul et al., 2004;Lukhtanov et al., 2009;Wiemers et al., 2009) have largely confirmed a monophyletic genus Polyommatus sensu Hesselbarth et al. (1995) to the exclusion of Glabroculus, Elviria and Kretania.The latter appears to be closely related to Plebejus (Plebejides).The position of Cyaniris and Lysandra, however, differs depending on the molecular markers employed.In this paper, we mainly follow the systematic treatment of Hesselbarth et al. (1995), with the exception of Lysandra.This genus was synonymized with Meleageria by Hesselbarth et al. (1995), but none of the previous molecular studies have suggested a close relationship between these two genera.According to this delineation, the genus Polyommatus consists of approximately 200 species, of which more than half belong to the subgenus Agrodiaetus sensu Eckweiler & Häuser (1997).
Several molecular phylogenetic studies have been conducted on Polyommatus s. str., a group of taxa closely related to the type species P. icarus (Vodolazhsky & Stradomsky, 2008a, b;Vodolazhsky et al., 2009).Other studies focused on the subgenus Agrodiaetus, a large radiation of about 130 species, which is peculiar for its extraordinary variation in chromosome numbers (Lesse, 1960;Wiemers, 2003;Kandul et al., 2004;Lukhtanov et al., 2005;Kandul et al., 2007;Wiemers et al., 2009).Of these latter studies only Wiemers (2003) and Wiemers et al. (2009) include a near comprehensive selection of other Polyommatus taxa.Several species that might be important for understanding the phylogeny of Polyommatus were not included in these analyses, e.g.Central Asian taxa of the stoliczkanus-group (Bryna) and several species of the subgenus Plebicula, including the NW African P. atlanticus (Elwes, 1905), which has the highest chromosome number known in Metazoa (n = 221-223, De Lesse, 1970).
Such large deviations from the modal chromosome number of n = 24 are only known from three subgenera (Agrodiaetus, Lysandra and Plebicula) and therefore some authors assume them to be closely related (Lorkovi , 1990).Resolving the relationships between these subgenera might also increase the level of understanding of chromosomal evolution in Polyommatus.
Also of interest is the evolution of larval food plant associations in Polyommatus.Although the food plants of several Asian species are unknown, all Polyommatus species appear to be restricted to plants of the family Fabaceae.Only a Spanish population of Cyaniris semiargus is also known to utilize another family, Plumbaginaceae (Rodríguez et al., 1991(Rodríguez et al., , 1993)).By contrast, closely related genera (like Plebejus s. l.) also use several other plant families (e.g.Geraniaceae, Cistaceae, Labiatae, Ericaceae, Primulaceae, Polygoniaceae and Chenopodiaceae).While larvae of the ubiquitous P. icarus accept many different genera, most other species appear to be oligophagous on a single or two closely related genera of Fabaceae (Hesselbarth et al., 1995;Tuzov et al., 2000;Tolman & Lewington, 2008).
This paper presents the first comprehensive study using both mitochondrial and nuclear DNA sequences to infer the phylogeny of Polyommatus.
In this paper, we also re-assess the taxonomic status and distribution of P. celina (Austaut, 1879).This taxon (type locality: Sidi-Bel-Abbès in Algeria) used to be regarded as a subspecies, synonym or form of P. icarus until molecular studies (Wiemers, 2003;Wiemers & Fiedler, 2007) discovered strong genetic differentiation (p-distance: 5.9-6.8% in COI and 1.5-2.5% in ITS2) between a population from the High Atlas mountains in Morocco and Eurasian populations.Vodolazhsky & Stradomsky (2008b) confirmed this result using another specimen from the Moroccan Anti-Atlas and raised P. celina to species level based not only on molecular differences but also on a morphological feature, the presence of a broad marginal darkening on the upper surface of the male forewing, which is always absent in P. icarus but typical of P. eros.However, the extent of the marginal darkening in Moroccan specimens appears to be quite variable and can be nearly absent in some cases.The question therefore arises whether P. icarus also occurs on the African continent, e.g. in coastal districts, or is completely replaced by P. celina.assigned to ssp.celina (Wiemers, 1995), likewise needs clarification.

MATERIAL AND METHODS
60 specimens of 23 species of the subgenus Polyommatus and outgroups from all other subgenera of the genus Polyommatus were included in this phylogenetic analysis.Included were taxa from all species groups of Polyommatus s. l. according to Bálint & Johnson (1997) apart from the cyane-and eurypilus-group (subgenera Glabroculus, Elviria & Kretania), which do not seem to belong to Polyommatus (see above).In Agrodiaetus, representatives of each of the clades found by Wiemers (2003) were chosen, plus Agrodiaetus actinides (Staudinger, 1886).Only a single representative of Lysandra was selected, because earlier studies indicate that this subgenus is monophyletic.Previous taxonomic assignments of taxa selected for this study are indicated in Table 2. Chilades trochylus (Freyer, 1845) was used as an outgroup to root the tree, because Chilades is the only Holarctic genus in the subtribe Polyommatina that does not belong to the genus Polyommatus, i.e. does not belong to Polyommatus s.l. or Plebejus s. l.No molecular data are currently available for the other three Polyommatina genera (Pseudolucia Nabokov, 1945;Madeleinea Bálint, 1993;Paralycaeides Nabokov, 1945), which are all restricted to the Neotropical region (mostly confined to the high Andes).
Selected specimens with voucher codes, locality data and GenBank accession numbers are listed in Table 3.The names of the subgenera used follow Hesselbarth et al. (1995) with the exception of Lysandra, which is listed as a distinct subgenus.
An additional 8 specimens of presumed Polyommatus celina from different areas in Morocco, Tunisia and the Canary Islands were used in a COI haplotype analysis (Table 3, Fig. 1).
DNA extraction, PCR and sequencing were carried out according to the protocols described in Vodolazhsky & Stradomsky (2008a) for the specimens with "ILL" voucher codes and in Wiemers (2003) for the remaining specimens.
The 5' (barcode) section of the mitochondrial gene Cytochrome c Oxidase I (COI) and the (mostly complete) nuclear Internal Transcribed Spacer 2 (ITS2) were sequenced because these two sections have proved most successful for resolving the phylogeny of young radiations (Wiemers et al., 2009).
The following cycling protocols were used: For the primer pairs PolF/PolR, PolF/PolRR and PiF/PiR an initial 4 min denaturation at 94°C and 36 cycles of 40 s denaturation at 94°C, 40 s annealing at 58°C and 40 s extension at 72°C; for k698/Nancy an initial 4 min denaturation at 94°C and 35 cycles of 30 s dena-turation at 94°C, 30 s annealing at 55°C, 1 min extension at 72°C and a final extension at 72°C for 7 min; for ITS3/ITS4 an initial 1 min denaturation at 94°C and 40 cycles of 1 min denaturation at 94°C, 1 min annealing at 48°C, 1 minute extension at 72°C and a final extension at 72°C for 1 min.
The alignment of ITS2 sequences was based on the alignment by Wiemers et al. (2009) taking secondary structure information into account.Additional ITS2 sequences were subsequently aligned with ClustalW (Thompson et al., 1994) and some minor manual corrections were carried out.COI sequences were aligned manually.
A Bayesian approach for estimating phylogeny using MrBayes 3.1.1was used (Ronquist & Huelsenbeck, 2003).The data set was divided into 4 partitions, one for each COI codon position and one for ITS2.Model parameters were estimated separately for each partition using MrModeltest 2.2 (Nylander, 2004) and unlinked across partitions.Standard model parameters were applied for the partition containing the 2 nd COI codon position.For the other partitions, a General Time Reversible model was applied, with a proportion of invariable sites for the partitions containing the 1 st and 3 rd COI codon positions and a gamma-shaped distribution of rates across sites.The overall evolutionary rate was allowed to differ for the different partitions.The standard 4by4 substitution model and a flat Dirichlet prior were used for this analysis.Four completely independent analyses, each with four MCMC chains, were run for 10,000,000 generations and sampled every 100 th generation.The heating temperature was set to 0.05 to improve chain mixing (compared to the standard setting of 0.20).The first 200 trees were discarded as burn in.PAUP 4.0 beta 10 (Swofford, 1998) was used to calculate Maximum Parsimony bootstrap values (1000 replicates).Incongruence between the mitochondrial and the nuclear data partitions was determined using the Incongruence Length Difference (ILD) test, implemented as "Partition Homogeneity Test" in PAUP (Farris et al., 1995; but see Barker & Lutzoni, 2002;Planet, 2006).Invariant characters were excluded from the data set (Cunningham, 1997).Separate analyses for the mitochondrial and nuclear partitions were also conducted with the same parameter settings.MEGA 4.1 (Tamura et al., 2007) was used to calculate sequence statistics and pairwise distances (Kimura 2 Parameter model).We used the same programme to date the age of several major nodes by calculating the mean pairwise uncorrected distance of the descendant species of this node.In order to avoid sampling bias, only one sequence per species of the nominate subspecies was used.Standard error estimates were obtained by a bootstrap procedure (1000 replicates).The mean distance was divided by a substitution rate of 1.5% per million years, which appears to be a conserved rate for the COI of arthropods (Quek et al., 2004).
Statistical parsimony networks of COI sequences were calculated with TCS 1.21 (Clement et al., 2000) using a statistical connection limit of 95%.
Mesquite 2.72 was used for ancestral character state reconstruction of karyological and food plant traits using unordered parsimony reconstruction (Maddison & Maddison, 2009).

RESULTS
Because there were no insertions or deletions in the dataset the alignment of COI sequences was straightforward.Of 690 positions, 219 were variable and 153 parsimony-informative.Variability was mainly confined to silent positions.The translation had only 35 (= 15.2%) variable amino acid positions of which 19 (= 8.3%) were parsimony-informative.The aligned ITS2 dataset had 719 positions of which 199 were variable and 81 parsimony-informative (with gaps treated as missing data).
The resulting trees from the separate analyses of the mitochondrial and nuclear partitions (not shown) revealed that most of the phylogenetic signal was concentrated in the ITS2 character set.Resolution of the COI tree was mainly confined to closely related species groups, e.g.among members of Polyommatus s. str.The ILD test indicated a significant difference between the mitochondrial and the nuclear partition (p < 0.001).However, the only major difference revealed by the two analyses is the position of Lysandra, which is sister to Neolysandra coelestina in the COI tree, with a Bayesian support of 0.99, whereas it has the most basal position in the trees inferred from the ITS2 dataset.Another minor difference is in the Central Asian Bryna subclade.COI data suggest a sister relationship between P. stoliczkanus and P. ariana, whereas the ITS2 data set favours P. stoliczkanus and P. venus as sister species.As the difference is small the two data sets were pooled and reanalysed.
Three clades contain species with high chromosome numbers, i.e. much higher than the modal value of 24: the clades containing the subgenera Plebicula, Agrodiaetus and Lysandra (Figs 2-3).None of these clades appear to be sister clades.Furthermore, the sister species to the Plebicula clade is a species with n = 24 (P.thersites).Therefore, marked chromosomal diversification appears to have taken place at least three times in the genus Polyommatus.These results are in agreement with Kandul et al. (2004).
The COI haplotype analysis (Fig. 4) recovers all specimens of Polyommatus celina from Africa and the Canary islands in a single network, which is not connected to the other Polyommatus networks.The lowest interspecific pairwise distance of 3.9% was found between the celina specimens BA09010 and MW02035 and P. eros erotides (ILL068).
With a minimum of only two nucleotide differences and a mean pairwise distance of 0.5% ± 0.002, the populations from Fuerteventura and Tunisia are hardly differentiated from the Moroccan populations.No geographical structure can be detected within the Moroccan populations, which are from different mountain ranges (High Atlas, Anti-Atlas and Middle Atlas).

DISCUSSION
Our molecular phylogeny produces monophyletic groupings, which are largely congruent with morphologybased taxonomic units.An exception appears to be the subgenus Neolysandra, which is not recovered by our analysis.However, this might be due to the strong differences among members of Neolysandra in the COI gene (Wiemers, 2003), while all Neolysandra appear to be closely related in the ITS2 analysis (Wiemers et al., 2009).We also do not find evidence for synonymizing Lysandra under Meleageria, and therefore suggest that Lysandra is kept as a distinct subgenus.
Of special interest are those taxa whose taxonomic relationships have been debated for decades.These are P. escheri, P. amandus and P. thersites, which have been attributed to the (sub-)genera Plebicula (Higgins, 1969;Higgins & Riley, 1978), Agrodiaetus (Higgins, 1975;Higgins & Riley, 1983;Tolman & Lewington, 1997, 2008) or Polyommatus (Tolman & Lewington, 1998).Bálint & Johnson (1997) attached P. thersites to Agrodiaetus in their "actinides-group", whereas they combined P. amandus and P. escheri into a separate "icarius-group".This "icarius-group" (P.icarius Esper, [1793] is a synonym of P. amandus) was equated with "Plebicula s. str.", even though the authors did not include the type species of Plebicula in this group, which is P. dorylas ([Denis & Schiffermüller], 1775) (= Papilio argester Bergsträsser, [1779]).One reason for such shifts could simply be the inability of the Linnean taxonomic system to accurately represent highly nested phylogenetic relationships.However, the analysis of molecular data appears to result in ambiguous relationships between some of these taxa.While ITS2 data suggest a sister relationship between P. escheri and P. amandus if information about secondary structure is taken into account (Wiemers et al., 2009), this is not supported by the COI data (Wiemers, 2003).This conflict might explain why the position of P. escheri is not resolved by the combined analysis.The sister relationship between P. amandus and (a part of) Neolysandra in our tree is surprising and has not previously been suggested.Although the considerable differences in wing pattern do not seem to support this result, the larval host plant associations do: P. amandus and Neolysandra (as far as currently known) share their specificity for Vicia, an otherwise unusual food plant of Polyommatus.On the other hand, a close relationship of P. thersites to Agrodiaetus, which is often suggested partly due to their shared use of Onobrychis as a host plant, is not confirmed by our analysis.The position of P. thersites as a sister to Plebicula, which is mainly based on ITS2 data, appears well supported.It is therefore suggested that P. thersites be included in Plebicula.
Although the subgenus Neolysandra does not appear as a monophyletic unit in this study, in the light of other evidence it is suggested this subgenus is retained pending further study, and that P. amandus and provisionally also P. escheri are attached to this subgenus.It should be mentioned that Bálint & Johnson (1997) split both P. amandus and P. escheri into several allopatric morphospecies.This has not been generally adopted and is not supported by our analysis of molecular data, even though there are large COI distances between African and Asian populations of P. amandus (Wiemers & Fiedler, 2007).
Within the subgenus Polyommatus the monophyly of Polyommatus s. str.sensu Vodolazhsky et al. (2009) is clearly established.This subclade includes the "icarusgroup" sensu Bálint & Johnson (1997) and also P. eros (Ochsenheimer, 1808) from their "eros-group", but not the Central Asian taxa P. amor (Staudinger, 1886) and P. venus (Staudinger, 1886).P. ciloicus de Freina & Witt, 1983, which Bálint & Johnson (1997) placed in the "stoliczkanus-group", also belongs to the Polyommatus subclade and appears as a sister species to P. icadius (Grum-Grshimailo, 1890).The position of P. buzulmavi as a sister to P. icarus is also of interest, because the species status and relationship of P. buzulmavi is debatable.Carbonell (1991) suggested that P. icarus, P. ciloicus and P. stoliczkanus are close relatives, whereas Hesselbarth et al. (1995) mention the possibility of conspecifity with P. icadius.In contrast to P. icarus and P. ciloicus, which have sympatric distributions with P. buzulmavi in Hakkari Province (SE Turkey), P. icadius and P. buzulmavi have allopatric distributions and their genetic distance (COI) is only 1.2-1.7%.Bálint & Johnson (1997) attach P. buzulmavi to the daphnis species-group (Meleageria), but the molecular data clearly preclude such a treatment.A second monophyletic subclade , which consists of only Central Asian species (Bryna), is included by Bálint & Johnson (1997) either in the "stoliczkanus-group" [P.stoliczkanus, P. ariana Moore, 1865 andP. erigone (Grum-Grshimailo, 1890)] or in the "eros-group" (P.amor, P. venus).According to our analysis, the latter species is not especially closely related to P. amor, but a sister to P. stoliczkanus, while the former is closely related to P. erigone.In recent decades, Bryna is hardly ever recognized as a distinct (sub-)genus and usually synonymized with Polyommatus.The molecular results also confirm the close relationship of Bryna with the remaining species of Polyommatus s. str.including P. celina, a taxon which replaces the morphologically and ecologically similar P. icarus in North Africa and the Canary Islands.Its sister relationship remains unresolved and it forms a tritomy with the Polyommatus and Bryna subclades.The basal position of P. celina in the tree suggests an early divergence between African P. celina and its Eurasian sister taxon, which subsequently radiated on the Eurasian continent into the taxa now found in the Polyommatus subclade .
According to our dating estimates, this group is the result of a recent radiation, which occurred during the late Pliocene and Pleistocene, well after the last connection between Northwest Africa and Europe during the Messinian salinity crisis 5.3-5.6 MYA (Garcia-Castellanos et al., 2009).One of the oldest major clades is Agrodiaetus, which evolved 2.9-3.4MYA.This time estimate matches well with that of Kandul et al. (2004) who postulated an age of 2.51-3.85MYA for this subgenus, using a slightly different dating method and selection of genes (COI + COII).The fact that this species-rich radiation with ca 130 species is not represented on the African continent confirms its recent evolution occurred in the Palaearctic region (Wiemers et al., 2009).Most other major clades have representatives in Northwest Africa, but these are mainly species that are also found on the Eurasian continent.They must have reached Northwest Africa either via the Asian land bridge or via dispersal across the strait of Gibraltar.The latter colonization pathway seems most probable for the Moroccan endemic P. atlanticus, whose closest relatives are the Spanish endemics P. nivescens Keferstein, 1851 and P. golgus (Hübner, [1813]).Like almost all Northwest African Polyommatus species, P. atlanticus is confined to higher altitudes, which could explain its genetic differentiation.Nevertheless, the strait of Gibraltar, which is a mere 14 km wide at its narrowest point, appears to be a strong barrier to genetic exchange.This is supported by the strong genetic differentiation between the Eurasian P. icarus and the Northwest African P. celina.Both taxa are among the most ubiquitous butterflies within their range, occurring in almost all open habitats from sea level to high altitudes (up to 3200 m in the High Atlas Mountains).Because these two taxa hardly differ morphologically, their strong genetic differentiation, both in mitochondrial and nuclear DNA sequences, came as a big surprise (Wiemers, 2003).Apparently, P. celina is not even the sister species of P. icarus, and both taxa must have split from their common ancestor about 1.1-2.9MYA, i.e. during the early Pleistocene.This is surprising because both species should have easily been able to cross the strait of Gibraltar.P. icarus is found on almost all Mediterranean islands, even those which have never been in contact with the continent (Dennis et al., 2000), and P. celina is the only Polyommatus species that has reached the Canary Islands.Their colonization involves crossing a minimum of almost 100 km of open sea.The low mean genetic distance between Canarian and Northwest African populations of P. celina indicates a relatively recent colonization, but not in historical times and thus not aided by man.On the Eurasian continent, the last common ancestor of P. celina and P. icarus gave rise to two distinct radiations.One of them, the icarus-subclade , includes the widespread Palaearctic species P. icarus and P. eros, a monophyletic Southeast Palaearctic radiation (P.amorata), the Central Asian P. icadius, as well as some local Kurdistanian and Iranian endemics, i.e.P. buzulmavi, P. ciloicus and P. forsteri.The low genetic divergences among species of this subclade, hybridization between P. icarus and P. icadius indicated by mtDNA introgression (Lukhtanov et al., 2009) and missing synapomorphic molecular characters for P. eros (Vodolazhsky et al., 2009), support its recent evolution.The other radiation, the stoliczkanus-group (subclade ), is confined to the mountains of Central Asia, especially the Pamirs and Himalayas.

CONCLUSIONS
Based on our molecular analysis and additional evidence the following systematic treatment of the genus Polyommatus is suggested: , two further species, which were originally included with Plebicula by Higgins (1969), are located at other positions on the tree.P. amandus forms clade together with Neolysandra fatima (Eckweiler & Schurian, 1980) and N. corona (Verity, 1936), whereas the position of P. escheri is not well resolved.The subgenus Polyommatus s. str.(clade ) forms a monophyletic sister clade to clade .The resulting clade clusters together with Sublysandra, Agrodiaetus and Polyommatus escheri, as a sister to clade .Meleageria and Lysandra form two clearly separate clusters at the base of the tree, and Cyaniris branches off outside Polyommatus (clade ).

Fig. 1 .
Fig. 1.Map showing the distribution of Polyommatus celina.Dots indicate populations that were subjected to DNA analysis.

Fig. 3 .
Fig. 3. Ancestral parsimony reconstruction of food plant use and evolution of chromosome number in Polyommatus.Actual and inferred larval host plant use and the incidence of chromosome numbers strongly deviating from the modal value are mapped onto the combined Bayesian phylogenetic tree (see Fig. 2 for further information).

Fig. 4 .
Fig. 4. Statistical Parsimony Network of Polyommatus celina COI haplotypes.Open nodes indicate unsampled haplotypes.Numbers along branches indicate the positions of nucleotide substitutions in the alignment.

TABLE 1 .
The identity of the P. icarus populations in the Canary Islands, which have also been Previous systematic arrangement of the genus Polyommatus.

TABLE 3 .
List of material with voucher codes and GenBank accession number.