Introduction

The Caucasus and neighboring Iran and Turkey are situated between the Levant and Europe, and can be considered as one of the potential pathways during the colonization of Europe by modern humans about 40 000 kya ago.1 A number of Early Upper Paleolithic sites, dating back to about 30–40 kya, have been excavated in this region suggesting an early presence of modern humans in the region.1, 2, 3, 4, 5 Importantly, the Upper Palaeolithic sites were found only in the mountain ridges along the passage between the Caucasus Mountains and the Black Sea that connects the South Caucasus with the southern part of East Europe.1 Thus, genetic study of this region can provide some insights into ancient migrations and can facilitate a reconstruction of modern human migration routes from the Levant to Europe.

Several genetic studies of a number of different groups from this region have been carried out in the last decade.6, 7, 8, 9, 10, 11, 12, 13, 14, 15 These studies were mostly based on sequences of the first hypervariable segment of the noncoding mtDNA control region (HV1) and Y chromosome binary marker variation in various groups from this region. mtDNA studies reveal a high level of diversity, exceeding that within all of Europe and only slightly lower than West Asian mtDNA diversity, which might indicate an old age of human populations from this region.10 Overall, the Caucasus groups showed greater similarity with West Asian than with European groups for both genetic systems, although this similarity was much more pronounced for the Y chromosome than for mtDNA, suggesting that recent male-mediated migrations from West Asia have influenced the genetic structure of Caucasus populations.11

Previous studies of complete mtDNA genome sequences have in general first obtained HV1 sequences (and sometimes, genotyped some haplogroup-defining SNPs), and then selected particular individuals of interest for complete mtDNA genome sequencing. However, such biased sampling renders the resulting sequences unsuitable for many demographic analyses. We have therefore used a high-throughput sequencing approach to generate 147 complete mtDNA genome sequences from random samples of individuals from three groups from the Caucasus (Armenians, Azeri and Georgians), and one group each from Iran and Turkey. We used these sequences to investigate the genetic structure of these groups and more accurately infer their demographic history. In particular, we have examined population size changes through time via Bayesian Skyline Plots (BSPs), and we have used an approximate Bayesian computation (ABC) approach to estimate divergence times between groups; neither of these analyses can be applied to HV1 sequences alone because the HV1 sequences lack sufficient information for such analyses. Our results amply show the value of random sampling of complete mtDNA genome sequences that is afforded by this high-throughput sequencing approach.

Materials and methods

Samples and DNA extraction

Samples were obtained from unrelated individuals, representing three populations from the Caucasus: Armenians (30 cheek cell swabs from Erevan), Azeri (30 cheek cell swabs from Baku) and Georgians (28 saliva samples from Batumi). Samples were also obtained from Turks (29 saliva samples from Ankara) and Iranians (30 blood samples from Tehran). A map of the sampling locations is shown in Figure 1. Genomic DNA from cheek cell swabs was extracted using a standard salting-out procedure.16 DNA from saliva samples was extracted using a previously described method.17 Blood samples were processed using the QIAamp DNA Blood Mini Kit (Qiagen, GmbH, Germany), following the instructions of the manufacturer. Informed consent and information about birthplace, parents and grandparents was obtained from all donors.

Figure 1
figure 1

Geographic location of sampling sites.

Sequencing complete mtDNA genomes

The entire mtDNA genome was amplified in two overlapping products of about 9.7 and 7.3 kb, using primer pairs L12279/H2986 and L2603/H12314.18 Long-range PCR was carried out using the Expand Long-Range dNTP pack (Roche, Berlin, Germany) and 3 ng of template DNA in a 50 μl volume, using the protocol provided by the manufacturer. The annealing temperature was 57 °C. PCR products were purified with SPRI beads (Agencourt, Denver, MA, USA) using the manufactureŕs instructions. The two PCR products for each individual were mixed in equimolar ratios and nebulized using a Bioruptor Sonicator UCD-200 (Diagenode, Liege, Belgium). The size range of the sonicated DNA fragments was between 150 and 450 bp. MinElute spin columns (Qiagen) were used to remove small DNA fragments, and the purified, nebulized DNA was eluted in 20 μl of elution buffer. About 400 ng of DNA was used for tagging the nebulized PCR product for each individual with a specific tag sequence, as described elsewhere,19 and Illumina Genome Analyzer II (Illumina Inc., San Diego, CA, USA) libraries were prepared according to the protocol described elsewhere.20 Fifty individuals with unique tags were pooled in each library in equimolar ratios. Subsequently, libraries were sequenced with 36 and/or 76-bp reads in one lane of an Illumina flow cell (Cluster Generation kit V2, FC-103–300 × sequencing chemistry, San Diego, CA, USA) according to the manufacturer's instructions.

mtDNA sequence assembly

Each run was processed with RTA 1.5 (Illumina Inc.). Afterward, the PhiX174 control reads of a dedicated control lane were aligned to the corresponding reference sequence to obtain a training data set for the alternative base caller Ibis.21 Raw sequences called from Ibis were separated by sample using their index read, allowing one mismatch and the loss of the first base.20 Reads were then filtered for sequence quality and complexity. In this step, reads having >5 bases with a quality score below 10 (PHRED scale) and reads with sequence entropy (calculated by summing –p*log2(p) for each of the four nucleotides) below 1.0 were removed. Sequence reads lacking an expected tag or with more than one tag were removed. After untagging, reads with identical start and end positions were also removed, since they may represent clones of a single sequence. Sequence reads that passed these filters were sorted according to unique sequence tags. Consensus mtDNA sequences were assembled using the software MIA22 by mapping reads to the revised Cambridge reference sequence (rCRS).23 A multiple alignment of the consensus sequences was carried out using the software MAFFT v6.708b.24 The mtDNA genome sequences have been deposited into Genbank (accession numbers HM852756–HM852902). The list of polymorphic sites and undetermined two positions is shown in the Supplementary Table 1.

Haplogroup assignment

Sequences were assigned to haplogroups according to Phylotree.org Build 6,25 using a custom PERL script. Sequences were assigned to the closest matching haplogroup. As in Phylotree, positions 309.1C(C), 16182C, 16183C, 16193.1C(C) and 16519 were not used to assign haplogroups, because these are highly mutable positions.

Data analysis

DnaSP 4.026 was used to calculate basic parameters of genetic diversity. Analyses of molecular variance (AMOVA) were carried out using Arlequin.27 The statistical significance of Fst values was estimated by permutation analysis, using 10 000 permutations. The software package GENESYN28 was used to calculate the number of polymorphic sites, substitutions, and nonsynonymous and synonymous substitutions. The STATISTICA package (StatSoft Inc., Tulsa, OK, USA) was used for multi-dimensional scaling (MDS) analysis.29

BSPs were produced from the coding region sequences (positions 577–16023) using MCMC sampling in version 5.1 of the program BEAST.30 The plots were obtained with a piecewise linear model and ancestral gene trees were based on a general time-reversible substitution model with invariant sites (GTR+I). A Bayes factor computed via importance sampling indicated that the strict molecular clock could not be rejected and was therefore used for the analysis. We allowed 20 discrete changes in the population history, using a coalescent-based tree prior with a linear model in which population size grows and declines between changing points. Each MCMC sample was based on a run of 40 000 000 generations sampled every 4000 steps, with the first 4 000 000 generations regarded as burn-in. Three independent runs were made for each population, and a mutation rate of 1.691 × 10−8 was used.31

To estimate the divergence time between each pair of populations, we implemented the ABC approach described previously,32 with the exception that all data were simulated with the coalescent simulator ms. Details of the rejection-regression algorithm are as described previously.32 This basically involves fitting a local-linear regression of all simulated parameter values to simulated summary statistics. Furthermore, the observed summary statistics are then substituted into a regression equation. All parameters were transformed with logtan before the actual regression analysis.33 Distances between observed and simulated summary statistics were calculated as Euclidean distances. We simulated 1 million sequence data sets for each estimation procedure. Each data set was simulated as a 15 kb region with a mutation rate of 1.69 × 10−8 substitutions per site per year.31 The recombination rate was set to 0. The demographic history underlying each population was a consensus history from all five populations, obtained from the BSP results: an ancestral population of an effective size of 400 experienced a sudden growth 45 000 years ago to an Ne of 4000, and after a long period of constant size a second sudden growth 15 000 years ago increased the size of the population to 40 000. A simple 1– parameter model of a single population split between two populations with no migration was assumed, and the time of divergence TDIV was the parameter to be estimated. The prior for this parameter was TDIVU[1002000] in generations. Each estimation was done as a pairwise population comparison, resulting in 10 ABC estimations in total. Out of the 1 million simulations for each ABC estimation, the 10 000 simulations with the smallest Euclidean distances were kept and used to estimate the posterior distribution. The Euclidean distances for the retained best simulation ranged from 0.2 to 1.2. We used six summary statistics in total: S (the number of segregating sites for each population in every pairwise comparison); π (the average number of pairwise differences for each population in every pairwise comparison); the percentage of shared polymorphisms per pairwise population comparison; and the pairwise Fst value.

Results

After quality filtering and removal of potential duplicate reads, there were an average of 33 219 reads per individual. The average coverage was 69.3 for 90 samples sequenced with 36-bp reads (Supplementary Figure 2A); 211.1 for 21 samples sequenced with 76-bp reads (Supplementary Figure 2B) and 102.5 for 36 samples sequenced twice with 36-bp reads, because of poor coverage in the first sequencing (Supplementary Figure 2C). Overall, the average coverage of each position was 87-fold (Supplementary Table 2).

Summary statistics describing genetic diversity are shown in Table 1. The number of mean pairwise differences (MPDs) ranged from 31.8 to 33.5 for all groups except the Azeri, for which the MPD is 39.3. All groups exhibit similar nucleotide diversity values, as well as an excess of low-frequency variants that is characteristic of a recent population expansion, as shown by significantly negative values for Tajima's D (Table 1).

Table 1 Summary statistics for Armenians, Azeri, Georgians, Iranians and Turks based on complete mtDNA genome sequences

A total of 855 polymorphic sites were detected (Table 2). The highest number of variable sites was found in the control region (192 sites), which is more than three times higher than expected (57 sites) based on the length of the control region. The number of transitions significantly exceeds the number of transversions for all parts of the mtDNA genome (Table 2). The ratio of transitions vs transversions was higher for protein-coding genes (16.71) than for the control region (4.49), but not significantly so (P=0.387). The number of synonymous substitutions was higher than nonsynonymous substitutions for all protein-coding genes except ATP6 and ATP8 (Table 2). Among the 13 protein-coding genes, there were 164 sites with nonsynonymous changes and 385 sites with synonymous changes (Table 2). The ratio of the number of polymorphic nonsynonymous sites per total nonsynonymous sites to the number of polymorphic synonymous sites per total synonymous sites (pn/ps values), is less than one in all cases, but elevated for ATP6 and ATP8 (Table 2).

Table 2 Number of variable sites, transitions, transversions, synonymous and nonsynonymous differences, and pa/ps ratio

A plot of the number of differences in the HV1 sequences vs the number of differences in the coding region sequences for each pair of individuals showed that the pairwise comparisons with no differences in the HV1 sequences goes up to a maximum of 34 differences in the coding region (Figure 2). Thus, there can be appreciable coding region variation among individuals with identical HV1 sequences.

Figure 2
figure 2

Plot of the number of differences in the HV1 sequences vs the number of differences in the coding region sequences for each pair of individuals. The best-fit regression line is indicated.

The haplogroup assignment for each individual, according to the nomenclature of Phylotree.org,23 is given in Supplementary Table 1. As none of the sequences in this study matched any previous haplogroup exactly, a sequence was assigned to a haplogroup if it contained all the mutations that define that haplogroup. All of our sequences therefore contain all of the mutations that define the haplogroup plus additional mutations.

A total of 97 different haplogroups were identified (Supplementary Table 1), which fall into 16 major haplogroups (Table 3). Haplogroups H, HV, J, T and U were found in all groups, and are generally among the most frequent West Eurasian mtDNA haplogroups. Several haplogroups typical for Central and East Asia (A, C, D and F) were found only in Azeri and Turks (Table 3). One haplogroup (X) was restricted to the three Caucasus groups.

Table 3 Haplogroup frequencies

In order to visualize the relationships of these five groups based on complete mtDNA genome sequences, an MDS plot was constructed from the pairwise Fst values (Supplementary Table 3). The results (Figure 3a) showed a cluster of groups from the Caucasus in the first two dimensions, but a similar placement of Azeri and Turks in the third dimension For comparison, we constructed an MDS plot for just the HV1 region (Figure 3b), which showed a similar cluster of the Caucasian groups in the first two dimensions, but did not associate Azeri and Turks in the third dimension. Thus, the complete mtDNA sequences reveal an association between the two Turkic-speaking groups that is not seen in the HV1 sequences.

Figure 3
figure 3

MDS plot for five groups from the South Caucasus, Iran and Turkey. (a) Based on complete mtDNA genome sequences. (b) Based on HV1 sequences only.

The geographic and linguistic structure of the Caucasian, Iranian and Turkic groups was investigated by an AMOVA (Supplementary Table 4). Approximately 99% of the variance was due to the within-population component. A geographic classification of populations gave a slightly better fit to the genetic data than did a linguistic classification, although permutation tests showed that the higher among-group components are not significantly different from zero.

Many of the above descriptive analyses can be carried out on HV1 sequences as well as complete mtDNA genome sequences. However, HV1 sequences are not amenable to demographic inference, and we utilized our random samples of complete mtDNA genome sequences to estimate population size change over time and pairwise population divergence times. To investigate changes in population size over time, BSPs were constructed using the coding region (positions 577–16022). The BSPs for the three Caucasus groups as well as the Iranian group (Figure 4) are generally similar, showing a population expansion around 45–50 kya, followed by a constant population size, and then another expansion around 15–18 kya. However, the BSP for the Turkic group differs, with an increase from 35 to 50 kya followed by a prolonged period of constant population size, and no indication of a second growth period.

Figure 4
figure 4

Bayesian skyline plots. The Y axis for each plot is the effective population size and the X axis is time in years. (a) Georgians; (b) Azeri; (c) Armenians; (d) Iranians and (e) Turks. The gray lines represent the 95% posterior density intervals.

An ABC approach was used to estimate the divergence time between each pair of populations. For each estimation, 1 million simulations were carried out and compared with empirical data based on six summary statistics. The posterior distributions based on the best-fitting 10 000 simulations show pronounced peaks (Supplementary Figure 1), and the Euclidean distances between the simulated and observed parameters vary between 0 and 1.2, indicating good support for the divergence time estimates (Table 4). In general, the oldest divergence times are between the group from Turkey and other groups from the South Caucasus and Iran (400–600 generations). The divergence time for the South Caucasus groups from each other is about the same as their divergence from Iran (average Iran-Caucasus is 361 generations, average within Caucasus is 365 generations).

Table 4 Divergence time (Tdiv) (in generations) between pairs of populatiuons

Discussion

Overall, the randomly sampled complete mtDNA genome sequences indicated extraordinarily high genetic diversity in the groups from the South Caucasus, Iran and Turkey. For the Georgians, Armenians and Iranians, all individuals had different sequences, while for the Azeri two individuals shared the same sequence, and for the Turks there were 27 haplotypes among 29 individuals (Table 1). Overall, there were 144 different sequences among the 147 individuals. By contrast, when only the HV1 sequences are considered, the level of haplotype sharing is much higher (Table 1), indicating that individuals with identical HV1 sequences often had different complete mtDNA genome sequences.

The assignment of sequences to haplogroups further reinforces the extraordinary diversity in the complete mtDNA genome sequences, as none of the 147 sequences in this study was identical to a sequence in the Phylotree.org database.23 This is all the more remarkable when one considers that the complete mtDNA genome sequences present in the public databases are heavily biased toward Eurasia; for example, nearly 40% of >5000 sequences analyzed by Pereira et al34 are from Europe and the Middle East. Thus, even in well-studied areas of the world, there is much mtDNA diversity yet to be discovered.

The high level of diversity extends to the level of major haplogroups, as the individual haplogroups fall into 16 major haplogroups (Table 3). Most of these haplogroups (H, HV, I, J, K, M, N, R, T, U and Z) are typical for West Eurasia,34 while others (A, C, D and F) occur mostly in Central and East Asia.34 Haplogroup X is widespread, albeit at low frequencies, across Eurasia.34 Among the West Eurasian haplogroups, the Caucasian groups are characterized by relatively high frequencies of haplogroups U and X. The Central/East Asian haplogroups were notable in that they were found only in a few individuals from the Azeri and Turkish groups (the two Turkish-speaking groups in the study), suggesting some Central Asian influence specifically on these groups. Indeed, there is historical documentation of such contact via the Oguz migrations from Central Asia to Anatolia and the South Caucasus in the eleventh century.35 These groups brought the Turkic language into the Azeri and Turkish populations, and presumably left some genetic footprint along with their language. Although the specific contact(s) between Azeri and Turks with Central Asian groups that brought in these Central Asian mtDNA lineages is unknown, overall the low frequency of these mtDNA lineages is in good agreement with previous estimates of low levels of gene flow from Asia into Anatolia.36

Overall, the complete mtDNA genome sequences indicate that the three South Caucasus groups are genetically similar, although they represent three different language families. This can be seen in the MDS plots (Figures 3a and b), where Indo European-speaking Armenians and Turkic-speaking Azeri cluster with the Caucasian-speaking Georgians. The clustering of groups on the basis of geographic rather than linguistic relationships is in keeping with previous studies of genetic diversity among Caucasian populations.10, 11 However, the complete mtDNA genome sequences do reveal some additional genetic similarity between the two Turkic-speaking groups (Azeri and Turks) that was not evident in previous studies. Presumably, the sharing of some Central Asian mtDNA haplogroups by Azeri and Turks (Table 3) accounts for this signal of genetic similarity between the Azeri and the Turks. Thus, the greater genetic resolution afforded by the complete mtDNA genome sequences both confirms and extends previous studies.

We also investigated patterns of variation in the mtDNA genome itself (Table 2). Of the 855 variable positions, about 22% occurred in the control region, which is significantly more than expected if variable positions occur randomly across the mtDNA genome (P<0.0001). The excess in polymorphisms and transversions for the control region most likely reflects weaker functional constraints on this major noncoding region of the mtDNA genome. Among the 13 protein-coding genes, there was a significant excess of nonsynonymous polymorphisms in the overlapping ATP6 and ATP8 genes, relative to the other mtDNA protein-coding genes (Table 2). Previously, an excess of nonsynonymous polymorphisms in the ATP6 gene was found in Siberian populations and hypothesized to reflect positive selection for cold adaptation.37 However, subsequent studies suggested that relaxation of functional constraints, rather than positive selection, can explain the higher pa/ps ratio for ATP6.38 The finding of a significantly elevated pa/ps ratio for ATP6 in this study, as well as in a previous study of Filipino groups,39 argue against the cold adaptation hypothesis.

In the interval of ca 30–40 kya, Europe underwent numerous changes known as the Upper Paleolithic revolution, which involved dramatic changes in technologies, hunting techniques, human burials and an artistic traditions revealed in the archeological records.40 These changes presumably also involved population size increases, and indeed the BSPs for all five groups indicate a strong expansion around this time (Figure 4). Following this initial expansion, the BSPs indicate constant population size during the period of the Last Glacial Maximum (LGM), about 18–30 kya. The dates for the second expansion event evident in the BSPs from the Caucaus and Iranian groups are in good agreement with the start of the continuous deglaciation of the ice shield 18 kya ago, when environmental conditions were almost fully glacial.41 Curiously, the BSP for the mtDNA sequences from Turkey do not show any evidence of this second expansion (Figure 4). The different BSP for Turkey cannot be explained by Central Asian mtDNA sequences in this population, as a BSP with the Central Asian sequences removed is identical to the BSP for all of the sequences (data not shown). This suggests that the ancestors of the group from Turkey have a different history than the ancestors of the Caucasian and Iranian groups in this study. Specifically, these results suggest that the ancestors of the group from Turkey did not expand after the LGM. This could be explained by the fact that Turkey was not influenced heavily by glaciation during the LGM. The most extensive LGM glaciers descended only to an altitude of 2150 m above sea level in central Turkey,42 suggesting that the end of the LGM may not have caused as dramatic changes in the lowland environment as occurred in the Caucasus.

In addition, we estimated the divergence time between pairs of populations using an ABC–MCMC approach with uniform priors.32 Assuming a generation time of 28 years for mtDNA,43 the earliest divergence occurred between the group from Turkey and the other four groups about 11.2–16.8 kya. This event coincides with the second expansion event for the South Caucasus and Iranian groups (Figure 4). This was a time of post-glacial recolonization; in the wake of climatic amelioration, temperate species expanded their distribution range to the north, following the expansion of favorable habitats (‘habitat tracking’).44 Most likely, modern humans followed the expansion of temperate species and recolonized the South Caucasus, which was mostly glaciated during the LGM.45 Along with the colonization of new lands and territories, human groups settled in different parts of the region, giving rise to new populations. This scenario is supported by divergence time estimates among South Caucasus groups, ranging from 5.6 to 11 kya, as well as divergence time estimates between Iran and South Caucasian groups of 8.9–10 kya. Thus, in this part of West Eurasia there appear to have been two major periods of population expansion and divergence after the LGM.

Admittedly, some caveats are in order. Divergence time estimates were obtained assuming no subsequent migration between groups, and hence should be viewed as minimum estimates of the actual divergence time. We did attempt to incorporate migration into the model, however, we were unable to obtain reliable estimates of both migration rate and divergence time with the ABC–MCMC approach (results not shown). Nevertheless, the pronounced peaks in the posterior distributions (Supplementary Figure 1), small Euclidean distances, and concordance between the BSPs and the divergence time estimates indicate that the demographic inferences are probably reliable. Moreover, the relatively low level of sequence sharing between groups further indicates that recent migration among these groups has been low.

In conclusion, the finding of extraordinarily high mtDNA diversity and resulting demographic inferences for the South Caucasian and West Asian groups studied here was enabled by random sampling of individuals, which in turn was made possible by the high-throughput sequencing approach implemented here. The speed, low cost and reliability of the resulting sequences demonstrated by this and similar studies39 further indicate that this is the approach of choice for generating complete mtDNA genome sequences.