Divergent patterns in the mitochondrial and nuclear diversity of the specialized butterfly Plebejus argus ( Lepidoptera : Lycaenidae )

Plebejus argus is a model species for studying the biology, population ecology and genetics of butterflies. It is patchily distributed throughout most of its European range and considered to be sedentary. Habitats of the butterfly have to encompass two vital larval-resources, i.e. specific food plants and ants, since caterpillars are obligatorily myrmecophilous. The genetic structure of nine P. argus populations (85 individuals) was studied at an intermediate geographical scale (Eastern Poland, diameter of about 400 km) using two kinds of molecular markers i.e. COI (mtDNA) and EF-1 (nuclear gene). Both markers were highly variable with as many as 16 haplotypes and 39 alleles, respectively. Great genetic differentiation in the COI gene was detected (overall FST = 0.411, P < 0.001) but little genetic differentiation in the EF-1 gene (FST = 0.021, P < 0.001). The number of COI haplotypes (ranging from one to seven) and their distribution varied considerably among P. argus populations. The possibility that this heterogeneity was related to Wolbachia was excluded as this endoparasitic bacterium was not detected in samples from any of the populations studied. PCA and SAMOVA analyses divided the sampled populations into two or three groups, which could indicate different colonization routes. Moreover, the differences in genetic differentiation with respect to mtDNA and nuclear markers may suggest male-biased dispersal of P. argus at a larger scale. The hypothesis that females are philopatric is consistent with direct observations of the restricted colonization abilities of the butterfly, while the relatively homogeneous genetic structure revealed by previous allozyme studies in some areas might be explained by the possible higher mobility of males.


INTRODUCTION
The specialized Palaearctic silver-studded blue butterfly Plebejus argus (Lepidoptera: Lycaenidae) may be considered as a model for studying the biology and ecology of insects.In Europe it is still widespread (Kudrna, 2002) and in some areas in the south where appropriate biotopes are continuous, e.g. in SW Spain, it may even be very abundant (Péténian & Nève, 2003).However, in Central and Northern Europe P. argus is patchily distributed and its habitats, as in the case of many other specialized species, are becoming more and more fragmented and isolated.As a result this butterfly is declining in some countries e.g. in the UK (Asher et al., 2001) and Czech Republic (Beneš et al., 2002).
The habitat of P. argus is determined by, amongst other factors, two different resources vital for successful larval development.Caterpillars feed on some species of Ericaceae, Fabaceae and Cistaceae, depending on locality, and are also tended by workers of the ants Lasius niger, L. alienus or L. platythorax (Thomas, 1985;Jordano & Thomas, 1992;Thomas et al., 1999;Péténian & Nève, 2003).Myrmecophily of P. argus is considered as mutu-alistic and obligatory, which is a very rare phenomenon in European butterflies, and reported only for another member of the same genus, P. idas (Fiedler, 2006).Caterpillars produce secretions that are imbibed by workers, which protect them against natural enemies.Pupation occurs within or close to ant nests, and pupae, as well as recently emerged adults, are attended by worker ants.Jordano et al. (1992) report that even oviposition is biased towards plants with appropriate ants at the base, although it is not clear if females respond directly to the presence of ant cues.The habitats of P. argus also feature important structural elements (non consumable resources) e.g.patches of taller vegetation preferred for roosting and mating (Dennis, 2004).
P. argus is given as an example of an extremely sedentary species (Hovestadt & Nieminen, 2009).Imagoes often occur in high densities, which is typical for a closed population structure, most of them move less than 20 m per day and only a very few move more than 50 m.Observations on the movements of imagoes between patches of habitat in metapopulation systems indicate that only a small fraction disperses up to 1.5 km.However, the maximum recorded distance inferred from coloniza-tion events is 4 km (Lewis et al., 1997;Asher et al., 2001).In addition, individuals from isolated localities are generally smaller and have smaller thoracic muscles (i.e. they show selection against mobility) (Thomas et al., 1998).
P. argus has also been the subject of molecular analyses, both in ecological and phylogenetic studies.In Britain allozymes were used to assess metapopulation systems by indirectly evaluating gene flow (Lewis et al., 1997) and to investigate the effects of bottlenecks (Brookes et al., 1997).The population structure in contrasting landscapes in Finland and Spain have been compared using these markers (Péténian & Nève, 2003).Most other analyses involving P. argus have used mitochondrial DNA, particularly the cytochrome oxidase subunit I gene (COI) (Gompert et al., 2008;Wiemers & Fiedler, 2007;Vila & Bjorklund, 2004;Lukhtanov et al., 2009;Dinca et al., 2011).However, none of these studies focused on P. argus and therefore there are only sequences for a very small number of specimens from very different parts of the range.
The aims of the present contribution were to obtain an insight into the genetic differentiation of P. argus at an intermediate geographical scale and determine whether the pattern revealed reflected ecological variation in some way, or resulted from other factors, e.g.isolation by distance.The sampling was carried out in an area with a diameter of about 400 km (Eastern and Southern Poland) and the distances between the sampled localities surpassed those recorded in metapopulation systems.P. argus there inhabits two different types of biotopes and uses different larval food plants.
The COI mitochondrial DNA gene and the nuclear gene elongation factor 1 (EF-1 ), a gene that is commonly used in phylogenetically oriented butterfly surveys to complement mtDNA markers, were sequenced (e.g.Wahlberg et al., 2003Wahlberg et al., , 2005;;Kodandaramaiah & Wahlberg, 2009).Simultaneous studying two genes that differed in their mode of inheritance made it possible to assess whether genetic differentiation might result from bias in dispersal of one of the genders.Moreover the samples were screened for the presence of Wolbachia, which is reported in several lycaenid butterflies (Russell et al., 2009), including close relatives of P. argus (Nice et al., 2009).This endoparasitic bacterium is transmitted in the cytoplasm of the egg and may result in an increase mtDNA differentiation among populations due to heterogeneity of infections among populations (Hurst & Jiggins, 2005).Therefore inspection of P. argus samples for Wolbachia was important for the interpretation of the pattern revealed by this study.

Sampling
The samples were collected in southern and eastern Poland where P. argus is widespread but occurs rather locally (Buszko, 1997).A total of 85 adults were collected in 2008 and 2009 from nine sites (Fig. 1, Table 1).Butterflies with worn wings (from 8 to 11 individuals at every locality) were netted and put into 95% ethanol.Pairwise distances among sites varied from 17 to over 400 km.Localities differed in terms of the biotopes and 538 Fig. 1.Location of the sites from which P. argus were collected in Poland (for full names and details see Table 1) and distribution of COI haplotypes.Haplotypes shared by at least two populations are in depicted in the same colour and unique haplotypes present only in single populations are not coloured (see larval food plants.Most of them were heathlands where Calluna vulgaris was used as the host plant.However, at two sites in the south (CHE and HUT) caterpillars fed on Lotus corniculatus and Coronilla varia, growing in xerothermal meadows.At one site (JAS) it was difficult to determine the larval food plant, as both heather and leguminous plants were present and no feeding caterpillars were found.However, imagoes concentrated in patches of Calluna, which could indicate that this was the preferred plant.

DNA extraction and sequencing
Genomic DNA was extracted from the thorax using Genomic Mini Kit (A&A Biotechnology).The primers used for the COI mtDNA fragment amplification were Ron and Hobbes, for the EF-1 gene ef51.9and efrcM4 were used (Monteiro & Pierce, 2000).The PCRs were performed on GeneAmp PCR System 9700 (Applied Biosystems) and conducted in a 5µL reaction volume containing: Qiagen Multiplex PCR Master Mix (1×), 0.2 µM of each primer, about 10-20 ng of genomic DNA and RNase-free water.Each PCR started with an initial activation step at 95°C for 15 min, followed by 40-45 cycles, with denaturation at 94°C for 30 s, annealing at 57°C for 90 s, extension at 72°C for 60 s and final extension at 60°C for 30 min.Direct sequencing was performed using the BigDye TM Terminator Cycle Sequencing Ready Reaction Kit 3.1 (Applied Biosystems).Both primers were used for sequencing both strands.The sequencing products were run on an ABI PRISM 3130 capillary automated sequencer (Applied Biosystems).
In addition, PCRs that amplify 16S rDNA were used to test for the presence of Wolbachia in P. argus, using the abovementioned protocols and the W-Specf and W-Specr primers (Werren & Windsor, 2000).About 50% of the individuals from all of the populations sampled were chosen at random and screened.As positive controls, samples of Phengaris alcon, in which this bacterium occurs, were used (Sielezniew et al., in press).

Data analysis
DNA sequences were aligned in BioEdit v 7.0.4(Hall, 1999) and revised manually.The haplotype reconstruction of the nuclear EF-1 gene from unphase/genotype data was conducted using the algorithms provided in PHASE as implemented in DNASP 5.0 (Librado & Rozas, 2009).The number of haplotypes, haplotype diversity (h) and nucleotide diversity ( ), and number of segregating sites, were calculated using ARLEQUIN 3.11 (Excoffier et al., 2005).Additionally, a network of the COI haplotypes was constructed using the 95% parsimony criterion implemented in the program TCS (Clement et al., 2000).Population divergence estimates were obtained as FST values in ARLEQUIN.The significance of pairwise FST values for the COI and the EF-1 genes was ascertained by 1000 permutations.The isolation by distance (IBD) pattern was tested by comparing genetic differentiation between populations, as measured by pairwise FST / (1 -FST), to the logarithm of geographical distance using IBD software (Bohonak, 2002).Principal component analysis (PCA) was performed on pairwise FST data and the PC1 scores for each population were obtained using GENALEX 6.0 (Peakall & Smouse, 2006).Analysis of molecular variance (AMOVA; Excoffier et al., 1992) using ARLEQUIN (with 10,000 permutations) was performed to assess structuring within the data.A spatial AMOVA using SAMOVA ver.1.0 (Dupanloup et al., 2002) was used to identify the grouping of sampling sites that maximized the FCT value based on 10,000 simulated annealing steps for K = 2 to K = 8 partitions of sampling sites.To compare the dispersal of the sexes, the migration parameters Nm and M were estimated from FST values (FST = 1 / (4Nm + 1), and M is 4Nem for autosomal and 2Nfm for mtDNA.Nm can be estimated when neutral markers are used, thus Fu and Li's (1993) D* and F* tests were used to check for the assumption of neutrality of mutations and detect the presence of background selection using DNASP 5.0.

RESULTS
Mitochondrial DNA sequences (COI) were obtained for the 802-bp fragment for 84 individuals.Altogether 16 haplotypes were found (GenBank accession nos.JN581045-JN581060).Fourteen sites were variable as defined by twelve transitions and two transversions.A minimum-spanning network of haplotypes is presented in Fig. 2. The estimated nucleotide ( ) and haplotype (h) diversity values for the COI gene were 0.0026 ± 0.16 SE and 0.814 ± 0.026 SE, respectively.For the populations studied, ranged from zero to 0.42% and h from 0.000 to 0.821 (Table 2).The number of COI haplotypes varied considerably among P. argus populations.Two populations were fixed for a single COI haplotype and in the HUT population as many as seven haplotypes were found.None of the haplotypes were present in all populations and the most common were Hap4 and Hap1, which were detected in six and four populations, respectively.On the other hand, eleven haplotypes were populationspecific.Hap4 was the most common haplotype, detected in as many as 29 individuals (Table 2, Fig. 1 and 2).A fragment of Wolbachia 16S rDNA was not amplified in any sample.
Nuclear DNA sequences (EF-1a) were obtained for the 488-bp fragment for 85 individuals.Thirty nine alleles were found (GenBank accession nos.JN581061-JN581099) as defined by 29 heterozygous segregating sites (31 substitutions, 25 transitions and 6 transversions).A minimum-spanning network of haplotypes is presented in Fig. 2. Number of alleles ranged from eight in the PIA, ORC and CHE populations to eleven in the SOW and JAS populations.Only one allele was present in all populations, while 23 alleles were population-specific and eight alleles were singletons.For the EF-1 gene, the average was 0.0056 ± 0.0033 SE (range: 0.0035-0.0066)and h was 0.863 ± 0.024 SE (range: 0.742-0.942)(Table 3).
Great genetic differentiation among P. argus populations in the COI gene was detected (FST = 0.411, P < 0.001).The pairwise FST values for COI were highly vari-able and ranged from 0.000 to 1.000, and the majority of comparisons (29 out of 36) were significant (Table 4).On the other hand, most of the FST values between P. argus populations for the EF-1 gene were small to moderate (0.00-0.063) and only four of 36 pairwise comparisons were significant (Table 4).The overall genetic differentiation in the EF-1 gene was low, albeit significant (FST = 0.021, P < 0.001).
Fu & Li's (1993) D* and F* tests were non-significant for both COI (D* = -1.734,P > 0.1, and F* = -1.618,P > 0.1) and EF-1 (D* = -1.676,P > 0.1, and F* = -1.892,P > 0.05) genes, thus leading to the rejection of background selection occurring in the dataset analyzed.Therefore, the difference between FST values for nuclear EF-1 genes and mtDNA COI may suggest that for P. argus the effective number of migrants (Nem = 11.65)exchanged among the populations studied was about 60 times larger than the effective number of migrant females (Nfm = 0.36).Even between populations, which were geographically close (e.g.PIA-SOW, CHE-SUK, HUT-SUK and JAS-TYS), there was still some evidence for sex-biased dispersal (male dispersal: 3.2, 13, 3.5, and 12 times larger than for females, respectively).
An IBD pattern was not detected in any population of P. argus, with respect to either the COI gene (r 2 = 0.  P > 0.05) or the EF-1 gene (r 2 = 0.0001, P > 0.05).The first and second axes of the PCA explained 41.42% and 35.93% of the total variability in P. argus in the COI gene, respectively.The north eastern populations (SOW, ORC, HOR and PIA) and SUK had positive PC1 scores and remaining southern populations: CHE, HUT, JAS and TYS had negative PC1 scores.As far as PC2 is concerned, values were positive for SOW, PIA, SUK and HUT and negative for CHE, TYS, ORC and HOR, while the PC1 value for JAS was close to zero.For the EF-1 gene, PC1 and PC2 explained 76.01% and 15.82% of the total variability in P. argus, respectively.Two groups were identified with respect to PC1 scores: SOW, ORC and HOR populations had positive PC1 scores (likewise for the COI gene), while the remaining six populations had negative PC1 values (Fig. 3).
Genetic structuring among populations of P. argus was supported by the AMOVA results, where all sampling sites were treated as a single group (FST = 0.450, P < 0.001 for the COI gene and FST = 0.019, P < 0.05 for the EF-1 gene).SAMOVA was then used to identify the subdivision that is most likely to explain the genetic structure observed in P. argus in Poland.For the COI gene the data were best explained by assuming three groups of populations: group I consisted of PIA, HUT CHE, SUK JAS and TYS; group II included ORC and HOR populations, while the SOW population formed group III.The percentage of variation among these three groups was the highest: 47. was 36.78% (FCT = 0.470, P < 0.01, FSC = 0.306, P < 0.001, FST = 0.632, P < 0.001).For the EF-1 gene, the data were best explained by assuming two groups of populations: group I consisted of HOR, while the remaining populations formed group II.The percentage of variation among these two groups was low: 6.29%, while among populations within groups it was 0.35% and within populations it was 93.39% (FCT = 0.063, P < 0.001, FSC = 0.004, ns, FST = 0.066, ns).
It was not possible to amplify a 16S rDNA fragment in any sample of P. argus from any of the populations studied.On the other hand, using the same reaction mix and PCR conditions an approximately 430 bp fragment of DNA was always obtained from the P. alcon samples studied (positive controls).Those fragments, after DNA sequencing reaction and BLAST in the GenBank database turned out to be Wolbachia 16S rDNA (Sielezniew et al., in press).

DISCUSSION
This study indicates that P. argus is a relatively variable species both for mtDNA and nuclear DNA markers.However, it is unknown whether there is a constant and high level of polymorphism throughout its range.Interestingly, results of phylogeographical studies on the nymphalid Melitea cinxia, whose variation in COI surpasses any other Eurasian butterfly studied so far, suggest that in Poland there is a contact zone between Central and Eastern lineages (Wahlberg & Saccheri, 2007).The highest variation was found in the same region in which one of the studied populations of P. argus occurred (e.g.HUT).However, the authors did not sample from other Polish localities of M. cinxia and, in contrast, only P. argus from Poland were included in this study.
The PCA analyses of the samples suggested some geographic structuring in P. argus.The pattern may indicate two different post-Pleistocene colonization routes and the mixed origin of some populations.Unfortunately, the sequences of COI of P. argus deposited in the NCBI database are shorter and moreover most of them do not fully overlap the fragments obtained in this study, and only one sequence of EF-1 is recorded there.Hence, at present it is difficult to draw any inferences about the phylogeographical relationships between specimens collected in Poland and other localities.
Differences in the level of mtDNA polymorphism among localities may also be caused by spatial isolation of populations.In the tropical nymphalid Mycalesis orseis the number of haplotypes decreases with the degree of isolation of habitat patches, but not with population size or patch size (Benedick et al., 2007).In this study the locality of P. argus with the highest polymorphism in the COI gene (HUT) was in part of a large area with many patches inhabited by the butterfly, so it was definitely not isolated.In addition, the HOR population isolated by clear barriers (lands with intensive agriculture and forests) yielded only a single haplotype.On the other hand, the same lack of variation was found in a large population (JAS), which inhabited a mountainous area on the Polish-Slovakian border, where a higher polymorphism was expected.The observed fixation of mtDNA haplotypes in some populations may be explained by genetic drift and/or pattern of colonization.The possibility that the pattern of differentiation in COI was related to Wolbachia infection was excluded, i.e. heterogeneity in its presence/absence, as this bacterium was not found the samples.
There was no evidence of genetic structure that could be explained by ecological variation.High pairwise FST values in mtDNA were found among populations using the same larval food plants.PCA and SAMOVA analyses for both markers also did not support the existence of any connection between larval-plant use and genetic variation.The studies of another lycaenid Jalmenus evagoras in Australia also suggest that ecological variation is not related to genetic variation in the COI gene, and major influences appear due to fragmentation (Eastwood et al., 2006).However, it is also important to mention, that the exact level of larval food plant specialization of populations of P. argus in Poland is unknown as there are no studies on females' or caterpillars' preferences.
The level of genetic differentiation in mtDNA observed in P. argus did not correlate with nuclear DNA differentiation, which was much lower than that recorded for mtDNA.Results obtained might be due to several factors.The effective population size of mtDNA is a fourth of the nuclear autosomal sequences.As a consequence mitochondrial lineages present a much faster lineage sorting rate, and random elimination of haplotypes is much faster than the nuclear allele extinction rate.Therefore, evolutionary relationships could be oversimplified and his- torical events within and among populations (e.g.splitting into isolated units, mixing of different genetic pools, geographical expansions, changes in effective population size) may not be detected correctly using mtDNA data (Zhang & Hewitt, 2003).The present genetic structure based on nucDNA may reflect the past status when the butterfly was less patchily distributed in Poland.This kind of relationship is found at the metapopulation level in M. cinxia (Orsini et al., 2008).

542
A clear discrepancy between mtDNA and nucDNA inferred population structure is often interpreted as evidence for sex-biased dispersal (Goudet et al., 2002).
There are examples of such analyses for many vertebrates (Lawson Handley & Perrin, 2007 and references therein).Usually, neutral mtDNA and autosomal markers (e.g.control region and microsatellites, respectively) are used to check for a sex bias in dispersal.As D* and F* tests were non-significant for both the COI and EF-1 genes, the occurrence of background selection in P. argus was rejected.Moreover, pairwise FST values calculated for microsatellite loci and for the EF-1 gene are similar and positively correlated in P. alcon (Sielezniew et al., in press), indicating that the EF-1 gene can be successfully used to estimate neutral genetic differentiation at nuclear loci, at a similar level to microsatellites.
The possibility that males of P. argus are much more mobile and that gene flow is generated by this sex seems plausible, as the majority of FST values for the EF-1 gene among the populations studied are not significant, while for the COI the opposite pattern is observed.However, butterflies have not been studied in this way and only Zakharov & Hellman (2007) suggest that malebiased dispersal may be one explanation of the discrepancy between mtDNA phylogeography and nuclear microsatellite data for hesperiid Erynnis propertius and Papilio zelicaon in North America.Male-based dispersal is generally poorly documented for insects and the only direct studies are for damselflies (Beirinckx et al., 2006) and grasshoppers (Maes et al., 2006), and molecular methods have revealed this phenomena in ants (Sundström et al., 2003).
Among butterflies some ecological data report differences in propensity to move both within and among habitat patches (Hovestadt & Nieminen, 2009 and references therein).However direct measurements of long distance flights of marked specimens are hardly practicable because of the low recapture probability, and this parameter is probably usually underestimated (Nève, 2009).All direct observations of P. argus dispersal (mark-release-recapture studies) also concern small spatial scales, which do not enable researchers to follow long distance movements.These data do not indicate whether the dispersal capabilities of the sexes differ (Lewis et al., 1997;Seymour et al., 2003).
If male-biased dispersal in P. argus really exists, it would not be surprising taking into consideration theories explaining the evolution of sex-biased dispersal (Gros et al., 2008(Gros et al., , 2009)).P. argus shows distinct protandry and in the first part of the flight period males fly considerable distances in search of females, and they strongly compete for virgin females.Moreover, males are polygynous and this mating system promotes male-biased dispersal in animal populations (Lawson Handley & Perrin, 2007).Another important factor may be related to the costs of emigration.Females of P. argus show specific microhabitat preferences (Jordano & Thomas, 1992) and it is possible that this sex sacrifices long distance dispersal for the sake of fecundity.
The little overall population genetic differentiation of the EF-1 gene is very similar to the values obtained based on studies of allozymes of P. argus in Spain and Finland but lower than those observed in Britain.However, all these studies were performed at smaller spatial scales i.e. distances among populations were 1-16 km in Finland, 6-47 km in Spain and up to about 2.5 km in Britain (Péténian & Nève, 2003;Lewis et al., 1997), and in the latter case some populations were also potentially influenced by the founder effect (Brookes et al., 1997).Therefore, one should be cautious when comparing different systems using different markers.However, the results do not contradict the suggestion that P. argus has low colonization abilities, as recorded in Britain (Asher et al., 2001).Simultaneously a possible male-biased gene flow may explain the relatively homogeneous population structure of the butterfly in some areas with a patchy distribution, as detected in Finland using allozymes.Péténian & Nève (2003) found that only populations isolated by a water barrier were significantly differentiated, and even forest did not impede gene flow.Therefore, P. argus in continental Europe is probably not threatened by genetic deterioration, providing that existing populations are not too isolated by distance and habitat discontinuity.
In the future it would also be interesting to evaluate the dispersal of individuals from Polish localities using markrelease-recapture and compare the results with available data.On the other hand, further indirect studies, using other markers (e.g.microsatellites), are needed to test the hypothesis about male-promoted gene flow in P. argus and if correct, to check whether this is a widespread phenomenon among butterflies and other insects.The data also encourage more extensive genetic studies of P. argus, not only to shed light on the origin of Polish populations, but also to reveal the phylogeographical structure of this interesting butterfly at a wider scale.

Fig. 2 .
Fig. 2. Minimum spanning network of P. argus haplotypes of the COI gene (mtDNA) and EF-1 nuclear gene recorded in Poland.Sizes of circles are proportional to the number of individuals with that haplotype.Black dots along the lines indicate the number of nucleotide substitutions separating connected haplotypes and dotted lines indicate alternative connections.For COI, haplotypes shared by at least two populations are depicted in a particular colour.

Fig. 3 .
Fig. 3. PCA analyses of the variation in the COI (mtDNA) gene and EF-1 nuclear gene among P. argus populations.

Table 3 )
. Sizes of circles are proportional to the number of individuals with that haplotype.

TABLE 1 .
Information on sampling locations and sizes of the samples of P. argus collected in Poland.

TABLE 2 .
Haplotype frequencies of the COI gene (mt DNA) in P. argus populations in Poland with haplotype diversity -h and nucleotide diversity -(%).Total denotes pooled sample of all the P. argus collected in Poland.

TABLE 3 .
Allele frequencies of the EF-1 nuclear gene in P. argus populations with gene diversity -h and nucleotide diversity -(%).

TABLE 4 .
Genetic differentiation (pairwise FST) of the COI gene (above diagonal) and EF-1 (below diagonal) among 10 populations of P. argus in Poland (for full names of the locations see Table1).Significant FST values are given in bold.