Phylogenetic Branch Lengths Are Variable among Y-Chromosome Haplogroups

To explore the reported pattern of branch length heterogeneity (Hallast et al. 2015), we analyzed the high-coverage Y-chromosome sequencing data of the 1000 Genomes Project (“the 1KG data set”). This data set contains 1,195 males from global populations (1000 Genomes Project Consortium 2015; Poznik et al. 2016), and has sufficient representation (more than ten samples) of 20 Y-chromosome haplogroups. We reconstructed a Y-chromosome phylogeny (fig. 1A), and rejected the hypothesis of evolutionary rate homogeneity throughout the tree (Kumar et al. 2012, 2018; Stecher et al. 2020). We subsampled data with one sample per haplogroup and uncovered significant interhaplogroup heterogeneity in evolutionary rate (likelihood-ratio test, P = 7.91 × 10⁻⁴¹), consistent with branch length variation (fig. 1A). These discoveries were robust to subsampling of individuals and nucleotide sites (supplementary note 1, Supplementary Material online). As a surrogate for branch length, we calculated normalized branch-specific evolutionary rate of terminal branches in the phylogeny using the RelTime method (Tamura et al. 2012). This method showed that the evolutionary rate varied greatly among haplogroups (fig. 1B), and the pattern of variation was consistent with figure 1A and previous reports (Hallast et al. 2015) (e.g., haplogroup R1b had shorter branch lengths than average).

Fig. 1.

Y-chromosome phylogenetic branch length heterogeneity. (A and C) Y-chromosome phylogeny reconstructed using the 1KG (A) and HGDP (C) data sets. The phylogenies were reconstructed using the maximum-likelihood method in MEGA and visualized with FigTree. Haplogroups with more than ten samples were included, along with a sample from haplogroup A. Although all samples were collected at the present time, their branch lengths vary substantially, for example, haplogroup R1b samples have much shorter branch lengths (dashed circle) than samples from most other haplogroups. (B and D) Distribution of branch length (branch-specific evolution rate) by haplogroup in the 1KG (B) and HGDP (D) data sets. Evolutionary rates were estimated using the RelTime method and normalized to a mean of 1. Outlier samples (three in 1KG and two in HGDP) were removed. Substantial interhaplogroup evolutionary rate variation was observed, which was consistent between data sets (ρ = 0.81) and with (A) and (C).

To critically assess these findings, we analyzed the Y-chromosome sequencing data from the Human Genome Diversity Project (“the HGDP data set”) (Bergström et al. 2019). This data set contains 554 male samples and has sufficient representation of 16 Y-chromosome haplogroups. A phylogeny was reconstructed (fig. 1C), and the hypothesis of evolutionary rate homogeneity was rejected, even when using one sample per haplogroup (P = 3.93 × 10⁻¹⁴⁸) and regardless of the subsampling (supplementary note 1, Supplementary Material online). We calculated branch-specific evolutionary rate as described above. Reassuringly, the pattern of variation observed in the HGDP data set (fig. 1D) was very similar to that of the 1KG data set. Taken together, we found that interhaplogroup branch length heterogeneity in human Y-chromosome trees does exist, whereas its cause requires further investigation.

Within-Cell-Line Y-Chromosome Mutations Are Informative for Germline Mutational Processes

Previous studies suggested that branch length heterogeneity might be caused by nongenetic factors, for example, paternal age variation across populations, acting over many generations (Jobling and Tyler-Smith 2017). Another possibility is variation in mutation rate among Y-chromosome haplogroups. We propose to use “within-cell-line mutations” on the Y chromosome to distinguish the two possible causes. Within-cell-line mutations denote sites polymorphic in a given cell line found in the nonpseudoautosomal regions of the sex chromosomes in males, which indicates the presence of an alternative allele in a subset of cells. These mutations can only occur somatically along the hematopoietic lineage or in vitro during culturing of LCLs, and therefore are unlikely to be affected by nongenetic (e.g., cultural, demographic, and environmental) factors, including paternal age.

We identified within-cell-line mutations in the 1KG and HGDP data sets, requiring at least three reads supporting the alternative allele. We applied filters to remove false positives (supplementary Materials and Methods, Supplementary Material online). We resequenced the Y chromosome of five 1KG samples to ∼8-fold coverage, and successfully validated 92.9% (13/14) of within-cell-line mutations for which the alternative alleles were supported by five or more reads in the original samples. These mutations were selected for validation because the fraction of cells carrying the alternative alleles was likely higher, and thus more likely to also be present in the separately obtained samples used for resequencing.

DNA replication timing is correlated with the spatial distribution of mutations in humans (Stamatoyannopoulos et al. 2009; Koren et al. 2012), suggesting that replication timing can be used as a proxy for mutational processes. However, replication timing of the human Y chromosome has not been previously explored. We inferred a consensus Y-chromosome replication profile in LCLs using a method we described previously (Koren et al. 2014; Ding et al. 2020) (fig. 2A). Briefly, we inferred Y-chromosome replication timing for each cell line based on fluctuation of read depth across the chromosome. We averaged replication timing profiles for all samples to generate the consensus (supplementary Materials and Methods, Supplementary Material online). We uncovered a negative relationship between replication timing and within-cell-line mutation density (fig. 2A [lower panel] and B). Note that rigorously establishing statistical significance of this trend entails compensation for two sequential data sets that have distinctly different patterns of autocorrelation along the genome sequence. It is probably more convincing to note that the negative trend is robust to different approaches (supplementary note 2, Supplementary Material online).

Fig. 2.

Y-chromosome replication timing and its relationship with mutations. (A) Top panel: replication timing profile of 1KG samples (green). Red: consensus. Bottom panel: mutation density tracks. Red: smoothed. (B and C) Y-chromosome replication timing shows a negative relationship with within-cell-line (B) and germline mutations (C). Replication timing windows were assigned to one of 25 bins based on timing. Mutation counts and timing of windows within each bin were averaged and plotted. Red, linear trend line.

Germline mutations (using SNPs in the 1KG data set) also showed negative relationship with LCL replication timing (fig. 2A [lower panel] and C). Germline and within-cell-line mutations have similar mutational spectrum (supplementary note 2, Supplementary Material online). These observations indicate that although mutational processes are known to differ between somatic cells (LCLs) and the germline, mutations in LCLs could nevertheless provide valuable insights into germline mutational processes.

Substantial Variation in Y-Chromosome Somatic Mutation Rate across Haplogroups

As reasoned above, within-cell-line Y-chromosome mutations can be a valuable tool in distinguishing hypotheses regarding branch length heterogeneity. For a given sample, we define “relative Y-chromosome somatic mutation rate” as the counts of within-cell-line mutations on the Y chromosome (calculated as genetic distance using the Kimura two-parameter model) divided by that of the X chromosome. X-chromosome within-cell-line mutations were utilized to correct for non-Y-chromosome-specific factors affecting mutation accumulation, for example, variation in donor age and cell culturing.

We observed significant variation in the Y-chromosome somatic mutation rate among haplogroups using the 1KG data set (fig. 3A, P = 5.30 × 10⁻¹⁰, Kruskal–Wallis rank-sum test). Strikingly, the interhaplogroup variation in somatic mutation rate and phylogenetic branch length were positively correlated (Spearman ρ = 0.54, P = 1.54 × 10⁻², fig. 3B), suggesting that similar variation in mutation rate is likely also present in the germline. For example, haplogroups E and R had the shortest phylogenetic branch lengths and the lowest somatic mutation rates. Consistent with the findings in the 1KG data set, significant interhaplogroup variation in the somatic mutation rate (P = 3.52 × 10⁻⁹, fig. 3C) was observed in the HGDP data set, which was positively correlated with branch length heterogeneity (ρ = 0.78, P = 2.62 × 10⁻³, fig. 3D). These findings were robust to analytical approaches (supplementary note 3, Supplementary Material online).

Fig. 3.

Substantial Y-chromosome somatic mutation rate variation across haplogroups. (A and C) Distribution of the Y-chromosome somatic (within-cell-line) mutation rate among haplogroups in the 1KG (A) and HGDP (C) data sets. In both data sets, the mutation rate was significantly variable across haplogroups (see also, supplementary fig. S10, Supplementary Material online). (B and D) Variation in somatic mutation rate is correlated with branch length heterogeneity in the 1KG (B) and HGDP (D) data sets, suggesting that interhaplogroup mutation rate variation is a parsimonious explanation for branch length heterogeneity.

The extent of interhaplogroup variation in somatic mutation rate was substantial. The haplogroups with the highest and lowest median somatic mutation rate have an 83.3% difference in the HGDP data set (48.6% in the 1KG data set). There was no clear pattern that the mutation rates cluster by phylogeny, suggesting multiple events, at different time scales, affected Y-chromosome mutation rate. For example, assuming that higher mutation rate is ancestral, there were likely multiple slowdown events which occurred independently in the ancestors of haplogroups E and R. Our conclusions were unlikely driven by batch effects (supplementary note 4, Supplementary Material online). In summary, our findings indicate that there is substantial interhaplogroup variation in Y-chromosome mutation rate, and that such variation is a parsimonious explanation for phylogenetic branch length heterogeneity.

This all begs the question—why do Y haplogroups have different rates of mutation? What is the mechanism driving this variation? One plausible mechanism is variation in replication timing across haplogroups. We previously reported interindividual variation in replication timing, and showed that it was associated with genetic variants (rtQTLs) (Koren et al. 2014; Ding et al. 2020). Since replication timing is negatively linked to mutation rate, haplogroups associated with earlier replication may have a lower mutation rate. Y chromosomes have highly variable amounts of heterochromatin (Repping et al. 2006), which is usually late-replicating. One hypothesis would be that Y chromosomes with the greatest abundance of heterochromatin are also the latest replicating. But heterochromatin abundance changes rapidly, and it is not well correlated with haplogroup (Repping et al. 2006), so there are probably other factors, as yet unknown, that account for the among-haplogroup heterogeneity in replication timing and/or mutation rate.

In this study, we generated the first human Y-chromosome replication timing profile (supplementary file 1, Supplementary Material online), which provides quantification of regional mutation rate variation. We also revealed that the Y-chromosome mutation rate varies across haplogroups. This has important implications for evolutionary genetics, since inferences of divergence times could possibly be distorted by variation in mutation rate across haplogroups (supplementary note 5, Supplementary Material online).