Deep Phylogenetic Analysis of Haplogroup G1 Provides Estimates of SNP and STR Mutation Rates on the Human Y-Chromosome and Reveals Migrations of Iranic Speakers

Y-chromosomal haplogroup G1 is a minor component of the overall gene pool of South-West and Central Asia but reaches up to 80% frequency in some populations scattered within this area. We have genotyped the G1-defining marker M285 in 27 Eurasian populations (n= 5,346), analyzed 367 M285-positive samples using 17 Y-STRs, and sequenced ~11 Mb of the Y-chromosome in 20 of these samples to an average coverage of 67X. This allowed detailed phylogenetic reconstruction. We identified five branches, all with high geographical specificity: G1-L1323 in Kazakhs, the closely related G1-GG1 in Mongols, G1-GG265 in Armenians and its distant brother clade G1-GG162 in Bashkirs, and G1-GG362 in West Indians. The haplotype diversity, which decreased from West Iran to Central Asia, allows us to hypothesize that this rare haplogroup could have been carried by the expansion of Iranic speakers northwards to the Eurasian steppe and via founder effects became a predominant genetic component of some populations, including the Argyn tribe of the Kazakhs. The remarkable agreement between genetic and genealogical trees of Argyns allowed us to calibrate the molecular clock using a historical date (1405 AD) of the most recent common genealogical ancestor. The mutation rate for Y-chromosomal sequence data obtained was 0.78×10-9 per bp per year, falling within the range of published rates. The mutation rate for Y-chromosomal STRs was 0.0022 per locus per generation, very close to the so-called genealogical rate. The “clan-based” approach to estimating the mutation rate provides a third, middle way between direct farther-to-son comparisons and using archeologically known migrations, whose dates are subject to revision and of uncertain relationship to genetic events.


Introduction
Despite multiple studies of the phylogeography of individual Y-chromosomal haplogroups, haplogroup G1-M285 has not received attention so far. This is partly explained by its relatively low frequency in its main area of distribution in South-West Asia [10,42], and partly by its uneven geographic distribution with a maximum frequency in the Madjar population in Kazakhstan [5]. For this reason, study of the phylogeography of haplogroup G [44] dealt mainly with the G2 sub-branch, and the only statement about G1 is an estimate of its age from Y-STR markers (19,000 ± 6,000 years). However, newly accumulated data indicate that G1 is present over a wider area in the Eurasian steppe than in Madjars only [10], and it also reaches very high frequencies in geographically distant populations of the Armenian plateau (Table 1). Thus, haplogroup G1 might mark an ancient genetic link between Iranic speakers of South-West Asia and populations of the Central Asian steppes where Iranian speech predominated in the second and first millennia BC ( Fig 1A). However, the place of origin of this haplogroup remains unclear, and it is unknown whether South-West Asians and Madjars have the same or different subbranches of haplogroup G1, what the age of the branch(es) are, and which ancient migrations contributed to the contemporary distribution and diversity of this haplogroup.
These details of haplogroup G1 phylogeography have been hard to answer, because existing methods allowed only slow progress in discovering phylogenetically informative SNPs. Fortunately, during recent years the possibility for full resequencing of the Y-chromosome [17,41,43,49,50], and more particularly the Y-capture technologies which became commercially available in the year 2013, stimulated intensive discovery of phylogenetically informative SNPs. For example, during the last decade (from the first extensive papers in 2000 till 2011) only 485 SNPs were placed on the global Y-chromosomal phylogenetic tree, while in the three following years the number of SNPs has exceed 9,000 (www.isogg.org).
Within the last decade, there has been significant uncertainty in dating Y-chromosomal haplogroups due to a three-fold difference between so-called "genealogical" and "evolutionary" mutation rates of Y-STRs. The former rates were repeatedly obtained in a set of studies [18,22,46] comparing father-son pairs, while the latter was obtained in single study [54] where calibration was done using population events with known historical dates. Increasing datasets of complete Y-chromosomal sequences allowed new calculations of the mutation rates, this time focused on SNPs. Four mutation rates have been suggested so far, ranging from 0.6 to 1.0 ×10 -9 per bp per year: the pedigree-based rate [50], calibrations based on peopling of the Americas [41] and Sardinia [17], and the rate adopted from the pedigree rate for autosomal SNPs [37]. The two-fold difference between these rates makes further estimations necessary. In the current study we had the chance to calibrate the Y-chromosomal molecular clock using a historically reliable date of the most recent common genealogical ancestor of carriers of haplogroup G1 in Kazakh clans. Migration of Iranic-speaking populations between the Central Asian steppes and South-West Asian uplands is an important issue in human population history, directly related to the much-debated problem of the homeland and early migrations of Indo-Europeans. Followers of the Kurgan theory propose that the carriers of Iranic languages expanded from the Eurasian steppe southward to present-day Iran, from which region these languages received their name ( Fig 1B). The competing theory locating the Indo-European homeland in Eastern Anatolia proposes that the Iranic branch migrated from the Iranian plateau northward to the steppes ( Fig 1B). Thus, both theories agree on the area populated by ancient Iranic-speakers (both the Iranian-Armenian plateau and Central Asia steppes) and later replacement of Iranic languages in the steppes by the Turkic ones. But these theories suggested opposite directions of the population movements between the steppes and uplands [34]. This study presents a deep phylogeographic analysis of haplogroup G1 by combining traditional approaches with the new powerful options emerging from complete sequencing of the Y-chromosome. We set out to provide a new independent estimate of the mutation rate using the tight links between haplogroups and clans typical in patrilineal nomadic societies. In addition, we aimed to find which direction of the ancient migration of Iranic speakers better fits the haplogroup G1 phylogenetic pattern.

Methods Genotyping
We genotyped the commonly-used SNP M285 which defines haplogroup G1 (YCC, 2002) in multiple Eurasian populations using the TaqMan technique (Applied Biosystems) and identified 367 M285-derived samples in 27 populations. All these samples were then genotyped at 17 Y-chromosomal STRs using the Y-filer genotyping kit (Applied Biosystems). All sample donors gave their written informed consent (the study was approved by the Ethics Committee of the Research Centre for Medical Genetics, Russian Academy of Medical Sciences). Data available from the literature were also incorporated ( Table 1, S1 Table).
Then we selected 19 samples for high-throughput sequencing of the Y-chromosome. To capture maximum phylogenetic diversity and thus increase the cost-effectiveness of the analyses, we applied three criteria for selecting samples. The geographic criterion led to samples from both steppe and mountain parts of the haplogroup's area being included, particularly from populations where G1 frequency is high. The phylogenetic criterion led to samples from all clusters revealed on the STR network being included and represented by at least two samples for full sequencing, because STR-clusters might reflect real phylogenetic branches and a single sample would not allow us to distinguish phylogenetically-informative SNPs from private ones. The third criterion could be applied only to those populations where paternal clan structure is present: it led to representatives from different clans being included because members of the same clan have a high probability of sharing almost identical paternal lineages. As an outgroup for the 19 G1 samples, we also sequenced one sample from its brother haplogroup, G2.
Y-chromosomal genotyping was performed using a custom enrichment design created for the commercially available "BigY" product offered by Gene By Gene, Ltd. In total, the target regions attempt to sequence around 20 million base pairs with 67,000 capture probes, on the Illumina HiSeq platform. This design captured 11,383,697 bp within the non-recombining malespecific Y-chromosome, consistent with regions genotyped by previous Y sequencing studies [41] and the Y positions placed on the phylogenetic tree by the Y Chromosome Consortium [28]. Following BigY sequencing, and also as part of the product, downstream software analysis was performed using the Arpeggi Engine (AEngine) pipeline. This includes short read mapping, alignment post-processing, and variant calling. For quality control purposes, BigY samples are monitored for read totals, average coverage and average base quality, and should a sample fall below BigY standard thresholds, the sample is re-run. A regions file listing the genomic build 37 capture targets of BigY can be found at https://www.familytreedna.com/ documents/bigy_targets.txt. Variants found across the samples are classified as any deviation from the reference genome, and in addition, we reported genotypes for about 37,000 phylogenetically informative SNPs in the FamilyTreeDNA database (www.familytreedna.com).
In addition to the genotyping per sample, we wanted to ensure for this study that SNP positions examined were adequately covered across all samples. This is a concern, because many variant calling methods in high-throughput sequencing are ambiguous when not reporting a variant as to whether there was not enough coverage to genotype, or if there was a legitimate homozygous reference genotype. To discern such cases, each BigY sample was given a "confidence" region list determined by genotype quality scores for each base. The genotype quality is computed as the probability that the genotype is correct, according to a phred score. This probability is derived from AEngine's proprietary statistical model considering characteristics of read coverage, individual read mapping qualities, and base sequencing quality scored by the HiSeq. A base position is appended to the confidence regions for that sample if its genotype quality score is above 3.02. Thus, if there is no variant occurring at a base within confidence intervals for a sample, it can be assumed that the sample is reference genotype at that position. Variant calls were produced and handled as Variant Call Format (VCF) files, according to the established field standards (http://samtools.github.io/hts-specs/VCFv4.1.pdf). To this effect, the intersection of confident regions covered by the 20 samples studied was also recorded, and can be found within the Supplementary Data. More details on the BigY capture sequencing method are available at https://www.familytreedna.com/learn/wp-content/uploads/2014/08/ BIG_Y_WhitePager.pdf, which features the methods and the capture probe list and target regions.
To estimate the potential sequencing error rate, we applied a phylogenetic approach. We checked whether we found all SNPs in the BigY captured region which are known to be phylogenetically located between haplogroups A0 and G (www.isogg.org) and thus should be present in our samples. The proportion of missed SNPs was the false negative rate. We also checked whether we see SNPs known to define other haplogroups, which are therefore not expected to be present in our haplogroup G samples. The proportion of these unexpected SNPs was considered as the false positive rate. We note that this approach overestimates the error rate, because it considers parallel mutations as errors and ignores potential inaccuracies in identifying the SNP ancestral states in ISOGG database. The (over)estimated rates were 0.008 for false negatives and 0.005 for false positives. (See details in S2 Table).

Analyses based on G1 frequencies and STR-haplotypes
The frequency distribution map of haplogroup G1-M285 was created using data reported here for the first time (27 populations, Table 1), data from the literature (33 populations, Table 1) and published data on other 266 Eurasian populations where G1 frequency was zero. The map was created in the GeneGeo software as described previously [2,32] setting the weight function to 2 and radius of influence to 2500 km. The map is presented at two scales: a specific scale highlighting the distribution pattern of this haplogroup (Fig 2) and the "universal" scale routinely used in GeneGeo for mapping haplogroup frequencies (S1 Fig).
An analysis of molecular variance (AMOVA) was performed using Arlequin [13] on two groups of populations: those from the ancient area of Iranic speakers, compared to the group of all other Eurasian populations. We calculated variation among these two groups of populations using data on each haplogroup separately and identified haplogroups demonstrating highest differentiation between "Iranic" and "non-Iranic" populations (S3 Table).
Reduced median networks [4] of STR haplotypes (S1 Table) were created in the programs Network 4.6.0.0 and Network Publisher (Fluxus-Engineering, http://www.fluxus-engineering. com). We applied the preprocessing star contraction algorithm [16] and postprocessing Steiner maximum parsimony algorithm [40]. The allele sizes for locus DYS389II were determined with the subtraction of DYS389I. Loci DYS385a and DYS385b were excluded from network analyses. The main network was based on 15 STRs genotyped in 386 samples (S1 Table). To include data from the Madjar subclan of the Kazakhs [5], the second network restricted to the Argyn Kazakh population was based on 10 STRs (S1 Table). Haplotype diversity was calculated according to [38] as HD ¼ 1 À X p 2 i , where p i is the frequency of the i th haplotype. Data on 17 Y-STR haplotypes (DYS389II alleles were subtracted) belonging to G1-haplogroup from S1 Table were used. Neighbouring populations with small numbers of G1 haplotypes were pooled to reach a minimum sample size of six; the average sample size was 30. The values of haplotype diversity were plotted on a map using the Gen-eGeo software with the weight function set to 4.

Phylogenetic analysis of the Y-chromosomal sequence data
The BigY output VCF files (S1 Data) contained around 33,700-35,900 SNP calls in each of 20 samples. We combined these datasets into single table and filtered out (i) indels, (ii) SNPs with a call rate below 95% (i.e. not called in at least one out of 20 samples; BED files indicating called ranges for each sample are present in the S1 Data) and (iii) SNPs which demonstrated no polymorphism in our samples (i.e. all samples were either identical to or all were different from the reference at these positions). The resulted filtered dataset (S4 Table) consisted of 19 G1 samples, one G2 outgroup sample and 636 SNP positions with very little homoplasy.
The parsimony trees were constructed from this dataset using TNT [20] and Phylomurka (http://phylomurka.sourceforge.net) software. Only one optimal topology was obtained although the states of internal nodes could be marked in different ways. S5 Table presents the ages of branches estimated according to [45].
The same dataset was also subjected to analysis with BEAST software [11] which can reconstruct phylogeny and estimate divergence time by a number of Markov chain Monte Carlo methods. For the test we chose the GTR nucleotide substitution model and Gamma-distributed site heterogeneity with default parameters. We tested both strict and lognormal relaxed clock models and finally preferred the latter due to a positive posterior value of the rate variance between tree branches (S5 Table). For the tree prior we chose the Expansion Growth model assuming that the population grew exponentially since a relatively recent time. The prior for the mutation rate was set as a uniformly distributed value, initially equal to 1.2×10 -5 per SNP (S5 Table) while the age of Kazakh cluster was forced to be normally distributed with the mean of 627 years (the value obtained using genealogical records, see below) and standard error of 50 years. Sufficient ESS values were achieved with the MCMC chain size of 20,000,000 and higher. The consensus of 10,000 trees produced by BEAST is the same as our parsimony tree, and Bayesian age estimates show less than 20% difference from those obtained with Rho statistics (S5 Table). Bayesian methods such as BEAST assume random sampling from a population and interpretation of their output can be less straightforward when lineage-based sampling is used [25], and our dataset was restricted to a single haplogroup. However, the coinciding topologies of the trees generated by the different methods in our study shows that the phylogenetic structure is robust to this concern. Note that in our analyses, choosing the model of population growth had the major influence on the results-probably larger than the influence of sampling from a particular lineage-but we report the general consensus of the results across all settings used.
In an additional analysis we included two G1 samples from the 1000 Genomes Project (NA20858 and NA20870, Gujarati Indians sampled in Houston, Texas (GIH), 2-4X average coverage). Data were handled in the same way, although the lower coverage of the 1000 Genomes samples halved the number of SNP calls and the filtered dataset consisted of 22 samples and 393 SNPs (S4 Table). The parsimony method yielded two optimal topologies, and the one supporting the monophyly of all non-Indian lineages was preferred as a more likely reconstruction.

Results
We genotyped the haplogroup G1-specific marker M285 [28] in 5,346 individuals from 27 Eurasian populations (Table 1) and identified 367 M285-positive samples, which were further genotyped by 17 Y-STRs (S1 Table). For 19 haplogroup G1 samples and one outgroup G2 sample we performed complete sequencing of the "extended gold standard" regions of the Y-chromosome.

The frequency distribution of the haplogroup G1
The frequency distribution of haplogroup G1 in Eurasia is presented in Fig 2, which is based on the dataset from Table 1. This haplogroup is distributed over a large area from Italy in the west to Mongolia in the east, but is present at high frequencies only in an uninterrupted area including the Central Asian steppes and Iranian-Armenian plateau. Two frequency peaks can be seen at the opposite sides of this area, namely in North Kazakhstan (up to 80%) within the steppe part and in Armenia (up to 42%) within its mountainous part. In Kazakhs, haplogroup G1 is typical of the Argyn tribe: among 291 G1 samples with known tribal affiliation in Kazakhs, 262 (90%) belong to the Argyn tribe. In Armenians, this haplogroup is particularly frequent in Hemsheni Armenians (42%). Both populations are not small: according to a census performed in the beginning of the 20 th century-tribal affiliation was not recorded in later censuses-there were around 500,000 Argyns [36] and now the population is expected to be larger; the estimated present-day number of Hemsheni Armenians is 150,000 [35]. Thus, the increased frequency of G1 cannot be explained by recent genetic drift and likely indicates drift during the formation of these populations many centuries ago.
It is notable that the area of haplogroup G1, including the Eurasian steppes from the North Black Sea region to the Mongolian Altai and South-Western Asian uplands (Iran and historical Great Armenia), corresponds well with the area populated by Iranic speakers in the second and first millennia BC (Fig 1A). This correspondence was statistically confirmed by AMOVA (S3 Table).

STR-variation within haplogroup G1
On the network (Fig 3) four clusters are visible: two include only Armenian samples, while other two are specific to Kazakhs and Bashkirs, respectively. Samples from other populations are spread all around the network and do not form clusters. The Kazakh cluster is highly specific to the Argyn tribe within the Kazakhs; the Armenian cluster includes more Hemsheni Armenians than other Armenian populations; all but one members of the Bashkir cluster belong to the Kangly tribe of Bashkirs. The Y-STR pattern thus shows that haplogroup G1 is genetically diverse and widespread, with some sub-branches at high frequencies due to founder effects, while others remain at very low frequencies in occasional locations within the area of the haplogroup.
The haplotype diversity of haplogroup G1 varies drastically from 92% in Iran to zero in Mongolia ( Table 2). The map (Fig 4) reveals a cline of decreasing diversity from West Iran to the eastern parts of South-West Asia and further northward to the Eurasian steppes.

A sequence-based phylogenetic tree of haplogroup G1
We sequenced~11 Mb of the Y-chromosome in 19 samples selected using three criteria to cover the maximum diversity within haplogroup G1. The average coverage was 67x, ranging from 48x to 88x. Among the 766 SNPs in the filtered dataset (see Methods for details) 281 have already been listed by ISOGG (www.isogg.org) and/or YFULL (www.yfull.com), most of these SNPs defined our outgroup G2 sample. We gave the remaining SNPs (S4 Table) names from GG001 to GG388 (GG is the abbreviation of Gene Geography).   cluster and an Indian cluster (described below) form independent branches. The Mongol sample forms a branch on its own, although this Mongolian branch then joins the Kazakh cluster in agreement with common geographical and historical background of the two groups. This tree corresponds in general with the pattern revealed by the STR-based network (Fig 3). Kazakh, Armenian and Bashkir clusters are clearly visible on both plots. However, Armenians, which seemed to have two clusters and multiple single-haplotype mini-branches from the STR data, all turned out to belong to one single and compact cluster when complete Y-chromosomal resequencing was performed. Similarly, the Mongolian sample, which seemed to form a separate branch on the STR-based plot, actually joins the Kazakh cluster. We conclude that haplogroup G1 lineages actually form a restricted number of clusters, in contrast to the impression one can get from STR-data, with the caveat that the number of sequences examined thus far is limited.
The presence of additional clusters was confirmed when we included two GIH (Gujarat Indians from Houston) samples from the 1000 Genomes Project, which are the only publicly available data on haplogroup G1. Including the low coverage sequences halved the number of SNPs called in all samples (S4 Table), but tree revealed the same topology, and the Indian G1s formed their own cluster (Fig 5). One technical point is notable: the lengths of all the branches on the tree are similar, as they should be if the mutation rate is constant. The only exception is the very long branches of the samples from the 1000 Genomes Project, which is likely caused by the filtering criteria not being optimized for low coverage datasets. However, 26 SNPs were independently called in both samples, thus confirming the reality of the Indian-specific branch of haplogroup G1.
The Kazakh cluster fits the previously described G-L1323 branch (www.isogg.org), while the Bashkir, Armenian, Mongolian and Indian branches were not previously reported. Fig 5  approximates the phylogenetic relations between five branches found in our study and three previously known ones.

Estimating the mutation rate
The Argyn tribe in which haplogroup G1 predominates is believed to descend from a single male common ancestor (Argyn) and is divided into 12 clans (Fig 6B). Though there is no historical evidence for the existence of Argyn, who is known only from genealogical tradition, his great-grandson Karakhoja is a historical personality and is mentioned, among other sources, as ambassador of the Toshtamish khan-ruler of the Golden Horde-to Tamerlane in 1405. Most of the Argyn clans are believed to originate from Karakhoja (S4 Fig) and other clans are believed to originate from his brother Somdyk.
The genetic tree based on high-throughput sequencing of the Kazakh G1 chromosomes (Fig 6A) perfectly fits the genealogical tree: representatives of Argyn clans who originated from Karakhoja (from Kazakh1 to Kazakh6) form a single and young subcluster and all are equidistant from the MRCA, as predicted by the genealogy (Fig 6B). Thus, the de jure ancestor known from genealogical tradition and historical records was likely to be also the de facto biological ancestor of most present-day male members of the Argyn tribe. Considering the time span of 606 years between Karakhoja (who was around 50 years old in 1405 and then likely fathered his sons on average around 1385) and average date of birth of the 6 present-day Kazakhs sampled (1991), the total length of Y-chromosomal segments sequenced in each of these 6 Kazakh samples (10,005,352 bp), and the average number of accumulated mutations 4.67 (S2 Fig), we obtained the mutation rate for Y-chromosomal sequences 0.77×10 -9 per bp per year.
While this paper was under review, we obtained experimental data from three additional samples, representing clans claiming their origin from Karakhoja's brother Somdyk (Kazakh 7, Kazakh8, and Kazakh9, Fig 6B). These samples formed a subcluster of their own (Fig 6A), a genetic finding which fully agrees with the genealogical tradition. This allows us to reconstruct even minor details, like the observations that mutations 23081087 C->T and 23526483 A->G occurred between Akzhol and his great-great-grandson Ermen, and that no mutations occurred between Akzhol and his son Karakhoja within the sequenced regions. Considering the time span of 627 years between Akzhol (we added 30 years-the typical male generation time-to the time estimated above for Karakhoja) and the average date of birth of the 9 present-day Kazakhs sampled (1982), the total length of Y-chromosomal segments sequenced in each of these 9 Kazakh samples (9972660 bp), and the average number of accumulated mutations (4.89: S2 Fig), we obtained a SNP mutation rate for Y-chromosomal sequences of 0.78×10 -9 per bp per year. This value agrees well with the aforementioned initial estimation.
We applied the same approach to Y-STR data also. Including data on 15-Y-STR haplotypes in Argyns (S1 Table, Fig 3), we counted 21 generations (using a male generation time of around 30 years [14]) and 0.68 mutations on average. Thus we obtained the mutation rate for Y-STRs 0.0022 per locus per generation. It is very close to the "genealogical" rate [18,22,46], despite the time span used (600 years) better fitting the population events used for estimating the "evolutionary" rate [54].
Discussion generation gives half a million descendants in the 19 th generation and Kazakh families had 3.5 children on average [36]. Also, an even more impressive expansion up to 16 million descendants was suggested for the same medieval steppe societies [53]. Note that the traditional genealogical partitioning of Argyns into three clans corresponds well with the Y-STR data (S4 Fig). This finding also questions the hypothesis [5] about the relationship between the Argyn subclan Madjars and Magyars (Hungarians), because haplogroup G1 (comprising 82% of the Madjar gene pool) finds its place within other Argyn Kazakh (S4 Fig) while no G1 samples have been reported in Hungarians so far.
The expansion in the Hemsheni Armenian is genetically dated to 1150 YBP using our rate (S5 Table). It corresponds well with the historical evidence [30,48] that the Hemsheni originated from relatives and servants of Prince Shapuh Amatuni, who migrated in 791 from the Abbasid Persian state.
The expansion in the Kangly tribe of Bashkirs is genetically dated to the 15 th century AD (S5 Table). This tribe originated from the Pechenegs around the 8 th century AD, then joined the Bashkirs, and later expansion in a restricted part of the tribe might have been caused by demographic changes when it became part of the Golden Horde in the 14 th century and part of the Russian state in 16 th century.
We note that despite geographic proximity, the ancestor of the G1 cluster in Bashkirs had no close genetic relationship to the corresponding ancestor in Kazakhs. These branches (and the third branch detected in Mongolians) have survived in the Eurasian steppe perhaps since the Scythian epoch.
The remarkable coincidence between the genealogical tree of the Argyn Kazakh clan (Fig 6B) and the genetic tree obtained from full Y-chromosomal sequences (Fig 6A) allowed us to suggest an independent calibration of the mutation rate of Y-chromosomal SNPs. This "clan" rate has been tested only within the time frame in which it was obtained (a few centuries), and in cases when it is reasonable to suppose expansion of a single paternal line rather than multiple lineages in the founding population, and by applying the "gold standard" portion of the Y-chromosome (that included in the BigY technology used in our study). Provided these limitations are taken into account, this "clan-based" calibration might be at least as reliable as calibrations based on archeological evidence, because archeological dates are subject to revision and of uncertain relationship to genetic events. For example, the calibration of the Y-chromosomal mutation rate in [41] is based on "archeological evidence that humans first colonized America around 15 kya" while the study that provides the commonly-used calibration of the mitochondrial DNA control region [15] relies on a "major wave of migration [which] brought one population ancestral to Amerinds from north-eastern Siberia to America 20,000-25,000 years ago". Fortunately, despite differences in approaches, all mutation rates suggested for the "full" sequences of the Y-chromosome fall within the interval 0.6-1.0 ×10 -9 per bp per year, and this uncertainty may be further narrowed, as we demonstrated for the haplogroup G1.
Supporting Information S1 Data. The archive includes 20 VCF files for 20 sequenced samples, 20 BED files showing Y-chromosomal ranges where a given sample was sequenced reliably, and the BED files listing ranges intersecting between all 20 samples (used for tree construction) and between 6 Kazakh samples (used for calibrating the mutation rate).