Host tracking or cryptic adaptation? Phylogeography of Pediobius saulius (Hymenoptera, Eulophidae), a parasitoid of the highly invasive horse-chestnut leafminer.

Classical biological control is often advocated as a tool for managing invasive species. However, accurate evaluations of parasitoid species complexes and assessment of host specificity are impeded by the lack of morphological variation. Here, we study the possibility of host races/species within the eulophid wasp Pediobius saulius, a pupal generalist parasitoid that parasitize the highly invasive horse-chestnut leaf-mining moth Cameraria ohridella. We analysed the population genetic structure, host associations and phylogeographic patterns of P. saulius in Europe using the COI mitochondrial gene. This marker strongly supports a division into at least five highly differentiated parasitoid complexes, within two of which clades with differing degrees of host specialization were found: a Balkan clade that mainly (but not only) attacks C. ohridella and a more generalist European group that attacks many hosts, including C. ohridella. The divergence in COI (up to 7.6%) suggests the existence of cryptic species, although this is neither confirmed by nuclear divergence nor morphology. We do not find evidence of host tracking. The higher parasitism rates observed in the Balkans and the scarcity of the Balkan-Cameraria haplotypes out of the Balkans open the possibility of using these Balkan haplotypes as biological control agents of C. ohridella elsewhere in Europe.

Introduction

One of the most important biological interactions that shapes biodiversity is parasitism, with up to half of known species globally being parasites or parasitoids (Windsor 1998; Poulin and Morand 2004). This overwhelming diversity has been attributed to the intimate and specialized relationship parasites have with a limited number of host species, as a general feature of parasite life histories is host specificity (Fox and Morrow 1981). This entails the potential for disruptive selection associated with adaptation to different hosts (Price 1980; Ackermann and Doebeli 2004).

Based on the observed predominance of specialists over generalists in phytophagous and parasitoid insects (Futuyma and Moreno 1988; Mitter et al. 1988; Berenbaum 1990; Jaenike 1990; Groman and Pellmyr 2000; Weiblen and Bush 2002; Henry et al. 2008; Cook and Segar 2010), it has been postulated that specialization must be a widespread evolutionary trend, with generalists becoming rarer as biological radiations age (Simpson 1953; Schluter 2000). Parasitoids are mainly found in the orders Hymenoptera and Diptera and have an extremely intimate relationship with hosts because these harbour the parasitoid's offspring until maturity (Godfray 1994). At the family level, insect parasitoids have undergone extensive adaptive radiations as evidenced by the large numbers of species in the major parasitoid clades (Quicke 1997; Godfray and Shimada 1999; Irwin et al. 2003). In particular, Hymenoptera are the most species-rich group of parasitoid insects and might account for up to 20% of all insect species (La Salle 1993). However, this diversity is likely to be overlooked owing to the identification of morphologically similar, but genetically and ecologically isolated lineages as single, generalist species (Herre 2006; Bickford et al. 2007; Smith et al. 2008).

In this paper, we study the case of Pediobius saulius Walker, 1893 (Chalcidoidea, Eulophidae), a widespread, generalist parasitoid of leaf-mining insects throughout the Palearctic (Viggiani 1964). It is known as a primary ecto-, endo- and hyper-parasitoid of larvae and pupae of over 107 different host species of three insect orders: Coleoptera (11 species), Lepidoptera (76 species) and Hymenoptera (20 species) (Noyes 2003). Its main lepidopteran hosts are Gracillariidae leaf-mining moths (59 species). Pediobius saulius is the main parasitoid of the invasive horse-chestnut leaf-mining moth Cameraria ohridella Deschka & Dimic, 1984. Ever since its discovery in Macedonia in 1984 (with evidence of outbreaks since at least 1961 and a history in the Balkans dating back to at least 1879; Lees et al. 2011), this micromoth has become a pest of horse-chestnut trees (Aesculus hippocastanum) and has experienced a rapid range expansion, progressively colonizing almost all of Europe (Augustin et al. 2010). Both mitochondrial and nuclear markers have identified the southern Balkans as the most likely area of origin of C. ohridella (Valade et al. 2009; Lees et al. 2011) and additional data from historical herbarium collections from Albania and Greece now render any other hypothesis for its origin highly improbable (Lees et al. 2011). In all neocolonized regions where the moth is present, C. ohridella maintains permanent outbreak densities, causing severe aesthetic damage to the horse-chestnut, a highly valued ornamental and amenity tree (Freise and Heitland 2004). Until now, surveys in the Balkans have failed to identify specific natural enemies that could be used as biological control agents, in particular parasitoids, and parasitism rates are lower than those usually observed in native leaf miners (Grabenweger et al. 2005a). Research on parasitoid complexes of C. ohridella conducted during the past 10 years shows significant differences between the Balkans and the rest of Europe in the prevalence of P. saulius. This wasp is the dominant species in both artificial (gardens and parks) and natural (moist aspects of mountains and canyons) horse-chestnut populations where C. ohridella occurs in the Balkans (Grabenweger et al. 2010), while in central and western Europe, P. saulius is a common parasitoid of other leaf miners (Noyes 2002; Girardoz et al. 2007b) but rarely attacks C. ohridella (Freise et al. 2002; Stojanovic and Markovic 2004; Grabenweger et al. 2005b; Lupi 2005; Volter and Kenis 2006; Girardoz et al. 2007a; Grabenweger et al. 2010). Girardoz et al. (2007b) suggested the possible existence of cryptic host-differentiated races or sibling species of the wasp, based on the observation of high parasitism of P. saulius on the plane leaf miner Phyllonorycter platani, while P. saulius was totally absent from sympatric populations of C. ohridella in Switzerland. The possibility of undetected species was also suggested by Grabenweger et al. (2010), as a consequence of successful biological, behavioural or phenological adaptations to the new host.

In this study, we use mtDNA (COI), nuclear sequence data (28S-D2, ITS2) and morphology to test whether (i) P. saulius is a complex of specialized host races/cryptic species and (ii) whether a host-specific Balkan race/species of P. saulius has tracked the expansion of its host C. ohridella into central and Western Europe. The potential discovery of a Balkan host-specific race or cryptic species of P. saulius would open the possibility for its use as a biological control agent against the horse-chestnut leaf miner across Europe.

Materials and methods

Sampling and morphology analysis

Samples of P. saulius were collected from 38 localities in 10 European countries (n = 152) (Fig. 1, Appendix S1 in Supplementary Material). A total of 63 individuals were obtained from Balkan populations and 83 from central and Western Europe (Table 1). One hundred and four P. saulius adult wasps were reared from C. ohridella mines on A. hippocastanum, while 42 specimens were reared from other hosts: 12 adults were reared from P. platani on Platanus sp., nine from Phyllonorycter milleriella on Celtis australis, four from Phyllonorycter froelichiella on Alnus glutinosa, two from Phyllonorycter quercifoliella on Quercus spp., one from Phyllonorycter abrasella, one hyper-parasitoid from Pediobius sp., 11 from Orchestes quercus on Quercus spp., one from Orchestes fagi on Fagus sylvatica and one from Phyllonorycter coryli on Corylus avellana (Table 1, Supplementary Material Appendix S1). Samples from mined leaves were reared as described previously (Lopez-Vaamonde et al. 2005). Reared specimens were kept in 70% ethanol. The gaster of each specimen was removed for DNA analysis, and the remaining body was mounted on a rectangular card following Noyes (1982).

Figure 1

Map showing the 38 localities sampled in our study (for details about localities see Supplementary Material, Appendix S1).

Table 1

Number of individuals, haplotype designation, host associations and genetic diversity for sampled populations grouped according to geographical origin. Haplotypes present in the two regions are highlighted in bold characters

Region	Number of individuals sampled	Number of localities	Host	Number of haplotypes	Distribution of haplotypes	Haplotype diversity (±SD)	Nucleotide diversity (±SD)
Balkans	63	16	Cameraria ohridella (47) Phyllonorycter sp. (16)	28	H3(16), H6, H7, H9, H11, H12, H13, H15(2), H16, H18(17), H20, H25, H30, H32(2), H34(2), H36, H39(2), H40, H41, H43, H44, H56, H60, H61, H64, H66, H67, H70	0.85070 ± 0.059	0.01088 ± 0.00219
Rest of Europe	83	21	C. ohridella (57) Phyllonorycter sp. (15) Orchestes sp. (12)	58	H1, H2, H4(3), H5, H8(2), H10(10), H14, H17(2), H18, H19(9), H21(3), H22, H23, H24, H26, H27, H28, H29, H31, H33, H34(2), H35, H37, H38, H42, H45, H46, H47, H48, H49, H50, H51, H52, H53, H54(2), H55, H56, H57, H58, H59, H62, H63, H65, H68, H69, H71, H72, H73, H74, H75, H76, H77, H78, H79, H80, H81, H82, H83	0.97424 ± 0.011	0.02194 ± 0.00106
Total	146	37	C. ohridella (104) Phyllonorycter sp. (31) Orchestes quercus (11) Orchestes fagi (1)	83		0.80511 ± 0.020	0.01963 ± 0.00106

The DNA-barcoded and card-mounted specimens were analysed for any kind of variation in the external morphology, structural as well as colour-based. Specimens used for scanning electronic microscopy (SEM) and stereomicroscope photographs were selected from the barcoded card-mounted specimens (Supplementary Material Appendix S2). SEM micrographs were taken from uncoated specimens on their original card in a low-vacuum mode on a JEOL®JSM 5600LV SEM microscope (JEOL Skandinaviska, AB, Hammarbacken, Sollentuna, Sweden). The same specimens were used for the stereomicroscope photographs. In the stereomicroscope (Nikon SMZ 1500; Nikon Nordic AB, Solna, Sweden), several photographs with different focus of depth (i.e. different layers) were taken (with camera Nikon DS-Fi1; Nikon Nordic AB, Solna, Sweden) and were subsequently merged into one single picture with the software Helicon Focus 4.46. A morphometric analysis was performed based on measurements taken from 16 Balkan and 27 central and western European randomly selected specimens (Supplementary Material Appendix S3). To evaluate whether the material could be consistently sorted based on morphology, ratios from the measurements were calculated and plotted against each other (Supplementary Material Appendix S3.1).

To add to our morphological analysis, the newly discovered wing interference patterns (WIPs) (Shevtsova et al. 2011) were analysed (Supplementary Material Appendix S4). The WIP photographs (i.e. colour photographs of wings) were taken with another stereomicroscope [Nikon SMZ 1000 (Nikon Nordic AB, Solna, Sweden) with camera Nikon DS-Fi1]. To get comparable pictures of the wings, these were removed from the specimens and placed on a slide with a cover glass on top of them (to keep them flat). The slide was then placed against a black background, and 3–4 photographs per wing, in different layers, were taken. These layers were subsequently merged as described above. After this procedure, the wings were glued to the card with the specimen to which they belong.

As discussed by Shevtsova et al., it is the pattern in the wing that is taxon specific, not the hue. The colour varies with the thickness of the wing, and the thickness in turn depends on the size of the specimen, so large specimens and small specimens of same species have different colours, but the pattern remains the same. The morphometric and WIP analyses were supplemented with an ocular examination through a stereomicroscope of all barcoded specimens for colour differences.

DNA sequencing

DNA was extracted from adult abdomens, using both the DNAeasy tissue kit (Qiagen, Courtaboeuf, France) and using a silica-based 96-well automated extraction protocol (Ivanova et al. 2006). A total of 130 individuals were sequenced with the primer set LepF1/LepR1 (Hajibabaei et al. 2006, as ‘LepF/LepR’). The same fragment was amplified using a slightly different set of primers (LCO1490 and HCO2198; Folmer et al. 1994) for a subset of 16 samples. PCR amplifications were performed according to the standard protocol used in CCDB (Hajibabaei et al. 2006). A total of 25 individuals (from the main CO1 haplotype groups) were successfully sequenced for a 1005-bp fragment of the nuclear 28S large ribosomal subunit. PCR and sequencing were performed as described previously (Lopez-Vaamonde et al. 2001). Additionally, a set of 10 individuals collected from P. platani, O. quercus and C. ohridella (five of which were already barcoded) were sequenced for a 533-bp fragment of ITS2, using the primers CAS5p8sFc and CAS28sB1d (Ji et al. 2003). PCR conditions for ITS2 were denaturation at 95°C for 3 min, 35 cycles at 94°C for 45 s, 72°C for 1 min and a final extension step at 72°C for 7 min.

Genetic analyses

DNA sequences were edited using CodonCode Aligner version 3.0.1 (CodonCode Corporation, Dedham, MA, USA). No indels and stop codons were present in the alignments. COI sequences were trimmed to the same length (556 bp) to eliminate missing/bad data. A 452-bp fragment of the D2 expansion region of the nuclear 28S was used because the D1 and D3 expansion regions showed no variability. The ITS2 amplicons ranged in size from 465 to 534 bp. Records for all 146 barcoded specimens used in our analyses are gathered within the project ‘Phylogeography of Pediobius saulius’ (code PEDIO) in the Published Projects section of BOLD (Ratnasingham and Hebert 2007). Information on specimen vouchers (field data and GPS coordinates) and sequences (nucleotide composition and trace files) is found in this project by following the ‘view all records’ link and clicking on the ‘specimen page’ or ‘sequence page’ links for each individual record.

Sequences were aligned unambiguously by eye. Polymorphism analyses (haplotype and nucleotide diversity; (Nei 1987)) and Tajima's D (Tajima 1989) were performed using DnaSP version 5.0.3 (Rozas et al. 2003).

To test the effect of geography and host associations on the distribution of genetic diversity, samples were grouped according to origin (Balkans/rest of Europe) and to host (C. ohridella– other hosts). The Balkans/rest of Europe division takes into account the difference in parasitism rates on C. ohridella in these regions and assumes that P. saulius is a single species. We compared the partition of genetic variability among populations and among groups of populations by an analysis of molecular variance (Excoffier et al. 1992), estimated by computing conventional F-statistics for COI haplotypes using Arlequin version 3.1 software with 10 000 permutations. Detailed information about the history and relationships between haplotypes was obtained through Bayesian phylogenetic reconstructions. We analysed both COI and 28S-D2 data sets to compare the level of phylogenetic resolution and to compare the level of congruence between the markers. We used jModelTest version 0.0.1 (Posada 2008) to choose a model of nucleotide substitution for both data sets. The best fit chosen was a general time-reversible model with gamma-distributed rate variation and a proportion of invariant sites (GTR + I + χ). Bayesian phylogenetic reconstructions were performed using mrbayes version 3.1.2 (Ronquist and Huelsenbeck 2003) using a cold chain and three heated chains with T = 0.1. Starting trees for each chain were random, and we used mrbayes default settings (including prior distributions). Two independent metropolis coupled Markov chain Monte Carlo (MCMC) were run for 5 million generations, with trees sampled every 1000 generations. Burn-in was set at 1.5 million generations, at which point chains had converged to stable likelihood values <0.01. Posterior probabilities were used to assess clade support.

Quantitative species delimitation was performed with the general mixed Yule-coalescent (GMYC) method (Pons et al. 2006) to evaluate the possible presence of cryptic species. The method is implemented in an r package called GMYC (available from http://r-forge.r-project.org/projects/splits/). The GMYC model attempts to identify a species ‘boundary’ as a shift in branching rates on a tree that contains multiple species and populations. Branching patterns within the genetic clusters reflect neutral coalescent processes occurring within species (Kingman 1982), whereas branching among genetic clusters reflects the timing of speciation events (Yule 1924). The GMYC exploits the predicted difference in branching rate under the two modes of lineage evolution, assessing the point of highest likelihood of the transition (Pons et al. 2006). Independently evolving lineages are recovered as putative species, more distinct than predicted if the entire sample derived from a single species without genetic isolation (Acinas et al. 2004; Barraclough et al. 2003). We applied the GMYC method on an ultrametric ‘all-compatible’ consensus tree including only unique haplotypes, generated from FigTree version 1.3.1 after running the COI data set in Beast version 1.4.8. from an xml file created under Beauti version 1.4.8, both programs from (http://beast.bio.ed.ac.uk/). We used a GTR model with (1 + 2) and 3rd positions analysed separately, with a random starting tree under a coalescent process and an enforced molecular clock (relaxed lognormal), and after running an MCMC chain for 10 000 000 iterations. We checked for convergence of all parameters to an adequate effective sample size using this number of iterations using Tracer version 1.5 (http://beast.bio.ed.ac.uk/). A single threshold was applied in the GMYC script.

Results

Morphological variation

The analysis of SEMs does not reveal any significant morphological differences (Supplementary Material Appendix S3). The morphometric analysis did not show any separation into clusters. On the contrary, specimens from the two clades (Balkans versus Central and north-west Europe) were completely mixed in all graphs (Supplementary Material Appendix S3.1). The WIPs showed two things (Supplementary Material Appendix S3): (i) there were no notable differences in the colour patterns between specimens of the same sex from the two clades; (ii) there was a small but distinct sexual dimorphism in the WIPs, which was consistent in both clades. The morphometric and WIP analyses were supplemented with an ocular examination through a stereomicroscope of all barcoded specimens for colour differences. We could not find any difference between specimens belonging to the two clades, also in this analysis, specimens from the two clades were similar.

Mitochondrial DNA variation

A total of 83 haplotypes were found among the 146 individuals analysed in all localities, with 69 represented by a single individual (Fig. 2). Haplotype 18 (H18) is the commonest with 18 individuals, representing 21.69% of the total number of individuals, followed by H3 (n = 16; 19.28%), H10 (n = 10; 12.05%) and H19 (n = 9; 10.84%, Table 1). There is a clear divide between Balkan and European haplotypes, as only three haplotypes (H18, H34, H56) are present in both geographic regions. H34 and H56 have the same frequencies in both regions (two and one individuals, respectively), and H18 is predominant in the Balkans (17 individuals compared to one in Europe). Fifty-five haplotypes are restricted to central and northern European populations (Hd = 0.85070 ± 0.059), and 25 are exclusive to the Balkans (Hd = 0.97424 ± 0.011; Table 1, Fig. 2). Furthermore, only seven haplotypes (H3, H10, H17, H18, H19, H34 and H56) are present in more than one country (Fig. 2).

Figure 2

Geographic distribution of the 83 haplotypes among the 10 sampled countries. Each pie represents a country, and colours represent haplotypes present in more than one country. Noncoloured triangles represent unique haplotypes. The number of individuals sampled per country is proportional to the size of the pie charts.

The amova revealed that most of the variation is explained by differences between individuals within populations (62.83%), with a much smaller, but still significant, amount of variation (20.56%) found among populations within groups. When populations were clustered in Balkans versus the rest of Europe, variation among groups was found to be not significant (Table 2). When samples were grouped according to host, results were similar with variation within populations (57.96%) greater than that observed for populations within groups (36.42%). Variation among host groups was not significant (5.62%Table 2).

Table 2

Results of amova test on COI data. Comparisons were made for populations grouped according to (1) geographic origin (Balkans versus Rest of Europe) and (2) host association (Cameraria ohridella versus Phyllonorycter spp. + Orchestes spp.)

	Source of Variation	Variance components	Variation (%)
(1) Balkans versus Rest of Europe	Among groups	1.63602	16.61NS
	Among populations within groups	2.02610	20.56*
	Within populations	6.19013	62.83*
(2) Cameraria versus Other hosts	Among groups	0.53157	5.62NS
	Among populations within groups	3.44394	36.42*
	Within populations	5.48138	57.96*

P < 0.001; NS, not significant.

Geographic and host-associated patterns of mtDNA variation

Most of the Balkan individuals (45 out of 63) are grouped in a monophyletic clade (C1) in the COI Bayesian tree (Fig. 3A). Furthermore, 33 individuals of this clade belong to Balkan-only haplotypes (H3 and H18), including an H18 individual from the Czech Republic. Only one of the individuals in C1 was collected on P. platani, and all the others were reared from C. ohridella (Fig. 3B). The posterior probability value for the C1 node, however, is <0.90. Of the 18 Balkan individuals outside C1, 14 were reared from, and are clustered with, Phyllonorycter clades (Fig. 3B). There are eight monophyletic groups that contain almost exclusively central and western European individuals (Figs 3 and 4). Interestingly, there is a well-supported monophyletic clade (C2) with a mixture of Balkan and European individuals that seem to be specialized on Phyllonorycter (Figs 3B and 4).

Figure 3

Bayesian phylogenetic relationships for the 146 COI samples. Trees are presented in circular form with mid-point rooting for better visualization (outgroups not shown). Black-and-white circles represent the area of origin (A), and coloured bars represent host associations (B). C1 is the Balkan-Cameraria clade; C2 is the Balkn/rest of Europe-Phyllonorycter clade (see text for details).

Figure 4

Bayesian phylogenetic relationships for COI haplotypes. Coloured branches represent the independently evolving units estimated by GMYC analyses. Significant posterior probability values are represented by asterisks (*≥90%; **≥95%). Red bars on nodes highlight clades 1 and 2 (C1, C2) and correspond to those in Fig. 3.

Overall, genetic diversity of P. saulius populations was higher in central and north-west Europe than in the Balkans (Table 1), with the median number of haplotypes being 10.16 (1–23; n = 21 localities) in central and north-west Europe and 6.75 (range 6–13; n = 17 localities) in the Balkans. Similarly, both nucleotide and haplotype diversity were higher for central and western European populations (Table 1). The highest genetic diversity was found in France, where we found 23 haplotypes out of 25 individuals analysed (Supplementary Material Appendix S1, Fig. 1). Tajima's D was negative (D = −1.03907) and was not significantly different from zero (P > 0.10), further supporting the lack of important demographic structure (i.e. no major founder effects or sudden expansion).

Phylogenetic relationships and GMYC analysis

The Beast analysis of the COI haplotypes data set resulted in a number of well-supported groups (PP > 95%, Fig. 4). Deep relationships were poorly resolved. The maximum K2P genetic divergence between haplotypes was of 7.6% (between H61 and H69). Species delimitation for COI using the GMYC model revealed a total of five independent coalescent ML clusters (six ML entities; nodes A-E plus H45; Fig. 4), all but one with PP support >0.95. These clusters do not match the two clades C1 and C2 as previously highlighted (Fig. 3) based on their consistency in terms of geography and host specificity, but both are grouped within different coalescent clusters. Average divergence between GMYC clades was 2.8%. Lineage-through-time plots (LTT) for the Bayesian CO1 tree show a single recent shift to a higher branching rate that occurred near the branch tips (not shown). This matches the sharp peak in the likelihood surface for the GMYC model predicted for independently evolving species.

The 28S data set for 29 of the haplotypes resulted in sequences of 1005 bp long, of which 452 bp corresponded to the 28S-D2 ribosomal region. The alignment was indel-free, and genetic distances were lower than those for COI with the largest being of 1.04% between H37 and H35. Node support for the Bayesian tree was high only for one clade (Supplementary Material Appendix S5). Overall, 28S sequences showed little variability, and very divergent haplotypes for CO1 show no genetic divergence for 28S. GMYC species delimitation for 28S resulted in only two coalescent groups, one consisting of the outgroup and the other of all the P. saulius individuals. It should be noted that our 28S data set included haplotypes representing all five GMYC clusters, but no representative of clade 1 (Fig. 3, C1). The ITS2 sequences were 465–534 bp long, with no indels in the alignment. Differentiation was even smaller than for 28S, with the largest genetic distance of 0.4% between an individual collected in France on O. quercus and seven individuals from the Balkans collected on both P. platani and C. ohridella. The Bayesian tree resulted in a comb with virtually no differences between individuals, regardless of host or geographic origin (Supplementary Material Appendix S6).

Discussion

Is Pediobius saulius a complex of cryptic species?

The intraspecific genetic diversity shown by the COI barcode region was larger for central and north-western Europe than for the Balkans, and overall, it is much higher than typical values obtained for other Hymenoptera (Stone et al. 2001; Rokas et al. 2003; Stone et al. 2007; Smith et al. 2008). Considering recent findings based on DNA analysis of parasitoid wasp communities (e.g. Smith et al. 2008, 2009), revealing high levels of cryptic diversity, these results question the identification of P. saulius as a broadly distributed highly generalist species.

We found evidence of host-associated genetic differentiation shown by the Balkan haplotypes H3 and H18, which are clustered in a Balkan-only clade associated with C. ohridella (only two non-Balkan C. ohridella-specific haplotypes were clustered in this clade; C1, Figs 3 and 4). All but four of the Balkan individuals outside this clade were found on Phyllonorycter hosts (Fig. 3). Further evidence for host-associated differentiation is the fact that the sister group to the Balkan–Cameraria group (clade 1) consists of 15 central and western European individuals, of which all but one were parasitoids of C. ohridella. The fact that the Balkan–Cameraria clade gathers 71.4% of the individuals from the Balkans and that only one of those came from hosts other than C. ohridella is consistent with significantly higher rates of parasitism of C. ohridella by P. saulius in Balkan populations (5.4%) than in north-western and central European populations (0.2%) described in recent Pan-European studies (Grabenweger 2004; Grabenweger et al. 2005b; Girardoz et al. 2007a; Grabenweger et al. 2010). However, most clades are composed of haplotypes from individuals collected on more than one host. Consequently, the amova results for host-associated variation (Table 2) should be interpreted with caution, as sampling effort was much greater for C. ohridella, owing to the low frequencies of P. saulius found in other hosts. To further test the host specificity of the Balkan C. ohridella clade, more sampling is needed from hosts other than Cameraria, in particular from Phyllonorycter, to made comparisons with Central and Western Europe more meaningful. As the low frequency and abundance of P. saulius make extensive sampling difficult, indirect DNA-based techniques can also be used. In particular, molecular analysis of parasitoids’ gut content (MAPL; Rougerie et al. 2010) and barcoding of host tissue remnants (Hrcek et al. 2011) are low-cost and effective alternatives.

The reported existence of a host-specific Balkan clade in P. saulius, supported by higher parasitism rates in the Balkans, raises the question of how this association was formed and where it originated, as there is no evidence of a range expansion from the Balkans into northern Europe or vice versa. A possible explanation is that the Balkan populations constitute a different species. A number of recent studies have uncovered substantial cryptic variation and ecological specificity in a number of insect taxa (Molbo et al. 2003; Hebert et al. 2004; Smith et al. 2005; Challis et al. 2007; Smith et al. 2007; Stone et al. 2008). A general trend in all these studies is that COI is an effective marker for pointing at host- and eco-specific taxa in allegedly morphologically undistinguishable groups. There are, however, some exceptions to this general trend: the tortricid moth Homona mermeroides shows no COI differentiation, confirming it as a polyphagous, widespread species (Hulcr et al. 2007). Likewise, two morphospecies of generalist tachinid parasitoid flies were confirmed as a single generalist species by COI, while a single generalist morphospecies was recovered as two generalist species (Smith et al. 2007). In addition, two species of parasitoid microgastrine wasps were found to be generalists (Smith et al. 2008). Here, our results are particularly intriguing because the COI GMYC analysis identified five independently evolving units, suggesting the potential existence of five generalist (with a high degree of haplotype host specificity within clusters) cryptic species (Fig. 4). Average genetic distances between the units are slightly greater than those for microgastrine wasps (2.8% compared to c. 2%; Smith et al. 2007, 2008).

Lack of differentiation for nuclear markers

In spite of the deep mtDNA divergences, we found little differentiation for 28S-D2 and virtually none for ITS2 within P. saulius. There are two possible alternative explanations for the lack of nuclear divergence with respect to mitochondrial genetic differentiation. One possibility is that only mtDNA is truly divergent in the absence of any strong nuclear DNA differentiation. This would represent a situation of high asymmetrical gene flow owing to a number of causes (symbiont-related genetic sweeps, differences in molecular evolution, hybridization) and might not correspond to a situation of incipient species. Similar incongruence between nuclear and mitochondrial markers has also been reported for other hymenopterans with high mitochondrial and low nuclear divergences (Lopez-Vaamonde et al. 2005; Haine et al. 2006; Smith et al. 2008).

Maternally inherited endosymbionts like Wolbachia pipientis are in linkage disequilibrium with the mtDNA of their hosts and can therefore induce selective sweeps through hitch-hiking, decreasing genetic diversity over many generations (Narita et al. 2006; Galtier et al. 2009; Rodriguero et al. 2010; Sun et al. 2011). Among the many possible effects of bacterial symbionts on mtDNA diversity and evolution, by far the most common is a reduction in mtDNA diversity (Ballard et al. 1996; Jiggins 2003; Dyer and Jaenike 2004; Shoemaker et al. 2004; Hurst and Jiggins 2005; Haine et al. 2006; Rodriguero et al. 2010). In addition to producing selective sweeps, inherited symbiont presence may also promote balancing selection on different mtDNA haplotypes (if different strains in a host species are under balancing selection themselves). When this occurs, the observation is of greater than expected intraspecific (or interpopulation) mtDNA diversity and relatively deep nodes in mtDNA divergence within a species, with particular mtDNA haplotypes associated with particular infection strains (Schulenburg et al. 2002; Shoemaker et al. 2003; Charlat et al. 2009). The latter scenario is consistent with our observations of high mitochondrial diversity within (62.83%) and among populations (20.56%). In both cases, population structure would be expected to correspond to that of the endosymbiont, neither to geography nor to host.

It is also possible that the inconsistency between nuclear and mitochondrial differentiation is attributable to a higher rate of molecular evolution for the latter. It has been observed that the rate of mitochondrial to nuclear molecular evolution in parasitic Hymenoptera is accelerated with respect to other insects. Synonymous positions in mitochondrial protein-coding genes evolve 17–35.5 times faster than nuclear genes in Nasonia. Furthermore, divergence between Nasonia giraulti and Nasonia vitripennis is 3.1% for nuclear genes and 52% for mitochondrial. Also, the relative mitochondrial to nuclear synonymous substitution ratio for N. giraulti–Nasonia longicornis is 35.5, compared with only 2.4 for the Drosophila melanogaster–Drosophila simulans species pair, in spite of the enormous difference in divergence times (much higher for the Drosophila pair; Oliveira et al. 2008). Similarly, molecular evolution comparisons indicate that the parasitic lifestyle is not associated with an increased rate of mtDNA genetic divergence in parasitic Diptera but re-affirm that it is in parasitic Hymenoptera (Castro et al. 2002). Owing to the smaller effective population size of mitochondrial DNA (as haploidy and uniparental inheritance reduce the effective number of mitochondrial genes per locus in a diploid population to about one-quarter that of nuclear; Palumbi et al. 2001; Piganeau and Eyre-Walker 2009), divergence can be accomplished earlier than for nuclear genes after a recent speciation event.

Other explanations are also plausible, like the recent mixing of formerly separated incipient species (i.e. retention of mtDNA ancestral polymorphism; Navajas et al. 1998) or even ancient species (Dyer et al. 2011), naturally large intraspecific variation, or high immigration and hybridization from different populations (Smith et al. 2007; Linnen and Farrell 2007). In addition, the absence of nuclear ribosomal DNA variation has frequently been associated with strong concerted evolution (Elder and Turner 1995; Graur and Li 2000).

Alternatively, the failure to detect nuclear DNA variation could be an artefact owing to the use of inappropriate markers. Both 28S-D2 and ITS2 rRNA nuclear markers have proved useful for species diagnosis in various metazoan taxa, including hymenopteran species (Babcok et al. 2001; Sonnenberg et al. 2007; Smith et al. 2008; Gebiola et al. 2010; Nicholls et al. 2010; Li et al. 2010; Xiao et al. 2010; Zhang et al. in press; but see Heraty et al. 2007) and particularly in closely related eulophid species (Gumovsky 2002; Gebiola et al. 2009, 2010). A nuclear marker that could be considered in future analyses is elongation factor–1α (EF-1α), a conserved nuclear coding gene that can also be used to investigate recent divergences owing to the presence of rapidly evolving introns (Sanchis et al. 2001; Kawakita et al. 2003). Other nuclear markers that have been successfully used to discern species complexes include Rhodopsin and Wingless (Griffiths et al. 2011). An ideal marker to confirm the lack of congruence between nuclear and mitochondrial differentiation is microsatellites, which are far more sensitive in detecting recent or ongoing speciation events (Michel et al. 2010; Kobmoo et al. 2010). All evidence taken together, P. saulius seems to be a single species, composed of highly structured mitochondrial clades that could represent either incipient species/host races or an extreme example of asymmetric nuclear–mitochondrial gene flow.

Has Pediobius saulius tracked the invasion of Cameraria ohridella?

Given that P. saulius was already present in most of Europe prior to the time of C. ohridella's invasion, the only way to test with certainty the origin and possible host tracking of invasive populations by P. saulius would be by analysis of haplotypes through time (i.e. from sample series collected throughout the invasion process, maybe assisted by sampling from herbarium sheets: e.g. by examining abundant C. ohridella samples preserved from Albania in the early 1960s; Lees et al. 2011). This host-tracking scenario implies the spread of C. ohridella-specialized haplotypes/races from the Balkans into central and northern Europe. There is no evidence of such a spread in our study, as only one individual from a predominantly Balkan haplotype (H18) was found in Czech Republic, and only three haplotypes are present in both geographic regions at similar, low abundances. Accordingly, there is no sign of sudden demographic expansions (i.e. colonization events), as shown by the neutral Tajima's D result. Furthermore, amova results show that variation is not related to geographic origin, and it depends more on within-groups and within-population variation (i.e. no geographic structure).

Although the parasitoid indeed appears to have tracked as far as the first Macedonian outbreaks in planted horse-chestnut trees around Ohrid town (Grabenweger and Grill 2000), the fact that there is no signature of host tracking from the Balkans raises the question of why Cameraria-associated haplotypes from the Balkans (clade 1: Fig. 3) have not accompanied the introduction of C. ohridella to central and northern Europe. Because P. saulius is a pupal parasitoid that can be very abundant in overwintering pupae in dead leaves, the fact that there is little sign of host tracking from the Balkans suggests that dead leaves accidentally transported by humans is not a major means of long-distance spread, as previously hypothesized (Gilbert et al. 2005; Augustin et al. 2009). Instead, the spread of C. ohridella could be due rather to flying females carried in vehicles (Heitland et al. 1999; Buszko 2006) or, possibly, eggs and young larvae on seedlings (Gilbert et al. 2005).

It has been shown that the also invasive leafminers Macrosaccus robiniella and P. platani tend to be more heavily parasitized than C. ohridella in sympatric European populations (Girardoz et al. 2007a). Indeed, P. platani and Phyllonorycter leucographella have attained parasitism rates of 37.5% and 56.6%, respectively, at least the former also with tracking of specialists, with <10 years after their introduction into the UK (Godfray et al. 1995). Historical events show that C. ohridella was introduced from the Balkans to Austria by 1989 (Gilbert et al. 2004), from where it spread to Central and Western Europe. Our study suggests that this jump was most probably made without its associated Balkan P. saulius haplotypes. Grabenweger et al. (2010) noted that, in Austria, P. saulius was not found on C. ohridella until 1996, despite intensive samplings. Since then, parasitism rates by this wasp have gradually increased in Eastern Austria, which Grabenweger et al. (2010) attributed either to behavioural, phenological or biological adjustments of a local strain of P. saulius to a novel host, or to the arrival of a Balkan strain or cryptic species of the wasp in the mid-1990s together with the spread of C. ohridella from the Balkans. Our results suggest that both scenarios (lack of host tracking and local P. saulius adaptation to C. ohridella) may have occurred as we found only two ‘Balkan’Cameraria-specific haplotypes (C1, Fig. 3) in Austria and Czech Republic among a majority of Central European haplotypes.

Biological control potential and future prospects

The scarcity of the Balkan–Cameraria haplotypes of P. saulius and the very low parasitism rates in Central and Western Europe compared to the Balkans suggest that the likely spread of the Balkan–Cameraria haplotypes to other regions could be accelerated by intentional introductions. Comparative tests could be carried out between Balkan populations and others to identify bio-ecological and behavioural variations that may explain differences in efficiency (Klug et al. 2008). Laboratory tests could also involve cross-mating tests to check for reproductive isolation (Heraty et al. 2007) between Eastern and Western European haplotypes.

Comprehensive DNA barcode studies have revealed several cases of deep intraspecific divergences in different groups of insects (DeWaard et al. 2011; Hausmann et al. 2011; Park et al. 2011). Detailed integrative studies of those deep mitochondrial splits including morphology, nuclear markers and ecology may reveal a complex of cryptic diversity (Dinça et al. 2011). In our case, the lack of nuclear and morphological divergence calls for deeper studies on the mechanism causing high mitochondrial divergence and whether it is the result of host-associated divergence, geography, selective sweeps or a combination of factors. Nevertheless, the evidence presented here for the existence of a C. ohridella-specialized Balkans-‘race’ of P. saulius opens the possibility of investigating these Balkan haplotypes as potential biological control agents of C. ohridella in the rest of Europe.