Modeling the Multiple Facets of Speciation-with-Gene-Flow toward Inferring the Divergence History of Lake Whitefish Species Pairs (Coregonus clupeaformis)

Abstract

Parallel divergence across replicated species pairs occurring in similar environmental contrasts may arise through distinct evolutionary scenarios. Deciphering whether such parallelism actually reflects repeated parallel divergence driven by divergent selection or a single divergence event with subsequent gene flow needs to be ascertained. Reconstructing historical gene flow is therefore of fundamental interest to understand how demography and selection jointly shaped genomic divergence during speciation. Here, we use an extended modeling framework to explore the multiple facets of speciation-with-gene-flow with demo-genetic divergence models that capture both temporal and genomic variation in effective population size and migration rate. We investigate the divergence history of replicate sympatric species pairs of Lake Whitefish (normal benthic and dwarf limnetic) characterized by variable degrees of ecological divergence and reproductive isolation. Genome-wide SNPs were used to document the extent of genetic differentiation in each species pair, and 26 divergence models were fitted and compared with the unfolded joint allele frequency spectrum of each pair. We found evidence that a recent (circa 3,000-4,000 generations) asymmetrical secondary contact between expanding postglacial populations has accompanied Whitefish diversification. Our results suggest that heterogeneous genomic differentiation has emerged through the combined effects of linked selection generating variable rates of lineage sorting across the genome during geographical isolation, and heterogeneous introgression eroding divergence at different rates across the genome upon secondary contact. This study thus provides a new retrospective insight into the historical demographic and selective processes that shaped a continuum of divergence associated with ecological speciation.

Keywords: Coregonus; JAFS; demographic inference; ecological speciation; population genomics; speciation-with-gene-flow.

Fig. 1.

—Geographic locations of the lakes where sympatric whitefish species pairs were sampled, and overview of the extent of shared versus private polymorphism. Pie charts illustrate the amount of shared and private SNPs among lakes as well as between species among or within lakes.

Fig. 2.

—The 26 models implemented in this study. The models implemented here represent extensions of the four classical models of divergence: “Strict Isolation” (-SI), “Isolation with Migration” (-IM), “Ancient Migration” (-AM) and “Secondary Contact” (-SC). Briefly, T_S corresponds to the duration of complete isolation between diverging populations and T_AM and T_SC correspond to the duration of gene-flow in AM and SC models, respectively. The first extension of these four models accounts for temporal variation in effective populations size (-G models), allowing independent expansion/contraction of the diverging populations. The last categories correspond to “Heterogeneous gene flow” models, which integrate parameters allowing genomic variations in effective migration rate (-2 m), effective population size (-2 N) or both (-2m2N) to account for genetic barriers and selection at linked sites.

Fig. 3.

—Historical demography of the Lake Whitefish species pairs. (A) Observed joint allele frequency spectrum (JAFS) for normal (-N; x axis) and dwarf (-D; y axis) populations for each of four lakes (T = Témiscouata, E = East, I = Indian, and C = Cliff), obtained by projection of empirical data to 13 diploid individuals per population. For each JAFS, the color scale indicates the number of SNPs falling in each bin defined by a unique combination of the number of derived allele observed in normal and dwarf populations. (B) Predicted JAFS of the fittest model for each lake. (C) Representation of the fittest model for each lake.

Fig. 4.

—Model comparisons. Barplots showing the effect of taking into account a particular demographic or selective aspect in the models, assessed using model scores, with (A) the effect of including temporal variation in population effective size (-G), (B) heterogeneous migration rates among loci (-2 m) and (C) heterogeneous effective population sizes among loci (-2 N). The vertical bars indicate the variance of model scores within a given category of models and asterisks represent significant differences in average model scores between the compared categories of models. (D) Heat-map of the weighted AIC ($w_{AIC}$) showing the relative weights of the 26 models for each lake. The color scale corresponds to $w_{AIC}$ values ranging from 0 to 1. Warmer colors indicate the best models.

Fig. 5.

—Asymmetrical effective gene flow between normal and dwarf whitefish within lakes. Bar plot of the effective number of migrants per generation in both directions from dwarf to normal (black) and reciprocally (gray), obtained from estimated parameters of the fittest model for each lake, using the average gene flow formula (average gene flow=N×b × (P×m_e + (1−P)×m_e’) in each direction.

Fig. 6.

—Genetic structure and relationships among lakes and species. (A) Discriminant analysis of principal components (dAPC) of the different lakes (Cliff, Indian, Webster, East, Témiscouata) for either Dwarf or Normal whitefish (D or N), representing 3D relationships among populations. The first axis (LD1, 39.5%) captures the geographical signal of differentiation among lakes. The third dAPC axis (LD3, 16%) separates species pairs according to their residual genetic promixity to ancient glacial lineages represented by the least introgressed populations from Cliff Lake. Positive coordinates represent populations with high proportions of Atlantic ancestry whereas negative coordinates reflect increased proportions of Acadian ancestry. The fourth axis (LD4, 11%) tends to separate species pairs within each lake. The second dAPC axis LD2 is not shown here to avoid partial redundancy with LD1, but is provided as a Supplementary Material online. (B) Boxplot of individual coordinates along the fourth dAPC axis (LD4), highlighting the divergence parallelism between species among lakes.

Fig. 7.

—Shared ancestral genetic variation between allopatric populations and admixture events between sympatric species pairs. The least introgressed normal whitefish population from Cliff Lake was used to root the tree. Horizontal branch lengths are proportional to the amount of genetic drift in each branch, and the scale bar indicates 10 times the average standard error (SE) of the entries in the covariance matrix between pairs of populations. Color scale indicates the weight of inferred migration events.