Microsatellite and Wolbachia analysis in Rhagoletis cerasi natural populations: population structuring and multiple infections

Rhagoletis cerasi (Diptera: Tephritidae) is a major pest of sweet and sour cherries in Europe and parts of Asia. Despite its economic significance, there is a lack of studies on the genetic structure of R. cerasi populations. Elucidating the genetic structure of insects of economic importance is crucial for developing phenological-predictive models and environmental friendly control methods. All natural populations of R. cerasi have been found to harbor the endosymbiont Wolbachia pipientis, which widely affects multiple biological traits contributing to the evolution of its hosts, and has been suggested as a tool for the biological control of insect pests and disease vectors. In the current study, the analysis of 18 R. cerasi populations collected in Greece, Germany, and Russia using 13 microsatellite markers revealed structuring of R. cerasi natural populations, even at close geographic range. We also analyzed the Wolbachia infection status of these populations using 16S rRNA-, MLST- and wsp-based approaches. All 244 individuals screened were positive for Wolbachia. Our results suggest the fixation of the wCer1 strain in Greece while wCer2, wCer4, wCer5, and probably other uncharacterized strains were also detected in multiply infected individuals. The role of Wolbachia and its potential extended phenotypes needs a thorough investigation in R. cerasi. Our data suggest an involvement of this symbiont in the observed restriction in the gene flow in addition to a number of different ecological factors.


Introduction
The European cherry fruit fly, Rhagoletis cerasi L. (Diptera: Tephritidae) (Fig. 1), is a pest of major agricultural importance infesting mainly sweet (Prunus avium) and sour (P. cerasus) cherries as well as honeysuckle (mainly Lonicera tatarica and L. xylosteum) and wild growing prunus. Rhagoletis cerasi has dispersed and become established in almost all European countries, spreading from the original habitat, the Caucasian area of western Asia (Fimiani 1989;White and Elson-Harris 1992). Rhagoletis cerasi follows a univoltine life cycle with obligatory pupa dormancy that lasts over 9 months. Adults emerge in spring and oviposit usually a single egg on ripe or ripening fruits (Fletcher 1989;Daniel and Grunder 2012). Larvae feed on mesocarp, destroying fruits and causing considerable economic loss exacerbated by secondary fungal and bacterial infections (Fimiani 1989). Because of the economic importance, the control of the European cherry fruit fly received much attention in the framework of the sterile insect technique (SIT) in the early 1970s (Boller et al. 1971. A number of studies have focused on the ecology of R. cerasi (Boller and Prokopy 1976;Boller et al. 1998;Kovanci and Kovanci 2006;Papanastasiou et al. 2011;Moraiti et al. 2012a,b) or the development of environmental friendly control methods, such as mass trapping, application of the oviposition deterring pheromone, bait sprays using plant-derived insecticides and biological control (Katsoyannos et al. 2000;Daniel and Grunder 2012). However, very little is known regarding the biology of this pest at the molecular, genetic, and population levels. The cytogenetics (including the development of polytene chromosome maps) have been recently studied (Kounatidis et al. 2008), though the genetic structure of the European cherry fruit fly populations has received much less attention. Allozyme polymorphisms (Schwarz et al. 2003) raise the question of host race formation in this species, due to host shift, and a small-scale analysis based on microsatellite markers (Augustinos et al. 2011a) points toward elevated genetic differences even in a short geographical range.
In fruit flies (other than Rhagoletis species), microsatellites have been used to shed light on: recent invasion events and dispersion routes of several invasive species (Meixner et al. 2002;Bonizzoni et al. 2004;Augustinos et al. 2005;Nardi et al. 2005;Aketarawong et al. 2007; Khamis et al. 2009;Zygouridis et al. 2009;Virgilio et al. 2010); the degree of genetic structuring in natural populations that can be associated with incipient speciation (Michel et al. 2007); and pattern of remating in nature (Bonizzoni et al. 2002;Song et al. 2007). Lately, microsatellite markers have been also used to address phylogenetic issues in the Bactrocera musae complex of species (Drew et al. 2011).
Microsatellite markers have been recently developed for R. pomonella (Velez et al. 2006), R. indifferens (Maxwell et al. 2009), and R. cerasi (Arthofer et al. 2009a) and used to: (a) address issues regarding speciation events in R. pomonella (Michel et al. 2007(Michel et al. , 2010Forbes et al. 2009), and (b) analyze native and introduced populations of R. completa (Chen et al. 2010). The recent development of 13 microsatellite markers for R. cerasi (Arthofer et al. 2009a), along with the evaluation of cross-amplified microsatellite markers from other Rhagoletis species (Augustinos et al. 2011a), provide a solid basis for exploring the population genetics of this species.
Wolbachia pipientis (Wolbachia for the purpose of this manuscript) is perhaps the most widespread and abundant symbiont among insect species. Extensive genetic diversity has been recorded in the Wolbachia group with the different strains assigned in 13 supergroups, namely A to F and H to N (Augustinos et al. 2011b). The genotyping of the Wolbachia strains is based on either single gene (such as wsp, ftsZ, gltA, groEL, dnaA) or Multi Locus Sequence Typing (MLST) approaches (Casiraghi et al. 2005;Baldo et al. 2006;Paraskevopoulos et al. 2006).
The widespread distribution of Wolbachia, its ability to manipulate reproductive properties, as well as other important physiological functions of its arthropod hosts, have attracted the attention of investigators regarding the role of this symbiont in host biology, ecology, and evolution (Serbus et al. 2008;Werren et al. 2008;Saridaki and Bourtzis 2010). In addition, Wolbachia symbiosis is currently being explored and harnessed toward population control of insect pests and disease vectors (Zabalou et al. 2004(Zabalou et al. , 2009Brelsfoard and Dobson 2009;Apostolaki et al. 2011;Atyame et al. 2011;Walker et al. 2011).
Revealing the close association between R. cerasi and Wolbachia, Riegler and Stauffer shed important light on the reproductive incompatibility reported in the mid-1970s between northwest and southeast European populations of this pest Matolin 1976;Riegler and Stauffer 2002). So far, five Wolbachia strains of R. cerasi have been described (wCer1 to wCer5). The identification of the first two strains, wCer1 and wCer2, was based mainly on wsp gene divergence. Wild R. cerasi populations were found to be infected either with wCer1, or super infected with wCer1 and wCer2. Moreover, wCer2 was proposed to be the causal factor for the reproductive incompatibility of the aforementioned populations, now known as Wolbachia-induced cytoplasmic incompatibility (CI) (Riegler and Stauffer 2002). Based on a polymerase chain reaction (PCR) screen with Wolbachia-specific primers for 16S rRNA and wsp, Kounatidis et al. (2008) reported that the four Greek populations tested were also 100% infected. The Wolbachia strains wCer3, wCer4, and wCer5 were only recently isolated from R. cerasi populations (Arthofer et al. 2009b. Interestingly, wCer4 was identified only after Ceratitis capitata transinfection with Wolbachia, using R. cerasi as donor (Zabalou et al. 2004). Arthofer et al. (2011) proposed Allelic Intersection Analysis as a novel tool for the MLST assignment of multiply infected hosts and described the MLST profile of R. cerasi Wolbachia strains wCer1 to wCer5 (except wCer3), indicating that wCer1, wCer2, and wCer4 belong to Wolbachia supergroup A, while wCer5 belongs to supergroup B. The Wolbachia strain wCer3 is supposed to be a recombinant A/B strain, according to the available wsp gene sequence data (Arthofer et al. 2009b).
Recent studies showed that R. pomonella ) and R. cingulata (Drosopoulou et al. 2011) are infected by Wolbachia, while R. completa seems to be free of Wolbachia (Drosopoulou et al. 2010). Rhagoletis pomonella is infected with the wPom strain, which seems to be very similar (or identical) to wCer2, while other sequence variants for different MLST genes (not attributed to specific strains) have also been reported . Rhagoletis cingulata is infected with at least two strains very similar (or identical) to wCer1 and wCer2 (Schuler et al. 2009;Drosopoulou et al. 2011), named wCin1 and 2 (Schuler et al. 2009).
Wolbachia strains derived from R. cerasi have been recently used for the development of an alternative and environment-friendly strategy to control two major agricultural pests: the Mediterranean fruit fly, C. capitata (Zabalou et al. 2004(Zabalou et al. , 2009 and the olive fly, Bactrocera oleae (Apostolaki et al. 2011). This technique, the Incompatible Insect Technique (IIT), is in principle analogous to the SIT. The main difference between them is the sterilization factor; the IIT is based on Wolbachia-induced CI while the SIT is based on irradiation (Bourtzis and Robinson 2006). IIT holds great potential for the control of agricultural pests and disease vectors (Zabalou et al. 2004;Brelsfoard and Dobson 2009;Saridaki and Bourtzis 2010;Apostolaki et al. 2011;Atyame et al. 2011).
Because of its wide geographical distribution, univoltine life history, patchy distribution of its hosts, low adult dispersion abilities, and infection with different Wolbachia strains (Phillips and Dirks 1933;Jones and Wallace 1955;Boller and Prokopy 1976;Fletcher 1989;Kneifl et al. 1997), R. cerasi bears several properties of becoming an important model for addressing issues regarding effect of ecological factors and symbionts on the differentiation of genetic, demographic, and behavioral traits (Augustinos et al. 2011a;Moraiti et al. 2012b).

Collection of fly samples and DNA extraction
A total of 465 adults from 18 different R. cerasi popualtions were genotyped using 13 microsatellite markers. Fifteen of the populations came from Greece (12 mainland and 3 island) and additionally three from Germany (two populations) and Russia. In all cases, field-infested sweet cherries were sampled and kept in the laboratory until the collection of pupae. Adults emerging from pupae were stored at À80°C in 95% ethanol. Genomic DNA was extracted as described in Augustinos et al. (2011a). Collection sites and the number of flies used in the study are shown in Figure 2 and Table 1.

Microsatellite genotyping and data analysis
PCRs, genotyping, and specific allele sequencing of microsatellite markers were performed as described in Augustinos et al. (2011a). PCR conditions for markers Rce76-1 and 83-44 were similar to those described in Arthofer et al. (2009a).
Genetic variability was measured as the mean number of alleles per locus (n a ), effective number of alleles (n e ), observed (H o ) and expected heterozygosity (H e ), using POPGENE, version 1.31 (Yeh et al. 1999). Allelic diversity was determined after correction for sample size, using FSTAT (Goudet 2001). Deviations from Hardy-Weinberg Equilibrium (HWE) were tested with the G 2 likelihood ratio test in POPGENE. Genotypic disequilibrium was tested with Genepop , using Fisher's exact test, for all pairs of loci in all populations and across them. Presence of null alleles was tested with ML-Null software. Hardy-Weinberg tests for heterozygote deficiency were performed , after performing 10,000 randomization tests using Monte-Carlo randomization as described by Guo and Thomson (1992). Genetic distances were measured according to Nei (1972), using POPGENE (Yeh et al. 1999). Population differentiation was estimated using the FSTAT software, using a significance level of 0.05, adjusted for multiple comparisons. Software GENALEX 6.5 (Peakall and Smouse 2006) was used for the estimation of the pairwise population PhiPT values. This PhiPT matrix was used in Genalex 6.5 to perform Principal Components Analysis (PCA). PHYLIP 3.6c (Felsenstein 1994) was used for the construction of an unrooted UPGMA dendrogram after 1000 bootstrap resamples, using allele frequencies, as described in Augustinos et al. (2005). STRUCTURE v2.2 software was used to determine the number of possible genetic clusters (Pritchard et al. 2000). We used all four different models implemented in this program, assuming the presence of admixture and without admixture, for independent or correlated allele frequencies, with a burnin period of 50,000 and 500,000 Marcov Chain Monte Carlo (MCMC) repetitions after the initial burn-in. We tested for K = 1 to K = 12 (where K stands for the assumed number of populations). We also used the modification described by Evanno et al. (2005) to predict the true number of clusters in complex situations more accurately. For this purpose, we used the "no admixture" model with correlated frequencies, with a burn-in period of 50,000 and 500,000 MCMC repetitions after the initial burn-in. We ran 10 repetitions of the above model, assuming K = 1 to K = 12. We chose the specific model because, as it is suggested by the authors of the program, it can better resolve "complicated" or "weak" cases of structuring. Software GeneClass 2.0 (Piry et al. 2004) was used to perform population assignment and exclusion test and to calculate the probability of origin for each individual and each population, under a Bayesian model using the Rannala and Mountain (1997) criterion with a 0.05 threshold. Isolation by distance was tested using the "Mantelize it" option of the FSTAT software (Goudet  2001) for plotting genetic distances against geographic distances, and Pearson's r coefficient was measured using Excel. Finally, the BOTTLENECK software was used (Cornuet and Luikart 1996) to detect any recent bottleneck phenomena. This software runs under the assumption that recent bottlenecks lead to a shift from the Lshaped distribution of allele frequencies along with a quicker loss of rare alleles than heterozygosity. However, the low number of microsatellite loci used, especially if not highly polymorphic, can lead in reduced "ability" to detect bottleneck phenomena. In the frame of this study, the Mode shift and the Wilcoxon sign rank tests implemented in this software were used.

Wolbachia genotyping and data analysis
Wolbachia infections from six individual flies belonging to separate genetic clusters (according to microsatellite analysis) and/or of different geographic origin were genotyped by sequencing of the 16S rRNA, wsp, and MLST genes. A 16S rRNA gene fragment (438 base pairs) was amplified with the Wolbachia-specific primers wspecF and wspecR (Werren and Windsor 2000). Gene fragments of wsp and the MLST genes (gatB, coxA, hcpA, fbpA, and ftsZ) were amplified using the respective primers and conditions reported previously (Casiraghi et al. 2005;Baldo et al. 2006;Ros et al. 2009).

Cloning and sequencing
After detecting multiple peaks in preliminary sequencing reactions, PCR products of all genes were cloned, sequenced, and assembled as described in Augustinos et al. (2011b). Three to six colonies were directly subjected to PCR, using the primers T7 and SP6 and double-strand sequenced.

Phylogenetic analysis
All Wolbachia gene sequences generated were aligned using the ClustalW algorithm (Thompson et al. 1994), implemented in MEGA (Tamura et al. 2011). Sequences obtained from GenBank representing all currently known supergroups of Wolbachia were included in the analysis. Phylogenetic reconstruction was performed using the Maximum Likelihood (ML) method, after 500 bootstrap resamples of the original data.

Recombination analysis
Sequences for all genes (MLST, wsp, 16S rRNA) were tested for recombination events using the RDP3 package (Martin et al. 2010). All the available modules (RDP3, BOOTSCAN, GENECONV, MAXCHI, CHIMAERA, SISCAN, 3SEQ, VisRD, LARD, and LDHAT) were implemented under the default options. As the authors of the software state, the "default options" are sufficient to detect recombination events for most of the datasets even without the need of reference nonrecombinant sequences.
Screening for wCer specific strains Ten to twenty flies per population were screened for the presence of wCer1 to wCer5 Wolbachia strains using a wsp gene-based PCR assay as previously reported (Arthofer et al. 2009b).

Nucleotide sequence accession numbers
All microsatellite and Wolbachia gene sequences generated in this study have been deposited into GenBank under accession numbers JX870457-JX870517 and as well as to the Wolbachia MLST database.

Microsatellite marker variability and HWE
The markers used have on average 6.23 n a with the mean n e being 2.16 ( Table 2). The mean observed heterozygosity was lower than the expected (0.39 vs. 0.47). Out of the 234 locus-sample tests performed for HWE, 25 showed deviation according to the G 2 criterion at a significance level of 0.05. Most of these deviations are concentrated in marker Rp11 (six) and can be partially attributed to the presence of null alleles, as this locus exhibits an excess of different classes of homozygotes. Analysis performed with the ML-Null software points to the presence of null alleles for this locus in some of the populations. No evidence for extensive presence of null alleles for the other loci was found.

Population polymorphism and HWE
Populations studied exhibit a varying degree of polymorphism with 1.9 to 2.51 n e per locus (Dafni1-GR and Krasnodar-RU, respectively). According to all measures, Agia Larissa-GR and Krasnodar-RU are the most polymorphic populations while Kernitsa-GR and the two neighboring islands of Chios-GR and Lesvos-GR present the lowest variability. Sporadic HWE deviations (up to three loci) were observed in few of the populations studied (Table 1).

Genetic distances and population differentiation
Pairwise genetic distances (Nei 1972) among populations were determined ranging from 0.0140 to 0.1630 (Table 3).   (Table 4). There are no significant differences among populations from Thessaly, except the one observed between Kamari Pilio-GR/Pertouli-GR. The four populations from Macedonia are also quite homogeneous, with the exception of the difference between Kastoria and Dafni2. The single population from Epirus (Konitsa-GR) shows differentiation with some Thessaly/ Macedonia populations and at the same time, it is not significantly differentiated from the Dossenheim population. The populations from Peloponnesus (Kernitsa-GR), Crete (Chania-GR), Aegean islands (Lesvos-GR and Chios-GR), and Russia (Krasnodar-RU) are all differentiated from all previously reported and among them. The German populations seem to form a homogeneous cluster.

Analysis of the structure of R. cerasi natural populations
To further unravel the genetic structure of R. cerasi natural populations, different programs and assumptions were used. A UPGMA dendrogram was constructed in Phylip 3.6c (Fig. 3). Bayesian clustering analysis was performed using STRUCTURE 2.2 (Fig. 4), along with the modification in Evanno et al. (2005), which is supposed to give a better resolution of the true number of clusters in complex cases of population structuring. PCA was also performed using Genalex 6.5 (Fig. 5). Although these approaches do not produce exactly the same results, the main conclusion is that there is some structuring with the populations from Krasnodar-RUS and Aegean Sea (Chios-GR and Lesvos-GR) to be strongly differentiated from all others. Differentiation is also evident for Kernitsa-GR (Peloponnesus) and the two German populations. Although there is not strong statistical support, two additional loose clusters may be present which include populations from Thessaly-GR and Macedonia/Epirus-GR.
We used software GeneClass v2.0 to assign individuals to respective groups. Eight groups, based on geographic criteria, were formed: Thessaly-GR, Epirus-GR, Macedonia-GR, Peloponnesus-GR, Crete-GR, Aegean Islands-GR, Germany, and Russia. Some populations are well separated; the vast majority of their individuals are correctly assigned to the respective groups (Fig. 2). This is true for populations from Kernitsa-Peloponnesus (29/30 correctly assigned), the Aegean islands (56/60), Russia (28/32), and Stecklenberg (28/32). There was a high number of false assigned individuals for the remaining populations (ranging from 26% to more than 50% per population), especially for Macedonia-GR, Thessaly-GR, and Dossenheim-GER.

Isolation by distance
The hypothesis whether the observed genetic distances can be correlated with geographic distances between sampling locations was tested in pairwise comparison of F ST values and the logarithm of the geographic distances. There was only a weak correlation (r = 0,33, P < 0.1) indicating that other factors may also play a significant role in the development of the observed gene flow restriction between R. cerasi populations.

Detection of bottlenecks
The program BOTTLENECK was used for the detection of recent bottlenecks. Results have been handled with care    Genetic distances according to Nei (1972).   as only few of the populations consist of at least 30 individuals and the loci used in this analysis are <20 and, in addition, they are not highly polymorphic. Given these limitations, no evidence for such phenomena was found with the exception of the Kernitsa population where a shifted mode of allele frequencies from the L-shaped distribution was observed.  PCR products were obtained from these six flies for the 16S rRNA, wsp as well as the five MLST genes (coxA, fbpA, ftsZ, gatB, and hcpA). All PCRs produced the expected amplicon and these products were cloned. Two to six clones were sequenced per product and were found to be Wolbachia specific: sixteen 16S rRNA, twenty-one coxA, twenty fbpA, twenty-two ftsZ, twenty gatB, 15 hcpA, and eighteen wsp clones were sequenced revealing three to five alleles per locus (Table 5). The major findings of the sequencing analysis (Table 5) can be summarized as follows: (a) the majority of the 16S rRNA sequences (14/16) form a cluster in supergroup A ( Figure S2). There are two exceptions: first, one sequence derived from the fly from Chania clusters in another clade, but still with sequences belonging to supergroup A and second, one sequence derived from the fly from Stecklenberg clusters with sequences belonging to supergroup B; (b) the 21 coxA sequences revealed the presence of four alleles (Figure S3). Most of the sequences (16) correspond to allele coxA-84, which has been reported in wCer1 and wCer4 strains, and were found in Greek populations (Agia, Chios, Kernitsa) as well as in Krasnodar and in Stecklenberg. The allele coxA-1, which is present in wCer2, was detected in Stecklenberg too. Interestingly, two new alleles were detected in the individual from Chania-GR. These new alleles are named coxA-188 and coxA-189. CoxA-188 differs only in 1 bp from allele coxA-130 and coxA-189 differs in 1 bp from allele coxA-58. However, both cluster near the allele coxA-84, which is common in wCer1 and wCer4 strains. None of our coxA sequences group with wCer5, a Wolbachia strain belonging to supergroup B; (c) the 20 fbpA sequences suggested the presence of four alleles (fbpA-1, -4, -79, and -160) which correspond to alleles characterized so far to four different wCer strains. Most of the clones (15) correspond to allele fbpA-160, present in strain wCer1, and originate from Greek populations (Agia, Chania, Chios, Kernitsa) as well as from Stecklenberg and Krasnodar. Two clones correspond to fbpA-1 of strain wCer2 and were detected in Stecklenberg, two clones correspond to fbpA-79, which is assigned to wCer4, and were present in Chania and in Kernitsa while one clone (from Chios-GR) corresponds to fbpA-4, which is found in wCer5, the only strain belonging to supergroup B; (d) the 22 ftsZ sequences revealed the presence of four alleles (ftsZ-3, -22, -70, and -79) described for wCer strains. All sequences cluster with wCer strains belonging to supergroup A, with the exception of the sole sequence from Stecklenberg, which clusters with wCer5 of supergroup B (ftsZ-22). Thirteen of the clones correspond to the allele ftsZ-79 assigned to wCer1 and originate from Greek populations (Agia, Chios, Kernitsa) and from Krasnodar. Six clones correspond to ftsZ-70 of wCer4 (originating from Chania, Kernitsa, and Krasnodar) while two clones correspond to ftsZ-3 assigned to wCer2; (e) the 20 gatB sequences suggested the presence of three wCer alleles belonging to supergroup A (gatB-1, -8, and -53). Eleven of them correspond to allele gatB-8 of strain wCer1 originating from Greek populations (Agia, Chios, Kernitsa) as well as from Stecklenberg and Krasnodar. Seven clones are identical to gatB-53 (wCer4) and were detected in Chania and Krasnodar while two clones are identical to the allele gatB-1, which is assigned to strain wCer2, and were found in Stecklenberg; (f) the 15 hcpA  sequences revealed the presence of five wCer alleles (hcpA-40, -85, -103, -212, and -213). All of them cluster with wCer sequences of supergroup A, with the exception of two of the three clones from Chios, which are identical to allele hcpA-40, which is assigned to wCer5 ( Figure S4). Eight clones correspond to allele hcpA-103 of strain wCer1 originating from Greek populations (Agia, Chios, Kernitsa) and from Krasnodar. Only two sequences correspond to hcpA-85 of wCer4, both detected in Kernitsa. Interestingly, in the individual from Chania, two new alleles were present (-212, -213). HcpA-213 is very similar to hcpA-103 of strain wCer1 (different in the first 3 bp), while hcpA-212 is very similar to hcpA-85 of strain wCer4 (different in the first 3 bp also); and (g) the 18 wsp sequences suggested the presence of three known wCer alleles (wsp-113, -335, and -581). All of them cluster in supergroup A, with the exception of one of the five clones from Chios, which is identical to allele wsp-581, which is found in strain wCer5. Fifteen of the sequences were identical to wsp-335 of wCer1 originating from Greek populations Agia, Kernitsa, Chania, Chios) and from Krasnodar, while two clones from Chania were identical to wsp-113 of strain wCer4. Taken together these data strongly suggest the presence of multiple Wolbachia infections in all R. cerasi populations studied. The presence of multiple Wolbachia infections was also confirmed by a wsp gene-based diagnostic PCR assay previously developed (Arthofer et al. 2009b). The primers used in this assay selectively amplify different alleles of the wsp gene that have been attributed to the different wCer strains. Ten to twenty individuals were screened per population (Table 6). In accordance with previous studies, all individuals were found infected with Wolbachia and all of them harbored at least the wCer1 strain (Riegler and Stauffer 2002; Arthofer et al. 2009bArthofer et al. , 2011. The majority of the Greek individuals were found doubleinfected with wCer1 and wCer4 strains while a few individuals were found triple-infected with different combinations of wCer1, wCer 2, wCer4, and wCer5 strains or infected with four Wolbachia strains (wCer1, wCer2, wCer4, and wCer5). It has to be noted, however, that wCer1 single-infected individuals were detected in several populations while one of them, Lesvos-GR, was found to be 100% single-infected with this strain. The majority of Krasnodar-RUS individuals were double-infected with wCer1 and wCer4 strains while triple-infected (wCer1, wCer2 and wCer4) flies were also present. The German R. cerasi populations were almost fixed for the presence of four different Wolbachia strains: wCer1, wCer2, wCer4, and wCer5. None of the individuals tested was found infected with the wCer3 strain.

Discussion
Rhagoletis cerasi present some unique properties when compared with the medfly and the olive fly. Medfly is a multivoltine and polyphagous species (Papadopoulos et al. 2001), the olive fly is a multivoltine and monophagous species (Tzanakakis 2003), while the cherry fruit fly is a univoltine and stenophagous species (Moraiti et al. 2012a,b). These three species also differ in their Wolbachia infection status: natural populations of olive fly and medfly are Wolbachia free (Bourtzis et al. 1994;Zabalou et al. 2004;Apostolaki et al. 2011; but see Rocha et al. 2005 for an infected Brazilian population) while Wolbachia has been detected in all natural populations of the European cherry fruit fly tested so far. Our results indicate the presence of structuring in natural populations of R. cerasi, which, especially for Greece, is not in accordance with data for the olive fly (Augustinos et al. 2005;Zygouridis et al. 2009), and the medfly (K. Economou et al. unpubl. data), which point to rather homogeneous groups for both. Is then the genetic structuring observed in natural cherry fruit fly populations due to ecological factors or rather due to the presence of different/multiple Wolbachia strains?

Variability of microsatellite markers and HWE
We used 13 microsatellite markers; 10 cross-amplified from other Rhagoletis species (see Augustinos et al. 2011a), one from the olive fly (Augustinos et al. 2008), and two of the set that has been developed de novo for R. cerasi (Arthofer et al. 2009a). Some of the cross-amplified markers identified low polymorphism compared with species in which they were developed. This is not surprising as it is known that cross-amplification of microsatellite markers includes three major risks: (a) amplification of nonhomologous regions, (b) amplification of homologous regions with fewer or disrupted repeats, leading therefore to lower mutation rates, and (c) poor amplification in some or most of the individuals, due to changes in primer binding sites. Sequencing analysis of crossamplified products verified that we amplified homologous regions (data not shown), with the exception of RcMic (Boms3b) from the olive fly and RcRp1 from R. pomonella. In most of the cases, microsatellites with fewer and/or interrupted repeats or with new motifs (data not shown) were present, explaining the reduced polymorphism. There was also evidence of size variation not attributed to motif changes, even within R. cerasi (data not shown), implying the participation of other events, such as unequal crossing over and/or indels affecting the microsatellite marker size. Although 13 de novo developed microsatellite markers were available for R. cerasi (Arthofer et al. 2009a), we managed to use only two of them.
In preliminary trials, we found that scoring of the remaining 11 markers was very difficult and unreliable, also due to secondary bands and poor quality of main bands (faint bands or not sharp enough to be scored). This is probably because these markers harbor multiple Table 6. PCR screening of Rhagoletis cerasi populations with the primer pairs developed for the detection of Wolbachia wCer 1-5 strains (Arthofer et al. 2009b).
motifs, one of which being one base pair repeat motif, making their analysis problematic, especially in polyacrylamide gels. Deviations from HWE were found in few of the populations for some of the cross-amplified loci but not for the two species-specific microsatellite markers (Table 2). This could partially be attributed to the small sample size (in some cases) or the presence of null alleles, as mentioned above. Evidence for a relative extended presence of null alleles in some of the samples was found only for Rcmic(Rp11).

Rhagoletis cerasi natural populations' polymorphism
All populations exhibited a comparable degree of polymorphism. The most polymorphic population is the one from Krasnodar-RUS. This can be explained considering that R. cerasi moved to Europe from western Asia through the Caucasus. This population should therefore be closest to the original expansion area of the species. Populations near the expansion center of a species are expected to be more polymorphic, if no drastic bottlenecks have occurred and populations have maintained a relative large size through time. On the other hand, populations from the Eastern Aegean islands (Chios and Lesvos) along with Kernitsa from Peloponnesus present the lowest polymorphism. This can be attributed to an island effect, as isolated populations are expected to be less polymorphic, due to restriction of gene flow and probably small effective population size. As far as the sample from Kernitsa (Peloponnesus) is concerned, Peloponnesus is not a true island and communicates with the rest of mainland Greece; however, cherry orchards in the area are few and dispersed. Another factor to be taken into account is that the sample from Kernitsa consists of individuals that have undergone prolonged dormancy, while all other samples consist of flies that emerged after annual dormancy. Although the prolonged dormancy in R. cerasi pupae seems to be a plastic response to environmental cues (C. A. Moraiti and N. T. Papadopoulos unpubl. data), its possible impact on population structuring including bottleneck phenomena deserves further investigation.

Structuring of populations
Despite the low overall polymorphism observed, the markers used here identified well and poorly separated clusters (Figs. 2,3,4,5). The well separated groups included populations from Aegean islands-GR and Krasnodar-RUS followed by Kernitsa-GR and German samples. The two loose groups were those from Macedonia/Epirus-GR and Thessaly-GR. Despite its geographical isolation, the single population from Chania was not well differentiated probably due to trade and airport connections between Crete and many mainland regions in Greece and other European countries including Germany.
The genetic difference of the samples collected from islands or "island"-like sites (Lesvos, Chios, Kernitsa) is accompanied by a slightly lower degree of polymorphism (Table 1). This is in accordance with the hypothesis of partially isolated populations. No evidence of recent bottlenecks was found (with the exception of Kernitsa population discussed above) suggesting that we are probably dealing with long-standing populations that manage to keep a relative high population density through time.
Another finding, although with low support (Figs. 3,5), is the topology of the Macedonian-GR group of samples (especially the sample from Kastoria): as said before, they form a loose group, placed near the Thessaly-GR cluster, having, however, a tendency to link with the German cluster. If we take into account that Thessaly is in Central Greece and Macedonia-GR is mainly in Northern Greece and is the region that neighbors with Albania, Former Yugoslav Republic of Macedonia and Bulgaria, we can speculate in favor of the existence of gene flow from Northern Europe in this area. The lack of samples from this area limits our analysis.
Samples from the two Eastern Aegean islands are strongly differentiated from all other Greek and German samples. They seem to form a "loose" cluster with the sample from Russia, pointing toward gene flow from Asia Minor, most probably through Turkey, a major cherry producer. Analysis of samples from the coastal area of Turkey all the way to the coast of Black Sea would definitely shed light on gene flow from Asian to European R. cerasi populations. It is worth noting that these island populations do not form one panmictic population, despite the fact that the distance between them is only 90 km.
Despite differences in genetic distance within the same geographic scale between R. cerasi, C. capitata, and B. oleae, the degree of polymorphism was comparable or even smaller in R. cerasi (though the same markers couldn't be used). Therefore, these differences are not marker biased. They are also not population biased, as relatively large samples were used in all studies (20-30 individuals per population). It seems that the fragmented landscape of mainland Greece, together with a patchy distribution of cherry orchards, the low dispersion ability of the cherry fly and local adaptation from close synchronization of adult flight with the sweet cherry ripening period may restrict gene flow to relatively short distances (Phillips and Dirks 1933;Jones and Wallace 1955;Boller and Prokopy 1976;Fletcher 1989;Kneifl et al. 1997 diapause intensity between coastal and highland R. cerasi populations (Papanastasiou et al. 2011) and extended difference in adult life history traits between populations (Moraiti et al. 2012b). Low reproduction rates of R. cerasi because of univoltinism and stenophagy may also partially explain differences in genetic structure against both C. capitata and B. oleae. On the other hand, infection of all natural populations of R. cerasi with different strains of Wolbachia may be yet another factor contributing toward structuring and genetic isolation.
Wolbachia is assumed to invade natural populations rapidly and is implicated in speciation phenomena and restriction of gene flow between natural populations (Shoemaker et al. 1999;Bordenstein et al. 2001;Jaenike et al. 2006;Koukou et al. 2006;Miller et al. 2010;Branca et al. 2011). This can be either through CI and/or rapid fixation of specific genotypes during invasion. Natural European R. cerasi populations are infected with multiple Wolbachia strains and show different infection patterns. For example, in Poland, Italy, and Austria, there are populations infected with all five wCer strains, in the Czech Republic and Portugal there are populations infected with four of the five strains (except wCer2), while in Switzerland there was an infection with all wCer strains except wCer3 (Arthofer et al. 2009b. Moreover, it has been shown that the presence of different Wolbachia strains in R. cerasi could be the causal factor (Riegler and Stauffer 2002) of the incompatibility observed among R. cerasi natural populations in previous studies Matolin 1976).
All individuals in this study (10-20 per population) were infected with Wolbachia. The analysis of six individuals from different samples (genetically or geographically distant) through cloning and sequencing analysis, as well as the PCR screening of more individuals with the diagnostic PCR primers, pointed toward (a) the presence of complex patterns of infections with the four of the five known wCer strains (1, 2, 4, and 5) and (b) the possible existence of new, still uncharacterized strains (see Chania population). However, not all genes resulted into the same infection pattern (Tables 5, 6). This can be attributed to the preferential recovery of specific sequences after PCR and/or cloning. Overall, these two approaches (PCR screening and cloning plus sequencing) produced similar results and supported the presence of multiple infections. However, both PCR and cloning can lead in bias in favor of specific sequences and to the misinterpretation of the "actual" infection status. On the other hand, PCR screening allows the analysis at population level but currently is limited by the development of strain-specific primer pairs only for the wsp gene and by the inability to detect "hidden" strains. Arthofer et al. (2011) highlighted the difficulties in Wolbachia strain characterization through MLST in multiple-infected hosts. They also provided a tool to analyze complex infection types (Allelic Intersection Analysis) and evaluated this tool for wCer Wolbachia strains of R. cerasi by assigning specific MLST alleles to four of the five wCer strains (Table 5). It is evident that the MLST-and wspbased typing approaches do not have the power to resolve complex infection patterns such as those observed in R. cerasi. There is probably the need for new tools and approaches to address such questions. Current use and evaluation of different types of molecular markers, including Variable Number Tandem Repeat markers (VNTRs) (Riegler et al. 2005(Riegler et al. , 2012Schneider et al. 2013) can probably contribute toward this direction.
Multiple and/or low titer strains have been shown to contribute to genetic differentiation of their hosts (Arthofer et al. 2009bSchneider et al. 2013). Our data confirm the presence of multiple infections, different types of infections in different samples and polymorphic infection status within the samples. They also raise the possibility of the presence of still unknown strains. Given also the ability of at least two strains (wCer2 and wCer4) to induce CI Riegler et al. 2004;Zabalou et al. 2004Zabalou et al. , 2009Apostolaki et al. 2011), Wolbachia is not a factor that can be overlooked when addressing the question of population structuring in R. cerasi.
Recombination has perhaps been the most significant contributing factor to the evolution of Wolbachia and may, at least partly, explain the diversity of the wCer strains present in European populations of R. cerasi (Jiggins 2002;Baldo and Werren 2007;Arthofer et al. 2009b;Klasson et al. 2009). Low titer infections are an additional, hidden source of Wolbachia diversity, as documented in different hosts (Arthofer et al. 2009b;Miller et al. 2010;Augustinos et al. 2011b). There is also a number of studies showing that the transfer of Wolbachia genes (in some cases even the entire genome) into host chromosomes might be a rather common phenomenon (Kondo et al. 2002;Hotopp et al. 2007;Nikoh et al. 2008;Aikawa et al. 2009;Klasson et al. 2009;Nikoh and Nakabachi 2009;Woolfit et al. 2009;Doudoumis et al. 2012). If this is the case, then nucleusencoded Wolbachia genes may also contribute to the observed genetic diversity. However, our data do not provide evidence for either (intragene) recombination or horizontal gene transfer events.
Do all wCer Wolbachia strains induce CI? The fact that R. cerasi is a univoltine species and an efficient and robust lab rearing system is lacking makes any experimentation to address questions regarding the CI properties of wCer strains very difficult. Interpretation of historical data in European populations of R. cerasi as well as transinfection experiments with medfly and olivefly suggest that wCer2 and wCer4 induce strong CI into their hosts and are bidirectionally incompatible Riegler and Stauffer 2002;Riegler et al. 2004;Zabalou et al. 2004Zabalou et al. , 2009Apostolaki et al. 2011). On the other hand, Drosophila simulans was not shown to be a suitable host for the establishment and/or the expression of CI for either wCer1 or wCer2 strains ). Thus, transinfection experiments in suitable hosts like medfly can probably unravel the CI properties of wCer3, wCer5, and the two potentially new wCer strains.

Concluding Remarks
Our results reveal some structuring of natural European cherry fruit fly populations, which can probably be generalized to its whole geographic range. A number of factors may have contributed to this structuring, such as habitat fragmentation as well as Wolbachia infection. Microsatellite analysis of more populations obtained from geographically distant areas, along with further genetic characterization of Wolbachia strains and their CI properties are required to understand the population dynamics of this economically important species and its associated Wolbachia strains more completely. This knowledge will be needed before any novel population control method is applied on a large scale for this important pest species.

Acknowledgments
We thank Antigone Zacharopoulou (Department of Biology, University of Patras, Greece) for her generous support and encouragement throughout this study. Thanks are also extended to Ludmilla Vassilieva, Kirsten Koeppler, Byron Katsoyannos, Ilias Kounatidis, Alex Diamantidis, Stella Papanastasiou, and Charalampos Ioannou for either submitting R. cerasi samples or participating in collecting the different populations. We also thank Wolfgang Miller and Daniela Schneider for providing control DNA samples used in the wCer strain-specific PCR screen presented in Table 6. This research was supported in part by the European Community's Seventh Framework Programme CSA- SA_REGPROT-2007-1 (under grant agreement no 203590 (K. B.). This study was partially supported by IKYDA and the EPEAEK II -Pythagoras (cofounded by the European Social and National Resources and the Greek ministry of Education) awarded to N. T. P. P. M., N. T. P., and K. B. benefitted from travel grants in the frame of the EU COST Action FA0701 "Arthropod Symbiosis: from fundamental studies to pest and disease management". The authors are grateful to the FAO/IAEA Coordinated Research Program "Use of symbiotic bacteria to reduce mass-rearing costs and increase mating success in selected fruit pests in support of SIT application" for the overall support of this study. The authors also thank Stefan Oehler for his comments on an earlier version of the manuscript. The funding agencies had no role in study design, data collection or analysis, decision to publish, or preparation of the manuscript. Figure S1. Posterior probability of possible number of populations in STRUCTURE. Note the potential presence of six genetic clusters indicated by the arrow in diagram (D). Figure S2. Maximum Likelihood dendrogram based on the 16S rRNA gene sequences. Gray shade: sequences generated in this study. Wolbachia 16S rRNA gene sequences representing current known Supergroups were retrieved from Genebank. For these sequences, the host scientific name, the Genebank Accession Number and the respective supergroup are shown. Analyses were conducted in Mega 5 (Tamura et al. 2011). The evolutionary history was inferred using the Maximum Likelihood method based on the Tamura-Nei model. The tree with the highest log likelihood is shown and percentages of trees in which the associated taxa clustered together are indicated next to the branches. Figure S3. Maximum Likelihood dendrogram based on the coxA gene sequences. Gray shade: sequences generated in this study. Darker gray shade and bold: sequences of known wCer strains. Different Wolbachia coxA alleles were retrieved from the Wolbachia MLST database. Alleles representing all available supergroups were chosen. New alleles coxA-188 and coxA-189 are marked with asterisks (*). Analyses were conducted in Mega 5 (Tamura et al. 2011). The evolutionary history was inferred using the Maximum Likelihood method based on the Tamura-Nei model. The tree with the highest log likelihood is shown and percentages of trees in which the associated taxa clustered together are indicated next to the branches. Figure S4. Maximum Likelihood dendrogram based on hcpA gene sequences. Gray shade: sequences generated in this study. Darker gray shade and bold: sequences of known wCer strains. Different Wolbachia hcpA alleles were retrieved from the Wolbachia MLST database. Alleles representing all available supergroups were chosen. New alleles hpA-212 and hcpA-213 are marked with asterisks (*). Analyses were conducted in Mega 5 (Tamura et al. 2011). The evolutionary history was inferred using the Maximum Likelihood method based on the Tamura-Nei model. The tree with the highest log likelihood is shown and percentages of trees in which the associated taxa clustered together are shown next to the branches.