A discrete-beta model for testing gene flow after speciation

Due to the polymorphism in the ancestral species and the stochastic fluctuation of the coalescent process, the divergence time between sequences from two closely related species varies throughout the genome. This variation is described by an exponential distribution. Gene flow between the species causes additional variation in the sequence divergence time, beyond the expected variation from the ancestral coalescent process. This prediction can be used to test for gene flow between species using genomic sequence data. A previous implementation of the model used an empirical beta distribution to describe the variation in species divergence time across the genome, with numerical integration used to calculate the 3D integrals involved in the likelihood function. However, with a huge number of loci, the numerical integration is not accurate enough, so the false-positive rate of the likelihood ratio test may exceed the significance level. In this study, we replace the beta model with a discretized version, so that the likelihood calculation involves 2D integrals. We implement the new discrete-beta model in the program 3s and also extend our previous implementation to accommodate loci with different data configurations. We use computer simulation to examine the false-positive rate, the power and the robustness of the likelihood ratio test. The results show that the test had low false-positive rates and had even higher power in detecting gene flow than a likelihood ratio test based on a proper implementation of the isolation-with-migration model. It is, however, not robust to high levels of recombination within locus. Application of the test to two empirical data sets (from the hominoid genomes and from the Drosophila fruit flies) suggests the presence of gene flow for both data sets.