A new isolation with migration model along complete genomes infers very different divergence processes among closely related great ape species

Abstract

We present a hidden Markov model (HMM) for inferring gradual isolation between two populations during speciation, modelled as a time interval with restricted gene flow. The HMM describes the history of adjacent nucleotides in two genomic sequences, such that the nucleotides can be separated by recombination, can migrate between populations, or can coalesce at variable time points, all dependent on the parameters of the model, which are the effective population sizes, splitting times, recombination rate, and migration rate. We show by extensive simulations that the HMM can accurately infer all parameters except the recombination rate, which is biased downwards. Inference is robust to variation in the mutation rate and the recombination rate over the sequence and also robust to unknown phase of genomes unless they are very closely related. We provide a test for whether divergence is gradual or instantaneous, and we apply the model to three key divergence processes in great apes: (a) the bonobo and common chimpanzee, (b) the eastern and western gorilla, and (c) the Sumatran and Bornean orang-utan. We find that the bonobo and chimpanzee appear to have undergone a clear split, whereas the divergence processes of the gorilla and orang-utan species occurred over several hundred thousands years with gene flow stopping quite recently. We also apply the model to the Homo/Pan speciation event and find that the most likely scenario involves an extended period of gene flow during speciation.

Figure 1. Isolation-with-migration model.

Our isolation-with-migration model considers two separated populations (sub-species or species) derived from a shared ancestral population in the recent past. The model assumes that the ancestral population split into two populations in the past, at time , and that these two populations exchanged genes with migration rate until a later time, , where gene flow stopped. The coalescence process in this model is parameterized with a coalescence rate (inverse of the effective population size), , and a recombination rate, . The model is translated into a finite-state hidden Markov model by discretizing time into time intervals with break points .

Figure 2. Estimation accuracy.

The box-plot shows the distribution of parameter estimates for six different simulation scenarios. In all scenarios the coalescence rate and the recombination rate parameters are kept fixed, while the end of gene flow, , the initial population split, , and the migration rate, , varies between scenarios. For each simulation scenario, 10 independent data sets were generated and analyzed. The dashed horizontal lines indicate the simulated values for the five parameters. The recombination rate is consistently under-estimated while the remaining parameters are well recovered.

Figure 3. The effect of mutation rate variation.

The figure shows the effect on parameter estimation when the mutation rate is varied along the genome alignment. We split the alignment into segments geometrically distributed with mean length 500 bp and 2 kbp, and the mutation rate is then scaled by a random value chosen uniformly in the range 0.75 to 1.25 or 0.5 to 1.5. The dashed lines show the simulated values. The largest effect on varying the mutation rate is seen in the top-most parameters, the coalescence rate and the mutation rate. Varying the mutation rate increases the variance in coalescence times scaled with mutation rate, which is interpreted by the model as a decreased coalescence rate, while segments with low mutation rates are seen as more recent coalescence rates which the model interprets as evidence for migration. Consequently, variation in mutation rate decreases our estimates of the coalescence rate and increases our estimates of migration rates.

Figure 4. The effect of recombination rate variation.

The figure shows the effect on parameter estimation when the recombination rate is varied along the genome alignment. To simulate variation in the recombination rate, we sampled random 10 Mbp segments of the human genome, extracted the DeCODE recombination map for these segments, and scaled the recombination rate in the simulations according to the variation in this map. The dashed lines show the simulated values of the parameters. For most parameters, the effect of varying the recombination rate is seen as an increased variance in the estimates, while they do not appear to be biased. The exception is the recombination rate that becomes even more underestimated than for a constant recombination rate.

Figure 5. The effect of using a random genotype phase.

We simulated the situation where the genotype phase is unknown by simulating two genomes and selecting a random allele for all heterozygotic sites. The plot shows the effect on parameter estimates of not knowing the phase.

Figure 6. Split times estimates for the three great ape genera.

The box plot shows the estimated split times using either the isolation model or the isolation-with-migration model for the three great ape comparisons. The box plots on the left shows the split time estimate in the isolation model while the box plots on the right shows both the initial population divergence and the end of gene flow. The variation in estimates is from each 10 Mbp segment of the genome.

Figure 7. Chromosome wise split time estimates.

The box plots show the estimates of the initial split time and the end of gene flow in the isolation-with-migration model for each 10 Mbp segment for each chromosome.

Figure 8. Model comparison between the isolation and the isolation-with-migration model.

The box plots show the Akaike Information Criteria (AIC) for the isolation model against the isolation-with-migration model. For each 10 Mbp genomic segment we have plotted the AIC for the model including migration minus the model without. The model with the smallest AIC should be preferred, so values below zero prefers the isolation model while values above zero prefers the migration model.

Figure 9. Parameter estimates for the human/chimpanzee split with the isolation model.

The histograms show the distribution of parameter estimates for the human/bonobo speciation (blue) and the human/chimpanzee speciation (red) using the isolation model.

Figure 10. Parameter estimates for the human/chimpanzee split with the isolation-with-migration model.

The histograms show the distribution of parameter estimates for the human/bonobo speciation (blue) and the human/chimpanzee speciation (red) using the isolation-with-migration model.

Figure 11. Model comparison for the human/chimpanzee and the human/bonobo split.

The histograms show the distribution of AIC differences for the isolation and isolation-with-migration model for the human/chimpanzee comparison and the human/bonobo comparison. Negative values indicate a preference for the isolation model while positive values indicate a preference for the isolation-with-migration model. The overall result points toward a preference for a prolonged speciation for the Homo/Pan split.

Figure 12. Split times scaled in years.

The figure shows the inferred split times when scaled with a mutation rate, , ranging from to . The solid lines show mean estimates while the dashed lines the 95% confidence interval (SEM). The Homo/Pan slit is annotated with key fossils, the chimpanzee/bonobo split with the formation of the Congo River, and the orang-utan split with glacial period where sea level was low and migration between orang-utans possible.

Figure 13. Ancestral recombination graph and state space.

On the left is shown an ancestral recombination graph for two genomes with two nucleotides. Lineages, in the notation we use for constructing the CTMCs, are shown in red. On the right is shown the corresponding list of transitions int he CTMC with the type of transitions on the arrows: recombination (R), migration (M) and coalescence (C). The transition from the two separate populations to the ancestral population is a special transition – the projection matrix in the CTMC – shown in red.