Sequencing Y Chromosomes Resolves Discrepancy in Time to Common Ancestor of Males Versus Females

See allHide authors and affiliations

Science  02 Aug 2013:
Vol. 341, Issue 6145, pp. 562-565
DOI: 10.1126/science.1237619

Examining Y

The evolution of human populations has long been studied with unique sequences from the nonrecombining, male-specific Y chromosome (see the Perspective by Cann). Poznik et al. (p. 562) examined 9.9 Mb of the Y chromosome from 69 men from nine globally divergent populations—identifying population and individual specific sequence variants that elucidate the evolution of the Y chromosome. Sequencing of maternally inherited mitochondrial DNA allowed comparison between the relative rates of evolution, which suggested that the coalescence, or origin, of the human Y chromosome and mitochondria both occurred approximately 120 thousand years ago. Francalacci et al. (p. 565) investigated the sequence divergence of 1204 Y chromosomes that were sampled within the isolated and genetically informative Sardinian population. The sequence analyses, along with archaeological records, were used to calibrate and increase the resolution of the human phylogenetic tree.


The Y chromosome and the mitochondrial genome have been used to estimate when the common patrilineal and matrilineal ancestors of humans lived. We sequenced the genomes of 69 males from nine populations, including two in which we find basal branches of the Y-chromosome tree. We identify ancient phylogenetic structure within African haplogroups and resolve a long-standing ambiguity deep within the tree. Applying equivalent methodologies to the Y chromosome and the mitochondrial genome, we estimate the time to the most recent common ancestor (TMRCA) of the Y chromosome to be 120 to 156 thousand years and the mitochondrial genome TMRCA to be 99 to 148 thousand years. Our findings suggest that, contrary to previous claims, male lineages do not coalesce significantly more recently than female lineages.

The Y chromosome contains the longest stretch of nonrecombining DNA in the human genome and is therefore a powerful tool with which to study human history. Estimates of the time to the most recent common ancestor (TMRCA) of the Y chromosome have differed by a factor of about 2 from TMRCA estimates for the mitochondrial genome. Y-chromosome coalescence time has been estimated in the range of 50 to 115 thousand years (ky) (13), although larger values have been reported (4, 5), whereas estimates for mitochondrial DNA (mtDNA) range from 150 to 240 ky (3, 6, 7). However, the quality and quantity of data available for these two uniparental loci have differed substantially. Whereas the complete mitochondrial genome has been resequenced thousands of times (6, 8), fully sequenced diverse Y chromosomes have only recently become available. Previous estimates of the Y-chromosome TMRCA relied on short resequenced segments, rapidly mutating microsatellites, or single-nucleotide polymorphisms (SNPs) ascertained in a small panel of individuals and then genotyped in a global panel. These approaches likely underestimate genetic diversity and, consequently, TMRCA (9).

We sequenced the complete Y chromosomes of 69 males from seven globally diverse populations of the Human Genome Diversity Panel (HGDP) and two additional African populations: San (Bushmen) from Namibia, Mbuti Pygmies from the Democratic Republic of Congo, Baka Pygmies and Nzebi from Gabon, Mozabite Berbers from Algeria, Pashtuns (Pathan) from Pakistan, Cambodians, Yakut from Siberia, and Mayans from Mexico (fig. S1). Individuals were selected without regard to their Y-chromosome haplogroups.

The Y-chromosome reference sequence is 59.36 Mb, but this includes a 30-Mb stretch of constitutive heterochromatin on the q arm, a 3-Mb centromere, 2.65-Mb and 330-kb telomeric pseudoautosomal regions (PAR) that recombine with the X chromosome, and eight smaller gaps. We mapped reads to the remaining 22.98 Mb of assembled reference sequence, which consists of three sequence classes defined by their complexity and degree of homology to the X chromosome (10): X-degenerate, X-transposed, and ampliconic. Both the high degree of self-identity within the ampliconic tracts and the X-chromosome homology of the X-transposed region render portions of the Y chromosome ill suited for short-read sequencing. To address this, we constructed filters that reduced the data to 9.99 million sites (11) (Fig. 1 and fig. S2). We then implemented a haploid model expectation-maximization algorithm to call genotypes (11).

Fig. 1 Callability mask for the Y chromosome.

Exponentially weighted moving averages of read depth (blue line) and the proportion of reads mapping ambiguously (MQ0 ratio; violet line) versus physical position. Regions with values outside the envelopes defined by the dashed lines (depth) or dotted lines (MQ0) were flagged (blue and violet boxes) and merged for exclusion (gray boxes). The complement (black boxes) defines the regions within which reliable genotype calls can be made. Below, a scatter plot indicates the positions of all observed SNVs. Those incompatible with the inferred phylogenetic tree (red) are uniformly distributed. The X-degenerate regions yield quality sequence data, ampliconic sequences tend to fail both filters, and mapping quality is poor in the X-transposed region.

We identified 11,640 single-nucleotide variants (SNVs) (fig. S3). A total of 2293 (19.7%) are present in dbSNP (v135), and we assigned haplogroups on the basis of the 390 (3.4%) present in the International Society of Genetic Genealogy (ISOGG) database (12) (fig. S4). At SNVs, median haploid coverage was 3.1x (interquartile range 2.6 to 3.8x) (table S1 and fig. S5), and sequence validation suggests a genotype calling error rate on the order of 0.1% (11).

Because mutations accumulate over time along a single lengthy haplotype (13), the male-specific region of the Y chromosome provides power for phylogenetic inference. We constructed a maximum likelihood tree from 11,640 SNVs using the Tamura-Nei nucleotide substitution model (Fig. 2) and, in agreement with (14), observe strong bootstrap support (500 replicates) for the major haplogroup branching points. The tree both recapitulates and adds resolution to the previously inferred Y-chromosome phylogeny (fig. S6), and it characterizes branch lengths free of ascertainment bias. We identify extraordinary depth within Africa, including lineages sampled from the San hunter-gatherers that coalesce just short of the root of the entire tree. This stands in contrast to a tree from autosomal SNP genotypes (15), wherein African branches were considerably shorter than others; genotyping arrays primarily rely on SNPs ascertained in European populations and therefore undersample diversity within Africa. Two regions of reduced branch length in our tree correspond to rapid expansions: the out-of-Africa event (downstream of F-M89) and the agriculture-catalyzed Bantu expansions (downstream of E-M2). Among the three hunter-gatherer populations, we find a relatively high number of B2 lineages. Within this haplogroup, six Baka B-M192 individuals form a distinct clade that does not correspond to extant definitions (11) (fig. S7). We estimate this previously uncharacterized structure to have arisen ~35 thousand years ago (kya).

Fig. 2 Y-chromosome phylogeny inferred from genomic sequencing.

This tree recapitulates the previously known topology of the Y-chromosome phylogeny; however, branch lengths are now free of ascertainment bias. Branches are drawn proportional to the number of derived SNVs. Internal branches are labeled with defining ISOGG variants inferred to have arisen on the branch. Leaves are colored by major haplogroup cluster and labeled with the most derived mutation observed and the population from which the individual was drawn. Previously uncharacterized structure within African hgB2 is indicated in orange. (Inset) Resolution of a polytomy was possible through the identification of a variant for which hgG retains the ancestral allele, whereas hgH and hgIJK share the derived allele.

We resolve the polytomy of the Y macro-haplogroup F (16) by determining the branching order of haplogroups G, H, and IJK (Fig. 2 and fig. S6). We identified a single variant (rs73614810, a C→T transition dubbed “M578”) for which haplogroup G retains the ancestral allele, whereas its brother clades (H and IJK) share the derived allele. Genotyping M578 in a diverse panel confirmed the finding (table S2). We thereby infer more recent common ancestry between hgH and hgIJK than between either and hgG. M578 defines an early diversification episode of the Y phylogeny in Eurasia (11).

To account for missing genotypes, we assigned each SNV to the root of the smallest subtree containing all carriers of one allele or the other and inferred that the allele specific to the subtree was derived (fig. S8). We used the chimpanzee Y-chromosome sequence to polarize 398 variants assigned to the deepest split—a task complicated by substantial structural divergence (11, 17).

We estimated the coalescence time of all Y chromosomes using both a molecular clock–based frequentist estimator and an empirical Bayes approach that uses a prior distribution of TMRCA from coalescent theory and conducts Markov chain simulation to estimate the likelihood of parameters given a set of DNA sequences (GENETREE) (11, 18) (Table 1). To directly compare the TMRCA of the Y chromosome to that of the mtDNA, we estimated their respective mutation rates by calibrating phylogeographic patterns from the initial peopling of the Americas, a recent human event with high-confidence archaeological dating.

Table 1 TMRCA and Ne estimates for the Y chromosome and mtDNA.

Pop., population.

View this table:

Archaeological evidence indicates that humans first colonized the Americas ~15 kya via a rapid coastal migration that reached Monte Verde II in southern Chile by 14.6 kya (19). The two Native American Mayans represent Y-chromosome hgQ lineages, Q-M3 and Q-L54*(xM3), that likely diverged at about the same time as the initial peopling of the continents. Q is defined by the M242 mutation that arose in Asia. A descendent haplogroup, Q-L54, emerged in Siberia and is ancestral to Q-M3. Because the M3 mutation appears to be specific to the Americas (20), it likely occurred after the initial entry, and the prevalence of M3 in South America suggests that it emerged before the southward migratory wave. Consequently, the divergence between these two lineages provides an appropriate calibration point for the Y mutation rate. The large number of variants that have accumulated since divergence, 120 and 126, contrasts with the pedigree-based estimate of the Y-chromosome mutation rate, which is based on just 4 mutations (21). Using entry to the Americas as a calibration point, we estimate a mutation rate of 0.82 × 10−9 per base pair (bp) per year [95% confidence interval (CI): 0.72 × 10−9 to 0.92 × 10−9/bp/year] (table S3). False negatives have minimal effect on this estimate due to the low probability, at 5.7x and 8.5x coverage, of observing fewer than two reads at a site (observed proportions: 3.1% and 0.6%) and due to the fact that the number of unobserved singletons possessed by one individual is offset by a similar number of Q doubletons unobserved in the same individual and thereby misclassified as singletons possessed by the other (11) (figs. S9 and S10). This calibration approach assumes approximate coincidence between the expansion throughout the Americas and the divergence of Q-M3 and Q-L54*(xM3), but we consider deviation from this assumption and identify a strict lower bound on the point of divergence using sequences from the 1000 Genomes Project (11). As a comparison point, we consider the out-of-Africa expansion of modern humans, which dates to approximately 50 kya (22) and yields a similar mutation rate of 0.79 × 10−9/bp/year.

We constructed an analogous pipeline for high coverage (>250x) mtDNA sequences from the 69 male samples and an additional 24 females from the seven HGDP populations (11) (fig. S11). As in the Y-chromosome analysis, we calibrated the mtDNA mutation rate using divergence within the Americas. We selected the pan-American hgA2, one of several initial founding haplogroups among Native Americans. The star-shaped phylogeny of hgA2 subclades suggests that its divergence was coincident with the rapid dispersal upon the initial colonization of the continents (23). Calibration on 108 previously analyzed hgA2 sequences (11) (fig. S12) yields a point estimate equivalent to that from our seven Mayan mtDNAs, but within a narrower confidence interval. From this within-human calibration, we estimate a mutation rate of 2.3 × 10−8/bp/year (95% CI: 2.0 × 10−8 to 2.5 × 10−8/bp/year), higher than that from human-chimpanzee divergence but similar to other estimates using within-human calibration points (24, 25).

The global TMRCA estimate for any locus constitutes an upper bound for the time of human population divergence under models without gene flow. We estimate the Y-chromosome TMRCA to be 138 ky (120 to 156 ky) and the mtDNA TMRCA to be 124 ky (99 to 148 ky) (Table 1) (11). Our mtDNA estimate is more recent than many previous studies, the majority of which used mutation rates extrapolated from between-species divergence. However, mtDNA mutation rates are subject to a time-dependent decline, with pedigree-based estimates on the faster end of the spectrum and species-based estimates on the slower. Because of this time dependency and the need to calibrate the Y and mtDNA in a comparable manner, it is more appropriate here to use within-human clade estimates of the mutation rate.

Rather than assume the mutation rate to be a known constant, we explicitly account for the uncertainty in its estimation by modeling each TMRCA as the ratio of two random variables. We estimate the ratio of the mtDNA TMRCA to that of the Y chromosome to be 0.90 (95% CI: 0.68 to 1.11) (fig. S13). If, as argued above, the divergence of the Y-chromosome Q lineages occurred at approximately the same time as that of the mtDNA A2 lineages, then the TMRCA ratio is invariant to the specific calibration time used. Regardless, the conclusion of parity is robust to possible discrepancy between the divergence times within the Americas (11). Using comparable calibration approaches, the Y and mtDNA coalescence times are not significantly different. This conclusion would hold whether or not an alternative approach would yield more definitive TMRCA estimates.

Our observation that the TMRCA of the Y chromosome is similar to that of the mtDNA does not imply that the effective population sizes (Ne) of males and females are similar. In fact, we observe a larger Ne in females than in males (Table 1). Although, due to its larger Ne, the distribution from which the mitochondrial TMRCA has been drawn is right-shifted with respect to that of the Y-chromosome TMRCA, the two distributions have large variances and overlap (Fig. 3).

Fig. 3 Similarity of TMRCA does not imply equivalent Ne of males and females.

The TMRCA for a given locus is drawn from a predata (i.e., prior) distribution that is a function of Ne, generation time, sample size, and demographic history. Consider the distribution of possible TMRCAs for a set of 100 uniparental chromosomes. Although the Mbuti mtDNA Ne is twice as large as that of the Baka Y chromosome, the corresponding predata TMRCA distributions overlap considerably.

Dogma has held that the common ancestor of human patrilineal lineages, popularly referred to as the Y-chromosome “Adam,” lived considerably more recently than the common ancestor of female lineages, the so-called mitochondrial “Eve.” However, we conclude that the mitochondrial coalescence time is not substantially greater than that of the Y chromosome. Indeed, due to our moderate-coverage sequencing and the existence of additional rare divergent haplogroups, our analysis may yet underestimate the true Y-chromosome TMRCA.

Supplementary Materials

Materials and Methods

Supplementary Text

Figs. S1 to S13

Tables S1 to S3

Data File S1

References (2651)

References and Notes

  1. Materials and methods are available as supplementary materials on Science Online.
  2. Acknowledgments: We thank O. Cornejo, S. Gravel, D. Siegmund, and E. Tsang for helpful discussions; M. Sikora and H. Costa for mapping reads from Gabonese samples; and H. Cann for assistance with HGDP samples. This work was supported by National Library of Medicine training grant LM-07033 and NSF graduate research fellowship DGE-1147470 (G.D.P.); NIH grant 3R01HG003229 (B.M.H. and C.D.B.); NIH grant DP5OD009154 (J.M.K. and E.S.); and Institut Pasteur, a CNRS Maladies Infectieuses Émergentes Grant, and a Foundation Simone et Cino del Duca Research Grant (L.Q.M.). P.A.U. consulted for, P.A.U. and B.M.H. have stock in, and C.D.B. is on the advisory board of a project at 23andMe. C.D.B. is on the scientific advisory boards of Personalis, Inc.; InVitae (formerly Locus Development, Inc.); and M.S. is a scientific advisory member and founder of Personalis, a scientific advisory member for Genapsys Former, and a consultant for Illumina and Beckman Coulter Society for American Medical Pathology. B.M.H. formerly had a paid consulting relationship with Variants have been deposited to dbSNP (ss825679106–825690384). Individual level genetic data are available, through a data access agreement to respect the privacy of the participants for transfer of genetic data, by contacting C.D.B.
View Abstract

Navigate This Article