Biases in demographic modelling affect our understanding of recent divergence

Estimating patterns of gene flow during the early stages of divergence is central to understanding whether reproductive isolation arises through gradual erosion of gene flow or via successive stages of strict isolation and secondary contact. Such scenarios can be tested by comparing the joint allele frequency spectrum (jAFS) of a set of populations to jAFS simulated under scenarios of isolation with migration (IM) and secondary contact (SC). However, the potential effect of unaccounted demographic events (such as population expansions and bottlenecks) on model choice and parameter estimation remains largely unexplored. Using simulations, we demonstrate that under realistic divergence scenarios with constant gene flow, failure to account for population size (Ne) changes in daughter and ancestral populations leads to biases in divergence time estimates as well as model choice. On the other hand, when the simulations included long periods of strict isolation the correct gene flow scenario was usually retrieved. We illustrate these issues reconstructing the recent demographic history of North Sea and Baltic Sea turbots (Schopthalmus maximus) by testing 16 IM and 16 SC scenarios, modelling changes in Ne as well as the effects of linked selection and barrier loci. Failure to account for changes in Ne resulted in selecting SC models with long periods of isolation and divergence times preceding the formation of the Baltic Sea. In contrast, models accounting for Ne changes suggest that the Baltic Sea turbot population originated from a recent (