Host plant use drives genetic differentiation in syntopic populations of Maculinea alcon

András Tartally; Andreas Kelager; Matthias A. Fürst; David R. Nash

doi:10.7717/peerj.1865

Host plant use drives genetic differentiation in syntopic populations of Maculinea alcon

András Tartally^1,2, Andreas Kelager^2,3, Matthias A. Fürst^2,4, David R. Nash ²

1Department of Evolutionary Zoology and Human Biology, University of Debrecen, Debrecen, Hungary

2Centre for Social Evolution, Department of Biology, University of Copenhagen, Copenhagen, Denmark

3Centre for Macroecology, Evolution and Climate, Natural History Museum of Denmark, University of Copenhagen, Copenhagen, Denmark

4IST Austria (Institute of Science and Technology Austria), Klosterneuburg, Austria

DOI: 10.7717/peerj.1865

Published: 2016-03-31
Accepted: 2016-03-09
Received: 2015-09-30

Academic Editor: Robert Toonen

Subject Areas: Conservation Biology, Ecology, Entomology, Evolutionary Studies, Genetics
Keywords: Phenological separation, Immigration, Disruptive selection, Conservation units, Myrmica, Gentiana, Host specificity, Maculinea rebeli

Copyright: © 2016 Tartally et al.
Licence: This is an open access article distributed under the terms of the Creative Commons Attribution License, which permits unrestricted use, distribution, reproduction and adaptation in any medium and for any purpose provided that it is properly attributed. For attribution, the original author(s), title, publication source (PeerJ) and either DOI or URL of the article must be cited.

Cite this article: Tartally A, Kelager A, Fürst MA, Nash DR. 2016. Host plant use drives genetic differentiation in syntopic populations of Maculinea alcon. PeerJ 4:e1865 https://doi.org/10.7717/peerj.1865

The authors have chosen to make the review history of this article public.

Abstract

The rare socially parasitic butterfly Maculinea alcon occurs in two forms, which are characteristic of hygric or xeric habitats and which exploit different host plants and host ants. The status of these two forms has been the subject of considerable controversy. Populations of the two forms are usually spatially distinct, but at Răscruci in Romania both forms occur on the same site (syntopically). We examined the genetic differentiation between the two forms using eight microsatellite markers, and compared with a nearby hygric site, Şardu. Our results showed that while the two forms are strongly differentiated at Răscruci, it is the xeric form there that is most similar to the hygric form at Şardu, and Bayesian clustering algorithms suggest that these two populations have exchanged genes relatively recently. We found strong evidence for population substructuring, caused by high within host ant nest relatedness, indicating very limited dispersal of most ovipositing females, but not association with particular host ant species. Our results are consistent with the results of larger scale phylogeographic studies that suggest that the two forms represent local ecotypes specialising on different host plants, each with a distinct flowering phenology, providing a temporal rather than spatial barrier to gene flow.

Introduction

Larvae of Maculinea van Eecke (Lepidoptera: Lycaenidae) butterflies start their development on specific host plants. A few weeks later they are adopted into the nests of suitable Myrmica Latreille (Hymenoptera: Formicidae) colonies, where they act as social parasites of the ants (Thomas et al., 1989). This unusual life cycle has shaped their evolution, as different populations are strongly selected to adapt to different initial host plants and Myrmica species depending on their availability (Thomas et al., 1989; Witek et al., 2008).

Larvae of the Maculinea alcon Denis & Schiffermüller group follow a rather specialised development compared to other Maculinea species, as they are not simply predators of ant brood, but are fed by Myrmica workers in preference to their own brood–a behaviour that has been described as a “cuckoo” strategy (Thomas & Elmes, 1998). Because they are constantly interacting with worker ants, this means that they need to adapt precisely to the local host ant species, e.g., by mimicking the odour (Akino et al., 1999; Nash et al., 2008; Thomas & Settele, 2004) and sounds (Barbero et al., 2009) of the ants, in order to be accepted by a suitable Myrmica colony. While the initial host plants of this group are all species of gentian (Gentiana L. and Gentianella Mönch), they can occur in very different open habitats, such as lowland and mountain meadows or wet and dry swards (Munguira & Martín, 1999; Oostermeijer, Vantveer & Dennijs, 1994; Settele, Kühn & Thomas, 2005; Tartally, Koschuh & Varga, 2014). Based on these different habitat types, several forms or (sub)species of the M. alcon group have been described. The most widely accepted separation within this group is that the nominotypic M. alcon occurs on humid meadows and there is another xerophilous form which has usually been referred to as M. rebeli Hirschke (Thomas et al., 2005; Thomas & Settele, 2004; Wynhoff, 1998). Both forms are patchily distributed (Wynhoff, 1998) and have been considered as endangered in many European countries (Munguira & Martín, 1999), with the xerophilous form considered to be a European endemic (Munguira & Martín, 1999). However, several papers (Habeler, 2008; Kudrna & Belicek, 2005; Kudrna & Fric, 2013) have made the case that the xerophilous form is most likely not synonymous with the nominotypic M. rebeli, which is found at higher altitude, and has a unique host plant and host ant usage (Tartally, Koschuh & Varga, 2014). Furthermore, recent molecular phylogenetic studies (Als et al., 2004; Ugelvig et al., 2011b; Bereczki et al., 2015) suggest that the hygrophilous and xerophilous forms of M. alcon, while distinct from other congeners, are not two distinct lineages, and show very little variation in genes normally used for phylogenetic inference. This has been confirmed by several regional population genetic studies (Bereczki et al., 2005; Bereczki, Pecsenye & Varga, 2006; Sielezniew et al., 2012; Bereczki et al., 2015), where there is no consistent separation of the two forms. This has led to the current situation where xerophilous and hygrophilous M. alcon are not distinguished for conservation purposes, and the species is now considered as “of least concern” in Europe (Van Swaay et al., 2010). To avoid confusion, we will refer to the “typical” hygrophilous form of M. alcon as ‘M. alcon H’ and the xerophilous form as ‘M. alcon X’ throughout the rest of this paper, following Tartally, Koschuh & Varga (2014).

The host plant and host ant usage of the two M. alcon forms are different, because different gentian and Myrmica species are available on the hygric sites of M. alcon H and xeric sites of M. alcon X. While M. alcon H starts development typically on the marsh gentian (Gentiana pneumonanthe L.), M. alcon X typically uses the cross gentian (G. cruciata L.), and there is some evidence that enzyme systems related to host plant use may be diverging in the two forms (Bereczki et al., 2015). The development of M. alcon X typically continues in nests of Myrmica schencki Viereck and My. sabuleti Meinert but M. alcon H most often uses My. rubra L., My. ruginodis Nylander or My. scabrinodis Nylander as host ants. Furthermore, some other minor or locally important host plant and host ant species have been recorded for both forms (summarized in Witek, Barbero & Markó, 2014).

Despite these differences in the host plant and ant usage of M. alcon H and M. alcon X, phylogenetic reconstruction using morphological and ecological characters suggests that western Palaearctic M. alcon H are closer to European M. alcon X than Asian M. alcon H (Pech et al., 2004). In combination, all these results suggest local ecological but not genetic differentiation of the two forms between hygric and xeric sites. Until recently this could only be tested by comparing sites that were separated by tens of kilometres or more, but in the last decade a site has been recorded from Răscruci (Transylvanian basin, Romania) where patches supporting M. alcon H and M. alcon X occur in a mosaic separated by tens of meters. The two forms use different host plants and mostly different host ants on this site (Tartally et al., 2008), and their flying periods are largely separated based on the phenology of their host plants (Czekes et al., 2014; Timuş et al., 2013). In addition, most previous genetic studies have been based on collecting samples from either flying adults or caterpillars as they emerge from the host plant, which means that any separation by host ant species could not be examined directly. Our aim was therefore to investigate the genetic differentiation between the two forms of M. alcon at this unique syntopic site, to relate this to differences in both host plant and host ant use, and to make recommendations for the conservation of the forms based on their shared and predicted future histories (c.f. Bowen & Roman, 2005).

Materials and Methods

Field methods

Two sites in Transylvania (Fig. 1) were visited in the summers of 2007 and 2009 to record host plant and host ant usage and to collect genetic samples of M. alcon. Host ant specificity results from 2007 have already been published in Tartally et al. (2008). The first site is at Răscruci (46°54′N; 23°47′E; 485 m a.s.l.), which is predominantly an extensively grazed tall-grass meadow steppe with Gentiana cruciata (the host plant of M. alcon X), but also with numerous small marshy depressions with tall-forb vegetation in which G. pneumonanthe (the host plant of M. alcon H) is common (Czekes et al., 2014). This site gave the unique possibility to compare the host ant specificity and population genetics of M. alcon H and M. alcon X within the same site. To collect samples, two nearby patches were chosen within this mosaic site where G. pneumonanthe and G. cruciata were well separated from each other (there was a ca. 20 m wide zone without gentians). In other parts of this site border effects (because of the potential migration of Myrmica colonies) or the co-occurrence of the two gentians made it difficult to find M. alcon larvae originating clearly from G. pneumonanthe or G. cruciata. The patch with G. pneumonanthe will henceforth be referred to as ‘Răscruci wet’ (M. alcon H patch), while the patch with G. cruciata will be referred to as ‘Răscruci dry’ (M. alcon X patch). The nearest known M. alcon site (a M. alcon H site) to Răscruci is at Şardu (46°52′N; 23°24′E; 480 m a.s.l.), 29 km west of Răscruci, which was chosen as a control site. Şardu is a tall-grass, tall-sedge marshy meadow with locally dense stands of G. pneumonanthe. The two sites are separated by a range of hills without suitable M. alcon habitat (Fig. 1).

Figure 1: Site layout and Bayesian clustering of samples.
(A) Map showing the relative position of the two sample sites. Shaded areas in the detailed map correspond to 50 m contour lines. (B) Initial food plants of the hygric (*G. pneumonanthe*) and xeric (*G. cruciata*) forms of *Maculinea alcon*. Images modified from WikiMedia commons. (C) Posterior probabilities of the number of clusters (K) identified by the Bayesian population assignment program Structure. The solid line shows the mean posterior probability for each value of K, with error bars representing the standard deviation across simulations. The dashed line shows the ΔK values of the posterior probabilities from Structure using the method of Evanno, Regnaut & Goudet (2005) shown relative to the maximum value of ΔK. Peaks in the value of ΔK may represent different levels of population substructure. (D) Comparison of genetic clustering of samples into two to four and ten groups using the Bayesian clustering program Structure without a location prior. Each column represents an individual, and is divided according to its probability of membership of cluster 1 (orange), 2 (blue), 3 (dark purple) 4 (green), or 5–10 (other colours).

Download full-size image

DOI: 10.7717/peerj.1865/fig-1

To obtain data on the host ant specificity and to get samples for genetic analysis, Myrmica nests were searched for within 2 m of randomly selected Gentiana host plants, which is considered to be the approximate foraging zone of worker ants of the genus Myrmica (Elmes et al., 1998). Searches were made no earlier than four weeks before the flying period of M. alcon at both sites, so that any M. alcon caterpillars or pupae found must have survived the winter in the ant nest, and hence have become fully integrated (Thomas et al., 2005). Nests were excavated carefully but completely, after which the ground and vegetation were restored to as close to the original conditions as possible. All M. alcon caterpillars, pupae and exuviae were counted, placed in 98% ethanol, and stored at −20 °C until DNA could be extracted. Five to ten worker ants were also collected from each ant nest and preserved in 70% ethanol for later identification in the laboratory using keys by Seifert (1988) and Radchenko & Elmes (2010). For further details, see Tartally et al. (2008).

Host ant specificity

Host ant specificity (deviation from random occurrence in nests of different Myrmica species) was calculated based on the number of fully grown butterfly larvae, pupae and exuviae in two ways: P1 is the 2-tailed probability from a Fisher exact test of heterogeneity in infection of host ant nests (as implemented at http://www.quantitativeskills.com/sisa/), and P2 is the probability from a randomization test of ant nests between species, using the software MacSamp (Tartally et al., 2008). Published (Tartally et al., 2008) and newly-collected data on host ant specificity were combined for these analyses. In the case of Răscruci, host ant specificity results were calculated separately for Răscruci wet and Răscruci dry and also based on the combined data from both patches (‘Răscruci both’ below).

DNA extraction and microsatellite analysis

DNA was extracted from approximately 1–2 mm³ of tissue from caterpillars or pupae using a 10% Chelex-10 mM TRIS solution with 5 µl Proteinase K. Samples were incubated at 56 °C for minimum 3.5 h or overnight and boiled at 99.9 °C for 15 min. The supernatant was collected and stored at 5 °C or −20 °C for short or long term storage, respectively. For each sample, nine polymorphic nuclear microsatellite loci developed for Maculinea alcon were amplified: Macu20, Macu26, Macu28, Macu29, Macu30, Macu31, Macu40, Macu44, and Macu45 (Table 1; Ugelvig et al., 2011a; Ugelvig et al., 2012) using a REDTaq^® ReadyMix™ PCR r eaction m ix (Sigma-Aldrich). These primer pairs (see concentrations in Table 1) were amplified using standard PCR conditions: initial denaturation for 5 min at 95 °C followed by 30 cycles of 30 s at 95 °C, 30 s at locus-specific annealing temperature (see Table 1) and 30 s extension at 72 °C, finishing with elongation of 15 min at 72 °C run on a Thermo PCR PXE 0.2 Thermal Cycler. Total reaction volume was 10 µl of which 1 µl was template DNA. PCR products were run on a 3130xL Genetic Analyzer with GeneScan 500 LIZ (Life Technologies) as internal size standard and analyzed with GeneMapper ^® Software version 4.0 (Applied Biosystems). Locus Macu40 could not be scored consistently (excessive stutter bands) and was omitted from all further analysis. The overall proportion of alleles that could not be amplified was 4.6% (see Table S1).

Table 1:

Details of microsatellites used in this studyis.

Primer concentration in the PCR mix is given below the table.

Primer	SSR motif	Primer sequences 5′–3′	SR	T_A	Dye and Multiplex	N_a	Genbank accession	Ref.
Macu20	(CT)_n(AT)_n(CT)_n	F: TGGCCCGATTTCCTCTAAAC R: TGCGTGTTTATTTTCATTTTAACAG	92–122	57	Fam 1	9	HM535963	U12
Macu26	(CA)_n	F: CTCCCGGGATAGCATTGAC R: CATTGTCGCGGTCGTAATTC	92–128	57	Ned 2	7	HM535964	U12
Macu28	(CA)_n(CGCA)_n(CA)_n	F: TTTTAATCAAAATCGGTTCATCC R: TCAACCACAAAGCAAGTGAGTC	195–223	57	Fam2	12	KT851400	New
Macu29	(TC)_n	F: AAACGCGCTTATGGCTAAAC R: CGGTATGTCCCGTTACATCG	81–143	57	Vic 3	15	KT851401	New
Macu30	(TG)_n	F: GACGCGCTGTTATGTATTGC R: CGTCTAGCGTGACCGTAACA	93–109	57	Pet 4	5	HM586096	U11
Macu31	(GTA)_n(GTC)_n(GTA)_n	F: GTTCTGTCCCCCGAACTAGG R: AAACCTGGGATTGGTTAAAAAC	110–173	62	Ned 5	5	HM586097	U11
Macu40	(CA)_n(GA)_n(CA)_n(GA)_n (CA)_n(GA)_n(CA)_n	F: CCGTTTGGGAGATACGATGT R: CGCGTGTGCGTATATGTGAT	110–220	57	Pet 1	–	KT851402	New
Macu44	(AC)_n	F: ATAAGTCAGCACGTCAAAGCTG R: TGCAAATACTCCGAATAAATAACTG	170–220	57	Ned 3	10	HM535965	U12
Macu45	(AC)_n(GC)_n(AC)_n	F: TGTGTGACTGCGGTTCTTATC R: TGTAATCGCAGGAGAGATGTG	145–217	57	Vic 4	20	HM535966	U12

DOI: 10.7717/peerj.1865/table-1

Notes:

SR: product size range (base pairs)
T_A: Annealing Temperature (°C), Primer dye and Multiplex group
N_a: Number of alleles
Ref.: Reference source (U11, Ugelvig et al., 2011a, U12, Ugelvig et al., 2012; New, This study)

Primer concentration in PCR mix: 0.1 ng/μl: Macu20, Macu26, Macu 29, Macu30; 0.2 ng/μl: Macu 28, Macu 31, Macu 40, Macu 44, Macu 45.

Tests for Hardy-Weinberg and linkage disequilibrium

The eight microsatellite loci analysed were tested for linkage disequilibrium (genotypic disequilibrium) between all pairs of loci in each sample and for deviations from Hardy–Weinberg proportions using exact tests in FSTAT version 2.9.3.2 (Goudet, 1995) based on 480 and 1,260 permutations, respectively. The software package Micro-Checker version 2.2.3 (Van Oosterhout et al., 2004) using 1,000 iterations and a Bonferroni corrected 95% confidence interval, was employed to test for possible null-alleles.

Population structure, genetic differentiation and kinship

We studied the genetic clustering of individual genotypes using the Bayesian algorithm implemented in Structure version 2.3.4 (Falush, Stephens & Pritchard, 2003; Falush, Stephens & Pritchard, 2007; Hubisz et al., 2009; Pritchard, Stephens & Donnelly, 2000). The most likely number of genetically distinct clusters (K) was estimated for each K in the range 2–12, allowing for sub-structuring of samples. A burn-in length of 50,000 MCMCs was used to secure approximate statistical stationarity, followed by a simulation run of 500,000 MCMCs using an admixture model with correlated allele frequencies as recommended by Pritchard, Stephens & Donnelly (2000). No location prior was used, and LnP(D) values were averaged over 20 iterations. The most likely value of K (number of clusters) was estimated using the ΔK method of Evanno, Regnaut & Goudet (2005). To check whether the assumptions inherent in Structure were biasing our genetic clustering, we also used the Bayesian genetic clustering programs BAPS version 5.2 (Corander et al., 2008) and InStruct version 1.0 (Gao, Williamson & Bustamante, 2007), which gave essentially identical results (see Additional Analysis S1, Figs. S1 and S2).

For more detailed population differentiation, samples were explored individually as well as in four different partitions: (a) pre-defined populations (Pop: Răscruci dry, Răscruci wet and Şardu) which also relates to host plant use (Răscruci wet and Şardu: G. pneumonanthe. Răscruci dry: G. cruciata), (b) host ant use (Ant: Myrmica scabrinodis, My. sabuleti, My. schencki and My. vandeli), (c) host ant nests (Nest: specific nest ID within Pop), and (d) year of sampling (Year: 2007 and 2009), the latter to test for potential temporal differences.

We studied the overall population differentiation between pre-defined populations (Pop) calculating Weir & Cockerham’s (1984) estimate of F_ST (θ) using FSTAT version 2.9.3.2 based on 1,000 permutations. As the magnitude of the value of θ is related to the allelic diversity at the marker loci applied, we further calculated the standardized $G_{S T}^{″}$ (Meirmans & Hedrick, 2011) , and the estimator D_EST (Jost, 2008) as alternative quantifications of genetic differentiation, making comparisons with studies based on other marker loci possible (Meirmans & Hedrick, 2011). $G_{ST}^{″}$ and Jost’s D_EST for pairs of Pop samples were calculated using GenoDive version 2.0 b27 (Meirmans & Van Tienderen, 2004). Hierarchical AMOVA (Analysis of Molecular Variance: Excoffier, Smouse & Quattro, 1992) was calculated for Pop, Ant, Nest, and separately for Pop and Year using the R-package HierFStat (Goudet, 2005) with 9,999 permutations to estimate the variance components and their statistical significance. Individual-based Principal Coordinate Analysis (PCoA) with standardized covariances was employed to obtain a multivariate ordination of individual samples based on pairwise genetic distances, as implemented in the software GenAlEx version 6.502 (Peakall & Smouse, 2012). The PCoA were explored for Nest within Ant within Pop across Year using nested MANOVA based on the sum of the variances of the different coordinates, as implemented in JMP 12.02 (SAS Institute).

To examine whether the low dispersal ability of Maculinea alcon females could lead to high relatedness between samples of individuals collected in the same nest, pairwise measures of kinship (Loiselle et al., 1995) and relatedness (Queller & Goodnight, 1989) between samples were estimated using GenoDive and GenAlEx respectively. Both of these values estimate the probability that samples share alleles by descent, based on the distribution of alleles in the whole set of samples, with possible values ranging from −1 to +1. Negative values show that the two individuals compared are less similar in the alleles they share than two randomly picked individuals. Values of kinship and relatedness were then compared between the different partitions of the data using stratified Mantel tests as implemented in GenoDive. To further test the hypothesis that individuals found in the same nest were likely to derive from eggs laid by the same female, the program Colony (v 2.0.6.1; Jones & Wang, 2010) was used to give a maximum likelihood estimate of the probability that any two sampled individuals were likely to be either full or half siblings. Since males of M. alcon are much more mobile than females, who tend to oviposit in a limited area (Kőrösi et al., 2008; DR Nash, 2009, unpublished data), this was based on a mating system assuming female polygyny and male monogyny, with other parameters kept at their default values.

Results

Host ant specificity

A total of 135 Myrmica ant nests were found within 2 m of host gentian plants on the two sites, and 90 Maculinea larvae, pupae and exuviae were found in 26 infested nests (Table 2) including 87 nests and 56 Maculinea already published in Tartally et al. (2008). Altogether four Myrmica species were found. Only My. scabrinodis was present at all sites, and was the most abundant ant species (59% of all ant colonies found). This species was used as a host on all three sites. Only a single M. alcon X was found in a nest of My. scabrinodis at Răscruci dry despite the dominance of this ant there and its frequent usage by M. alcon H on Răscruci wet (Fisher’s exact test, P = 0.032). The much greater exploitation rates of My. sabuleti and My. schencki led to significant overall host ant specificity at Răscruci dry (Table 2).

Table 2:

Details of sampled Myrmica nests.

The number of nests found within 2 m of gentians at each site, their infection with M. alcon H or M. alcon X, the number of individual M. alcon used for genetic analysis (“Genetic samples”: listed in Table S1), and statistical tests of host ant specificity within each site: P1, probability from Fisher exact test and P2, probability from a randomization test of ant nests between species. Significant P-values (P < 0.05) are marked in bold.

Site	Maculinea	Myrmica	No. nests	No. with M. alcon	P1	No. of M. alcon	Range	P2	Genetic samples
Răscruci dry	alcon X	sabuleti	10	5	0.004	17	1–8	0.002	13
		schencki	6	2		18	1–15		5
		scabrinodis	23	1		1			1
Răscruci wet	alcon H	scabrinodis	31	9	–	30	1–7	–	28
Răscruci both	both	as above			0.078			0.021
Şardu	alcon H	vandeli	27	2	0.147	9	2–7	0.495	2
		scabrinodis	38	7		15	1–4		11

DOI: 10.7717/peerj.1865/table-2

Genetic diversity and inbreeding

Measures of genetic diversity and F-statistics generated by FSTAT for each locus are listed in Table S2 in the supporting information. Analysis with Micro-Checker revealed that Macu29 had a highly significant (P < 0.001) excess of homozygotes and cases of non-amplification consistent with the presence of a relatively high proportion (>20%) of null-alleles, and was therefore excluded from further analysis. All other loci showed no significant deviations from Hardy-Weinberg proportions. Tests for linkage disequilibrium revealed only a few sporadic significant results showing no overall pattern (Table S3), so all loci were retained in further analysis, which was thus based on seven polymorphic loci. All three of the pre-defined populations showed no evidence of inbreeding (Table 3), and in fact showed negative values for the inbreeding coefficient F_IS (meaning an excess of heterozygotes), although not significantly so (Table 3).

Table 3:

Pairwise differentiation between, and inbreeding and genetic diversity within predefined populations.

Values above the diagonal in the matrix (with blue background) are θ(F_ST), values along the diagonal (with green background) are F_IS, values below the diagonal (with yellow background) are

G_{ST}^{″} ∕ D_{EST}

. Values in bold differ significantly from zero (P < 0.001). Below the matrix are mean values (±SE) of four different measures of within-population genetic diversity. The effective number of alleles per locus (A_E), the observed heterozygosity (H_O), the expected heterozygosity (H_E) and the unbiased expected heterozygosity (uH_E). P-values for comparisons between pre-defined populations based on mixed model comparison across loci are shown on the right.

	Răscruci dry	Răscruci wet	Şardu
Răscruci dry	−0.050	0.093	0.059
Răscruci wet	0.302/0.255	−0.106	0.103
Şardu	0.207/0.151	0.330/0.221	−0.052
				P
A_E	4.000 ± 0.606	3.365 ± 0.418	4.594 ± 0.909	0.131
H_O	0.707± 0.067	0.717± 0.056	0.733± 0.090	0.944
H_E	0.707± 0.046	0.674± 0.034	0.709± 0.063	0.781
uH_E	0.729± 0.048	0.687± 0.035	0.738± 0.066	0.613

DOI: 10.7717/peerj.1865/table-3

Population structure

Structure analysis revealed rather invariable log-likelihood values for partitioning of the data into genetic clusters, but the highest change in log-probability value was for K = 2, with lower maxima at K = 4 and K = 10 (Fig. 1). There was a clear overall distinction between samples from Răscruci wet in one genetic group and Răscruci dry and Şardu in another group. Levels of admixture between genetic clusters were generally low, but four individuals from Răscruci wet (from two different My. scabrinodis nests) showed high affinity to the Răscruci dry-Şardu group, irrespective of the value of K. One individual found in a My. sabuleti nest at Răscruci dry (sample code: DA14) appeared genetically more similar to those from Răscruci wet. For values of K higher than 2 there was no additional partitioning between the pre-defined populations, but some substructure in Răscruci wet became apparent for K = 3, with two partitions that were relatively dissimilar, while for K ≥ 4 no additional grouping of individuals was apparent (Fig. 1).

Genetic differentiation

We found significant overall genetic differentiation between pre-defined populations (θ = 0.090, D_EST = 0.215; Table S2). Pairwise genetic differentiation measures θ, $G_{ST}^{″}$ and D_EST were significant for all population comparisons after Bonferroni adjustment (P < 0.003; Table 3). There was no evidence of inbreeding, either overall (F_IS = − 0.074, Table S2), or within pre-defined populations (Table 3).

Hierarchical AMOVA (Table 4) revealed that most genetic variance (93.4%) was within individuals, but that significant variation was also explained by Pop and Nest. The proportion of variation and inbreeding coefficients for individuals within Nest and Ant were both negative, indicating that there was greater heterozygosity between individuals in the same nest and between samples across different Myrmica species than between randomly chosen samples from the data set. Samples from different years explained only 0.12% of the genetic variance in a separate AMOVA (P = 0.354). The Principal Coordinate Analysis retained a total of six principal coordinates with eigenvalues greater than 1, which together explained 52% of the variance in genetic distance (13.5%, 10.3%, 8.4%, 8.0%, 7.2% and 4.5% for coordinates 1–6 respectively; see Fig. S3 for more details). These showed a similar result to the AMOVA where Year samples (2007 and 2009) overlapped completely in genetic ordination space (F_1,36 = 8.91 × 10⁻¹⁶, P = 0.999), while samples from Răscruci wet were separated from those from Răscruci dry and Şardu (Fig. 2; F_2,36 = 6.08, P = 0.005). We found a pronounced structuring of samples when examining nests within pre-defined populations (F_18,36 = 3.59, P < 0.001), with samples from the same nest clustering together, but there was no consistent clustering of samples from the nests of the same host ant species (F_3,36 = 0.887, P = 0.457). The single sample from Răscruci dry that was collected from a My. scabrinodis nest (sample code DB15) had a first principal component that was more characteristic of samples from Răscruci wet (which all used this host ant species; Fig. 2), but was not assigned to this population in the Bayesian analysis (Fig. 1).

Figure 2: Ordination of samples based on principal coordinate analysis.
Each symbol represents an individual, coloured according to its pre-defined population (blue, Răscruci wet, orange, Răscruci dry, purple, Şardu). Coloured lines are convex hulls enclosing all samples from each pre-defined population, while coloured regions are convex hulls enclosing samples collected from the same nest. The single individual (sample DB15) collected from a *My. scabrinodis* nest at Răscruci dry is shown with a larger symbol.

Download full-size image

DOI: 10.7717/peerj.1865/fig-2

Figure 3: Relatedness and parentage analysis of samples.
The pairwise matrix shows the estimated Queller & Goodnight (1989) relatedness of each pair of individuals (excluding those with negative relatedness) above the diagonal, and the probability that each pair are full siblings based on maximum likelihood estimates from Colony (Jones & Wang, 2010) below the diagonal. Comparisons between samples from the same pre-defined population are shaded according to the same colour scheme as Fig. 2 (blue, Răscruci wet, orange, Răscruci dry, purple, Şardu). Individuals sharing the same ant nest are outlined with lines in these same colours, and those sharing the same *Myrmica* species as host are outlined with black lines. The area and shade of each data point is proportional to the relatedness or probability of being full siblings for that pair of individuals.

Download full-size image

DOI: 10.7717/peerj.1865/fig-3

Table 4:

Hierarchical analysis of molecular variance.

Calculated using HierFStat (Goudet, 2005). The F-coefficient gives the estimated inbreeding coefficient (excess of homozygotes) at each hierarchical level. P-values are based on 1,000 re-samplings of the data.

Source	d.f.	Variance component	%variance	F-coefficient	P
Between Pop	2	0.577	10.6	0.106	0.027
Ant within Pop	3	−0.257	−4.7	−0.053	0.956
Nest within Ant	18	0.609	11.2	0.119	<0.001
Individuals within Nest	36	−0.569	−10.5	−0.126	>0.999
Within Individuals	60	5.067	93.4
Total	119	5.427

DOI: 10.7717/peerj.1865/table-4

Kinship and parentage analysis

Overall pairwise relatedness of individuals (Fig. 3) sampled from the same nest (0.311) was significantly higher than that of those sampled from different nests within the same site (−0.032; Stratified Mantel test: r² = 0.069, P ≤ 0.0001). Looking at individual sites, the same pattern was found at both M. alcon H sites (Răscruci wet: within-nest relatedness = 0.32, between nests = 0.10, r² = 0.130, P = 0.0002; Şardu: within-nest = 0.249, between nests = 0.065, r² = 0.126, P = 0.002), but relatedness was not significantly different within and between nests at Răscruci dry (within-nest = 0.003, between nests = 0.081, r² = 0.001, P = 0.448). Similar results were found when analysing pairwise kinship (see Fig. S4 and Additional Analysis S2). Maximum-likelihood analysis using Colony identified 38 pairs of individuals as potential full siblings (with probabilities ranging from 0.002 to 0.935; Fig. 3), and 161 as potential half siblings (with probabilities ranging from 0.002 to 0.742; Fig. S4). For the 21 pairs with high (>0.5) probability of being full siblings, 13 (62%) were from the same nest, six were from nests of the same ant species at the same site, and only two were from different sites, both including individual SCA86-2 from Răscruci wet, which appears closely related to two individuals (DA14 and SAB67-1) from different Myrmica nests from Răscruci dry (Fig. 3). This last result probably reflects the non-amplification of characteristic loci for these individuals (see Table S1), so that they share common alleles without being related. Within sites, a high proportion of individuals from both Răscruci wet (60.4%) and Şardu (37.5%) from multiply-infested nest had individuals estimated to be full siblings in the same nest (overall 52.9%, Generalized linear model with binomial errors and Firth corrected maximum likelihood, comparing sites: Likelihood-ratio χ² = 0.73 d.f. = 1, P = 0.391). However, none of the M. alcon X individuals sharing nests at Răscruci dry were estimated to be full siblings (M. alcon X vs. M. alcon H: likelihood-ratio χ² = 9.72, d.f. = 1, P = 0.002).

Although there was evidence for strong within-nest relatedness of individuals, the patterns of genetic diversity and differentiation were not strongly affected by this, and were unaffected when analyses were repeated with only a single individual from each nest (see Table S4).

Discussion

This study gives the first comparison of the host ant specificity and genetic composition of M. alcon H and M. alcon X within the same site.

The host ant specificity found in this study confirms the earlier results of Tartally et al. (2008) that these populations use the typical host ants found in other Central European studies (Höttinger, Schlick-Steiner & Steiner, 2003; Sielezniew & Stankiewicz, 2004; Steiner et al., 2003; Tartally et al., 2008; Witek et al., 2008). M. alcon H was found exclusively with My. scabrinodis at Răscruci wet and also with My. vandeli at Şardu, but M. alcon X was found mainly with My. sabuleti and My. schencki at Răscruci dry. Interestingly only one M. alcon X was found with My. scabrinodis, despite this Myrmica species being the most numerous at Răscruci dry (Table 2) and being the main host of M. alcon X in two other sites in the Carpathian-Basin (Tartally et al., 2008). My. scabrinodis usage could therefore be a potential link between the M. alcon H and M. alcon X populations at Răscruci (and probably in some other regions), but M. alcon X shows a clear separation from the M. alcon H in the proportional usage of this host ant. The background of this separation in the host ant specificity of M. alcon H and M. alcon X at Răscruci is not clear, but could reflect the dynamic arms race between the different genetic lineages of M. alcon and local host ants (Nash et al., 2008).

Our genetic results (Figs. 1 and 2) show strong genetic differentiation between M. alcon H and M. alcon X at Răscruci, indicating limited gene flow between these two groups, although it is interesting to note that a few individuals had genotypes more characteristic of the other population. This differentiation is likely due to separation in time rather than space because of the different phenology of the host plants, which results in largely non-overlapping flying seasons of M. alcon H and M. alcon X (Timuş et al., 2013). This may be reinforced by lowered fitness of any hybrid individuals that would emerge during the approximately 2-week gap when neither host plant is suitable for oviposition.

The lowest level of between-population differentiation, on the other hand, was between M. alcon X from Răscruci and M. alcon H from Şardu, and Bayesian population assignment suggests that these are so similar that they have almost certainly been part of a single population. This supports previous findings of no overall phylogenetic differentiation between the two forms of M. alcon (Als et al., 2004; Fric et al., 2007; Ugelvig et al., 2011b), and that the two forms tend to be more genetically similar regionally than either is to more distant populations that use the same host plant (Bereczki et al., 2005; Bereczki, Pecsenye & Varga, 2006; Pecsenye et al., 2007). Hence the two forms cannot be regarded as host races (Drès & Mallet, 2002), since they do not fulfill the criterion of spatial and temporal replicability. Genetic analysis of several Polish and Lithuanian M. alcon populations using microsatellite markers (Sielezniew et al., 2012) gave similar results to ours (Figs. 1 and 2) in that there was no clear pattern reflecting genetic division into two ecotypes. Sielezniew et al. also found that the M. alcon X ecotype was less polymorphic, and its populations more differentiated than those of the M. alcon H ecotype. Their data also suggest that M. alcon H populations form a single clade but M. alcon X can be split into more clades, suggesting that M. alcon H is an ancestral form and that M. alcon X represents a group of independently evolved M. alcon H populations that have switched to use dryer habitats with the locally available Gentiana and Myrmica species. They propose that the background of this pattern may be independent specialisations on different host ant species, since in their study clades of M. alcon X largely reflected host ant use. However, we find no evidence of genetic differentiation associated with host-ant usage at Răscruci or Şardu (Figs. 1– 3), and no difference in genetic diversity in populations of the two ecotypes (Table 3). Due to the relatively large distances and potential barriers between Răscruci and Şardu it is unlikely that there has been recent gene-flow between the two sites, which suggests that Răscruci was likely colonized at least twice from two different gene pools, and that the ancestors of the Răscruci wet population are no longer locally extant (or have evaded detection).

The higher within-nest than between-nest relatedness between individuals of M. alcon H is consistent with observations of limited dispersal of ovipositing females (Kőrösi et al., 2008) which is likely to lead to substantial within-population substructure between nests, as found here. The additional grouping of individuals from K = 2 to K = 4 in our population assignment analysis groups families of relatives within populations. Our parentage analysis confirms the relatively high probability that M. alcon H caterpillars infesting the same Myrmica nest are relatives, and in around 50% of cases may be full siblings. The lack of any association between relatedness and infestation of nests for M. alcon X from the Răscruci dry site is consistent with the difference in oviposition strategy and mobility of butterflies from this population compared with those from the Răscruci wet site (Czekes et al., 2014; Timuş et al., 2013).

Given the small size of these M. alcon populations (Timuş et al., 2013) and their low dispersal (Kőrösi et al., 2008), it is interesting to note that there is no evidence of inbreeding among the individuals that we examined, and that estimated inbreeding coefficients (F_IS and G_IS) were negative. This means that individuals were more heterozygous than expected if mating were random (albeit not significantly so), which is probably a result of difference in mating strategy of males and females. Males are highly mobile, and tend to patrol a large area while seeking females, whereas females are rather sedentary, and are often mated immediately after emerging from the Myrmica nest in which they developed, and then go on to lay eggs on host plants relatively close by (Kőrösi et al., 2008). This means that females are unlikely to be related to their mates, and may in fact be less related to them than to a randomly chosen male, because males pupate and emerge several days before females (Meyer-Hozak, 2000). This, together with the observation that caterpillars that develop in the same Myrmica nest may be offspring of the same female, can easily lead to the negative inbreeding indices observed both overall, and for individuals within Myrmica nests and host ant species within pre-defined populations within our hierarchical AMOVA.

Regardless of its origin, it is clear that the M. alcon X population at Răscruci is ecologically highly differentiated from the local M. alcon H populations in terms of its host plant and host ant use, as well as in its behaviour (Czekes et al., 2014; Timuş et al., 2013). The pattern of differentiation we see is not that typically associated with speciation via host race formation (Drès & Mallet, 2002), since the two forms do not fulfill the criterion of spatial and temporal replicability (i.e., “are more genetically differentiated from populations on another host in sympatry (and at the same time) than at least some geographically distant populations on the same host”; Drès & Mallet, 2002, p. 473–4). M. alcon X cannot, therefore, represent an evolutionarily significant unit in conservation terms as usually defined, but we would argue that it should still be regarded as a functional conservation unit (Casacci, Barbero & Balletto, 2013). The site at Răscruci represents the only known area where both forms of M. alcon occur syntopically, and so is of particular value to research on speciation, and has great potential for examining adaptation at non-neutral genetic markers. This is enhanced by the occurrence of two other Maculinea species on the same site; M. teleius Bergsträsser (Tartally & Varga, 2008) and M. nausithous kijevensis Sheljuzhko (Rákosy et al., 2010; Tartally & Varga, 2008), as well as the Myrmica parasites Microdon myrmicae Schönrogge et al. (Diptera: Syrphidae; Bonelli et al., 2011) and Rickia wasmannii Cavara (Ascomycota: Laboulbeniales; Tartally, Szűcs & Ebsen, 2007). The Maculinea spp. parasitoid Ichneumon eumerus Wesmael (Hymenoptera: Ichneumonidae) is also present (Tartally, 2008). Most of these species are also found at Şardu (except M. alcon X and M. nausithous: Tartally, 2008). It should be emphasized that all of these species can be found in the nests of, and ultimately depend on, My. scabrinodis (as well as other Myrmica species: see Witek, Barbero & Markó (2014) for a review), providing a unique opportunity to examine a complex set of parasitic interactions revolving around a single keystone ant species.

Conclusion

Our analysis of Maculinea alcon from a unique site where both the xerophilous and hygrophilous forms of this butterfly are found within tens of meters of each other has demonstrated strong genetic differentiation between the two forms. However, the xerophilous form was not significantly differentiated from the next closest known population of the hygrophilous form. This supports other recent work suggesting that the hygrophilous and xerophilous forms are not separate species or even subspecies, and that the name M. rebeli has frequently been applied to the xerophilous form incorrectly. There is some overlap in host ant species use between the two forms, so the most likely proximate reason for the local genetic differentiation found is differences in host plant phenology. We suggest that since the two different forms of M. alcon do not have a separate evolutionary history, they cannot be regarded as “evolutionarily significant units” for conservation as the term was originally used (Bowen & Roman, 2005; Casacci, Barbero & Balletto, 2013). However, since they represent current ecological diversification of the group, and may have different evolutionary potentials (e.g., through selection on different enzyme systems; Bereczki et al., 2015), they should continue to be treated as separate management units for long-term conservation (Bowen & Roman, 2005; Casacci, Barbero & Balletto, 2013). Hence, we support the continued separation of the two forms in future studies to further explore their evolutionary trajectories and conservation potential.

Supplemental Information

Combined supplementary information

Table S1. Multilocus genotypes of samples at the seven microsatellite loci used for analysis (raw data). Table S2. Summary diversity indices and F-statistics for loci analyzed in Maculinea alcon. Table S3. Linkage (genotypic) disequilibrium for Transylvanian Maculinea alcon populations. Additional Analysis S1. Bayesian clustering of samples using the population assignment programs Structure, BAPS and InStruct. Figure S1. Comparison of the posterior probabilities of the number of clusters (K) identified by the Bayesian population assignment programs Structure, BAPS and InStruct. Figure S2. Bayesian clustering of samples for selected values of K. Figure S3. Matrix plot of principal components. Figure S4. Kinship and half-sibship analysis of samples. Additional Analysis S2. Distribution of Relatedness and Kinship coefficients. Figure S5. Distribution of Relatedness and Kinship coefficients Table S4. Genetic diversity, inbreeding and differentiation when only one individual is selected from each nest. References

DOI: 10.7717/peerj.1865/supp-1

Download

[1] Akino T, Knapp JJ, Thomas JA, Elmes GW. 1999. Chemical mimicry and host specificity in the butterfly Maculinea rebeli, a social parasite of Myrmica ant colonies. Proceedings of the Royal Society of London, Series B: Biological Sciences 266:1419-1426

[2] Als TD, Vila R, Kandul NP, Nash DR, Yen SH, Hsu YF, Mignault AA, Boomsma JJ, Pierce NE. 2004. The evolution of alternative parasitic life histories in large blue butterflies. Nature 432:386-390

[3] Barbero F, Thomas JA, Bonelli S, Balletto E, Schönrogge K. 2009. Queen ants make distinctive sounds that are mimicked by a butterfly social parasite. Science 323:782-785

[4] Bereczki J, Pecsenye K, Peregovits L, Varga Z. 2005. Pattern of genetic differentiation in the Maculinea alcon species group (Lepidoptera, Lycaenidae) in Central Europe. Journal of Zoological Systematics and Evolutionary Research 43:157-165

[5] Bereczki J, Pecsenye K, Varga Z. 2006. Geographical versus food plant differentiation in populations of Maculinea alcon (Lepidoptera: Lycaenidae) in Northern Hungary. European Journal of Entomology 103:725-732

[6] Bereczki J, Rácz R, Varga Z, Tóth JP. 2015. Controversial patterns of Wolbachia infestation in the social parasitic Maculinea butterflies (Lepidoptera: Lycaenidae) Organisms Diversity and Evolution 15:591-607

[7] Bonelli S, Witek M, Canterino S, Sielezniew M, Stankiewicz-Fiedurek A, Tartally A, Balletto E, Schönrogge K. 2011. Distribution, host specificity, and the potential for cryptic speciation in the hoverfly Microdon myrmicae (Diptera: Syrphidae), a social parasite of Myrmica ants. Ecological Entomology 36:135-143

[8] Bowen BW, Roman JOE. 2005. Gaia’s handmaidens: the orlog model for conservation biology. Conservation Biology 19:1037-1043

[9] Casacci LP, Barbero F, Balletto E. 2013. The evolutionarily significant unit concept and its applicability in biological conservation. Italian Journal of Zoology 81:182-193

[10] Corander J, Marttinen P, Sirén J, Tang J. 2008. Enhanced Bayesian modelling in BAPS software for learning genetic structures of populations. BMC Bioinformatics 9:539

[11] Czekes Z, Markó B, Nash DR, Ferencz M, Lázár B, Rákosy L. 2014. Differences in oviposition strategies between two ecotypes of the endangered myrmecophilous butterfly Maculinea alcon (Lepidoptera: Lycaenidae) under unique syntopic conditions. Insect Conservation and Diversity 7:122-131

[12] Drès M, Mallet J. 2002. Host races in plant-feeding insects and their importance in sympatric speciation. Philosophical Transactions of the Royal Society of London, Series B: Biological Sciences 357:471-492

[13] Elmes GW, Thomas JA, Wardlaw JC, Hochberg ME, Clarke RT, Simcox DJ. 1998. The ecology of Myrmica ants in relation to the conservation of Maculinea butterflies. Journal of Insect Conservation 2:67-78

[14] Evanno G, Regnaut S, Goudet J. 2005. Detecting the number of clusters of individuals using the software STRUCTURE: a simulation study. Molecular Ecology 14:2611-2620

[15] Excoffier L, Smouse PE, Quattro JM. 1992. Analysis of molecular variance inferred from metric distances among DNA haplotypes—application to human mitochondrial-DNA restriction data. Genetics 131:479-491

[16] Falush D, Stephens M, Pritchard JK. 2003. Inference of population structure using multilocus genotype data: linked loci and correlated allele frequencies. Genetics 164:1567-1587

[17] Falush D, Stephens M, Pritchard JK. 2007. Inference of population structure using multilocus genotype data: dominant markers and null alleles. Molecular Ecology Notes 7:574-578

[18] Fric Z, Wahlberg N, Pech P, Zrzavy J. 2007. Phylogeny and classification of the Phengaris–Maculinea clade (Lepidoptera: Lycaenidae): total evidence and phylogenetic species concepts. Systematic Entomology 32:558-567

[19] Gao H, Williamson S, Bustamante CD. 2007. A Markov chain Monte Carlo approach for joint inference of population structure and inbreeding rates from multilocus genotype data. Genetics 176:1635-1651

[20] Goudet J. 1995. FSTAT (version 1.2): a computer program to calculate F-statistics. Journal of Heredity 86:485-486

[21] Goudet J. 2005. HIERFSTAT, a package for R to compute and test hierarchical F-statistics. Molecular Ecology Notes 5:184-186

[22] Habeler H. 2008. Die subalpin-alpinen Lebensräume des Bläulings Maculinea rebeli (Hirscke, 1904) in den Ostalpen (Lepidoptera, Lycaenidae) Joannea - Zoologie 10:143-164

[23] Höttinger H, Schlick-Steiner BC, Steiner FM. 2003. The Alcon blue Maculinea alcon (Lepidoptera: Lycaenidae) in eastern Austria: status and conservation measures. Ekologia 22:107-118

[24] Hubisz MJ, Falush D, Stephens M, Pritchard JK. 2009. Inferring weak population structure with the assistance of sample group information. Molecular Ecology Resources 9:1322-1332

[25] Jones OR, Wang J. 2010. COLONY: a program for parentage and sibship inference from multilocus genotype data. Molecular Ecology Resources 10:551-555

[26] Jost L. 2008. G_ST and its relatives do not measure differentiation. Molecular Ecology 17:4015-4026

[27] Kőrösi Á, Örvössy N, Batáry P, Kövér S, Peregovits L. 2008. Restricted within-habitat movement and time-constrained egg laying of female Maculinea rebeli butterflies. Oecologia 156:455-464

[28] Kudrna O, Belicek J. 2005. On the Wiener Verzeichnis, its authorship and the butterflies named therein. Oedippus 23:1-32

[29] Kudrna O, Fric ZF. 2013. On the identity and taxonomic status of Lycaena alcon rebeli HIRSCHKE, 1905—a long story of confusion and ignorance resulting in the fabrication of a ghost species (Lepidoptera: Lycaenidae) Nachrichten des Entomologischen Vereins Apollo 34:117-124

[30] Loiselle BA, Sork VL, Nason J, Graham C. 1995. Spatial genetic structure of a tropical understory shrub, Psychotria officinalis (Rubiaceae) American Journal of Botany 82:1420-1425

[31] Meirmans PG, Hedrick PW. 2011. Assessing population structure: F_ST and related measures. Molecular Ecology Resources 11:5-18

[32] Meirmans PG, Van Tienderen PH. 2004. GENOTYPE and GENODIVE: two programs for the analysis of genetic diversity of asexual organisms. Molecular Ecology Notes 4:792-794

[33] Meyer-Hozak C. 2000. Population biology of Maculinea rebeli (Lepidoptera: Lycaenidae) on the chalk grassland of Eastern Westphalia (Germany) and implications for conservation. Journal of Insect Conservation 4:63-72

[34] 1999. Action plan for the Maculinea butterflies in Europe. In: Munguira ML, Martín J, eds. Nature and environment. Strasbourg: Council of Europe Publishing. 97

[35] Nash DR, Als TD, Maile R, Jones GR, Boomsma JJ. 2008. A mosaic of chemical coevolution in a large blue butterfly. Science 319:88-90

[36] Oostermeijer JGB, Vantveer R, Dennijs JCM. 1994. Population structure of the rare, long lived perennial Gentiana pneumonanthe in relation to vegetation and management in the Netherlands. Journal of Applied Ecology 31:428-438

[37] Peakall R, Smouse PE. 2012. GenAlEx 6.5: Genetic analysis in Excel. Population genetic software for teaching and research—an update. Bioinformatics 28:2537-2539

[38] Pech P, Fric Z, Konvicka M, Zrzavy J. 2004. Phylogeny of Maculinea blues (Lepidoptera : Lycaenidae) based on morphological and ecological characters: evolution of parasitic myrmecophily. Cladistics 20:362-375

[39] Pecsenye K, Bereczki J, Tihanyi B, Toth A, Peregovits L, Varga Z. 2007. Genetic differentiation among the Maculinea species (Lepidoptera: Lycaenidae) in eastern Central Europe. Biological Journal of the Linnean Society 91:11-21

[40] Pritchard JK, Stephens M, Donnelly P. 2000. Inference of population structure using multilocus genotype data. Genetics 155:945-959

[41] Queller DC, Goodnight KF. 1989. Estimating relatedness using genetic markers. Evolution 43:258-275

[42] Radchenko AG, Elmes GW. 2010. Myrmica ants (Hymenoptera: Formicidae) of the old world. Warsaw: Natura Optima Dux Foundation.

[43] Rákosy L, Tartally A, Goia M, Mihali C, Varga Z. 2010. The Dusky Large Blue—Maculinea nausithous kijevensis (Sheljuzhko, 1928) in the Transylvanian basin: new data on taxonomy and ecology. Nota Lepidopterologica 33:31-37

[44] Seifert B. 1988. A taxonomic revision of the Myrmica species of Europe, Asia Minor, and Caucasus (Hymenoptera, Formicidae) Abhandlungen und Berichte des Naturkundemuseums Görlitz 62:1-75

[45] 2005. Settele J, Kühn E, Thomas JA, eds. Studies on the ecology and conservation of butterflies in Europe - Vol. 2: species ecology along a European gradient: Maculinea butterflies as a model. Sofia: Pensoft.

[46] Sielezniew M, Rutkowski R, Ponikwicka-Tyszko D, Ratkiewicz M, Dziekanska I, Svitra G. 2012. Differences in genetic variability between two ecotypes of the endangered myrmecophilous butterfly Phengaris (=Maculinea) alcon—the setting of conservation priorities. Insect Conservation and Diversity 5:223-236

[47] Sielezniew M, Stankiewicz AM. 2004. Simultaneous exploitation of Myrmica vandeli and M. scabrinodis (Hymenoptera: Formicidae) colonies by the endangered myrmecophilous butterfly Maculinea alcon (Lepidoptera: Lycaenidae) European Journal of Entomology 101:693-696

[48] Steiner FM, Sielezniew M, Schlick-Steiner BC, Höttinger H, Stankiewicz A, Górnicki A. 2003. Host specificity revisited: new data on Myrmica host ants of the lycaenid butterfly Maculinea rebeli. Journal of Insect Conservation 7:1-6

[49] Tartally A. 2008. Myrmecophily of Maculinea butterflies in the Carpathian Basin (Lepidoptera: Lycaenidae) PhD thesis, University of Debrecen. Available at http://ganymedes.lib.unideb.hu:8080/dea/bitstream/2437/78921/5/ertekezes.pdf thesis

[50] Tartally A, Koschuh A, Varga Z. 2014. The re-discovered Maculinea rebeli (Hirschke, 1904): host ant usage, parasitoid and initial food plant around the type locality with taxonomical aspects (Lepidoptera, Lycaenidae) ZooKeys 406:25-40

[51] Tartally A, Nash DR, Lengyel S, Varga Z. 2008. Patterns of host ant use by sympatric populations of Maculinea alcon and M. ’rebeli’ in the Carpathian Basin. Insectes Sociaux 55:370-381

[52] Tartally A, Szűcs B, Ebsen JR. 2007. The first records of Rickia wasmannii Cavara, 1899, a myrmecophilous fungus, and its Myrmica Latreille, 1804 host ants in Hungary and Romania (Ascomycetes: Laboulbeniales; Hymenoptera: Formicidae) [Abstract]. Myrmecological News 10:123

[53] Tartally A, Varga Z. 2008. Host ant use of Maculinea teleius in the Carpathian-Basin (Lepidoptera: Lycaenidae) Acta Zoologica Academiae Scientiarum Hungaricae 54:257-268

[54] Thomas JA, Elmes GW. 1998. Higher productivity at the cost of increased host-specificity when Maculinea butterfly larvae exploit ant colonies through trophallaxis rather than by predation. Ecological Entomology 23:457-464

[55] Thomas JA, Elmes GW, Schönrogge K, Simcox DJ, Settele J. 2005. Primary hosts, secondary hosts and ‘non-hosts’: common confusions in the interpretation of host specificity in Maculinea butterflies and other social parasites of ants. In: Settele J, Kuhn E, Thomas JA, eds. Studies on the ecology and conservation of butterflies in Europe - Vol. 2: species ecology along a European gradient: Maculinea butterflies as a model. Sofia: Pensoft. 99-104

[56] Thomas JA, Elmes GW, Wardlaw JC, Woyciechowski M. 1989. Host specificity among Maculinea butterflies in Myrmica ant nests. Oecologia 79:452-457

[57] Thomas JA, Settele J. 2004. Evolutionary biology: butterfly mimics of ants. Nature 432:283-284

[58] Timuş N, Craioveanu C, Sitaru C, Rus A, Rákosy L. 2013. Differences in adult phenology, demography, mobility and distribution in two syntopic ecotypes of Maculinea alcon (cruciata vs. pneumonanthe) (Lepidoptera: Lycaenidae) from Transilvania (Romania) Entomologica Romanica 18:21-30

[59] Ugelvig LV, Andersen A, Boomsma JJ, Nash DR. 2012. Dispersal and gene flow in the rare, parasitic Large Blue butterfly Maculinea arion. Molecular Ecology 21:3224-3236

[60] Ugelvig LV, Nielsen PS, Boomsma JJ, Nash DR. 2011a. Reconstructing eight decades of genetic variation in an isolated Danish population of the large blue butterfly Maculinea arion. BMC Evolutionary Biology 11:201

[61] Ugelvig LV, Vila R, Pierce NE, Nash DR. 2011b. A phylogenetic revision of the Glaucopsyche section (Lepidoptera: Lycaenidae), with special focus on the Phengaris-Maculinea clade. Molecular Phylogenetics and Evolution 61:237-243

[62] Van Oosterhout C, Hutchinson WF, Wills DPM, Shipley P. 2004. MICRO-CHECKER: software for identifying and correcting genotyping errors in microsatellite data. Molecular Ecology Notes 4:535-538

[63] Van Swaay C, Cuttelod A, Collins S, Maes D, Munguira ML, Šašic M, Settele J, Verovnik R, Verstrael T, Warren M, Wiemers M, Wynhof I. 2010. European red list of butterflies. Brussels: IUCN.