Founder effects drive the genetic structure of passively dispersed aquatic invertebrates

Genetic submodel

All individuals are considered to have n neutral loci and n loci under selection. All loci are biallelic and no mutation is assumed. Absence of linkage is assumed among neutral loci and among loci under selection. Contrarily, the model accounts for physical linkage between each neutral locus and a corresponding selected locus, assuming several recombination levels, including absence of linkage. Loci under selection act additively on growth rate. Consequently, no dominance and no epistatic effects are assumed. Local adaptation requires a genotype-environment interaction on fitness, which is modeled through δ_i,j,l, which is the effect on the intrinsic growth rate (see below) of allele i (i: 1, 2) at locus j (j: 1,…, n) in locality l (l: 1, 2). The assumptions are (1) δ_1,j,1 = δ2, j, 2 , and (2) δ_i,j,l =  − δ_j≠i,j,l; so, in the case of homozygotes for a given selected locus, they will experience an increase or decrease of their growth rate by —2 δ— depending on the locality. Hence, the growth rate for each genotype g in each locality l (r_g,l) can be decomposed into r (basal growth rate) and θ (deviation of each genotype), so that ${\displaystyle r_{g,l}=r+{\theta }_{g,l}}$

where g is the genotype, l is the locality, and θ_g,l is the summation of the fitness components (δ) in locality l of the alleles carried by a genotype g in the n loci under selection. Thus, in any given locality, the growth rate during the asexual reproduction will vary between the limits r ± 2nδ.

Sexual reproduction is assumed to be panmictic and, for simplicity, is considered to be synchronic and at the end of the growing season (t = τ). As linkage disequilibrium can occur due to selection and genetic drift, gametic frequencies are computed. Gametes are then drawn to produce the diapausing eggs.

Genetic distance between populations was estimated based on neutral loci as ${\displaystyle F_{\mathrm{ST}}=\frac{\overline{H_T}-\overline{H_S}}{\overline{H_T}},}$

where $\overline{H_T}$ is the average expected heterozygosity for the two populations considered as a single one for the neutral loci, and $\overline{H_S}$ the average of the mean expected heterozygosity within each populations for the neutral loci (Hedrick, 2011). Allelic frequencies for each locus were computed using the total number of alleles. Similarly, a genetic distance for loci under selection (F_STQ) was computed (Le Corre & Kremer, 2012). F_ST and F_STQ values were obtained just after hatching of diapausing eggs.

Population growth

The asexual phase spans from time t = 0 to τ, which is the moment when sexual reproduction takes place. During the asexual phase the population grows deterministically according to a logistic growth model: ${\displaystyle \frac{d{N}_{l,g}}{dt}=N_{l,g}{r}_{l,g}\left[1-\frac{{\Sigma }_g{N}_{l,g}}{k}\right]}$

where N_l,g is the density of the genotype g in the locality i, r_l,g is its intrinsic population growth rate during the asexual phase, and K the carrying capacity. Note that K is genotype-independent. At the onset of each asexual growth season (t = 0), N_l,g is the sum of the hatched diapausing eggs, a fraction of them having been locally produced H_l,g, and the rest being migrants M_l,g.

At t = τ of the sexual generation y, the number of diapausing eggs produced P_l,g (y) is stochastically computed from N_l,g(τ, y) assuming a sexual proportion m (fraction of the females that becomes sexual), a sex ratio sr and an effective fecundity e (number of diapausing eggs produced per sexual female).

Mortality of diapausing eggs in the sediment with egg bank was assumed to be age-independent (annual survival rate γ). Empirical information supporting this assumption for field populations is not available. However, our model can account for fast senescence when it assumes the absence of egg bank. When a new planktonic growing season starts (t = 0) a fraction λ of the diapausing eggs in the sediment hatches.

Source population and local population founding

The two populations are founded at time y = t = 0 by F diapausing eggs randomly drawn from a single source population. The source population is assumed to be in Hardy-Weinberg equilibrium and of infinite size, so that extraction of migrants does not change genotype frequencies. All loci are considered neutral in the source population, so no preadaptation to any of the populations exists.

Model implementation

The impact of carrying capacity (K), growth rate (r), migration (M), selection pressure (δ) and recombination rate on F_ST’s were analyzed by exploring a range of realistic values for zooplanktonic organisms. K was varied from 2 × 10² to 2 × 10⁷ individuals, which is equivalent to densities from 0.001 to 100 individuals/L in a small pond of 200 m² and 1 m depth, in good agreement with reported average densities of cladocerans and rotifers (Carmona, Gómez & Serra, 1995; Ortells, Gómez & Serra, 2003; Tavernini, 2008). r was explored from 0.05 to 1 days⁻¹. Cladocerans show maximum r of 0.2–0.6 days⁻¹ and rotifers 0.2–1.5 days⁻¹ (Allan, 1976). The number of population founders (F) was set to 1 diapausing egg across most simulations. That is, foundation is considered a rare event. Note that as the model assumes cyclical parthenogenesis, a single diapausing egg is enough for population establishment. Studies on passively dispersed aquatic invertebrates show that the number of founders is low (Haag et al., 2005; Louette et al., 2007; Ortells, Olmo & Armengol, 2011; Badosa et al., 2017). For example, Haag et al. (2005) reported and average of 1.7–1.8 colonizers for two species of Daphnia, and 57% of studied populations were likely founded by a single individual. The effect of the numbers of founders (F) was also explored (1, 2, 5, 50 diapausing eggs). Other parameter values used in the simulations are shown in Table 1.

Simulations considered two scenarios regarding diapausing egg banks: (1) an annual, age-independent, diapausing egg survival rate on the sediment (γ = 0.763) (i.e., a diapausing egg bank); and (2) γ = 0.763 for eggs of age = 1 and a γ = 0 for older eggs (i.e., absence of diapausing egg bank). Parameters for the diapausing egg bank (γ and λ, the annual hatching rate) were estimated from rotifer diapausing egg banks (García-Roger, Carmona & Serra, 2006b) by adjusting them to the model described by García-Roger, Carmona & Serra (2006a).

The simulation model was implemented in C++ and based on Monte-Carlo procedures (code available at https://github.com/monpau/founder_effects). The Mersenne twister algorithm (Matsumoto & Nishimura, 1998) was used as random number generator. The logistic model was iterated numerically. 50 replicates for each parameter combination (but 100 for values of δ and recombination rate) were performed. For each replicate, a source population was randomly created by drawing from a uniform distribution the allelic frequencies of the n and s loci. After foundation of the two populations, 1,000 sexual generations (4,000 generations for some scenarios) were simulated.

Sampling effects were taken into account for hatching and survival of diapausing eggs if the total number of eggs in the population was lower than 1,000. Selection of migrants and gametes for mating were performed randomly regardless of the number of eggs/individuals involved.

Paralleling the procedure in an empirical study, a statistical assessment was performed. Differences between F_ST values under a neutral scenario and scenarios with selective pressure and different recombination rates were analysed with an ANOVA and a priori contrasts. Correlations between F_ST and F_STQ at different combinations of population size, recombination rates and selective pressure were also tested using Kendall’s Tau and Sperman’s Rho. All statistical analyses were performed using SPSS v. 17 (SPSS Inc., Chicago, USA).

Effect of migration

The effect of the number of migrants on genetic differentiation of neutral loci strongly depends on K (i.e., population size; Fig. 2). In both the small and the large populations, F_ST decreases with increasing migration rates, as expected under a neutral scenario (Wright, 1931). For the lowest carrying capacity tested (K = 2 × 10² individuals; Fig. 2A), F_ST decreased rapidly down to very low levels with increasing migration. In contrast, for the highest K tested (K = 2 × 10⁷ individuals; Fig. 2B), F_ST was rather insensitive to the effect of migration, and populations remained highly differentiated (F_ST > 0.2) even at high levels of migration, as predicted by persistent founder effects. Note that the high level of genetic differentiation observed is due to the strong bottlenecks during colonisation (see below and Fig. S1). The number of migrants needed to cause a noticeable decrease of genetic differentiation on neutral loci is in the order of 100 and 1,000 individuals/sexual generation for the situation without and with diapausing eggs respectively.

Effect of population size

Carrying capacity (i.e population size) had strong effects on F_ST (Fig. 3; 1,000 sexual generations). In small populations (i.e., low K) populations did not differ genetically, while in large populations, F_ST remained as high as the values observed just after population foundation, reflecting the importance of migration and persistent founder effects respectively. In the case of the low carrying capacity scenario (low K), migration seems to be stronger than the effect of genetic drift and a reduction of genetic differentiation is observed. At higher K, genetic differentiation can only operate during the earliest generations before K is reached. When K is very high, high population densities are achieved very quickly and genetic differentiation remains equal to the genetic differentiation after founding (i.e., persistent founder effect). At intermediate K high population sizes are not achieved as soon, so genetic drift can operate during a few more generations. In other words, the highest F_ST values are found at intermediate population sizes. The F_ST-K pattern is qualitatively similar with and without diapausing egg bank, but in absence of a egg bank a lower maximum F_ST at a higher K was found.

These results are robust to changes in the maximum number of sexual generations explored (results for maximum y = 100, 500, 2,000 and 4,000 generations, data not shown). However, at 100 and, to a lesser extent, 500 sexual generations, the peak of F_ST at intermediate population sizes was less pronounced than at later sexual generations. The long-term (from 1st to the 4,000th sexual generation) time course is further explored in Fig. 4. In the absence of a diapausing egg bank (Fig. 4A), F_ST decreases with time at low population size, and the situation is reversed when K increases, to finally become virtually constant (i.e., determined by the initial condition) at the largest population size explored (K = 2 × 10⁷). A qualitatively similar pattern is found when a diapausing egg bank is present (Fig. 4B), although the shift to an increasing F_ST time course, and also to F_ST constancy, occurs at lower population sizes. Note that the small negative change found at K = 2 × 10² (regardless whether a bank is assumed) is associated to the very low initial F_ST values (Figs. 4C, 4D). Also note that F_ST values are calculated after hatching of residents and migrants; for instance, at y = 1, F_ST value is not the value after foundation but after migration.

In summary, population size and presence or absence of a diapausing egg bank are key to predict the main force shaping the genetic structure. Decreasing F_ST indicates that migration is becoming the dominant factor, while increasing values show that drift becomes dominant. Our finding of practically stationary F_ST from the first generations after foundation indicates the importance of persistent founder effects on the shaping of the genetic structure of populations.

Population growth rate interacts with population size in determining the level of genetic differentiation (Fig. 5). Low growth rates result in low genetic differentiation, regardless of population size, indicating a high impact of migration. However, for population growth rates above 0.1 d⁻¹, which are common for zooplanktonic organisms, genetic differentiation becomes sensitive to variations in population size.

Effects of the number of founders

As expected, increasing the number of population founders F results in a dramatic decrease of F_ST values just after foundation (Fig. 6); for instance, if compared to F = 1, F_ST is reduced by half for F = 2, and approaches 0 for F = 50. This effect is persistent, as after 4,000 sexual generations, the level of population differentiation still shows a negative relationship with the number of founders. Given the strong effect of the number of founders, we explored in further simulations how F affects the relationships between population differentiation and other factors. Our results suggest that the patterns outlined above are qualitatively maintained for F > 1 (data not shown and Fig. S1).

Effect of local adaptation

Above, a selectively neutral scenario was assumed. The effect of local adaptation was explored at two levels of K (2 × 10⁴ and 2 × 10⁷ individuals), which are realistic values for cladocerans and rotifers respectively. Two different selection scenarios (δ = 10⁻⁴ days⁻¹, weak selection, and 10⁻² days⁻¹, strong selection) in the presence/absence of diapausing egg bank, and six recombination rates—from complete linkage to unlinked genes—were tested (Fig. 7 summarizes the results for the scenario with diapausing egg bank; see Fig. S2, for the equivalent scenario without diapausing egg bank).

With strong selection, F_STQ reaches almost maximum values—i.e., populations are almost fixed for the locally adapted alleles—regardless of K (Fig. 7). In the case of low K, all F_ST values are statistically different from those obtained without selection (p-values <0.05 except at 0.5 recombination rate; p-value = 0.057). However, F_ST values are similar irrespective of the recombination rate. In contrast, for high K, only those values of F_ST with complete linkage (recombination rate = 0) are statistically different of those found without selection. This indicates that genetic hitchhiking in large populations acts only on neutral loci tightly linked to those under selection. Otherwise, linkage to the genes under selection does not reduce the persistence of founder effects.

With weak selection, F_STQ indicates the expected result that local adaptation becomes less important than with strong selection. In large populations (K = 2 × 10⁷), F_ST values do not statistically differ from the neutral scenario, showing the larger importance of founder effects over local adaptation when selection is weak. Moreover, F_STQ values also appear to be affected by persistent founder effects. In contrast to the situation with strong selection, genetic linkage does not alter differentiation at neutral loci. However, in small populations (K = 2 × 10⁴), local adaptation does play a role. Mean F_ST values statistically differ from the neutral scenario at all recombination rates (from 0.0 to 0.5), and the variance of the distribution of F_ST values is decreased (see Fig. 3 for comparison). Note that drift is the dominant factor in relatively small (K = 2 × 10⁴) populations with diapausing egg bank.

Effects of diapausing egg banks

In the presence of strong selection the effects of diapausing egg bank (see Supplementary Fig. S2) were minimal. In weak selection conditions: (1) at high population density (K = 2 × 10⁷) genes under selection are less affected by persistent founder effects that when no bank is present (Fig. 7), and populations show a trend to be locally adapted; (2) at low population density (K = 22 × 10⁴), F_ST values at recombination rates 0.0 and 0.1 are statistically different from the neutral scenario—unlike at higher recombination rates—, which indicates that genetic hitchhiking could be of some importance; (3) at K = 2 × 10⁴ F_ST and F_STQ had higher variance at all recombination rates than in the scenario with no diapausing egg bank (Fig. 7).

In the absence of a diapausing egg bank, populations reach maximum F_STQ values in about 40–50 sexual generations regardless of population size (data not shown). However, when a diapausing egg bank exists, advantageous alleles need a longer time to reach fixation (about 150 sexual generations for K = 2 × 10⁴, and about 300 generations for K = 2 × 10⁷).

We computed F_STQ vs. F_ST correlations within each tested parameter combination. Significant correlations were found only in the case of the low K (22 × 10⁴) without diapausing egg bank. Correlation coefficient is always positive, and the ranges are: Kendall’s tau = 0.66–0.53 and Spearman’s rho = 0.73–0.56 for strong selection; Kendall’s tau = 0.68–0.32 and Spearman’s rho = 0.80–0.38 for weak selection.

Concluding remarks

Molecular screening of natural population has uncovered an unexpectedly high genetic differentiation in taxa with high dispersal potential. These findings challenged classical views of the evolutionary processes in small multicellular organisms, and when focused on aquatic invertebrates, brought to postulate a combination of processes as causal factors for that genetic differentiation, the Monopolization Hypothesis (De Meester et al., 2002). Our analysis shows that a quantitative elaboration of this multifactorial hypothesis is able to dissect the relative weights of the different factors, and their interactions. Specifically, we found that founder effects drive the genetic structure of passively dispersed aquatic organisms. We conclude that although selective factors and migration have a role in explaining genetic structure of continental aquatic invertebrates, demographic processes are dominant. By studying which factors are important in what circumstances, our analysis can help understand relevant differences among the genetic structure of different species.