Open Access
Translator Disclaimer
1 July 2010 Phylogeography of Anterastes serbicus Species Group (Orthoptera, Tettigoniidae): Phylogroups Correlate with Mountain Belts, but not with the Morphospecies
Battal Çiplak, Sarp Kaya, İslam Gündüz
Author Affiliations +

Ten species of the genus Anterastes Brunner von Wattenwyl (Orthoptera, Tettigoniidae) show insular distribution in mountain meadows of Anatolia and the Balkans. Current understanding of the taxonomy and species relationships within the genus is based on morphological characters. However, the extent to which morphological characters are phylogenetically informative, when used to define taxonomic groups or to elucidate detailed evolutionary relationships within Anterastes, is in need of further examination. Moreover, because little is known about the historical biogeography and diversification factors in members of this genus, additional datasets are necessary to test the robustness of species, relationship hypotheses and associated biogeographic patterns, Here we specifically examined, using 16S rDNA sequences, the evolutionary relationships and species boundaries of three closely related species of Anterastes (i.e., the A. serbicus group, comprising A. serbicus, A. burri and A. antitauricus). Additionally A. tolunayi, a species not in the A. serbicus group, but morphologically very similar to members of the group, was included in the molecular analysis to locate the species of the A. serbicus complex within a phylogenetic frame. Hence, the phylogenetic relationships and taxonomic interpretation of the species complex appear more intricate than previously hypothesized, The current molecular data do not allow us to identify A. serbicus, A. burri and A. antitauricus as distinct phylogenetic species, but rather suggest that these morphospecies are themselves a complex of cryptic taxa. Despite the incongruencies among the phylogenetic trees and nonmonophyly of each the three morphospecies, the median joining network resulted in haplotype grouping consisting of four clusters that are definable by geography. Thus, based on the congruency between geography and gene clusters, and the molecular clock estimate, it can be interpreted that 1) a strong correlation between the radiation of the group and the topography of their ranges may exist, 2) the radiation of the group dates back to Late Pliocene or Early Pleistocene and 3) there is a break between the Anatolian and the European lineages, in respect to range change of cold-preferring forms, dating back prior to the last four glacial periods.


Discovering and describing species diversity is a fundamental part of evolution and biodiversity studies, and this task is far from complete, especially among poorly studied groups, i.e., almost all invertebrates. This is certainly the case for the bush cricket genus Anterastes Brunner von Wattenwyl, 1882 (Orthoptera, Tettigoniidae), particularly for the A. serbicus group distributed in Anatolia and the Balkans.

The complex includes three described species, Anterastes serbicus Brunner von Wattenwyl, 1882, A. burri Karabağ, 1951 and A. antitauricus Çıplak, 2004a. Despite the small number of species, and an effort to discover diagnostic morphological characters distinguishing members of this complex from each other and from others in the genus (Çiplak 2004a), the A. serbicus species complex appears to represent evolutionarily distinct lineages. Additional work is clearly needed to delimit the several members of these lineages, on morphological criteria alone, to test the morphospecies approach, to reveal the number, phylogenetic relationships and geographic distribution of evolutionary lineages, to propose historical biogeography and demography scenarios that accommodate the genetic variation observed among and within lineages and to detect, if any exists, hidden cryptic species diversity.

More accurate phylogenies have allowed authors to test previously suggested taxa (Elmer et al. 2007, Stephenes & Wiens 2009, Cicconardi et al. 2010) or to reconstruct their speciation in correlation with their historical ranges (Avis 2009). There is an overwhelming literature on phylogeography of West Eurasian lineages which suggests that present European taxa evolved mainly from the ancestral stocks inhabiting glacial refugia such as Iberia, Italy and the Balkans (Schmitt 2007). Although Anatolia (Asia Minor or Asian Turkey) has been considered to be a fourth glacial refugium in some studies (Hewitt 1996, 2004; Çıplak 2004a, 2004b, 2008), only a few regional or peninsula-wide phylogenetic studies showing different phylogeographic patterns in different taxa have been conducted (Bohlen et al. 2006; Dubey et al. 2007; Gündüz et al. 2005, 2007; Veith et al. 2003, 2008).

Any further phylogenetic study that may provide information on the phylogeographic history of the Anatolian fauna and flora is of interest. This is because 1) Anatolia has a very active tectonic history through the Neogene in the Mediterranean (Bozkurt 2001) which possibly has repeatedly led to dispersal/vicariance of the inhabiting populations (Kosswig 1955, Çıplak 2004b, Hberk 2002, Hberk et al. 2003), creating major genetic breaks within species. 2) It is an important crossroad at present (possibly also in the past) for the faunal/floral exchange between Europe, Asia and even Africa, via the Arabian Peninsula. 3) It was not covered by an ice sheet (the 52nd latitude is considered the border of permafrost; Hewitt 2000) during glacial periods; being located just south of the border of the permafrost, it may have played a refugial role. 4) It has a very heterogeneous topography/climate that leads to habitat fragmentation, geographic barriers and a biodiversity that is not uniformly distributed throughout the peninsula.

Large scale geographic definitions may indicate the importance of Anatolia as a general refugium, but these do not allow us to elucidate the importance of this region as a system of refugia and speciation areas. Firstly, Anatolia cannot be considered as a homogenous refugium, because most taxa could not spread into the whole Anatolian landmass; rather each taxon could extend its distribution only to particular districts, using suitable corridors during range shifts from the northern territories or vice versa (Hewitt 1996, Tarkhnishvili et al. 2000, Veith et al. 2003, Seddon et al. 2001, Çıplak 2004a). Secondly, large scale geographical definitions focus mainly on latitudinal range changes and these are trivially related to vertical changes. Anatolian topography might have played an important role in Pleistocene range changes both by constituting distribution corridors for latitudinal dispersals, and by providing suitable habitats during altitudinal shifts (Çıplak 2004a). Thirdly, the southern geographic location of the Anatolian refugium requires paying special attention to its role in the matter of rear-edge concept. According to the rear-edge concept (Hampe & Petit 2005), the cold-preferring lineages or populations should be fragmented from the south of their range. The expectation is that such forms disappeared in lowlands, and were confined to some highlands, as isolated populations, the fate of which is erosion in genetic diversity (Hewitt 2004) or drift (Knowles & Richards 2005). Thus, designation of a meaningful phylogeographic study or testing validity of present hypotheses requires considering this dual refugial role of Anatolia in combination with the influence of its extremely complex topography.

There are a few studies focusing on the horizontal distribution patterns around Anatolia (Rokas et al. 2003, Gündüz et al. 2007, Challis et al. 2007, Stone et al. 2007); however, the generality of these patterns needs to be tested using different lineages. Anatolia harbors a highly endemic biodiversity at different taxonomic levels (for Orthoptera see, Çıplak & Demirsoy 1995, 1996; Çıplak 2003) and radiation of the endemic mountainous forms is suggested to occur mainly due to climatic shifts of the Quaternary (see Çıplak 2003, 2004a, 2008). Possibly the vertical range changes were more important for the species or lineages endemic to Anatolia.

The present study, of a four-species lineage, mainly aims to address the following questions: 1) What is the portion of Anatolian forms in establishment of present populations of nonrefugial northern areas and how far could these forms have extended their range during the Holocene? 2) Are there lineages, reflecting past vertical range changes (in altitudinal heterogeneities) or fragmentations where this kind of range change could be defined by using the phylogeography of such lineages? Anatolia is one of the southernmost refugia and the Mediterranean Sea constitutes a strong barrier in the south for terrestrial species/populations. Therefore, Anatolia is a perfect area to test the rear-edge concept. 3) Is the rear-edge concept applicable to Anatolian lineages? The starting point for conservation activities is defining evolutionarily independent units and their genetic diversity. Speaking for Anatolian taxa, this is especially important for cold-preferring taxa isolated on its southern highlands or sky islands (Çıplak 2008). 4) Is each highland population isolated and is the isolation more prominent for those in the south?

Table 1.

Sampling localities and genetic diversity indices of four species. Shown from left to right are: the sample sizes (n), the number of segregating sites (s), the number of haplotypes without gaps (k), the number of haplotypes including gaps (K), haplotype diversity (h), nucleotide diversity (π); NA, not included in population genetics analyses.


The Anatolian-Balkanian genus Anterastes has recently been revised by Çıplak (2004a) who presented a morphology-based phylogeny and a biogeographic account. The genus includes nine species that are endemic to Anatolia, and one (A. serbicus) that also extends its range to the Balkanian highlands. As presently known, their populations prefer cold climates and are present on alpine/subalpine regions of highlands around/above 1500 m altitude. Some of the species are only known from a single upland, while others are known from more than one highland, but each also occurs as an isolated population. For example, A. serbicus that is present both in North Anatolia and the Balkans has not been recorded from lowlands of the Balkans and Anatolia (Çıplak 2004a, Çıplak unpublished observations). Since this is the distribution pattern/habitat preference observed in all species, the genus is a model candidate lineage to answer the questions just listed.

On the other hand, as mentioned, species of the genus differ in minor diagnostic characters (Çıplak 2004a). One of the species group in the genus is the A. serbicus group which consists of three species: A. serbicus (Anatolia+Balkans), A. burri (Northwest Anatolia) and A. antitauricus (South Anatolia). Each of these species is suitable for testing the above statements, since they are presently known from three or more different altitudes. On the other hand, A. tolunayi, a species not in the A. serbicus group, but very close through shared similarities (Çıplak 2004a), known from a single altitude, was also included in the study in order to locate the species of the A. serbicus complex within a phylogenetic frame. Using 16S rDNA sequences obtained from different populations of these four species we produced a phylogeny with the aim of testing the above addressed questions.

Materials and Methods

Sampling and DNA extraction.—Twelve populations of four ingroup species (A. serbicus, A. buni, A. antitauricus and A. tolunayi) and one population of an outgroup (A. uludaghensis), were sampled from Anatolia and the Balkans (Tables 1, 2; Fig. 1). Morphological species identification of the nominal species was based on differences given by Çıplak (2004). Total DNA was extracted from 153 adult specimens preserved in 96 % ethanol, kept at +4°C (25 belonging to A. serbicus, 56 to A. burri, 43 to A. antitauricus, 20 to A. tolunayi and 8 to outgroup A. uludaghensis) with the salt-extraction method (Aljanabi & Martinez 1997). The 16S rDNA gene was amplified using the universal primers LR-J-12887 and LR-N-13398 described by Simon et al. (1994). Amplifications were performed in a 50-µl volume containing 0.3 µl of each primer (100 µM), 1 µl dNTP mix (10 mM), 2 µl 50 mM MgCl2, 5 µl 10X Platinum PCR buffer (containing 200 mM Tris-HCl [pH 8.4], 500 mM KCl), 1.25 U Platinium Taq DNA polymerase (Invitrogen), and 0.5–1 µl of 50–70 ng template DNA. Temperature cycling was carried out in an Eppendorf Mastercycler Personal. Cycling conditions were 1 min denaturing at 95°C; (45 s at 95°C, 40 s at 49°C, and 50 s at 72°C) x 35. PCR products were purified and sequenced in both directions. The sequencing products were loaded onto an ABI 3730 XL automated sequencer.

Table 2.

Geographic information and indices of genetic diversity. Shown from left to right are; the sample sizes (n), the number of segregating sites (s), the number of haplotypes without gaps (k), the number of haplotypes including gaps (K), haplotype diversity (h), nucleotide diversity (π).


Phylogenetic analyses.—The sequences were aligned manually in BioEdit [version 5.09] (Hall 1999). DnaSP [version 5] (Librado & Rozas 2009) was used to determine unique haplotypes and to calculate their frequency. Sequences were deposited in GeneBank. The phylogenetic signal in the data partitions was estimated by maximum likelihood mapping method (Strimmer & von Haeseler 1997) using the program TREE-PUZZLE [version 5.2] (Schmidt et al. 2002). Initially, the evolutionary relationships between haplotypes were examined with the software Network version ( (Bandelt et al. 1999). Later, the aligned sequences were analyzed by maximum parsimony (MP) and maximum likelihood (ML), as implemented in PALIP [version 4.0b10] (Swofford 2000) as well as by using Bayesian inference of phylogeny as implemented in MrBayes [version 3.1.2] (Ronquist & Huelsenbeck 2003). Bayesian search was carried out using four simulations of Markov chains, five million generations and sampling every 100th generation. The software TRACER [version 1.5] (Rambout & Drummond 2003) was used to examine the parameters and determine the number of trees needed to reach stationarity (burn-in). Bayesian posterior branch probabilities were obtained by taking the majority rule consensus of the sampled trees, excluding the first 5000 trees as burn-in. Prior to calculations all ambiguously aligned positions and gap-loaded positions were excluded from the dataset. The MP analyses were carried out 10 times with the heuristic search approach using the TBR algorithm. The confidence of branching was assessedusing 100 nonparametric bootstrap resamplings (Felsenstein 1985). The parameters and the best-fit model were estimated using MODELTEST [version3.06] (Posada & Crandall 1998). The selected model was implemented in the ML analysis performed with PALIP. The ML tree search was conducted using the heuristic search approach, the ‘as is’ addition replicate. Nonparametric bootstrapping (Felsenstein 1985) was used to evaluate the support of nodes based on 1000 pseudoreplicates.

Fig. 1.

Localities sampled for one outgroup and four ingroup species (for abbreviations see Table 2).


Estimation of divergence times.—To estimate divergence times and to correlate them with the historical events, a generalized clock was applied. Divergence dates among the main mitochondrial clades were estimated using a Bayesian Markov Chain Monte Carlo (MCMC) approach, as implemented in BEAST [version 1.5.2] (Drummond & Rambaut 2007). Prior to running BEAST all ambiguously aligned positions and gap-loaded positions were excluded from the dataset. Several related analyses were run for the same alignment with the model selected by MODELTEST. The haplotype network in Fig. 3 and the phylogenetic tree in Fig. 4 were considered to constrain the analysis with regard to the topologies sampled. Dates of the divergence were inferred using a relaxed molecular clock, following the uncorrelated relaxed lognormal clock as implemented. The dataset was calibrated by general substitution rate (1.4%, per lineage, per 700 kyr) as suggested by Brower (1994) for 16S rDNA in Laupala (Orthoptera). All analyses were performed assuming a Yule process of diversification. BEAST was run for 20 million generations, sampling every 1000 and the convergence to stationarity and the effective sample size (ESS) of model parameters were checked by TRACER. The maximum clade-credibility tree built with TREEANNOTATOR, discarding the initial 10% of samples as burn-in [FIG-TREE version 1.3.1] (Rambaut 2008) was used to visualize the results, including the confidence intervals.

Analyses of genetic diversity.—The number of haplotypes, either without gaps (k) or with gaps (K), haplotype diversity (h), the number of polymorphic sites (S) and nucleotide diversity (π) were calculated for geographical populations and species separately. With unequal numbers of individuals, varying from 4 to 28 per population in the study, Pearson°s correlation tests were conducted to test value independence of the number of haplotypes (k), haplotype diversity (h) and nucleotide diversity (π), against the number of individuals (N). To evaluate the significance of differentiation between populations based on pairwise F ST, an exact test (Reymond & Rousset 1995) was performed with a Markov Chain length of 10,000 and 1,000 burn-in steps. Analysis of molecular variance (AMOVA) (Weir & Cockerham 1984) was conducted to measure the extent to which genetic variance is assigned to the three hierarchical levels of organization: 1) between morphospecies, 2) between geographic populations within morphospecies and 3) within geographic populations. The statistical significance of AMOVA was evaluated using 10,000 permutations. All above analyses were performed using Arlequin 3.01 software (Excoffier et al. 2005) under the TrN+Γ model, since it is suggested that it is the more general case of the HKY (Hasewaga et al. 1985).

Fig. 2.

Results of likelihood mapping analysis. The values in the panels indicate proportion of fully resolved (corners), partially resolved (along the sides) and fully unresolved quartets (in the center).



The outgroup A. uludaghensis is sampled from a single locality (it is so far known only from Bursa, Uludag Mt). Of the ingroup species, A. tolunayi was found from one locality (Izmir), A. serbicus (Sivas, Giresun and Kastamonu provinces of Turkey, and Bulgaria) and A. burri from four (Uşak, Bursa, Bolu and Çankırı provinces) and A. antitauricus from three (Karaman, Konya and Afyon provinces) (Tables 1, 2; Fig. 1). From the sequence analyses, 41 haplotypes were identified among the 153 samples belonging to one outgroup and four ingroup Anterastes species: single haplotype for A. uludaghensis (outgroup, n=8) and A. tolunayi (n = 20), 12 for A. serbicus (n = 25), 18 for A. burri (n = 56) and 9 for A. antitauricus (n = 43) (Table 1). Two of these 41 haplotypes (one from A. antitauricus and the other from A. burri) are definable by the gaps in three sites. The phylogenetic and time-estimation analyses were made using 523 bp of 39 haplotypes (excluding the gaps), but a dataset of all was used in comparative genetic analyses. Of the 523 sites, 443 are constant, 80 are variable (50 of which are parsimony informative). There is no haplotype shared either by two or more different morphospecies or by two or more geographic populations.

Phylogenetic analyses.—Results of likelihood mapping are presented in Fig. 2. Poor dichotomic phylogenetic signal was detected in the dataset.

The percentage of the quartets suggesting a star- or network- like phylogeny is 17.7%, an amount close to threshold, indicating that data are not reliable for a dichotomic phylogenetic analysis (Lemey et al. 2009). Thus, we initially performed a haplotype network analysis. This analysis suggested 1IZ and 2US to be distinctly diverged by 20 and 17 mutations from respectively most similar haplotypes. Remaining haplotypes constituted four clusters which are definable by geography: 1) Bulgarian cluster (A. serbicus), 2) Bursa-Uludag cluster (topotypical A. burri), 3) Black Sea Region cluster excluding Bolu (Çankırı-A. burri; Sivas, Giresun and Kastamonu - A. serbicus) and 4) Mediterranean + Aegean Anatolia plus Bolu province (from west part of Turkish Black Sea area) cluster (Bolu -A. burri, Karaman, Konya and Afyon- A. antitauricus). The haplotype 15KA occurs as ancestor of others in the last cluster.

MODELTEST, applied to the dataset, selected the HKY model (Hasegewa et al. 1985), according to hLRT with proportion of invariable characters I, of 0.6537, gamma correction Γ, of 0.6851 and ti/tv ratio of 3.9628. This model was implemented in BI and ML analyses. Each of BI, MP and ML analyses resulted in trees with different topologies; however, there are clades supported in all trees (Fig. 4).

The single haplotype from the Izmir population (1IZ, A. tolunayi) occurs as either an independent (BI, Fig. 4B) or a basal branch leading to all others (MP, ML; Fig. 4A, C). The latter clade accounts for all samples of A. serbicus, A. buni and A. antitauricus. The single haplotype (2US) from the Uşak-Murat Mt (Aegean Anatolia) population (A. burri) branches off basally from the remaining in all analyses (Fig. 4).

Fig. 3.

Haplotype network obtained from median-joining haplotype network analysis. Geographic area corresponding to haplotype clusters are indicated. Mutation numbers between haplotype pairs are indicated on the respective line; absence of the number indicates single mutation.


For the next step, MP and BI analyses suggested a polytomy of several branches, however, there are three internal phylogroups well supported in both trees (Fig. 4A, B). One of them accounts for the seven haplotypes from the Bulgarian population (A. serbicus), a second for the eight haplotypes from the Bursa-Uludag population (topotypical A. burri) and the last for the nine haplotypes from Sivas, Giresun, Kastamonu (A. serbicus) and Çankin (A. buni). Network analysis to the contrary (Fig. 3), the remaining 12 haplotypes (A. burri + A. antitauricus) occurred as independent branches of single or dual haplotypes in BI and MP analyses.

Table 3.

AMOVA of four species of Anterastes genus (*P<0.0001; df, degree of freedom).


ML analysis resulted in a different topology. After two successive branches (1IZ and 2US) there is a dichotomy of two haplotype clades: 1) seven haplotypes of A. serbicus from Bulgaria and 2) the remaining. On the next step the nine haplotypes from Sivas, Giresun, Kastamonu (A. serbicus) and Çankırı (A. burri) constitute a phylogroup, while the others make a second. The crown clade includes two phylogroups: 1) eight haplotypes from the Bursa-Uludag population (topotypical A. burri) and 2) 12 haplotypes from Mediterranean + Aegean Anatolia plus Bolu. All of the gene trees supported monophyly of their species clades (A. serbicus, A. burri and A. antitauricus), but not that of each morphospecies.

Fig. 4.

Results of MP (A), BI (B) and ML (C) analyses. Posterior probabilities (BI) and bootstrap supports (MP and ML) to branches are shown.




Estimation of divergence times.—Considering both haplotype network (Fig. 3) and phylogenetic analyses (Fig. 4), BEAST analysis was constrained with the following topologies: 1) outgroup and others, 2) the outgroup+1IZ and others, 3) the outgroup+1IZ+2US and others, and 4) three terminal groups, the first of which to include seven haplotypes from Bulgaria, the second seven haplotypes from Bursa-Uludag and the last 19 haplotypes originating from 15 KA (see Fig. 3). The BEAST analysis, under a HKY+I+Γ model, revealed that this segment of DNA showed a clock-like evolution (ucld.mean=0.52, ucld.stdev = 0.27; corresponding to a rate of 2% per million years; Fig. 5).

From BEAST, we estimated the time to most recent common ancestor (TMRCA) for the entire ingroup and internal clades as follows: 1) to ingroup origin 2.27 Myr (95% HPD interval = 1.534–3.059), 2) to three ingroup species clade (including A. serbicus, A. burri and A. antitauricus) 2.01 (95% HPD interval = 1.385–2.659) and 3) to the clade including Bulgarian plus its sister Anatolian phylogroup 1.56 (95% HPD interval = 1.077–2.094). The first corresponds to Late Pliocene/Early Pleistocene and the other two to Early Pleistocene. The diversification of terminal phylogroups (those corresponding to haplotype clusters suggested by haplotype network analysis) correspond to a time prior to Günz glacial (for TMRCA of each clade see Fig. 5).

Analyses of genetic diversity.—Strong geographic structuring of genetic variation was observed at the 16S rDNA marker. All of the haplotypes were unique to a single geographic population. Diversity indices in A. tolunayi were zero, since 20 samples, represent just a single haplotype. For the three ingroup species, the haplotype diversity (h) varies between 0.8239–0.8933 and the nucleotide diversity (π) between 0.007–0.028 (Table 1). Diversity indices for the geographic populations were calculated separately, since monophyly of morphospecies was not supported and a strong geographic structuring was found. These indices are zero for the Izmir (A. tolunayi, 20 samples), Uşak (A. burri, 16 samples), Sivas (A. serbicus; 4 samples) and Giresun (A. serbicus, 5 samples) populations, since each have a single haplotype. This is the case for the Karaman population (A. antitauricus) also, since only a single individual was available for study (from our observation, this population seems to have gone extinct very recently when its restricted habitat was destroyed).

In the other geographic populations haplotype diversity (h) varies between 0.0053–0.9048 and nucleotide diversity (π) between 0.0019–0.0136 (Table 2). The highest haplotype diversities were found in the Bursa, Kastamonu and Bulgarian populations, while nucleotide diversity was highest in the Bulgarian population. There is no evidence that heterogeneous levels of genetic variability were produced by the unequal sample sizes for each population (Pearson's correlation test between diversity indices and sample size: r= 0.196 and P= 0.564 for number of haplotypes, r= 0.108 and P= 0.752 for haplotype diversity, r= - 0.007 and P= 0.983 for nucleotide diversity).

A significantly high differentiation was observed for both interpopulation, within species, and interhaplotypes, within population comparisons. Based on an AMOVA of the entire 16S rDNA dataset, 22.03% of genetic variation was found to be partitioned between morphospecies, 67.28% among sampling sites inside each species, and 10.68% inside each sampling site, the last two being statistically significant (Table 3). Possibly this is because each geographic population has its unique haplotypes (or there are no haplotypes shared by two or more populations). Mean interspecies pairwise FST values were in the range of 0.1–0.52. A. tolunayi is significantly divergent from all others. The interpopulations differentiation is even higher. Mean interpopulation pairwise FST range between 0.30–1.0. FST values are below 0.5 in a few comparisons, indicating highly diverged geographic populations (Table 4).

Table 4.

Pairwise FST values among localities of four species (only significant values are represented).



The first and possibly most striking aspect of the data is that the gene tree does not support morphospecies (see Çıplak 2004a for morphological data). Of the four species, A. tolunayi seems to constitute a genetic clade, while the other three morphospecies, A. serbicus, A. burri and A. antitauricus, are para(poly)phyletic. Consistent with phylogenetic relationships, AMOVA analysis suggests that genetic divergence between species is not significant, while that between populations is obviously high. Furthermore, several haplotypes diverge from each other by several substitutions and nearly all are joined to each other by hypothetical ancestral haplotypes (see Fig. 3).

There are two possible explanations of this situation. First, the ingroup species, other than A. tolunayi, may be considered as a single species. The small morphological differences between the three species in the A. serbicus group (see Çıplak 2004a) may support this assumption. Second, this three-species clade may include several undescribed, morphologically cryptic, species. Present data support the second assumption, since 1) there are no clues in the data indicating gene flow between morphospecies, since there is no haplotype shared by two or more of them; 2) there are strongly diverged geographic populations (pairwise FST values are mostly higher than 0.5 and run up to 1.0 in most of the comparisons), each harboring only its unique haplotypes; 3) there are some internal phylogroups, such as Bulgarian (A. serbicus), Uşak-Murat Mt (A. burri), Bursa-Uludag Mt (topotypical A. burri), which seem to have been reproductively isolated since the Mid-Pleistocene (Figs 3, 5).

It is suggested that these animals are present only in highland habitats with restricted ecological conditions (Çıplak 2004a). If this habitat preference is the case in their history, it can be assumed that living in similar conditions may lead to a similar morphology because of similar selection pressures (Elmer et al. 2007, Stephenes & Wiens 2009, Cicconardi et al. 2010).

Although the Uşak-Murat Mt population is genetically the most divergent population of these three species, it is a typical A. burri in morphology. This is also the case for several other populations (e.g., Bulgarian, A. serbicus), and all these are consistent with the second assumption. However, it should be kept in mind that gene trees are not species trees (Edwards 2009), thus a nomenclatural decision at this stage will be premature,

A second aspect of the results to be evaluated is that analyses of the data by different procedures did not result in a robust phylogeny. For example, haplotypes of A. antitauricus from a single population (Konya) occur as independent branches in all phylogenetic analyses. Contrary to the phylogenetic analyses, these haplotypes differ from each other in a few base pairs and they constitute a single cluster originated from the haplotype found in the Karaman population in the network analysis (Fig. 3). Similarly, the four haplotypes from Çankırı population either constitute a phylogroup with those from Kastamonu, Giresun and Sivas (BI and ML trees) or independent branches (MP tree).

One reason underlying this inefficiency relates partially to poor phylogenetic signs in mtDNA sequences due to A-T bias (Simon et al. 1994). Results of the likelihood mapping analysis and the percentage of A-T (69.58%) and that of G-C (30.42%) support this assumption. However, there are signs in the data indicating another reason: if the local populations were independent evolutionary units with no or limited gene flow through their divergence from a polymorphic ancestral population, this could have led to loss of ancestral haplotyes. Such mosaic evolution may cause a limited number of synapomorphic sites (Simon et al. 1994).

Since such a process is suggested to have been more severe in recent radiations, the unresolved upper parts of the tree, particularly the polytomy of haplotypes found in A. antitauricus, are consistent with this statement. Additionally, the geographic populations exhibit strong genetic isolation, as suggested from the compelling finding that this phylogeny reflects relatively well, i.e., corresponds to, geography (i.e., haplotypes from the same or nearby localities are always genealogically exclusive and very high pairwise FST values exist among populations). This is as might be expected for species with limited dispersal capabilities (e.g., Bond et al. 2001, Hedin & Wood 2002). All these issues require discussing the phylogeography of the group in detail.

Phylogeography and divergence times.—We now address a few questions to be answered in a biogeographic perspective in the light of the data obtained. Anatolia is suggested as a glacial refugium, playing an important role in the establishment of fauna of non-refugial northern areas during the Holocene (Taberlet et al. 1998, Hewitt 2004).

The first question is: ‘what is the allotted portion of cold-preferring Anatolian forms extending to nonrefugial northern areas during the Holocene, or during the last four glacial periods?' The Balkanian population of A. serbicus is the only representative of this lineage extending beyond Anatolia. La Greca (1999) suggested that this dispersal occurred recently during the Holocene. Present data contradict this statement. Bulgarian haplotypes constitute a distinct and highly divergent phylogroup which share the last common ancestor with a sister Anatolian clade 1.56 My ago (see Fig. 5). This corresponds to a time much earlier than Günz glacial.

Fig. 5.

Chronogram for A. serbicus group. Divergence times inferred according to Bayesian inference with a relaxed molecular clock in BEAST. Bars at the nodes represent the 95% highest posterior density (HPD) credibility interval.


Species of Anterastes prefer cold mountainous habitats, hence a northward range expansion during warm interglacial periods, as suggested for such forms (Hewitt 1996; Çıplak 2004a, 2008), is an expectation. However, the range changes during glacial/interglacial periods are possibly related to dispersal ability. Anterastes are flightless insects and have strict habitat requirements (Çıplak 2004a): both may indicate reduced dispersal ability.

Although most of the hypotheses explaining the high biodiversity in the region refer to Anatolia as a glacial refugium (Taberlet et al. 1998; Hewitt 2004; Çıplak 2004a, b; Çıplak 2008), evidence from this study suggests that the Anatolian Peninsula is also home to old endemic biota. This situation is attested to by the high rate of specific endemism in several groups, such as Tettigoniinae with 83% endemism (Çıplak 2003), Pamphagidae with 75% endemism (Çıplak & Demirsoy 1996), and plants with 33% endemism (Davis 1965–1985, Médail & Quézel 1997, Médail & Diadema 2009). Evidence from phylogenetic studies and the presence of deep genetic divergence in populations (Cooper et al. 1995, Seddon et al. 2001, Rokas et al. 2003, Gunduz et al. 2007, Murienne et al. 2010) also suggest long biogeographic isolation. In this context, the long and complex palaeogeographic history, high habitat heterogeneity, and great topographic and climatic variations of the Anatolia are of fundamental importance in explaining its biodiversity.

The pattern discussed above considers the horizontal range changes in a latitudinal gradient. Although present data do not indicate such a pattern during the last four glacial periods, climatic shifts might have led to vertical range changes of cold-preferring Anatolian forms in an altitudinal gradient (Çıplak 2004a, 2008). Thus, our second question was ‘are there any clues in the data reflecting vertical range changes during warming/cooling cycles of the Pleistocene?’ Since the present populations were sampled during the Holocene warm period, an expectation for these mountainous forms is their restriction to uplands and disappearance from lowlands. Up to the present they have not been found in lowlands and each population is geographically isolated in its respective highland. So the distribution patterns fit the expectation.

The further consistent signs in the data are as follows. First, phylogroups or genetic clusters do not support monophyly of morphospecies, but they correlate well with altitudinal chains (see Fig. 3). Second, the populations from highlands within the same altitudinal chains, fall in the same phylogroups. Third, from our observations each local population is geographically isolated, as represented by unique haplotypes. Fourth, they differ from each other in very high values of pairwise FST; however, the lowest FST values were observed between populations of the same highland chains.

Another expectation related to restriction to uplands is reduction in population size and so loss of genetic diversity (Knowles & Richards 2005). Low haplotype diversity and high nucleotide diversity are suggested as clues to reduction in genetic diversity (Zhao et al. 2008). There are several populations exhibiting this pattern (Table 2). In the populations from Izmir, Uşak, Giresun and Sivas, diversity indices are zero. Also, in Afyon, Konya and Çankırı populations, when compared to others, the nucleotide diversity is higher and the haplotype diversity is lower. Thus, both distribution pattern and population genetic parameters indicate a recent reduction in population size consistent with the expectation.

Anatolia is one of the most southern of refugia, and the Mediterranean Sea constitutes a strong barrier in Anatolia's south, limiting ranges of terrestrial lineages, whose ranges were possibly affected by climatic changes of the Pleistocene. Thus, this peninsula constitutes an ideal area in which to test the rear-edge concept (Hampe & Petit 2005). Although the rear-edge concept applies to fragments of a single species, it may be applicable to lineages radiated in a similar way due to past fragmentations. Thus, our third question: ‘is the rear-edge concept applicable to Anatolian lineages?’

From the perspective of the concept (Hampe & Petit 2005), the fragmented populations on the southern edge are small, with reduced genetic diversity. In this case, the diversity indices should be lower for southern populations than for northern ones. From the latitudinal locations and the climate/vegetation characteristic of Anatolia (Çıplak 2003), of the 12 populations studied, the Konya, Karaman, Afyon, Izmir and Uşak populations may be considered as southern and the others as northern. A comparison of the diversity indices (k, K, h and π) supports the prediction. The diversity indices are zero in Izmir and Uşak and are lower in Konya, Afyon and Karaman, than in the north, although the highest sample size is from the Konya population. The highest diversity indices are higher for those populations from the north (Table 2).

Conservation activities require definition of the conservation units and the documenting of their genetic diversities to assess the plans for their future. Definition of conservation units can either be done on a taxonomic or evolutionarily independent unit basis (Frankham et al. 2006). For Anterastes, a taxon-based definition does not seem to work, since each species, whether valid or not when our results are considered, consists of populations isolated both geographically (on their sky island) and genetically (each has haplotypes not shared by any other). Thus, the conclusion is that each has its own evolutionary trend and needs to be considered separately in conservation plans.

Another issue relating to conservation is the genetic diversity per conservation unit, since the genetically homogenized populations are assumed to be more prone to extinction (Frankham et al. 2006). Haplotype and nucleotide diversities are zero in 5 of 12 populations, indicating a highly reduced genetic diversity. Combining the diversity values with the range width, suggests that some of these populations, especially that of Izmir (A. tolunayi) and Uşak (A. burri), are threatened. Çıplak (2008) suggested that cold-preferring taxa/populations of Anatolia are in general, isolated on highlands, particularly in its southern part. Thus, conservation plans for southern and northern populations require different considerations: an area conservation perspective for those in the north, but a conservation genetic perspective for those in the south. Present data are consistent with this suggestion and validate the statements. Moreover, this study provides a very good example of species-level “paraphyly” and adds to a growing body of literature suggesting that traditional methods for delineating species boundaries in different groups may underestimate actual species-level diversity.


We thank Dragan Chobanov (Bulgaria) for providing specimens from Bulgaria. Dr Deniz Sirin (Turkey) and M. Sait Taylan (Antalya) joined in one of the field studies. The young orthopterist Zehra Boztepe helped in preparation of the manuscript. Our research was supported by the Scientific and Technical Research Council of Turkey (Project No: 107T462). Studies were carried out in laboratories in the Department of Biology, Akdeniz University, and the paper was supported by the Akdeniz University Research Fund.



S.M. Aljanabi , I. Martinez 1997. Universal and rapid salt-extraction of high quality genomic DNA for PCR-based techniques. Nucleic Acids Research 25: 692–4693. Google Scholar


J.C. Avis 2009. Phylogeography: retrospect and prospect. Journal of Biogeography 36: 3–15. Google Scholar


H. J. Bandelt , F. Forster , A. Röhl 1999. Median-Joining Networks for Inferring Intraspecific Phylogenies. Molecular Biology and Evolution 16: 37–48. Google Scholar


J. Bohlen , A. Perdices , I. Doadrio , P.S. Economidis , 2006. Vicariance, colonisation, and fast local speciation in Asia Minor and the Balkans as revealed from the phylogeny of spined loaches (Osteichthyes; Cobitidae). Molecular Phylogenetics and Evolution, 39: 552–561. Google Scholar


J.E. Bond , M.C Hedin , M.G. Ramirez , B.D. Opell , 2001. Deep molecular divergence in the absence of morphological and ecological change in the Californian coastal dune endemic trapdoor spider Aptostichus simus. Molecular Ecology 10: 899–910. Google Scholar


E. Bozkurt 2001. Neotectonics of Turkey — a synthesis. Geodinamica Acta 14: 3–30. Google Scholar


A.V.Z. Brower 1994. Rapid morphological radiation and convergence among races of the butterfly Heliconius erato inferred from pattern of mitochondrial DNA evolution. Proceedings National Academy of Sciences United States of America 91: 6491–6495. Google Scholar


R. Challis , S. Mutun , J.L. Nieves-Aldrey , S Preuss , A. Rokas , A. Aebi , E Sadeghi , M Tavakoli , GN. Stone 2007. Longitudinal range expansion and cryptic eastern species in the western Palaearctic oak gallwasp, Andricus coriarius. Molecular Ecology 16: 2103–2114. Google Scholar


F. Cicconardi , F. Nardi , B.C. Emerson , F. Frati , P.P. Fanciulli 2010. Deep phylogeograpic and long-term persistence of forest invertebrates (Hexopoda: Collembola) in the northwestern Mediterranenan basin. Molecular Ecology 19: 386–400. Google Scholar


B. Çıplak , A. Demirsoy 1995. Türkiye'de Ensifera (Orthoptera, Insecta) alttakımının endemizm açısından değerlendirilmesi. Turkish Journal of Zoology 19: 213–220. Google Scholar


B. Çıplak , A. Demirsoy 1996. Caelifera (Orthoptera, Insecta) alt takımının Türkiye'de endemizm durumu. Turkish Journal of Zoology 20: 241–246. Google Scholar


B. Çıplak 2003. Distribution of Tettigoniinae (Orthoptera, Tettigoniidae) bush-crickets in Turkey: the importance of the Anatolian Taurus Mountains in biodiversity and implications for conservation. Biodiversity and Conservation 12: 47–64. Google Scholar


B. Çıplak 2004c. Biogeography of Anatolia: the marker group Orthoptera. Memorie della Societa Entomologica Italiana 82: 357–372. Google Scholar


B. Çıplak 2004a. Systematics, phylogeny and biogeography of Anterastes (Orthoptera, Tettigoniidae, Tettigoniinae): evolution within a refugium. Zoologica Scripta 33: 19–44. Google Scholar


B. Çıplak 2008. The analogy between glacial cycles and global warming for the glacial relicts in a refugium: a biogeographic perspective for conservation of Anatolian Orthoptera, pp. 135–163, In: Fattorini, S.: Insect Ecology and Conservation (Chapter 6). Research Signpost. Google Scholar


S.J. Cooper , K.M. Ibrahim , G.M. Hewitt 1995. Postglacial expansion and genome subdivision in the European grasshopper Chorthippus parallelus. Molecular Ecology 10: 2187–2198. Google Scholar


P.H. Davis 1965–1985. Flora of Turkey and the East Aegean Islands 1–9. 1965 (Vol. 1); 1967 (2); 1970 (3); 1972 (4); 1975 (5); 1978 (6); 1982 (7); 1984 (B); 1985 (9). Edinburgh. Google Scholar


A.J. Drummond , A. Rambaut 2007. BEAST: Bayesian evolutionary analyses by sampling trees. BMC Evolutionary Biology 7: 214–221. Google Scholar


S. Dubey , J.-F. Cosson , V. Vohralik , B. Kryštufek , E. Diker , P. Vogel 2007. Molecular evidence of Pleistocene bidirectional faunal exchange between Europe and the Near East: the case of the bicoloured shrew (Crocidura leucodon, Soricidae). Journal of Evolutionary Biology 20: 1799–1808. Google Scholar


K.R. Elmer , A.J. Dávila , S.C. Lougheed 2007. Cryptic diversity and deep divergence in an upper Amazonian leaf litter frog, Eleutherodactylus ockendeni. BMC Evolutionary Biology 7: 247–260. Google Scholar


L ExcoffIer , G Laval , S. Schneider 2005. Arlequin ver. 3.0: an integrated software package for population genetics data analysis. Evolutionary Bioinformatics Online 1: 47–50. Google Scholar


J. Felsenstein 1985. Confidence limits on phylogenies: an approach using the bootstrap. Evolution 39: 783–791. Google Scholar


R. Frankham , J.D. Ballou , D.A. Briscoe 2006. Introduction to conservation genetics. Cambridge University Press, London. Google Scholar


T. Garland , A.F. Bennett , E.L. Rezende 2005. Phylogenetic approaches in comparative physiology. Journal of Experimental Biology 208: 3015–3035. Google Scholar


I. Gündüz , R.V. Rambau , C. Tez , J.B. Searle , 2005. Mitochondrial DNA variation in the western house mouse (Mus musculus domesticus) close to its site of origin: studies in Turkey. Biological Journal of Linnean Society 84: 473–485. Google Scholar


I. Gündüz , M. Jaarola , C. Tez , C. Yeniyurt , P.D. Polly , J.B. Searle 2007. Multigenic and morphometric differentiation of ground squirrels (Spermophilus, Sciuridae, Rodentia) in Turkey. Molecular Phylogenetics and Evolution 43: 916–935. Google Scholar


T. Hall 1999. BioEdit: a user-friendly biological sequence alignment editor and analysis program for Windows 95/98/NT. Nucleic Acids Symposium Series 41: 95–98. Google Scholar


A. Hampe , R. Petit 2005. Conserving biodiversity under climate change: the rear edge matters. Ecology Letters 8: 461–467. Google Scholar


M. Hasegewa , H. Kishino , T. Yano , 1985. Dating of the human-ape splitting by a molecular cloack of mitochondrial DNA. Journal of Molecular Biology 22: 160–174. Google Scholar


M.C. Hedin , D.A. Wood , 2002. Genealogical exclusivity in geographically proximate populations of Hypochilus thorelli Marx (Araneae, Hypochilidae) on the Cumberland Plateau of North America. Molecular Ecology 11: 1975–1988. Google Scholar


G.M. Hewitt 1996. Some genetic consequences of ice ages, and their role in divergence and speciation. Biological Journal of Linnean Society 58: 247–276. Google Scholar


G.M. Hewitt 1999. Post-glacial re-colonization of European biota. Biological Journal of Linnean Society 68: 87–112. Google Scholar


G.M. Hewitt 2000. The genetic legacy of the Quaternary ice ages. Nature 405: 907–913. Google Scholar


G.M. Hewitt 2004. Genetic consequences of climatic oscillation in the Quaternary. Philosophical Transactions Royal Society London B 359: 183–195. Google Scholar


T. Hrbek , A. Meyer 2002. Closing of the Tethys Sea and the phylogeny of Eurasian killifishes (Cyprinodontiformes : Cyprinodontidae). Journal of Evolutionary Biology 16: 17–36. Google Scholar


T. Hrbek , F. Kucuk , T. Frickey , K.N. Stolting , R.H. Wildekamp , A. Meyer 2003. Molecular Phylogeny and historical biogeography of the Aphanius (Pisces, Cyprinodontiformes) species complex of central Anatolia, Turkey. Molecular Phylogeneticst an Evolution. 25: 125–137. Google Scholar


L.L. Knowles , C.L. Richards 2005. Importance of genetic drift during Pleistocene divergence as revealed by analyses of genomic variation. Molecular Ecology 14: 4023–4032. Google Scholar


C. Kosswig 1955. Zoogeography of the Near East. Systematic Zoology, 4: 50–74. Google Scholar


M. La Greca 1999. Il contributo degli Ortotteri (Insecta) alla conoscenza della biogeografia dell’ Anatolia: la componente gondwaniana. Biogeographia 20: 179–200. Google Scholar


P. Lemey , M. Salemi , A.M. Vandamme (Eds) 2009. The Phylogenetic Handbook: a Practical Approach to Phylogenetic Analysis and Hypothesis Testing, 2nd Edition, pp. 381–404, Cambridge University Press, Cambridge. Google Scholar


P. Librado , J. Rozas 2009. DnaSP v5: a software for comprehensive analysis of DNA polymorphism data. Bioinformatics 25: 1451–1452 Google Scholar


F. Médail , P. Quézel 1997. Hot-spots analysis for conservation of plant biodiversity in the Mediterranean basin. Annals of Missouri Botanical Garden 84: 112–127. Google Scholar


F. Médail , K. Diadema 2009. Glacial refugia influence plant diversity patterns in the Mediterranean Basin. Journal of Biogeography 36: 1333–1345. Google Scholar


J. Murienne , I. Karaman , G. Giribet 2010. Explosive evolution of an ancient group of Cyphophthalmi (Arachnida: Opiliones) in the Balkan Peninsula. Journal of Biogeography 37: 90–102. Google Scholar


D. Posada , K.A. Krandal 1998. MODELTEST: testing the model of DNA substitution. Bioinformatics 14: 817–818. Google Scholar


A. Rambout , A.J. Drummond 2003. Tracer v1.3. Available form: Scholar


A. Rambaut 2008. FigTree v1.2. Available from: <>. Google Scholar


M. Raymond , F. Rousset 1995. An exact test for population differentiation. Evolution 49: 1280–1283. Google Scholar


A. Rokas , R.J. Atkinson , M.I. Webster , G. Csoka , G.N. Stone 2003. Out of Anatolia: longitudinal gradients in genetic divesity support an eastern origin for a circum-Mediterranean oak gallwasp Andricus quercustozae. Molecular Ecology 12: 2153–2174. Google Scholar


F. Ronquist , J.P. Huelsenbeck 2003. MRBAYES 3: Bayesian phylogenetic inference under mixed models. Bioinformatics 19: 1572–1574. Google Scholar


H.A. Schmidt , K. Strimmer , M. Vingron , A. von Haeseler 2002. TREE-PUZZLE: maximum likelihood phylogenetic analysis using quartets and parallel computing. Bioinformatics 18: 502–504. Google Scholar


T. Schmitt 2007. Molecular biogeography of Europe: Pleistocene cycles and postglacial trends. Frontiers in Zoology 4: 11–23. Google Scholar


J.M. Seddon , F. Santucci , N.J. Reeve , G.M. Hewitt 2001. DNA footprints of European hedgehogs, Erinaceus europaeus and E. concolor: Pleistocene refugia, postglacial expansion and colonization routes. Molecular Ecology 10: 2187–2198. Google Scholar


C. Simon , F. Frati , A. Beckenbach , B. Crespi , H. Liu , P. Flook 1994. Evolution, weighting, and phylogenetic utility of mitochondrial gene sequences and a compilation of conserved polymerase chain reaction primers. Annals Entomological Society of America 87: 651–702. Google Scholar


P. R. Stephens , J.J. Wiens 2009. Bridging the gap between community ecology and historical biogeography: niche conservatism and community structure in emydid turtles. Molecular Ecology 18: 4664–4679. Google Scholar


G.N. Stone , R.J. Challis , R. J. Atkinson , G. Csoka , A. Hayward , G. Melika , S. Mutun , S. Preuss , A. Rokas , E. Sadeghi , K. Schorogge 2007. The phylogeographical clade trade: tracing the impact of human-mediated dispersal on the colonization of northern Europe by the oak gallwasp Andricus kollari. Molecular Ecology 16: 2768–2781. Google Scholar


K. Strimmer , A. von Haeseler 1997. Likelihood mapping: a simple method to visualize phylogenetic content of a sequence alignment. Proceedings National Academy Science, USA 94: 6815–6819. Google Scholar


D.L. Swofford 2002. PAUP* Phylogenetic Analysis Using Parsimony (* and other methods). v. 4.0 beta. M.A. Sinauer Associates, Sunderland. Google Scholar


G.A. Sword , L.B. Senior , J. Gaskin , A. Joern 2007. Double trouble for grasshopper molecular systematics: intra-individual heterogeneity of both mitochondrial 12S-valine–16S and nuclear internal transcribed spacer ribosomal DNA sequences in Hesperotettix viridis (Orthoptera: Acrididae). Systematic Entomology 32: 420–428. Google Scholar


P. Taberlet , L. Fumagalli , A.G. Wust-Saucy , J.F. Cosson 1998. Comparative phylogeography and postglacial colonization routes in Europe. Molecular Ecology 7: 453–464. Google Scholar


D. Tarkhnishvili R.S. Thorpe , J.W. Arntzen 2000. Pre-pleistocene refugia and differentiation between populations of the Caucasian salamander (Mertensiella caucasica). Molecular Phylogenetics and Evolution 14: 414–422. Google Scholar


B.S. Weir , C.C. Cockerham 1984. Estimating F-statistics for the analysis of population structure. Evolution 38: 1358–1370. Google Scholar


M. Veith , J.F. Schmidtler , J. Kosuch , I. Baran , A. Seitz 2003. Palaeoclimatic changes explain Anatolian mountain frog evolution: a test for alternating vicariance and dispersal events. Molecular Ecology 12: 185–199. Google Scholar


M. Veith , E. Lipscher , M. Oz , A. Kiefer , I. Baran , R.M. Polymeni , S. Steinfartz , 2008. Cracking the nut: geographical adjacency of sister taxa supports vicariance in a polytomic salamander clade in the absence of node support. Molecular Phylogenetics and Evolution, 47: 916–931. Google Scholar


L. Zhao , J. Zhang , Z.J. Liu , S.M. Funk , F.W. Wei , M.Q. Xu , M. Li 2008. Complex population genetic and demographic history of the salangid, Neosalanx taihuensis, based on cytochrome b sequences. BMC Evolutionary Biology 8: 201–228. Google Scholar
Battal Çiplak, Sarp Kaya, and İslam Gündüz "Phylogeography of Anterastes serbicus Species Group (Orthoptera, Tettigoniidae): Phylogroups Correlate with Mountain Belts, but not with the Morphospecies," Journal of Orthoptera Research 19(1), 89-100, (1 July 2010).
Received: 9 April 2010; Accepted: 1 May 2010; Published: 1 July 2010

Get copyright permission
Back to Top