Phylogeography of Trigonaspis synaspis ( Hymenoptera : Cynipidae ) from Anatolia based on mitochondrial and nuclear DNA sequences

Sequence data for the nuclear ITS2 region and part of the mitochondrial cytochrome b gene were used to reveal intraspecific phylogeography and population genetic structure of the oak gall wasp Trigonaspis synaspis from Anatolia. One hundred and sixty six individuals representing 21 populations yielded 50 unique cyt b haplotypes and 15 ITS2 alleles. Genetic diversity estimates indicated the existence of high molecular variation, with average gene diversities of 0.85541 and 0.3119 for cyt b and ITS2 respectively. Likewise, nucleotide diversity was 0.01346 for cyt b and 0.0015 for ITS2. Pairwise divergence ranged from 0.2 to 3.1% between cyt b haplotypes, and from 0.2 to 2.8% between ITS2 alleles. Phylogenetic analyses of the cyt b haplotypes conducted using three different approaches produced mostly similar topologies, and major clades generally grouped haplotypes with similar geographic origins. Contrary to the mtDNA results, ITS2 analyses produced largely similar but polytomous topologies lacking apparent geographic structure. Application of the widely used insect mitochondrial molecular clock linked the origin of major phylogeographic structure to the effects of the Quaternary climatic oscillations. The outgroup haplotypes of T. megaptera are inferred to have diverged from T. synaspis around the late Pliocene. Intraspecific node ages in T. synaspis link major clade structure to the Pleistocene glacial cycles. Geographical formations in Anatolia were also probably an important factor in shaping the phylogeographic structure of T. synaspis.


INTRODUCTION
Geographic distribution of genetic variation and lineages within a species range is shaped by current and historical factors (Avise, 2000).Areas with complex geological history, long-term temperature fluctuations and variation in climate even over short geographic distances are particularly likely to show significant phylogeographic patterning in associated organisms.Turkey is such an area, with a complex geological and paleogeomorphological history (Bozkurt, 2001).
Anatolia is formed by part of the Aegeid or the Irano-Anatolian plate, which is thought to have been formed by tectonic detachment from North Africa.Major formation of the Irano-Anatolian plate began in the Late Oligocene with a collision between the Indian and Asian Continental masses, followed by Africa and Eurasia.The last collision resulted in the formation of the major Anatolian mountain ranges (Rögl, 1999) which are thought to have shaped phylogeographic structure and species diversity in Turkey (Ekim & Güner, 1986;Çıplak, 2008).Such a major mountain line, the Anatolian Diagonal (Davis, 1971), has been shown to structure the Anatolian species and genetic diversity (Ansell et al., 2002;Mutun, 2010).A series of faunal exchanges between Anatolia, Asia, Africa and Europe through the Neogene and Quaternary (Koufos et al., 2005) have further shaped the contemporary fauna and the distribution of the genetic lineages in Anatolia.However, the greatest structuring of faunal composition and genetic variation of taxa can be linked to the Pleistocene (Rögl, 1999;Çıplak, 2004).It has been suggested that during the last glacial and interglacial cycles Anatolia served as a non-homogenous refugium (Çıplak, 2008).Taken together, these factors are thought to have influenced different taxa, resulting in a rich Anatolian fauna and flora that has acted as the genetic source for European populations (Dubey et al., 2007;Şekercioğlu et al., 2011).Further, recent studies have indicated that multiple oak gall wasp lineages diverged prior to the arrival of modern oaks in the western Palearctic, spreading westwards into Europe from an eastern, Anatolian origin (Challis et al., 2007;Stone et al., 2012).Indeed, recent studies have found genetic variation in Turkish and more easterly populations to be greatest among surveyed Western Palearctic populations, supporting an Anatolian origin for pre-Pleistocene colonization by European gall wasps (Rokas et al., 2003).
Trigonaspis synaspis (Hartig, 1841) is widely-distributed oak gall wasp species found from the Iberian Peninsula through Turkey to Iran (Melika, 2006).Here, we use Turkey and Anatolia interchangeably, because we found no samples from Thrace (Turkey in Europe), and all specimens are from the Asian part (Anatolia).The extent to which the genetic structure of Anatolian species has been shaped by the Pleistocene climatic cycles and Turkey's topography remains an issue to be resolved in different animal taxa, including oak gall wasps.Therefore, in this study we used sequence data for part of the mitochondrial cytochrome b gene and the complete nuclear ITS2 region for multiple gall wasp populations to: i) reveal intraspecific genetic variation across the distribution range of T. synaspis from Anatolia, and to determine whether eastern or central-western Anatolia shows higher variation, and ii) assess the signals of Pleistocene climatic changes in the pair (bp) fragment of the mitochondrial cyt b gene was amplified using the CB1 / CB2 primer pairs (Jermin & Crozier, 1994).Amplification protocols were performed as previously described in Stone et al. (2007).For amplifying ITS2 region of nuclear DNA ITS2-F-5´-ATGAACATCGACATTTCGAACGCATAT-3´ and ITS2-R-5´-TTCTTTTCCTCCGCTTAGTAATATGCTTAA-3´ primers were used (Ji et al., 2003).PCR reactions were carried out in a PTC-200 DNA Engine (MJ Research) thermal cycler.All reactions were carried out in a total volume of 40 µl including 10 µl 5X buffer, 2.5 µl MgCl 2 (25 mM), 3 µl dNTPs (10 mM), 1 µl of each primers (20 pmol), 0.5 µl template DNA and 1 U Taq polymerase (Vivantis); ddH 2 O was added to each reaction to a total volume of 40 µl.For the ITS2 region temperature rates, times and cycles were set as follows: 5 min at 95°C for first denaturation; 35 cycles of 40 s at 95°C for denaturation, 40 s at 52.2°C for annealing, 45 s at 72°C for elongation and 5 min at 72°C for last elongation.Successful PCR amplification was confirmed by running samples in a 1% agarose gel.Amplicons were sent to a company for clean-up and sequenced using Perkin-Elmer BigDye Terminator chemistry on an ABI automated sequencer (Applied Biosystems).To minimize PCR artifacts, ambiguities and base-calling errors both strands were sequenced.The resulting sequences are deposited in GenBank (Accession numbers KJ912832-KJ912880 for the cyt b gene and KJ912881-KJ912894 for the ITS2 region of T. synaspis, and KJ716686 and KJ716687 for the cyt b and ITS2 regions of T. megaptera, respectively).

Haplotype determination and population genetic estimates
Sequences were checked by eye to avoid misreading in both strands and then further analyzed using BioEdit 7.0.9software (Hall, 1999).Sequence alignment was performed separately for each sequenced fragment with default parameters and assembled sequences were collapsed into haplotypes.The open reading frame of each mitochondrial cyt b haplotype was confirmed with no internal stop codons, showing that sequenced haplotypes contemporary distribution of genetic variation among Anatolian populations of the species.

Sampling and laboratory protocols
Specimens of T. megaptera collected from Giresun and Edirne localities were used as outgroups in all phylogenetic analyses.Asexual generation galls of T. synaspis (Hartig, 1843) were collected from 21 distinct localities across Anatolia between 2011-2013 (Fig. 1 and Table 1).After rearing adults from the galls a total of 166 specimens were used for DNA extraction using the DNeasy Tissue Kit (QIAGEN cat.69504) following the manufacturer's protocol for insect DNA isolation.The use of outgroup specimens from the same locality with T. synaspis may be problematic due to the possibility of hybridization events and incomplete sorting of genetic variation between closely related gall wasp species within shared glacial refugia (Nicholls et al., 2012).However, such potential problems can be avoided with sampling and through comparison of phylogenetic patterns of unlinked loci, such as mitochondrial and nuclear genes.Thus, a 433 base Fig. 1.Sampling localities of T. synaspis.Names, abbreviations and coordinates of populations are provided in Table 1.Hewitt, 1996;Rokas et al., 2003).Population genetic statistics were estimated using both Arlequin 2001 (Schneider et al., 2000) and DnaSP 4.0 (Rozas et al., 2003).The number of polymorphic sites (S), nucleotide (π) and haplotype (h) diversity (Nei, 1987), the number of pairwise nucleotide differences (k) were also calculated (Tajima, 1983).Population demographic histories were characterized using the pairwise mismatch distribution as implemented in DnaSP 4.0 (Rozas et al., 2003).Observed pairwise number of differences between the eastern part of Anatolia (Bitlis, Muş, Bingöl, Erzincan, Elazığ, Tunceli and Adıyaman populations) and those from the Central-Western part of Anatolia (all other remaining populations) was used to compare against both constant population size and sudden population expansion models.Pairwise differences between haplotypes from expanding (but geographically unstructured) populations often generate a unimodal distribution, while samples drawn from populations at demographic equilibrium produce a multimodal pattern of numerous sharp peaks (Rogers & Harpending, 1992).The raggedness index (r) (Harpending, 1994) was calculated to check the deviations from neutrality and to quantify the smoothness of pairwise mismatch distributions.
To test for population expansion we applied two statistical tests often used to analyze demographic events.Tajima's D (Tajima, 1989) uses the frequency information of segregating sites, while Fu's Fs uses the distribution of haplotypes (Fu, 1997).Both are sensitive to factors such as population bottlenecks or expansion processes associated with a non-neutral process showing a shift in allele frequency compared to population expansion under neutral evolution.Statistical significance was assessed using 1000 simulations.
We analyzed the relationship between population pairwise genetic similarity and geographic distances using a spatial autocorrelation multivariate approach in GenAIEx 6.5 software (Peakal & Smouse, 2012).The 95% confidence intervals for two distance classes (50 km and 100 km) were estimated by random permutation of all individuals, and to establish the upper and lower confıdence intervals we recalculated correlation coefficient value (r) 1000 times.The probability was evaluated that a certain level of autocorrelation observed in each distance class may occur under the null hypothesis of no spatial structuring.Bootstrapping was used for the 95% confidence intervals around the r value, and the correlation coefficient was plotted as a function of distance to obtain a spatial genetic correlogram.Similarity in haplotypes sampled in neighboring sites, or positive autocorrelation signals, may be due to spatial patterns shaped by natural selection or neutral processes and are expected when gene flow is limited to short distances (Epperson, 1990).A monotonic decrease in the autocorrelation coefficient from significant positive values at short distances to significant negative values at large distances across multiple loci can be indicative of directional dispersal events (Stone & Sunnucks, 1993;Bertorelle & Barbujani, 1995).

Phylogenetic and phylogeographic analyses
The mean transition to transversion ratio and the nucleotide frequencies were calculated using PAUP*4.0(Swofford, 2002).The T. synaspis data were analyzed using maximum likelihood (ML) and maximum parsimony (MP) in PAUP*4.0, and using Bayesian inference in MRBAYES 3.0 (Huelsenbeck & Ronquist, 2001).MP analysis was performed using heuristic search options, the TBR branch-swapping algorithm and 100 replicates of random addition of taxa.1000 bootstrap replicates were used to assess branch support.The best-fit model of DNA substitution was chosen according to the Akaike information criterion (AIC) (Akaike, 1974;Posada & Buckley, 2004) as implemented in JModeltest (Posada, 2008).The HKY+I+G and GTR+I+G were used as the substitution models suggested for the mitochondrial cyt b gene and the nuclear ITS2 region, respectively.Bayesian phylogenetic analysis involved four simulations of Markov Chains, 10 7 generations with trees being sampled every 100 generations and a burn-in of 25%.The software tool TRACER version 1.5 (Rambaut & Drummond, 2003) was used to observe the parameters and to determine the number of generations needed to reach stationarity (burn-in).After discarding burn-in trees and evaluating convergence, the remaining samples were used to generate 50% majority rule consensus trees and for calculating mean, variance, 95% credibility intervals and posterior probabilities.
Parsimony and maximum-likelihood may be inappropriate for a population level study due to violation of their assumptions.Haplotype networks may be more useful to resolve population relationships within a species as gene flow may generate reticulate rather than hierarchical or treelike relationships between populations.Network analysis also provides information about the age of the haplotypes where interiorly located haplotypes with more mutational connections are assumed to be ancestral and older than tip haplotypes (Posada & Crandall, 2001).Reconstruction methods enforcing bifurcations may also not suitably show all conflicts between alternative topologies.We constructed a minimum spanning haplotype network using Arlequin (Excoffier et al., 1992), based on pairwise differences between haplotypes.The program HapStar Version 0.5 (C) (Teacher & Griffiths, 2011) was used to visualize the output.
To test if genetic differentiation is geographically and hierarchically structured, populations were clustered into different groupings and analyzed using Analysis of Molecular Variance (AMOVA) using Arlequin 2.0 (Excoffier et al., 1992).This method estimates the proportion of genetic variation at different hierarchical levels through using information drawn from the geographical distribution of haplotypes and the pairwise distances between haplotypes.In AMOVA analyses we tested several alternative grouping schemes including (i) eastern (BTL,MUS, BNG, ERC, ELZ, TuN, and ADY populations) versus centralwestern groups of populations (GRS,TOK, KAY, SVS, YOZ, COR, NEV, CKR, KIR, KAS, KAR, KON, and BOL), (ii) eastern (BTL,MUS, BNG, ERC, ELZ, TUN, and ADY), versus western (KAR, KON, AFY, BOL, CNK, KAS) versus central groups (GRS, SIV, KAY, NEV, YOZ, COR, TOK), and (iii) several small possible groups of populations based on their geographic proximity.AMOVA analysis was performed at three levels: among geographical groups of populations, among populations within each geographical group, and within populations, using 1000 permutations to estimate statistical significance.

Estimation of divergence times
Divergence times and their confidence intervals were estimated using a Bayesian Markov Chain Monte Carlo (MCMC) method implemented under BEAST version 1.5.2 (Drummond & Rambaut, 2007).The data set was calibrated by using a widely-applied substitution rate for insect mitochondrial DNA (2.3% sequence divergence per million years) (Brower, 1994).Dates of the divergence were inferred using the HKY+G+I model, a relaxed molecular clock, following the uncorrelated relaxed lognormal clock as implemented in BEAST.Both MRCAs (most recent common ancestors) and MACAs (most ancient common ancestors) were calculated and operators were optimized (Hayward & Stone, 2006) using Beauti.BEAST was run for 20 million generations sampling every 1000, with parameter convergence and effective sample sizes (ESS) checked using TRACER.The maximum clade credibility tree was built using TreeAnnotator, dicardin-ginitial 10% samples as burn-in, and FigTree version 1.3.1 (Rambaut, 2009) was used to visualize the results.

Mitochondrial DNA genetic diversity estimates of T. synaspis
The 433 bp cyt b gene of T. synaspis sampled from 166 individuals representing 21 Turkish localities revealed 50 haplotypes.No indels, stop codons or non-sense mutations were detected in the sequences.There were 391 constant and 42 polymorphic characters, of which 17 were parsimony uninformative and 25 were parsimony informative.Thirty three substitutions occurred at the third codon position (78%), one substitution was observed at a second codon position (2%), and eight were at first codon positions (19%).The average base composition (outgroup excluded) was A: 34.5, T: 42.2, C: 13.0 G: 10.3, indicating a strong anti-G bias which is typical of mitochondrial genes and not numt (Zhang & Hewitt, 1996).
There was only one transversion (at site 55) and 40 transitions (in all other substitution sites) which supports insect mitochondrial protein coding genes having higher transition than transversion rates (Brown et al., 1979).Transition/transversion bias calculated according to the equation provided in Tamura et al. ( 2004) was R = 5.19.Average number of nucleotide differences was k = 5.830.Tamura & Nei distances between haplotypes ranged from 0.2% to 3.1% (1 bp to 13 bp difference, respectively).A single basepair nucleotide difference was observed between 32 cases, with 13 bp as the highest difference detected between H37 (from Konya population) and H2-H3 (both from the Adıyaman population) (S1).Of the translated 144 amino acids, there were 6 replacement substitutions without any non-sense amino acids.
The results of the spatial genetic autocorrelation analysis across the study area for distance classes of 50 and 100 km with 95% confidence intervals are shown in Fig. 2. The correlogram for the 50 km-distance class size showed positive and significant values of r at 50 km (r = 0. 489, p < 0.01) and 100 km (r = 0.397, p < 0.01), with an x-axis intercept at 278.3 km (Fig. 2a).The correlogram for the 100 km distance class size detected positive and significant values of r at 100 km (r = 0.41, p < 0.01), but significantly negative values for distances beyond the x-axis intercept 308.3 km (Fig. 2b).These results imply that populations of T. synaspis separated by distances of less than 100 km are more genetically similar to each other are higher than expected by chance, while genetic differentiation between populations separated by more than ca.300 km is higher than expected by chance.

Nuclear DNA ITS2 Variation
After trimming both 5.8S and 28S rDNA nucleotides from both ends the amplified ITS2 fragment of T. synaspis is 433 bp in size without length variation.Fifteen alleles were detected with 414 constant and 19 variable characters, of which 5 were parsimony informative.There was only one multiple hit at site 309.Nucleotide uses were 29.64% for A, 34.84% for T, 20.08% for C and 15.44% for G.The transition and transversion ratios were k 1 = 2.917 (purines) and k 2 = 2.555 (pyrimidines).There were eight transversions (at sites 176, 219, 265, 284, 296, 305, 355 and 429) and 10 transitions (in other substitution sites) with the most frequently observed transition being C/T (18.8) followed by G/A (18.26).Transition/transversion bias was detected as R = 1.244.The average number of nucleotide differences for all alleles was k = 3.581.The mismatch distribution indicated that the Harpending's raggedness index was r = 0.0785, Tajima's D neutrality test conducted both on the Eastern and the Central-Western group of populations separately indicated that the Tajima's D = -1.092,p < 0.05 and Fu's Fs = -4.187,p < 0.009 for the Eastern group, and Of the 15 ITS2 haplotypes 8 were singletons and the remaining 7 were shared between two or more populations (S3).H1 (N H1 = 129) was the most frequent allele shared among all populations.Pairwise difference between alleles ranged from 0.2% to 2.8% (1 bp to 12 bp, respectively) (S4).The highest difference of 12 bp was found between H8 (from Elazığ) and H9 (from Giresun).Average allele and nucleotide diversities were 0.31 ± 0.10 and 0.0015 ± 0.0013, respectively (Table 1).Allelic diversity estimates ranged between 0.0 and 0.8 ± 0.1.Bingöl population has the highest diversity (0.8 ± 0.1) followed by Tokat population (0.69 ± 0.10).Further, nucleotide diversity ranged from 0.00 to 0.012 ± 0.010.The highest nucleotide diversity was found in Giresun (0.012 ± 0.010) followed by Elazığ (0.0030 ± 0.0023).
No spatial genetic autocorrelation was detected across the study area for distance classes of 50 or 100 km with 95% confidence intervals (results not shown).Thus, the nuclear ITS2 region did not resolve any spatial patterns for T. synaspis in Anatolia.

Phylogenetic and phylogeographic analyses
Three separate sets of analyses using maximum parsimony (MP), maximum likelihood (ML), and Bayesian analysis were conducted on both mitochondrial haplotypes and the non-coding ITS2 alleles separately.The MP and ML phylogenies for the cyt b data set were similar though with different bootstrap values, and were almost identical tree to the Bayesian consensus tree.Thus, only a single Bayesian consensus tree and the results of the Beast analysis are shown in Fig. 3. Five main clades are revealed.
At the most basal part of the tree, haplotypes H11 (from Bitlis) and H37 (from Konya) formed a monophyletic group referred as "Central-Eastern Populations".The most basally located haplotypes in the resulting tree indicate that the ingroup haplotypes and the outgroup haplotypes are from very distant localities (Konya and Bitlis), arguing against hybridization events between these two closely related species (Nicholls et al., 2012).On the basis of sequence divergence, the Central-Eastern Population haplotypes separated from the rest around 2.10 million years ago (MYA).All other haplotypes formed the remaining four major haplogroups; of these "Eastern Populations I" diverged from the others around 1.72 MYA.Likewise, the clade "Eastern Populations II" diverged about 1.53 MYA from a clade which bifurcated into a subclade which is a mix of haplotypes referred to as "Eastern, Central, Western Populations".It is also noteworthy that around 1.15 MYA a further bifurcation event produced two small haplogroups within this subclade.
A minimum spanning network analysis to resolve the possible reticulations and to show the evolutionary relationships among the cyt b haplotypes produced congruent results particularly to those in the Bayesian tree.Certain groupings among haplotypes are notable in the network where Bitlis and Konya haplotypes (H11 and H37, respectively) are connected to the small groups of haplotypes with some unsampled intermediate haplotypes.Besides, H49 found only in Tunceli as a singleton seems to be connected to most other eastern haplotypes (Elazığ, Erzincan, Bingöl, Bitlis) with two other haplotypes such as H43 and H44, both singletons found in Sivas population (Fig. 4a).
For revealing the hierarchical distribution of the genetic variation at different levels several grouping schemes co- sidering geographic origins of population were conducted (Table 2a).When populations were grouped into two major geographical groupings (Eastern and Central-Western) separated by the Anatolian Diagonal, among-group variation was 13% compared to among populations-within groups (45%), and within populations (40%).We tested several other grouping schemes; lower values were obtained in all analyses, and we do not consider these results further.
We also conducted three different sets of analyses on the ITS2 alleles of T. synaspis.As for the cyt b haplotype data set, both MP and ML analyses produced the same tree topology with different bootstrap values.The resulting ML tree was a consensus of 186 best found trees with CI = 0.57, RI = 0.128, however the MP tree was the consensus tree of 804 trees with CI = 0.78, RI = 0.32.In the ML/MP consensus tree, H9 (from Giresun) was the most basal allele.The remaining alleles formed a large polytomic group where H2, H15 and H5, from Adıyaman, Bingöl, Elazığ and Tunceli populations are grouped together.On the other hand, H7 (from Bitlis and Çorum) and H13 (from Tokat and Yozgat) formed a monophyletic group.The Bayesian consensus tree was highly similar to the MP/ML consensus tree topology with the exception of two monophyletic small clades grouping H2-H15, and H7-H13 respectively.As for the MP/ML trees, the Bayesian tree also supported H9 as the most basal allele, separated from the remaining alleles which form a large polytomy.Since none of the analysis provided enough resolution to resolve the phylogenetic relationships among the ITS2 alleles these trees are not presented here.
H9 was also found as the most distantly connected haplotype with seven hypothetical alleles connecting H9 to H1 in the ITS2 minimum spanning network (Fig. 4b).H1 (N = 129) being the most frequent and common allele is shared among all studied populations of T. synaspis, producing a star-shaped network.It is remarkable that H4 (N = 1, from  Bingöl) is connected to several haplotypes including H6 (from Bingöl) and two other unsampled haplotypes are connected to a singleton haplotype (H8) found in Elazığ.
In addition, H2 (N = 12 from Adıyaman, Bingöl, Elazığ and Tunceli) is connected to the two haplotypes H5 (N = 1) from Bingöl, and H15 (N = 1) from Tunceli.AMOVA analyses conducted for the ITS2 region are shown in Table 2b, in which grouped of populations into the East and Central-Western groups introduced above showed high and significant genetic variation (70%, p < 0.001) at the within-population level.Differences between east and the central-western regions accounted for 10.15% of the total variation, while 19.21% of the variance could be attributed to variation among populations within regions.

Genetic diversity estimates: mtDNA and nDNA
The large numbers of singleton haplotypes revealed in this study are quite striking.In fact, almost half of the cyt b haplotypes (22 out of 50) and more than half of the ITS2 (8 out of 15) alleles were singletons, and most were also private to a single population.The presence of many singletons that differ by small numbers of nucleotide differences from other haplotypes may indicate a recent origin for them (Crandal & Templeton, 1993).
Considering haplotype and nucleotide diversity we found that high haplotype and nucleotide diversity (0.855 and 0.013, respectively) was conspicuous for the cyt b gene segment of T. synaspis.High genetic diversity has been also reported for other animal species including oak gall wasps from Anatolia (Rokas et al., 2003;Challis et al., 2007;Mutun, 2010;Dinç & Mutun, 2011) and high genetic diversity in oak gall wasps between Europe and Iran (Stone et al., 2012).Çankırı and Yozgat populations showed the highest haplotype diversity while Karaman and Kırşehir displayed the highest nucleotide diversity; all are located in central or western Anatolia (see Table 1).Nonetheless, the nDNA ITS2 sequences also yielded significant allele and nucleotide diversity (0.312 and 0.0015, respectively).The highest allelic diversity was observed in Bingöl followed by Tokat population.On the other hand, the highest nucleotide diversity was in Giresun and Elazığ populations.Previous studies on oak gall wasps indicated that highest genetic variation was in the Anatolian populations.Indeed, in Andricus coriarius nucleotide diversity was 0.00529 and 0.00644 in Iranian and Lebanese populations, respectively; however nucleotide diversity was 0.01554 in Turkey (Challis et al., 2007).Likewise, in Andricus quercustozae the greatest genetic diversity was found in Turkey within a range of 0.2-4.2%,while in contrast in Italy it ranged between 0.12-0.17%(Rokas et al., 2003).Other similar investigations on Anatolian oak gall wasps showed that the genetic diversity across the Turkish populations was quite high.Average haplotype and nucleotide diversities in Andricus lucidus were 0.809 and 0.115, respectively (Mutun, 2011).Similarly, in Andricus caputmedusae haplotype and nucleotide diversities were 0.463 and 0.101 (Mutun, 2010), and in another study conducted on A. quercustozae using RFLP haplotype and nucleotide diversities were 0.45 and 0.05, respectively (Dinç & Mutun, 2011).
It is well-known that high genetic diversity and differentiation between populations within a species is usually coupled with the distinct phylogeographic structure (Avise, 2004).From the phylogeographic point of view, glacial refugia are assumed to harbor high diversity of genotypes and major lineages within a species (Hewitt, 2000).In T. synaspis, the existence of high genetic diversity in Turkey can be explained by concurrent results provided by this and previous studies on other oak gall wasps that identify Anatolia as a possible major refuge and genetic diversity hotspot area (Challis et al., 2007;Stone et al., 2007).In fact, Rokas et al. (2003) provided persuasive data that Anatolia was a cradle for oak gall wasp diversity.Conspicuously high amount of genetic diversity in species/lineages from Turkey has been also reported in plants (Ansell et al., 2002), and various other animal taxa such as the bicolored shrew (Dubey et al., 2007), the yellow-necked field mouse (Michaux et al., 2004), the Anatolian mountain frog (Veith et al., 2003), fishes (Hrbek et al., 2002), ground squirrels (Gündüz et al., 2007), and grasshoppers (Korkmaz et al., 2010).Thus, the current findings of the high genetic diversity in T. synaspis seems to provide more concrete data for this statement and shed more light on the existence of high molecular variation in this unique area and its significance for maintaining genetic richness.

Demographic analysis
Tajima's D values for cyt b were negative for both "East" and "Central-Western" groupings of T. synaspis populations, indicating a bias toward rare site variant haplotypes which may imply signals of recent population growth.Our results were, however, not statistically significant.Fu's F S test results were also significantly negative, indicating that an excess of rare haplotypes were present in each of the Eastern and the Central-Western population groups.It has been shown that Fu's F S test is more powerful than Tajima's D (Ramos-Onsins & Rozas, 2002), this may explain the differences in significance for T. synaspis analysis results.The higher value (Fu' F S = -14.31p < 0.001) of the Central-Western group compared to the Eastern group may imply recent population expansion, particularly in the Central-Western part of Anatolia.Similar analysis on ITS2 region for both groups of populations as Eastern and Central-Western produced significant negative Tajima's D and Fu's Fs implying historical population growth.Alternatively, significant negative values could indicate directional or balancing selection on the locus itself, or on a closely linked locus, although direct or indirect selection via hitchhiking on the mitochondrial genome is rare (Ruiz-Pessini et al., 2004).
Phylogenetic and phylogeographic structure of the Anatolian populations of T. synaspis use of both cytoplasmic and nuclear markers is preferred over use of a single gene to reveal the phylogeographic patterns (Presa et al., 2002;Mardulyn et al., 2011;Turmelle et al., 2011;Cornetti et al., 2014).In this study, we analyzed both mitochondrial and nuclear DNA markers with different resolving power through the application of three separate sets of phylogenetic analyses on each data set separately.All phylogenetic analyses conducted on the cyt b sequence data placed two haplotypes either as a monophyletic group or singly, into the basal part of the trees; an eastern haplotype (H11 from Bitlis) and a Central Anatolian haplotype (H37 from Konya).It is remarkable that both haplotypes have high numbers of nucleotide differences from other haplotypes in pairwise comparisons.Other studies on oak gall wasps have also found most eastern haplotypes to be basal compared to more western haplotypes (Mutun, 2010;Dinç & Mutun, 2011).
The Bayesian cyt b haplotype tree placed the haplotypes with respect to their geographic origin as East and Central-Western.Further support for the same division was provided by the minimum spanning network analysis.It is known that both distribution of species diversity and genetic lineages across Turkey point to an east-west separation in different plant and animal taxa (Ekim & Güner, 1986;Çıplak, 2004).Indeed, a well-defined east and west separation of the oak gall wasp haplotypes from A. caputmedusae (Mutun, 2010) and A. quercustozae (Dinç & Mutun, 2011) have been reported previously.Current findings of this study support previous studies with more superficial east and central-western distinction being observed in T. synaspis.
On the other hand, our AMOVA analysis based on F ST statistics revealed only a small amount of variation being partitioned between the two regions when the populations were clustered into East and Central-Western groups.Partitioning of genetic variation between the two regions is relatively low but significant (F ST = ~ 13%), while each region retains high genetic variation (~ 45%) with the highest partitioning of the variation being observed within populations.However, all our efforts to test different trial schemes produced lower among-group variation producing the highest and statistically the most significant differentiation being between East and Central-Western groupings (see Table 2a).Similar results were also obtained with the ITS2 region, with the same grouping considering the geographic locations of the sampled populations (see Table 2b).Overall, our AMOVA results imply that populations within each region located on the eastern and western part of the mountain lines had enough time to accumulate high genetic variation within each population.
Contrary to our cyt b analysis, ITS2 sequence analysis produced tree topologies with more polytomies.In recent years a number of studies have shown that both mitochondrial and nuclear genes can lead to markedly different phylogenies as a consequence of different inheritance pathways, selection pressures, and differences in response to evolutionary processes in both genomes (Shaw et al., 2002;Zhang & Hewitt, 2003).On one hand, nuclear DNA sequences in animals evolve more slowly than mtDNA (Brown et al., 1979), and mtDNA becomes monophyletic faster than nuclear genes because of its effective population size being one-fourth of that for nuclear genes (Moritz et al., 1987).On the other hand, ITS2 is a non-coding region and there is no known evolutionary pressure on this DNA region; our analyses of the ITS2 region yielded almost no resolution.Thus, the current results and differences in the power of resolution of both DNA regions may all be related to the slow rate of ITS2 evolution and may be considered to limit its use in intraspecific studies in oak gall wasps.Furthermore, the spatial genetic structure analyses produced more supporting results for the cyt b but not for the ITS2 region.An analysis including additional nuclear markers with more power of resolution could give a more complete intraspecific perspective for T. synaspis.

Implications for possible responses to the climatic changes
The Quaternary climatic oscillations resulted in repeated range contraction and expansion of species shaping genetic structure of populations (Hewitt, 2000).Such retreats/ expansions as well as the associated subdivision of species' populations among multiple glacial refugia has been reported to maintain high intraspecific genetic diversity (Weir & Schluter, 2007).The climatic oscillations started around 2.4 MYA and resulted in the growth of ice sheets in the Northern hemisphere, by which time the geological changes in Anatolia had already come to an end (Rögl, 1999).Rigorous statistical test gives confidence to the interpretation of the association of haplotype diversity and distribution, and allows us to infer scenarios that can be dated to the Pleistocene and the Quaternary with its successive climatic oscillations.Since there is no molecular clock for ITS2 our estimation of population divergence times was based only on the cyt b data set.Our estimations indicated that diversification of the outgroup haplotypes from the ingroup was around 2.3 MYA, which corresponds to the Late Pliocene.Considering the processes in the final formation and minor changes in the topographic structure of Anatolia in the late Pliocene it might be reasonable to imagine that both T. megaptera and T. synaspis separation might have taken place around the late Pliocene.
Application of the conventional insect mitochondrial DNA clock to the cyt b data set gives the following possible scenario for haplogroup formation in T. synaspis: in the beginning of the Pleistocene glacial times it seems that an mtDNA lineage gave rise to the small monophyletic clade including Bitlis and Konya haplotypes (H11 and H37) followed by a separation of the Konya-Karaman haplogroup.From Fig. 3 it is obvious that most of the deeper separation of the T. synaspis haplotypes occurred during different subglacial times.Around 1.91 MYA (Early Pleistocene) some of the eastern lineages (Bitlis, Muş, and Tunceli) formed a distinct subclade structure with high statistical support.Further separation between mitochondrial lineages seemed to occur around 1.72 MYA, giving rise to two clades: "Eastern Populations II" and a large clade comprising haplotypes from central, eastern, inner western and northern parts of Anatolia.The lineage giving rise to "Eastern Populations II" diversified from the rest around 1.53 MYA.The second group was also recurrently shaped from the middle to the Late Pleistocene suggesting that more recent diversification events also took place within the Turkish line-ages around this time period.For instance, a central/ west grouping was formed about 0.76 MYA which corresponds to the Günz glacial period (790,000-950,000 years BP) of the Pleistocene.Likewise, at the end of the Early Pleistocene it seems that dispersal may have occurred from the east towards the central and western distribution range of T. synaspis.Moreover, in all clades more recent and shallow diversification events are apparent that date to the Mindel (420,000 to 30,000 years BP) and Riss (140,000 to 80,000 years BP) periods.It seems that with the beginning of the Pleistocene glacial cycles some diversification events following isolation seem to have reoccurred several times among the T. synaspis lineages, which may support the Pleistocene history of Anatolia.We believe that more research in different taxa distributed in Turkey will be of value to establish the generality of this conclusion.

Fig. 2 .
Fig. 2. Correlograms showing spatial genetic correlation coefficient (r) as a function of geographic distances in the entire study area.a) 50 km distance class, b) 100 km distance class.Dashed lines indicate upper (U) and lower (L) boundaries of 95% confidence intervals under the null hypothesis of no geographic structure is present.95% bootstrapped confidence intervals are shown around each of the calculated r value.

Fig. 3 .
Fig. 3.The Bayesian consensus tree of the cyt b gene haplotypes with the node ages (tMRCA) are shown on the interior part of each related node.Posterior probability values are shown over the branches.

Fig. 4 .
Fig. 4. Results of minimum spanning network analysis of a) cyt b haplotypes and b) ITS2 alleles of T. synaspis.Color coded populations are shown in the squares.In both networks hypothetical haplotypes are shown by small black circles.Haplotype numbers are given inside the color coded circles where circle size is proportional to the frequency of the haplotypes.

Table 1 .
T. synaspis localities with their coordinates, abbreviations and sample size, haplotype (h) and nucleotide (π) diversity of each population.

Table 2 .
AMOVA analysis results of (a) the cyt b gene sequence data and (b) nuclear ITS2 region sequence data set.Only the results of two groups are provided here as first group includes eastern populations (BTL, MUS, BNG, ERC, ELZ, TUN, ADY), and the second group includes all other populations located in the western part with respect to the mountain line.