Different historical generation intervals in human populations inferred from Neanderthal fragment lengths and mutation signatures

After the main Out-of-Africa event, humans interbred with Neanderthals leaving 1–2% of Neanderthal DNA scattered in small fragments in all non-African genomes today. Here we investigate what can be learned about human demographic processes from the size distribution of these fragments. We observe differences in fragment length across Eurasia with 12% longer fragments in East Asians than West Eurasians. Comparisons between extant populations with ancient samples show that these differences are caused by different rates of decay in length by recombination since the Neanderthal admixture. In concordance, we observe a strong correlation between the average fragment length and the mutation accumulation, similar to what is expected by changing the ages at reproduction as estimated from trio studies. Altogether, our results suggest differences in the generation interval across Eurasia, by up 10–20%, over the past 40,000 years. We use sex-specific mutation signatures to infer whether these changes were driven by shifts in either male or female age at reproduction, or both. We also find that previously reported variation in the mutational spectrum may be largely explained by changes to the generation interval. We conclude that Neanderthal fragment lengths provide unique insight into differences among human populations over recent history.


S1 -Identification of archaic fragments in non-African individuals extant populations and ancient samples
We call archaic fragments in the samples of the Simons Genome Diversity Project (SGDP) 1 and ancient samples analysed in this study as described in Skov et al. 2020 and 2018 2,3 -a step by step tutorial is also available at https://github.com/LauritsSkov/Introgression-detection.
The method is described generally in the Methods section. In this section, we describe the specifics of the pipeline used in this study.
Outgroup variants set, window mutation rate and callability and derived allele polarization for the SGDP dataset To generate the set of variants seen in the outgroup, we merged all variants from the following populations: To generate the callability regions, we merged the following files:

S2 -Archaic fragment length gradient around the world is consistent with using other quantifications and filters
We studied the robustness of the difference in mean archaic fragment length among the 5 geographical groups studied applying multiple filters.

1) Median instead of mean
The mean is very sensitive to outliers. In our case, very long archaic fragments, for example in East Asians, could increase the mean and thus show an unrealistic pattern among regions.
To avoid that, we use median instead because it is more robust to outliers.

2) Vindija genome-like fragments
The method used in this study is able to find archaic fragments whose variation is not fully captured by the sequenced archaic individuals 2 . The difference in archaic fragment length can potentially be affected if there is a distinct archaic content among the extant populations studied here -for example, a greater and more recent Denisova component in Asia 11 .
It is known that the majority of the archaic component in Eurasia and America is from a Neanderthal population closely related to the Vindija genome 12 . Thus, we restrict fragments used in this analysis to share more variation with the Vindija Neanderthal genome than the Altai Neanderthal genome or the Denisovan genome.

3) High confidence archaic fragments
The method used in this study, returns the archaic fragments found in a genome with an associated mean posterior probability. We restricted archaic fragments compared to be of a high confidence (mean posterior probability >= 0.9).
When we study the archaic fragment difference among individuals in Eurasia and America applying the three different types of filters explained above, we can see that the pattern observed using all fragments holds (Supplementary Figure 2). We conclude that the difference in archaic fragment length is genuine and not depending on the factors exposed above. Median Archaic Fragment Length  High−confidence fragments Asians and West Eurasians.   Supplementary Table 4. Individual values are shown as dots, coloured depending on the region they belong to.

S5 -Differences on archaic fragments between West Eurasia and East Asia regions are replicated in the population level comparison
In our analysis, we divide the non-African individuals of the SGDP data into 5 main regions to compare them in terms of archaic fragment number, length and archaic sequence. In this section, we use complementary data to investigate whether similar differences are found if we use larger samples from more homogeneous populations to test for differences in archaic  Table 5). We used four populations that satisfy the following criteria:  Table 5, Supplementary Figure 6). This indicates that the regional variance observed in the SGDP data is likely to stem primarily from intra-population variance, rather than inter-population variance.
We then compare, for each SGDP region, if the two representative populations from the HGDP data set have distinctive distributions of archaic fragment length (Supplementary Figure 6). Eurasia and East Asia is replicated in an independent data set with more homogeneous populations.

Archaic fragments call in HGDP populations
To call archaic fragments in the HGDP data we follow the methodology described in the Methods section and in S1. However, since the HGDP data is mapped to the GRCh38 reference genome, we modify certain steps to create the source files.
First, we use the following individuals to generate our outgroup: To polarize alleles into ancestral and derived alleles we used the field "AA_ensembl" in the HGDP VCFs, which corresponds to the Ensembl's homo_sapiens_ancestor_GRCh38_e86 files.   There are also instances in which an archaic fragment does not share more SNPs with one of the archaic genomes but multiple, so we can't classify the affinity of the fragments; these fragments are called ambiguous fragments.

1) Denisova fragments
We only include archaic fragments which share more variants to Denisova genome than any of the two Neanderthal genomes.

2) nonDenisova fragments
In this analysis we exclude fragments used above from all the fragments. Thus, we include Vindija-like, Altai-like, ambiguous and unknown.

3) Neanderthal fragments
We only include archaic fragments that share more variants with either the Altai Neanderthal or the Vindija Neanderthal genomes than the Denisova genome. Neanderthal ambiguous fragments, fragments that share the same number of SNPs with Vindija or Altai but this number is higher than what is shared with the Denisova, are also included.
All results for the different filters are shown in Supplementary Table 8. The Denisova content is 3 times greater in East Asia than in West Eurasia (Denisova fragments filter). When this unequal component is removed (non-Denisova fragments filter), we can see that the collapsed archaic sequence is very similar between the two regions.
The analysis was repeated with fragments that share more variation with Neanderthal than with Denisova (Neanderthal fragments). In this case, we observe a 1.07 fold higher Neanderthal content in the East Asian group. We attribute this to the fact that since West Eurasia archaic fragments tend to be shorter, they do not contain enough SNPs to classify them to the category that they belong to. Thus, they are going to be more often classified as unknown compared to fragments in East Asia. Furthermore, the 2 method has higher false negative rate with short fragments, which will artificially decrease the total number of fragments in that region.

S7 -Simulations support a single Neanderthal pulse to the ancestors of East Asia and West Eurasia
In this study, we observe that East Asia individuals have longer and more archaic fragments compared to the rest of the world. This could be compatible with East Asians receiving archaic fragments from a much more recent Neanderthal admixture private to East Asia [19][20][21] . However, when East Asia and West Eurasia regions are compared by joining their fragments in each group, they seem to have similar amounts of total and shared archaic sequence.
In order to quantify the expected differences between populations in a scenario with a single When we join fragments of both groups, similar to what we do for the SGDP data (Methods) and described in the Methods section, we observe similar abundance of shared archaic sequences for both groups in the One Pulse scenario (Supplementary Figure 11). However, in the Two Pulses scenario, there is a large decrease in the shared archaic sequence for the EA group that received the second gene flow.
Overall, we observe that the Two Pulses scenario, compared to the null model, provides EA group with an excess of private sequences (~58%, 8 points more than WE) and slightly increases the mean fragment length by ~2.5% compared to WE group. In contrast, the SGDP samples from West Eurasia and East Asia show greater differences in terms of fragment length (~12%) and similar levels of total and shared sequence (53% and 56% respectively) as groups. Thus, we conclude that a second pulse scenario is not compatible with our observations.   Mean values for each distribution are shown as squares with their associated 95% CI as whiskers (computed as mean±(1.96*se); se = sd/sqrt(n)). Comparing the means of each distribution replicate, the ratio between EA and WE is shown (second subpanel rows) as black dots. a) Mean Archaic Fragment Length (bp) b) Number of archaic fragments c) Archaic sequence (bp).   Figure 11. Joined archaic fragments comparisons between WE and EA for the One Pulse and Two Pulses simulated scenarios. a) Joined archaic fragment length for each group (color coded) as bar plots in the 10 replicates (x-axis) of each simulated scenario (subpanels). Each bar is divided into shared (plain colour) and private sequence (transparent colour). b) For each replicate in a), the percentage of shared sequence between the EA and WE are shown as points per population (colour coded).

S8 -Derived alleles call outside regions with evidence of archaic introgression and acquired after the Out-of-Africa in SGDP samples
We retrieved the genotypes of all polymorphic loci for each individual in the 5 main regions and African samples as explained in the Methods section. We masked repetitive regions and regions of the genome in which there is some evidence of archaic introgression in the following way:

1) Neandertal introgressed regions
Neanderthals had a different mutation profile than modern humans 3 . Thus, differences in Neanderthal content per individual could influence those analyses that explore the mutation spectrum differences among populations. Also, by removing these regions, we will base the mutation analysis on regions of the genome that we haven't explored in the archaic fragment length part of the study. Thus, the tests are going to be independent of each other.
To do that, we disregarded any polymorphism localized in a region with evidence of archaic introgression in any of the individuals analyzed in this study (Methods, S1). For that, we joined all archaic fragments called in any individual included in this study using this command: bedtools merge -i ind1.bed ind2.bed … indN.bed > joined.bed where N denotes the total number of individuals.
In total, the joined archaic region adds up to 1,632,776,000 bp.

2) Repeats
We also excluded repetitive regions in which sequencing errors are expected to be more prevalent. For that, we downloaded the human reference genome by using the following command: Homozygous locus for the derived allele count as 2 mutations and heterozygous sites count as 1 for a given individual. The distribution of derived allele accumulation per region is shown in Figure 3 and the mean derived allele accumulation counts per region are provided in Supplementary Table 9.
Finally, we classified loci in different mutation types depending on the derived allele nucleotide, the ancestral allele nucleotide and their 5' and 3' nucleotide context. For example, as shown by the diagram below, a derived allele T that had an ancestral allele C with the context G and A (5' and 3' respectively) would be denoted as GCA>T. Because we do not make distinction of the strand in which the mutation occurred, we collapsed strand-symmetric mutations. This is the same as saying that GCA>T is equivalent to TGC>A. This way, we end up with 96 mutation types.

S9 -Estimation of the different parental generation time in West Eurasia and East Asia
As described in the main text, West Eurasia individuals have accumulated 1.09% more derived alleles than East Asians since the split with Africans (Out-of-Africa). Because we are only interested in the proportion of derived alleles accumulated after the split of West Eurasians and East Asians, we need to correct for the span of time since the Out-of-Africa event until the split of the two Eurasian populations (Supplementary Figure 13). Thus, we need to assume dates for the split between Africans and non-Africans and the split between Eurasians.
We note that in the literature dating the Out-of-Africa is widely discussed and controversial,

S10 -Mutation spectrum correlation with mean parental age and potential bias due to difference in mean paternal and maternal age in de deCODE dataset
The germline mutation spectrum is dependent on the parental sex and age at conception 30 . In this study, we observe differences in the abundance of derived alleles accumulated after the Out-of-Africa event when stratified by mutation type (Supplementary Figure 12, Supplementary Table 10). Here, we study to which extent these differences can be explained by changes in generation time in the 5 regions. For that, we compare the mutational patterns of de novo mutations (DNM) depending on parental age in trio studies 30,31 (deCODE data set) with the differences in mutation spectrum of extant populations with the mean archaic fragment length as a proxy of mean generation time (SGDP data set) as explained in the Methods section.
The estimates of the obtained linear models for each mutation type are given in Supplementary Table 11. Moreover, the correlations between the slopes of both data sets is shown in Supplementary Figure 14.
The probands of the deCODE data set have a bias towards fathers being older than mothers, with a mean of 2.77 years and the largest difference of more than 40 years (Supplementary Figure 15a). To study if the correlation of the mutation spectrum with the mean parental age is affected by the mentioned bias, we rerun the correlation test with the deCODE data set with only probands that have parents with an age difference of less than 4 years. This way, we retaining more than 50% of the data (Supplementary Figure 15b) and reduce the bias (mean = 0.94 differences in years, Supplementary Figure 15c). We then compared the slopes of the linear models calculated in the original deCODE data set and when we impose the parental age difference filter explained above (Supplementary Figure   16). We don't observe qualitative changes in the slopes when comparing the two and thus, we used all probands for our analysis.
Supplementary Figure 14. Slope coefficient correlation between SGDP data and deCODE data linear models. Dot plot graph illustrating the correlation between linear model slope coefficients derived from the SGDP data (x-axis) and deCODE data (y-axis) for each mutation type (color code   Figure 16. Slope coefficient correlation between linear models of deCODE data and deCODE data when only using probands with parental age difference less than 4 years. Dot plot graph illustrates the correlation between linear model slope coefficients derived from the deCODE data (y-axis) and the deCODE data when only using probands with parental age difference less than 4 (x-axis) for each mutation type (color code). 95%CI for each estimate are shown as error bars. The sample sizes of individuals for which summary statistics are derived from, together with other statistics, are indicated in Supplementary Table 11 and S10. The 1-to-1 correspondence is denoted by the gray dashed diagonal line.   Table 11. Linear models between mutation type fraction and mean generation time estimate in the SGDP and deCODE data sets. Two separate tables are given for the intercept and the slope of the linear models obtained using the R function lm().
For each mutation type and data set, the coefficients estimate, the SE, the t-test t value and the associated P value are provided.

X-to-A ratio
Due to the inheritance pattern of the X chromosome -2 copies transmitted in females while only 1 in males -compared to autosomes -2 copies in both females and males -, it is expected that the X chromosome has ¾ the diversity of the autosomes. However, this can be altered if the mutation rate changes disproportionately between females and males due to shifts in generation time between sexes. For example, an increase in the male mean generation time will decrease the yearly mutation rate in males and thus, proportionally less mutations are going to be accumulated in autosomes compared to the X chromosomes 32 . Therefore, the ratio of derived allele accumulation between the X chromosome and the autosomes will reflect variation on the generation time between males and females: higher values of the X-to-A ratio will be indicative of longer generation times in males compared to females and vice versa.
Although here we only consider generation time differences to affect the ratio, there are other factors that can perturbe this ratio such as reproductive variance between sexes 33 , demographic events 34 or differences in selection 35 .
We computed the X-to-A ratio as explained in the Methods section and we then correlated the ratio with the mean archaic fragment length for each individual (Figure 4a).

C>G maternaly enriched regions
As described in Jónsson et al. 2017 30 , there are regions of the genome in which DNM are clustered (cDNM). Those regions appear to be enriched in C>G mutations which originate in the maternal lineage. They also show that these clusters increase in number more rapidly with maternal than paternal age at conception.
Here we explore if there is a difference in the number of C>G segregating sites in cDNM genomic windows among the 5 regions.
To compute the C>G ratio between DNM cluster regions and the rest of the genome, we follow the procedure explained in the Methods section.
If the ratio = 1, it indicates that the C>G enrichment is similar in cDNM regions compared to the rest of the genome. If > 1, then there is an excess and if < 1, then there is a depletion.
Nonetheless, we are not interested in the actual ratio, but the comparison among geographical regions on this quantity. We then correlated the ratio with the mean archaic fragment length for each individual obtained in this study (Figure 4b).
Y chromosome Male individuals with shorter generation time are predicted to increase the mutation rate per year. Thus, Y chromosomes are expected to accumulate more derived alleles in individuals with a historically shorter mean generation time compared to others with longer ones.
To investigate that, we followed a similar procedure as explained in the Methods and in S8 to classify derived alleles into mutation types, changing certain steps and filters listed below: The accumulation of derived alleles in the Y chromosome per geographical region is shown in Supplementary Figure 17 (included in Data2_mutationspectrum.txt) and in Supplementary Table   13.