The Divergence of Neandertal and Modern Human Y Chromosomes.

Sequencing the genomes of extinct hominids has reshaped our understanding of modern human origins. Here, we analyze ∼120 kb of exome-captured Y-chromosome DNA from a Neandertal individual from El Sidrón, Spain. We investigate its divergence from orthologous chimpanzee and modern human sequences and find strong support for a model that places the Neandertal lineage as an outgroup to modern human Y chromosomes-including A00, the highly divergent basal haplogroup. We estimate that the time to the most recent common ancestor (TMRCA) of Neandertal and modern human Y chromosomes is ∼588 thousand years ago (kya) (95% confidence interval [CI]: 447-806 kya). This is ∼2.1 (95% CI: 1.7-2.9) times longer than the TMRCA of A00 and other extant modern human Y-chromosome lineages. This estimate suggests that the Y-chromosome divergence mirrors the population divergence of Neandertals and modern human ancestors, and it refutes alternative scenarios of a relatively recent or super-archaic origin of Neandertal Y chromosomes. The fact that the Neandertal Y we describe has never been observed in modern humans suggests that the lineage is most likely extinct. We identify protein-coding differences between Neandertal and modern human Y chromosomes, including potentially damaging changes to PCDH11Y, TMSB4Y, USP9Y, and KDM5D. Three of these changes are missense mutations in genes that produce male-specific minor histocompatibility (H-Y) antigens. Antigens derived from KDM5D, for example, are thought to elicit a maternal immune response during gestation. It is possible that incompatibilities at one or more of these genes played a role in the reproductive isolation of the two groups.

Introduction

A central goal of human population genetics and paleoanthropology is to elucidate the relationships among ancient populations. Before the emergence of anatomically modern humans in the Middle Pleistocene ∼200 thousand years ago (kya), archaic humans lived across Africa, Europe, and Asia in highly differentiated populations. Modern human populations that expanded out of Africa in the Upper Pleistocene received a modest genetic contribution from at least two archaic hominin groups, the Neandertals and Denisovans.2, 3, 4, 5 Especially in light of hypothesized genetic incompatibilities between Neandertals and modern humans, it is important to characterize differentiation between their ancestral populations and to investigate potential barriers to gene flow.

When populations diverge from one another, each retains a subset of the variation that existed in the ancestral population. Consequently, sequence divergence times usually exceed population divergence times, and this effect is more pronounced when the ancestral effective population size was large. In humans, a large fraction of genetic diversity is due to ancient polymorphisms that arose long before the emergence of anatomically modern traits. As a result, Neandertal and modern haplotypes are often no more diverged than modern human sequences are among themselves. This fact complicates the search for introgressed genomic segments, but two features facilitate their detection.6, 7 First, due to low levels of polymorphism among Neandertals, introgressed sequences are often quite similar to those of the Neandertal reference. Second, these regions have elevated linkage disequilibrium due to the relatively recent date of admixture, ∼50 kya.8, 9, 10 Although introgressed Neandertal sequences have been identified in modern human autosomes and X chromosomes, no mitochondrial genome (mtDNA) sequences of Neandertal origin have been reported in modern humans, and Neandertal Y-chromosome sequences have not yet been characterized.

Because uniparentally inherited loci have much smaller effective population sizes than autosomal or X-linked loci, the expected differences between sequence and population divergence times are smaller. Therefore, studying these loci can help to delineate an upper bound for the time at which populations last exchanged genetic material. To date, five Neandertal individuals have been whole-genome sequenced to 0.1× coverage or higher,2, 5 but all were female. Full mtDNA sequences are also available for eight individuals from Spain, Germany, Croatia, and Russia,11, 12 but the relationship between Neandertal and modern human Y chromosomes remains unknown.

In this work, we analyzed ∼120 kb of exome-captured Y-chromosome sequence from an ∼49,000-year-old (uncalibrated 14C) Neandertal male from El Sidrón, Spain. We compare it to the human and chimpanzee reference sequences and to the sequences of two Mbo individuals who carry the A00 haplogroup, the most deeply branching group known. We identify the relationship between the Neandertal and modern human Y chromosomes and estimate the time to their most recent common ancestor (TMRCA). We also examine coding differences and explore their potential significance for reproductive isolation.

Material and Methods

Sequence Data and Processing

We used the Y-chromosome sequences from the exome capture of a Neandertal from El Sidrón, Spain, and we downloaded the complete sequences of two A00 Y chromosomes. The Neandertal data included coding, non-coding, and off-target sequences, and all three sequences were mapped against the GRCh37 reference. Given that the A00 sequences were closely related,15, 16 we merged them to increase coverage. We called bases for both the Neandertal and A00 sequences by using SAMtools mpileup (v.1.1), specifying input options to count anomalous read pairs (-A), recalculate base qualities (-E), and filter out poor-quality bases (-Q 17) and poorly mapping reads (-q 20).

We then identified overlapping regions and excluded coordinates with unusually high coverage, filtering out sites with coverage greater than the mean plus five times its square root (Figure S1). Under a Poisson model, this cutoff would elicit the loss of less than one genuine site per 10,000. Finally, we removed sites with inconsistent base calls, discarding those with more than two reads differing from the consensus allele and those for which more than one third of the observed bases did not match the consensus. This filter should minimize the effects of postmortem DNA damage and of modern contamination.

Using the blastz file chrY.hg19.panTro4.net.axt.gz, we identified the subset of regions within which the human sequences align to the chimpanzee reference. This yielded a total of 118,643 base pairs (bp). In what follows, we refer to this set of sites as “filter 1.” We also identified a second, more restrictive, set of regions totaling 100,324 bp, “filter 2,” by further requiring that the alignment correspond to the chimpanzee Y chromosome rather than to another chimpanzee chromosome (Tables S1A and S1B).

For each position within these regions, we determined whether the Neandertal, A00, or both differed from the human reference sequence. We then used the corresponding chimpanzee allele as a proxy for the ancestral state in order to assign the mutation to the appropriate branch of the tree relating the four sequences (Figure 1A). In doing so, we discarded five sites: two at which the chimpanzee carries a third allele, one for which the chimpanzee carries a deletion, and two that were specific to A00 but only supported by a single read. Excluding these sites had little impact on our analyses.

Figure 1

Tree Inference

(A) A priori, three trees could feasibly have related the Y chromosomes of the chimpanzee (Chimp), the Neandertal (Neander), haplogroup A00, and the human reference (Ref). Mutations on branch a support topology i, with the Neandertal lineage as the outgroup to those of modern humans, whereas mutations on branches b and c support topologies ii and iii, respectively. Branches d, e, and f correspond to mutations private to individual lineages.

(B) Counts of SNVs consistent with each branch. Columns refer to sets of coordinates considered (see Materials and Methods). Incompatible sites are those that cannot be explained by a single mutation on any of the three trees.

Estimating TMRCA

To estimate the TMRCA of the Neandertal and modern human Y chromsomes (TNR), we decomposed this quantity (Figure 2) into the sum of the TMRCA of modern humans (TAR) and the time separating the most recent common ancestor of modern humans from its common ancestor with the Neandertal lineage (TNM):We then estimated TAR and used two methods to estimate α.

Figure 2

Estimating the TMRCA of Neandertal and Modern Y Chromosomes

The quantity of primary interest is TNR = TNM+ TAR. Branches are labeled as in Figure 1, and “M” denotes the most recent common ancestor of modern human lineages.

To estimate TAR, we used sequence data from the ancient Ust’-Ishim sample, first applying the filters described for the A00 sequences. To reduce the potential impact of postmortem DNA damage, we restricted this analysis to coordinates covered by at least three sequencing reads. We further restricted to the subset of Poznik et al. regions in which the human reference sequence is based on bacterial artificial chromosome clones derived from the RP-11 individual, a known carrier of haplogroup R1b. This left ∼7.83 Mb of sequence within which to assign variants to the appropriate branches (Figure S2, Appendix A). Using the known age of the Ust’-Ishim individual and the constrained optimization procedure described in Rasmussen et al., we obtained parametric bootstrap estimates for TAR as well as for the mutation rate and the TMRCA of haplogroup K-M526 (Appendix A). Briefly, we sampled from the process that generated the observed tree (Figure S2) by simulating the number of single nucleotide variants (SNVs) on each branch as a Poisson draw with mean equal to the observed number of mutations. To obtain bootstrap samples of the three parameters, we maximized their joint likelihood for each tree replicate.

In our first approach to estimate α, we used the relative numbers of mutations assigned to branches a, d, and e (Figure 1), assigning the four sites that did not fit the consensus topology to the A00 or reference lineages, as appropriate (Appendix A). The proportion of time represented by branch a is:

Therefore, assuming a time-homogeneous mutation rate, the number of branch-a mutations is binomially distributed with parameters p = (α – 1) ∕ (α + 1) and n equal to the total number of mutations. Estimating p from the data leads directly to a point estimate and confidence interval (CI) for α. This first method has the appealing property that it is independent of both the mutation rate and the absolute values of the times. However, the estimation error might be suboptimal due to uncertainty in both the numerator and the denominator.

In the second method, we estimated α via the ratio , making use of the fact that we can estimate TAR with greater certainty than we can TNM. To estimate TNM, we restricted our attention to sequences overlapping the ∼8.8 Mb of sequence analyzed by Karmin et al., leaving 80,420 bp (“filter 3”) or 75,596 bp (“filter 4”) when confining our analysis to those sites that passed filter 2 (Figure 1, Tables S2a and S2b). Let l equal the total length of sequence under consideration (e.g., 80.42 kb), let μ equal the mutation rate over the full 8.8 Mb, let r equal the ratio of the mutation rate within the smaller region to that of the larger, and let s equal the number of mutations shared by A00 and the reference sequence within the smaller region. With these, we constructed the estimator . Similarly, let L equal the subset of the 8.8 Mb for which the A00 sequence had 3× or greater coverage (also ∼8.8 Mb), and let S equal the number of mutations unique to either the reference sequence or to A00 over the entire 8.8 Mb. We can estimate TAR with and α with:

We estimated r by comparing the number of mutations unique to a single branch of the Y-chromosome tree of Karmin et al., both within the full 8.8-Mb region and within the ∼80-kb subset. These numbers, 32,853 and 279 (238 under filter 4), respectively, correspond to a relative mutation rate of 0.93 (95% CI: 0.82–1.04) (0.84 under filter 4 [95% CI: 0.74–0.95]). Because selection has the strongest effect on lower frequency mutations, we also estimated r by using only shared variants, and this yielded nearly identical point estimates.

Finally, to construct a CI for α, we sampled values of s and S from Poisson distributions with means equal to the observed numbers of mutations, and we sampled rl/L as the ratio of two Poisson random variables with means equal to 279 (or 238) and 32,853, respectively.

Functional Variation

We determined whether each mutation overlaps with annotated RefSeq genes and whether it overlaps with coding sequence (Figure 1, Data S5). For each coding SNV, we determined whether the mutation results in silent, missense, or nonsense mutations, but we did not consider frameshift mutations. For each coding non-synonymous mutation, we used the HumDiv model of PolyPhen-2 to evaluate ancestral-to-derived changes and MutationTaster to evaluate reference-to-alternative changes. We report findings from all sites for which these programs were able to make predictions.

Results

With the chimpanzee Y chromosome as the outgroup, three tree topologies could have related the lineages of the Neandertal, haplogroup A00, and the human reference (Figure 1A). To identify which of the three was consistent with the data, the key question was which of the three possible pairs of sequences is the most closely related. Of 118,643 sites (Figure 1B, filter 1) for which we had Neandertal data and human-chimpanzee reference alignments, we identified 24 biallelic SNVs for which the Neandertal sequence shared the chimpanzee allele and differed from both A00 and the human reference. In contrast, the chimpanzee and A00 sequences shared just four SNVs not present in the other sequences, and the chimpanzee and human reference sequences shared zero. Taken together, these data strongly support the tree that places the Neandertal Y as the most distantly related to the others (Figure 1A, tree i). Two of the four variants that are inconsistent with this topology are known to segregate within modern humans and are therefore the result of recurrent mutations or contamination (Appendix A).

Upon elucidating the topology of the tree relating the Neandertal Y chromosome to those of modern humans, our next goal was to estimate the divergence time. We decomposed the TMRCA, TNR, into the sum of two intervals (Figure 2): the TMRCA of A00 and the reference, TAR, and the time between their common ancestor and the common ancestor of the Neandertal lineage and that of modern humans, TNM. To estimate TNR, we estimated TAR and the ratio α ≡ TNR ∕ TAR, taking care to consider uncertainty both in the mutation rate and in the expected number of mutations to construct a CI. Because the numbers of mutations that accumulate on the branches of the tree are conditionally independent of one another and are nearly uncorrelated with the estimator of TAR, we estimated α and TAR independently (see Materials and Methods).

Leveraging data from an ∼45,000-year-old Siberian (Ust’-Ishim), we estimated that TAR = 275 kya (95% CI: 241–305 kya), and we estimated α by using two approaches that yielded similar results. In our first approach, we simply used the number of mutations shared by A00 and the reference (branch a of Figure 2) and the number of mutations unique to each (branches d and e) to estimate the relative times between splits. This method is insensitive to mutation-rate variability across the chromosome and led us to estimate α = 2.14 (95% CI: 1.64–2.89). In the second approach, we made use of the greater amount of data available for the denominator of the ratio and adjusted for mutation rate heterogeneity across the chromosome to estimate α = 1.82 (95% CI: 1.40–2.32). Because the main source of uncertainty is the limited sequence coverage for the Neandertal lineage, the CIs from the two approaches overlap substantially, but we prefer the first method, as it is simpler and potentially less biased. In both cases, we disregarded the number of variants unique to the Neandertal sequence (branch f) because this branch is enriched for false positives as a result of low coverage, DNA damage, and sequencing errors.

Combining the parametric bootstrap CIs of α and TAR, we estimated TNR = 588 kya (95% CI: 447–806 kya) with the first α estimate and TNR = 499 kya (95% CI: 375–656 kya) with the second.

Finally, we examined the potential functional relevance of the 146 mutations that differed among the Neandertal, A00, and reference sequences (Data S5). These included 11 non-synonymous changes and one nonsense mutation (Table 1). PolyPhen-2 predicted most missense mutations to have a benign effect, but it predicted possibly or probably damaging effects for Neandertal mutations in PCDH11Y (MIM: 400022) and USP9Y (MIM: 400005), for an A00 mutation in ZFY (MIM: 490000), and for a modern human mutation in KDM5D (MIM: 426000). The Neandertal nonsense mutation at codon 16 of TMSB4Y (MIM: 400017) might render its product non-functional, and MutationTaster predicts that it is probably deleterious.

Table 1

Protein-Changing Mutations

Coordinate	Gene	Lineage1	Substitution2	Effect	Tool	Function	MIM No.
2,844,774	ZFY	N	p.Val140Ala	B	P2	potential transcription factor	MIM: 490000
p.Val331Ala	B	P2
2,847,322	ZFY	A	p.Ile374Thr	B	P2	potential transcription factor	MIM: 490000
p.Ile488Thr	PrD	P2
p.Ile565Thr	B	P2
4,967,724	PCDH11Y3	N	p.Lys702Thr	B, B	P2, MT	protocadherin	MIM: 400022
5,605,569	PCDH11Y3	N	p.Ser1203Arg	PrD	P2	protocadherin	MIM: 400022
6,932,032	TBL1Y	N	p.Gly100Ala	B, B	P2, MT	–	MIM: 400033
14,832,610	USP9Y	N	p.Glu62Gly	PrD	P2	peptidase	MIM: 400005
14,832,620	USP9Y	R	p.Glu65Asp	B	P2	peptidase	MIM: 400005
14,838,553	USP9Y	N	p.Ala162Thr	B	P2	peptidase	MIM: 400005
15,816,262	TMSB4Y	N	p.Ser16∗	PrD	MT	actin sequestration	MIM: 400017
21,868,167	KDM5D	R, A	p.Arg1445Gln	B, B	P2, MT	demethylase	MIM: 426000
p.Arg1388Gln	B, B	P2, MT
p.Arg1476Gln	B, B	P2, MT
21,905,071	KDM5D	R, A	p.Ile69Val	PoD, B	P2, MT	demethylase	MIM: 426000
23,545,399	PRORY	A	p.Arg125Cys	–	–	–	–

Please see Data S5. Mutations variable among Neanderthal, A00, and the reference sequence, Data S6. Alignment (SAM) file with Neanderthal reads overlapping mutations in Table S3, Data S7. Neanderthal pileup lines for positions indicated in Table S3 for additional information on all mutations. Abbreviations are as follows: N, Neandertal; A, A00; R, reference; B, benign; PoD, possibly damaging; PrD, probably damaging; P2, PolyPhen-2 (ancestral to derived); MT, MutationTaster (reference to alternative).

Lineage(s) bearing the derived allele.

Multiple listings for a single coordinate reflect substitutions in different transcripts of the gene.

See Appendix B.

Discussion

We have estimated that the Neandertal Y chromosome from El Sidrón diverged from those of modern humans ∼590 kya, a value similar to TMRCA estimates for mtDNA sequences: 400 kya to 800 kya.11, 12 This time estimate and the genealogy we have inferred strongly support the notion that the most recent common ancestor of these Y chromosomes belonged to the population from which Neandertals and modern humans diverged, thereby refuting three alternative hypotheses. A priori, the Neandertal Y could have introgressed from a super-archaic population (Figure 3, scenario a), but this would have led to a far greater TMRCA estimate. Alternatively, it could have introgressed from the ancestors of modern humans after their divergence from Neandertals and prior to the most recent common ancestor of present-day Y chromosomes (scenario b) or from modern human populations subsequent to their migrations out of Africa (scenario c). We can also reject these hypotheses, as each requires a more recent split time.

Figure 3

Relationship of Neandertal Y Chromosome to Those of Modern Humans

The genealogy (red tree) can be parsimoniously explained as mirroring the population divergence (gray tree). We find no evidence for (a) a highly divergent super-archaic origin of the Neandertal Y chromosome, (b) ancient gene flow post-dating the population split, or (c) relatively recent introgression of a modern human Y chromosome into the Neandertal population.

The fact that the Neandertal Y-chromosome lineage we describe has never been observed in modern humans suggests that the lineage is most likely extinct. Although the Neandertal Y chromosome (and mtDNA) might have simply drifted out of the modern human gene pool, it is also possible that genetic incompatibilities contributed to their loss. In comparing the Neandertal lineage to those of modern humans, we identified four coding differences with predicted functional impacts, three missense and one nonsense (Table 1). Three mutations—within PCDH11Y, USP9Y, and TMSB4Y—are unique to the Neandertal lineage, and one, within KMD5D, is fixed in modern human sequences. The first gene, PCDH11Y, resides in the X-transposed region of the Y chromosome. Together with its X-chromosome homolog PCDH11X, it might play a role in brain lateralization and language development. The second gene, USP9Y, has been linked to ubiquitin-specific protease activity and might influence spermatogenesis. Expression of the third gene, TMSB4Y, might reduce cell proliferation in tumor cells, suggesting tumor suppressor function. Finally, the fourth gene, KDM5D, encodes a lysine-specific demethylase whose activity suppresses the invasiveness of some cancers.

Polypeptides from several Y-chromosome genes act as male-specific minor histocompatibility (H-Y) antigens that can elicit a maternal immune response during gestation. Such effects could be important drivers of secondary recurrent miscarriages and might play a role in the fraternal birth order effect of male sexual orientation. Interestingly, all three genes with potentially functional missense differences between the Neandertal and modern humans sequences are H-Y genes, including KDM5D, the first H-Y gene characterized. It is tempting to speculate that some of these mutations might have led to genetic incompatibilities between modern humans and Neandertals and to the consequent loss of Neandertal Y chromosomes in modern human populations. Indeed, reduced fertility or viability of hybrid offspring with Neandertal Y chromosomes is fully consistent with Haldane’s rule, which states that “when in the [first generation] offspring of two different animal races one sex is absent, rare, or sterile, that sex is the [heterogametic] sex.”