Biases in demographic modelling affect our understanding of the process of speciation

Estimating patterns of gene flow during the early stages of speciation is central to understanding whether reproductive isolation arises via the gradual erosion of gene flow or through successive stages of strict isolation and secondary contact. Such scenarios can be explicitly 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 scenarios of population divergence with constant gene flow, failure to account for population size changes in either a daughter population or the ancestral population leads to overestimated divergence time and to a bias towards the choice of SC models. 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 demographic history of North Sea and Baltic Sea turbots Schopthalmus maximus by testing 16 IM and 16 SC scenarios, modelling changes in effective population sizes as well as the effects of linked selection and heterogeneous migration rates across the genome. As in the simulated data, failure to account for changes in effective population sizes resulted in selecting SC models with a long period of strict isolation and divergence times preceding the formation of the Baltic Sea. In contrast, models accounting for population size changes suggest that the Baltic Sea turbot population originated from a very recent (