The Effect of Multiple Paternity on Genetic Diversity of Small Populations during and after Colonisation

Genetic variation within and among populations is influenced by the genetic content of the founders and the migrants following establishment. This is particularly true if populations are small, migration rate low and habitats arranged in a stepping-stone fashion. Under these circumstances the level of multiple paternity is critical since multiply mated females bring more genetic variation into founder groups than single mated females. One such example is the marine snail Littorina saxatilis that during postglacial times has invaded mainland refuge areas and thereafter small islands emerging due to isostatic uplift by occasional rafting of multiply mated females. We modelled effects of varying degrees of multiple paternity on the genetic variation of island populations colonised by the founders spreading from the mainland, by quantifying the population heterozygosity during both the transient colonisation process, and after a steady state (with migration) has been reached. During colonisation, multiple mating by males increased the heterozygosity by in comparison with single paternity, while in the steady state the increase was compared with single paternity. In the steady state the increase of heterozygosity due to multiple paternity is determined by a corresponding increase in effective population size. During colonisation, by contrast, the increase in heterozygosity is larger and it cannot be explained in terms of the effective population size alone. During the steady-state phase bursts of high genetic variation spread through the system, and far from the mainland this led to short periods of high diversity separated by long periods of low diversity. The size of these fluctuations was boosted by multiple paternity. We conclude that following glacial periods of extirpation, recolonization of isolated habitats by this species has been supported by its high level of multiple paternity.


I. INTRODUCTION
When new local populations are established in a metapopulation, genetic variation within the newly founded populations is initially governed by the genetic content of founders.At a later stage, during continued input of variation through migration, the genetic composition of migrants may potentially contribute new variation and hence counteract loss by drift and selection.In brooding and sexual species, empty sites are most likely colonised by single fertilised females that bring a brood of offspring along, while founders of virgin females, males, and juveniles fail to mate in an empty site [1].In brooding species, female promiscuity (the propensity to mate multiple males) may influence the genetic variation carried by the founders, and, if so, will have consequences on the effective population size of the new population.
Females mating multiple males have broods of offspring sired by more than one male, unless sperm competition or cryptic female choice prevent this.Female promiscuity, once believed to be rare in nature, is observed in a number of animal species including mammals, amphibians, fishes and invertebrates [2][3][4].Genotyping offspring of species with promiscuous females shows that a large proportion of the females releases offspring sired by 2 − 4 males [4].In some species of fish and invertebrates, levels of multiple paternity are even higher, since the number of males siring a brood regularly reaches 6 − 10 [5][6][7][8].An extreme example of high multiple paternity is the marine snail Littorina saxatilis with 15 − 23 males siring broods of single females [9].In this study we construct a mating model for this species and analyse how multiple paternity affects population genetic variation and structure in a metapopulation.We consider established populations in equilibrium, but also populations under establishment (during initial colonisation of a previously empty habitat).
Littorina saxatilis is strictly intertidal and most abundant in rocky shores in the north Atlantic, with population densities of tens to hundreds of snails per square metre [10].In contrast to many other marine snails, L. saxatilis does not have pelagic eggs or larvae, and therefore dispersal over a few metres range is infrequent.However, snails occasionally migrate among islands, probably by rafting.It has been estimated that within an archipelago of small and large islands, 3% of the small islands receive a migrant snail each year [11].In many areas, L. saxatilis forms metapopulations with local populations inhabiting discrete localities, such as islands of an archipelago, rocky outcrops and breakwaters intermingled by sandy substrates [12,13].During the retreat of the ice sheet 12000 − 15000 years ago, L. saxatilis spread from refuge areas both in the northern Atlantic and south of the ice-cap [9].Part of this postglacial expansion comprised colonisation of hundreds of islands in the archipelago along the Swedish west coast that successively became available by isostatic uplift, a process that is still ongoing.In this system, populations on the mainland and large islands are the oldest and largest, and these are likely to act as the ultimate sources of genetic variation during colonisation of emerging islands in a stepwise manner (Fig. 1A).We have re-analysed genetic data from L. saxatilis populations in the archipelago of west Sweden and found that the first principal component shows a largely linear relationship between population genetic variation and size/age of the islands/populations, with mainland populations at the one end, skerry populations (skerry sizes ≈ 10m 2 ) at the other end, and island populations (island sizes ≈ 10 5 m 2 ) in the middle (Fig. 1B).This suggests a simple linear stepping-stone model with the mainland population acting as a source for colonisation of islands at successively younger age, and at increasing distance from the mainland (Fig. 1C).
In this paper, we investigate how multiple paternity in L. saxatilis affects spatial and temporal structure of neutral genetic diversity in a metapopulation of this species.We analyse the effective genetic population size resulting from the mating behaviour observed in earlier studies, and derive simple approximations describing the genetic diversity of populations during colonisation, and in its equilibrium state with migration.We use simulations to assess the temporal effects of migration on genetic diversity as new mutations from the mainland spread to distant islands.

II. DESCRIPTION OF THE MODEL
We construct a stepping-stone colonisation model with the following basic assumptions mimicking how L. saxatilis colonises the post-glacial archipelago of western Sweden.(1) Colonisation is frequent and rapid, as rafted fertilised females release a few hundred offspring already in the first generation [11].(2) Small skerries are likely to be colonised within a few years after emergence, and hence all newly established populations are limited by the small size of the habitat, resulting in census sizes of ≈ 10 2 − 10 3 [11].(3) Colonisation is likely to take place in a stepping-stone manner with smaller and more distantly related islands being colonised from their closest already colonised islands.For simplicity, we consider a system consisting of a mainland and of linearly arranged islands of equal carrying capacities (substantially smaller than that of the mainland).

A. Stepping-stone colonisation model
In our model, islands are linearly arranged and numbered from 1 to k, with k being furthest away from the mainland (see Fig. 1C).We include high values of k in our model (such as k = 10) in order to be able to assess saturation effects.The mainland is labelled by 0. Generations are assumed to be nonoverlapping.At the generation when the process of colonisation starts, the mainland is the only populated habitat, and the population heterozygosity on the mainland is stationary (that is, the mainland population is assumed to be old).
Within our model, an empty island becomes fully colonised in a single generation after the arrival of one or more founder females from the nearest neighbouring island.This is motivated by the very large capacity for population growth of L. saxatilis in a suitable habitat [12].In our model the founder females give rise to 2N offspring in total, with equal sex ratio, where 2N denotes the carrying capacity of an island.Upon establishing a given island population, its population size remains constant over time.In our model, mating takes place before migration, which allows us to trace only the movement of females (males also migrate, but since they will not mate after migration, they do not contribute to the progeny on the island they migrated to).Individuals are equally likely to migrate to each of their closest neighbours (but the mainland and the last island have only a single neighbour, Fig. 1C).a neighbouring island, except for empty islands that only receive migrants.
In addition to the above, we assume that the population size on the mainland is much larger than the population size of a colonised island (unlike other models which assume that all habitats have the same carrying capacity, see, for example, [14]).This allows us to treat the mainland as the only source of genetic variation to the island populations.In our computer simulations (see Sections II-IV in the electronic supplementary material, ESM) we set the mainland heterozygosity to unity.This simplifies the simulations, since the dynamics of the mainland does not need to be simulated explicitly.

B. Mating model
In order to study the consequences of multiple paternity for genetic diversity, we introduce a mating model to describe different levels of multiple paternity in mating systems.
Based on known life-history traits [12], we assume that the duration of the reproductive cycle of females is the same for all females.Each female carries beneath her shell juveniles of varying degrees of maturity, and juveniles are released at an approximately steady rate.We also assume that after a successful mating, the mated female obtains a sperm package able to fertilise female eggs during its persistence time.Our observations show that sperm can be stored and used up to a year after mating.The number of eggs fertilised by a single sperm package is denoted by N eggs , and we assume that this number is the same for all sperm packages that a female receives during the reproductive cycle.The total number of sperm packages received by a female during her reproductive cycle is denoted by l.
The probability p that two eggs are fertilised by the same sperm package is p = l −1 assuming that all sperm packages a female received during her reproductive cycle persist until the end of the reproductive cycle, that all eggs are fertilised after all sperm packages have been collected, and that sperm packages are chosen with replacement to fertilise eggs.
The scheme presented above models the process of mating at an individual level.We assume that individuals belong to a well mixed diploid population of N m males and N f females, and we take N m ≫ 1 and N f ≫ 1.In our model, a female encounters s ≤ N m different males during her reproductive cycle.Having s < N m reflects the limited movement of snails during the reproductive cycle.To simplify the analysis, we assume that all males a female encounters are equally likely to be her mating partner in any of the matings she experiences.Moreover, we assume that all females are equally successful mothers.Within the model described, the effective size of a single local population is given by (see Section I of the ESM): Here, κ is the probability that two offspring having the same mother share a father Since κ decreases as the number of available mates s of a female increases (keeping the value of p < 1 fixed), we take κ −1 as a measure of the degree of multiple paternity.
Our model reduces to the model in [15] in the case a female encounters all males present in the population (upon substituting the number of matings in [15] by p −1 ).If each female mates with all males in the population and the probability that two eggs are fertilised by the same sperm package is p = 0, our model reduces to random mating.
We show in Fig. 2 that the effective population size increases as s increases.For the parameters set in Fig. 2, the maximum value of N e is equal to N m + N f , which corresponds to the effective population size under random mating.The increasing trend of N e saturates at s ≈ 10 for a given value of p.In summary, by increasing the degree of multiple paternity, the effective population size becomes larger (as found in [16] for a different mating model).
We compared the male family sizes (the number of offspring of a father involved in siring the brood of a given female) obtained under our model to those obtained under experimental conditions, as well as in natural populations.We conducted experiments such that 6 virgin females were placed in separate aquaria, and each was accompanied with s = 10 males.The sire of each offspring produced during a year was determined by genotyping 8 microsatellite loci.

III. RESULTS
In Fig. 3A we show the histograms of male family sizes obtained experimentally and in Fig. 3B we show similar data collected from females in natural habitats [9].For both sets of data we use computer simulations in order to find the parameters in our model resulting in male family sizes that are closest to those empirically observed (the brood sizes analysed in computer simulations correspond to those from the empirical data).For data obtained under experimental conditions, we vary the probability p in the computer simulations, and for data collected in natural habitats, we vary both s and p.We use a χ 2 -test to quantify the distance between the empirical and simulated data.The best fits obtained are shown in Fig. 3A-B by circles, and they correspond to p = 1/15 (Fig. 3A), and to s = 20, p = 0 (Fig. 3B).In summary, Fig. 3 suggests that our mating model describes empirical data obtained in the experimental setup well, whereas the agreement between the model and the empirical data taken from natural populations is less good.This is discussed in Section IV.
To address the question of how multiple paternity affects genetic variation and structure in our metapopulation, we analyse genetic diversity under the colonisation model described in Methods.We analyse separately two phases of population dynamics on each island: initial colonisation, and the equilibrium state that develops once the colonisation phase is over.For a given island, we compute the expected heterozygosity in the generation when the island is colonised (colonisation phase), as well as the heterozygosity in the equilibrium state.The corresponding analytical computations are described in detail in Sections II and III of the ESM.We also study the temporal changes of heterozygosity under our model by computer simulations.In the following two subsections we present separately the results for the colonisation phase and for the equilibrium state.

A. Colonisation phase
We compute the heterozygosity during colonisation analytically using a coalescent approach [17].We represent the taken from [9]; the data correspond to four broods of sizes 79, 77, 71, and 53.The width of the bins are chosen so that the expected number of counts in each bin is not smaller than 5.The probability assigned to each bar is proportional to the bar area.Symbols and error bars show the result of the best fit to the experimental data, together with their 95% confidence intervals: p = 1/15 in A, and p = 0, s = 20 in B. We simulated 10 3 independent realisations of the mating process.
population-size history of the population on island i as a sequence of i bottlenecks (i is the number of colonisation events that the population ancestral to that on island i went through before the island was colonised).Our analytical result is valid for small migration rate, M ≪ 1.Under this assumption, colonisation of empty islands occurs rarely, but when it does, an island typically receives a single founder female (see Section II of the ESM).The result is given in Eq. ( 14) of the ESM and in Fig. 4A.We see that the colonisation-phase heterozygosity decays as distance from the mainland increases.We note that the results of our computer simulations (see Fig. 1A-C in the ESM) agree well with the analytical results for low migration rates.For large migration rates (M = 0.5.i .e .0.5 females on average per generation) by contrast, the simulations assume somewhat higher values than the theory.The reason for this deviation is that for M = 0.5 it is probable that more than one founder female comes to an empty island to establish the population, and, consequently, will contribute with more genetic variation than just one founder female.For natural populations of L. saxatilis it has been estimated that 3% of empty islands receive a migrant each year [11].This estimate is close to the lower value of M used in our simulations (M = 0.05).An important result shown in Fig. 4A is that at any distance from the mainland, multiple paternity results in higher heterozygosity than single paternity.Mating two males (s = 2) increases the values of single-paternity heterozygosity by 10 − 100% for the parameters used in Fig. 4A, and mating ten males (s = 10) increases the values of single-paternity heterozygosity by 10 − 300% (Fig. 4A).The largest increase is observed at the island furthest from the mainland.We also note that mating more than about 10 males only marginally increases the heterozygosity (results not shown), as found in the case of freely mixing populations (see Subsection II B). les.

B. Equilibrium state
We show in Section III of the ESM how to compute the heterozygosity in the equilibrium state at distance i from the mainland using recursion relations (see, for example, [18]).Note that this derivation does not require M to be small.
The results are given in Section III of the ESM and in Fig. 4B.As in the colonisation phase, the equilibrium-state heterozygosity decreases as distance from the mainland increases.Also, by increasing the degree of multiple paternity, the heterozygosity at a given island increases (this effect saturates at s ≈ 10, results not shown).In contrast to the strong effect of multiple paternity during colonisation, the effect is substantially smaller in the equilibrium state.We find that the single-paternity heterozygosity in the equilibrium state increases by 10 − 20% for s = 2, and by 20 − 50% for s = 10 (Fig. 4B).As in the colonisation phase, the largest increase is observed at the island furthest from the mainland.
In addition, we examined the variation in heterozygosity over consecutive generations in a particular realisation of our model.We find that the heterozygosity shows strong temporal fluctuations (Fig. 5A).Notably, the fluctuations are strongest furthest from the mainland, with periods of high diversity separated by long periods of near or complete fixation.Hence the distribution of heterozygosity at large distance from the mainland is bimodal.The heterozygosity is expected to have a bimodal distribution in the case of a very small rate of income of new genetic material, as pointed out in [19] (see Fig. 1 in [19]).
In what follows, we analyse how the durations of the phases of low and high heterozygosity are affected by multiple paternity.Using the results of computer simulations, we compute the average durations of low-and high-heterozygosity phases at the island furthest from the mainland.We also derive corresponding analytical results under the assumption that the scaled migration rate M is small (see Section IV of the ESM).For island i = 10 we show in Fig. 5B the durations of lowand high-heterozygosity phases relative to their corresponding values for a single mate (s = 1).Fig. 5B shows that multiple paternity prolongs the duration of the high-heterozygosity phase, and decreases the duration of the low-heterozygosity phase.For example, the high-heterozygosity phase for the highest level of multiple paternity shown (s = 10) is prolonged by around 40% compared to its value under single paternity (s = 1).The low-heterozygosity phase is shortened by only around 10% for s = 10 (Fig. 5B).For comparison, Fig. 5C shows the equilibrium-state heterozygosity relative to its corresponding value for a single mate (s = 1).In conclusion, multiple paternity promotes heterozygosity by prolonging the duration of peaks of variation that reach the most distant islands.

IV. DISCUSSION AND CONCLUSIONS
In this study we analysed the effect of multiple paternity on genetic diversity over spatial and temporal scales in a metapopulation.We quantified the effect of multiple paternity during the colonisation of semi-isolated populations and in the equilibrium state developed after the colonisation phase.Our conclusions given below can be generalised to a metapopulation of patches that are partly isolated from each other, for example, by sandy beaches, harbors or other types of less suitable habitats.
We introduced a mating model which allows for different levels of multiple paternity in a population.In order to determine how realistic our mating model is, we compared the male family sizes of female broods obtained within our model, to empirically observed family sizes in populations of L. saxatilis from natural habitats and under experimental conditions.We found that male family size distributions from the experiments, where the true number of mates was known, are in very good agreement with our mating model.The best-fitting parameters indicate fewer matings than the brood size, suggesting that some of the eggs are fertilised by sperm packages retained between matings.By contrast, the corresponding empirical distribution in natural populations is best fitted by assuming an unlimited number of matings (i.e.no sperm retention).This is consistent with the high density of snails observed in the wild.However, the empirical distribution shows an excess of males with a single progeny compared to the mating model with the best-fitting parameters.This discrepancy could be explained by post-zygotic selection, where sperm from several males compete, resulting in uneven success among males.However, this effect should also be present in the experiments, but we did not find any evidence of it in the data.Another possibility is that this pattern is due to variation in individual snails movements in the wild, where some snails might move around more extensively and mate with a new partner each time while others stay within a small area and mostly re-mate with individuals in the close surroundings.This variation in mating behavior cannot happen in the experiments where the snails are confined to each other.
Within our mating model, by increasing the degree of multiple paternity, the effective population size increases, and thus the heterozygosity increases.However, since mating is considered to be costly [9], it is not clear whether or not mating multiple males is an evolutionary strategy to increase the heterozygosity.Recall that we estimated that in natural populations of L. saxatilis the probability that two eggs are fertilised by the same sperm package is likely to be very small.Under our model, this probability is equal to zero if each sperm package fertilises one egg, or if the actual number of sperm packages a female receives during her reproductive cycle is very large.If the latter applies, we find that it is unlikely that the heterozygosity increase is the main reason for the extreme multiple paternity in this species.As earlier suggested, it seems likely that the cost of rejecting an intercourse is higher than the cost of accepting it, and a consequence of convenience polyandry [9].Nevertheless, the consequences of multiple paternity for heterozygosity in relatively small and semi-isolated populations are substantial.This is summarised in the following.
At a given distance from the mainland, populations with high degree of multiple paternity establish higher heterozygosity than populations with low degree of multiple paternity.While this effect is substantial during colonisation, it is modest in the equilibrium state.This is explained in the following.Upon the arrival of founder females to an empty island, the carrying capacity of an island is reached within a single generation.Therefore, such a newly established population receives genetic material of most males that the founder females were inseminated by.By contrast, in the equilibrium state, all mothers present in an island contribute to the population in the next generation, and hence the impact of immigrant females to the next generation is rather small.From this reasoning, we find that it is possible that the effect of multiple paternity upon heterozygosity during colonisation might decrease if the growth rate of the island populations up to the carrying capacity were less than infinite (as assumed in our model).
The heterozygosity at distances far from the mainland fluctuates significantly.Long periods of almost complete loss of genetic variation are interrupted by bursts of high heterozygosity, and this effect is boosted by multiple paternity.The durations of high-and low-heterozygosity phases could be an important survival factor in natural populations.For example, the low-heterozygosity phase could be disadvantageous if a malignant disease appears in the population, assuming that only a particular mutation (not present in the population, or being rare) can survive the disease.
The wave-like nature of the spread of new alleles from the mainland population is also seen in the correlation of genetic diversity at neighbouring islands.We find that the correlation between heterozygosities at a pair of nearest-neighbouring islands increases as distance from the mainland increases (results not shown).These results suggest intermittent bursts of genetic diversity in remote islands, an effect which becomes stronger as the degree of multiple paternity increases.
The conclusions given above are confirmed by our computer simulations.In order to minimise computing time during simulations, we assumed that the mainland heterozygosity is equal to unity, which guarantees that whenever a migrant from the mainland comes to the first island, it carries genetic material that previously existed neither on the mainland nor on any of islands (and thus the population dynamics on the mainland does not need to be simulated explicitly).However, we emphasise that the conclusions given above qualitatively do not depend on the value of heterozygosity on the mainland.
We note that, unlike in our model, it is possible that the rate of successful colonisation in natural habitats is smaller than the rate of migration.For example, if an immigrant female carries a small number of progeny, her progeny alone may not be enough to colonise an empty island successfully.By allowing for the rate of successful colonisation to be smaller than the rate of migration in the colonisation model, the equilibrium-state values of heterozygosity remain equal to those obtained under the assumption that the colonisation and migration rates are the same (as in our model).The values of heterozygosity during colonisation, by contrast, are expected to differ from those found here.
In summary, this study can be used to quantify the gene flow between partly isolated natural populations using allelic frequencies at a number of neutral loci.Since our results show that the heterozygosity exhibits extreme fluctuations in populations founded through repeated founder events, we raise the question of whether similar fluctuations can be observed at any given time at neutral loci sampled genome-wide.In order to answer this and related questions, the effect of recombination needs to be analysed.Since island populations in our model experience at least one severe bottleneck, we expect that the degree of linkage disequilibrium in the colonisation phase is constant over a range of genetic distances, as shown in [20].However, how multiple paternity affects linkage disequilibrium during colonisation and in equilibrium is yet to be understood.It would also be interesting to analyse how selection combined with recombination affects genetic diversity in a metapopulation.Such results would provide an advance in the endeavor of identifying genes under selection, and especially, the genes underlying speciation [21][22][23].
Here ǫ τ is the inbreeding coefficient which stands for the probability that two alleles within a single randomly chosen individual are identical at time τ .The second term, χ τ , is the coancestry, that is the probability that two alleles sampled in generation τ from two different randomly chosen individuals are identical.In what follows, we compute F τ under our mating model, and thereafter we use Eq. ( 1) to determine the corresponding effective population size.We assume that mutation rate per generation per allele is µ ≪ 1, and we employ the infinite-alleles model [1].Our calculations are based on the approach used in [2].Under the assumption that all females are equally successful in producing offspring, it follows that the probability for two offspring to come from the same female, P f , is equal to P f = N −1 f .Let κ be the probability that two offspring coming from the same mother also share a father.The probability that two offspring come from the same male, P m , is thus: where N −1 m is the probability that two offspring having different mothers stem from the same father.The probability κ has two contributions: the probability that two offspring come from the same sperm package, p, and the probability that two offspring do not come from the same sperm package but they come from the packages coming from the same male, (1 − p)/s.Here, we assume that sperm packages are chosen with replacement to fertilise eggs.It follows that κ is given by where p is As mentioned in Section II of the main text, we take κ −1 as a measure of the degree of multiple paternity.The case of single paternity corresponds to κ = 1.
Using the expressions for P m and P f , we compute ǫ τ , and χ τ recursively.Under our model we find In order to compute the heterozygosity, we note that in our model, the mainland acts as the only source of genetic variation.This allows us to argue the following.First, if the MRCA of two alleles sampled randomly from the newly established population at island i was born on island j < i (j = 0), the two alleles sampled are identical.Second, if the MRCA was born on the mainland, the two alleles are expected to be identical with probability c .Therefore, in order to compute c , it suffices to determine the probability that the MRCA of two lines sampled from the newly established population at island i stem from an allele that was born on the mainland, P (0|i).
The probability P (0|i) has two contributions.The first contribution is the probability that the MRCA of two alleles sampled at island i is not found during a bottleneck.We find this to be equal to 1 − 1 8 (1 + κ).The second contribution is the probability that the MRCA of two alleles is not found between two successive bottlenecks.This term is equal to 2M N e (2M N e + 1) −1 .It follows that P (0|i) is given by Therefore, the colonisation-phase heterozygosity at island i is: Note that for the case described here, the population size switches between 2N (N males and N females during the waiting time before the colonisation of the next island) and unity (one inseminated mother during a bottleneck).Therefore, the probability κ does not enter Eq. ( 14) only through the expression for N e , Eq. ( 11).The heterozygosity in the colonisation phase for different parameters of our model is shown in Fig. 1A-C.The analytical result, Eq. ( 14), is shown as solid lines, and the results of our computer simulations are shown as symbols (the solid lines in Fig. 1B correspond to the solid lines in Fig. 4A in the main text).We see that the agreement between Eq. ( 14) and the results of computer simulations is good for M = 0.05, whereas for M = 0.5, Eq. ( 14) underestimates the results of computer simulations.This is discussed in Section III in the main text.

III. HETEROZYGOSITY IN THE EQUILIBRIUM STATE
In this section, the equilibrium state within the model introduced in Section II in the main text is analysed.Expressions for the equilibrium-state heterozygosity at distance i = 1, . . ., k from the mainland are derived under the assumption that all islands are populated (that is, the colonisation phase is over).As discussed in Section I of this supplementary material, the inbreeding coefficient in generation τ at island i, ǫ (i) τ , and the coancestry χ (i) τ contribute to the homozygosity F (i) τ at island i in generation τ .The coancestry between islands i, and j is equal to the inter-island homozygosity, and for this case we use the notation χ As mentioned in Section II of this supplementary material, we assume that the mainland is the only source of genetic variation.All habitats are assumed to have equal numbers of males and females.The population size on the islands i = 1, . . ., k is assumed to be large (and equal to 2N ), but much smaller than that of the mainland.As before, the heterozygosity on the mainland in generation τ = 0 is denoted by According to the spatial model introduced in Section II of the main text, the female-migration rate per island per generation is 2M for islands i = 1, . . ., k − 1, whereas for the mainland and for the island furthest from the mainland it is equal to M .Since the population size on the mainland is much larger than that of a populated island, the process of migration does not affect genetic variation on the mainland.Therefore, we have c .For the island populations, we consider separately the inter-island and the intra-island homozygosity.First, we treat sampling from two distinct islands, i = j.Second, we consider the case of sampling within a single island i = 1, . . ., k.Our calculations given below are based on the approach employed in [3].
The inter-island homozygosity F (i,j) τ +1 for i = 0, 0 < j ≤ k satisfies the following recursion: Here m = 2M/N ≪ 1 is the migration rate per island per female per generation, and δ j,k is equal to unity when j = k, and it is zero otherwise.The inter-island homozygosity for 0 < i < k, 0 < j < k, i = j, obeys: Lastly, when i = k, 0 < j < k, we find Finally, for the homozygosity at distance 0 < i < k from the mainland we find: For the island furthest from the mainland, i = k, the corresponding expression is: By setting in Eqs. ( 22)-( 27) the terms involving the time derivative to zero, one finds the expressions for the equilibriumstate homozygosity of the system.The equilibrium-state heterozygosity is obtained upon subtracting the equilibrium-state homozygosity from unity.Upon setting H (0) = 1, the results shown in Fig. 1D-F are obtained.We note that the lines in Fig. 1E correspond to the lines shown in Fig. 4B in the main text.

IV. DURATIONS OF LOW-AND HIGH-HETEROZYGOSITY PHASES
A scheme for computing the durations of low-and high-heterozygosity phases is shown in Fig. 2A.As indicated in this figure, we consider values of the heterozygosity smaller than 0.1 to be low, and values of the heterozygosity larger than 0.4 to be high.The maximum value for the low phase (0.1) is chosen because a locus is commonly considered monomorphic if the heterozygosity at this locus is ≤ 0.1 (see [4]).The minimum value for the high phase (0.4) is chosen since the typical maximum value that the heterozygosity has at the island furthest from the mainland is ≈ 0.5 for the parameters chosen in Fig. 2A.Note that the maximum value of the heterozygosity at a locus with only 2 allelic types is equal to 0.5.
The method used to record the times of transitions between low-and high-heterozygosity phases in computer simulations is explained next.Say the heterozygosity is in the high phase at time τ = 0. We record first the nearest point in time when the heterozygosity becomes less than 0.1.Say this occurs in generation τ 0 .Second, we search for the last generation before τ 0 in which the heterozygosity has a value larger than or equal to 0.25 (the middle point between 0.1 and 0.4).Say this happens in generation τ 1 < τ 0 .We take τ 1 to be the time of transition from the high-to the low-heterozygosity phase.The transitions from the low to the high phase are recorded using a similar scheme.The durations of the high-and low-heterozygosity phase, T high and T low , relative to their corresponding values for s = 1, computed as explained above, are shown as symbols in Fig. 2B,  D, and F. For comparison, we show in Fig. 2C, E, and G the corresponding equilibrium-state heterozygosities relative to their values for s = 1.
Next, we briefly explain how to estimate T low , and T high analytically.The heterozygosity can switch from the low to the high phase if new genetic material comes to the population, and if this material is not lost due to random genetic drift.Recall that, in our model, migration is the only process bringing new genetic material to islands.However, some migrations bring genetic material that already exists in a given island population (but note that when H (0) = 1, then the first island receives new genetic material with each migration).Yet, it is possible to estimate the effective successful migration rate per allele per generation, m (i) e , at a given distance i from the mainland using the analytical results derived for the equilibrium-state heterozygosity at island i, namely H (i) .Here, the term 'successful' implies that migrants bring new genetic material to the population.Note that under the assumption that a new allelic type is introduced to the population at distance i from the mainland at rate m (i) e per allele per generation, the equilibrium-state heterozygosity is computed as e N e 1 + 4m where 2m e N e is the total number of new allelic types introduced in the population per generation (this expression is analogous to the usual expression for the heterozygosity which involves the scaled mutation rate, see Eq. ( 9)).In the case m (i) e N e ≪ 1, it follows that typically every (2m e N e ) −1 generations, one mother comes to an island carrying a new allele.Furthermore, it is known [5] that a large haploid population with effective size 2N e spends on average 2/N e generations in the state with this allele having the frequency N e before fixation occurs.It follows that the typical waiting time for such a successful establishment of new genetic material at island i is equal to (4m (i) e ) −1 .This estimate agrees well with the results of our computer simulations (results not shown).
We now explain how to estimate the average duration of the high-heterozygosity phase, T high .In the case of rare income of new allelic types to island i, m (i) e N e ≪ 1, the island population at a given locus typically has at most 2 alleles.Thus, T high can be approximated by the time that a locus with two alleles in a Wright-Fisher population of N e diploid individuals needs to reach fixation.Note that under the condition m (i) e N e ≪ 1, fixation occurs typically much before new genetic material arrives to the population.The number of generations until fixation in a diploid population at a locus with two alleles, τ loss , is [4] τ loss (α) ≈ −4N e [α ln(α) + (1 − α) ln(1 − α)] . ( Here α is the initial frequency of a given allelic type.For a given value of the number of available mates, s, we compute T high upon integrating Eq. ( 29) from α = 0.1 to α = 0.9.The integral boundaries correspond to the value of heterozygosity ≈ 0.2, which is close to the value of 0.25 in the method used for determining the transition time between the low and the high phase.
We note that Eq. ( 29) results in underestimated time of fixation if the rate of income of new genetic material is not too small.This is confirmed by the results shown in Fig. 2B, D, and F. Note that the solid lines in Fig. 2B, and C correspond to the lines in Fig. 5B, and C (respectively) in the main text.

2 :FIG. 2 :
FIG.2: Effective population size (for the given number of available mates, s, and the probability that two eggs are fertilised by the same sperm package, p) relative to Nm + N f .The grey region depicts the average number of fathers in four broods of L. saxatilis[9].Parameters: Nm = N f = 10 3 .

FIG. 3 :
FIG.3: Histograms of male family sizes within broods of females.FIG.3: Histograms of male family sizes within broods of females.Bars in panel A show the empirical data obtained under experimental conditions for s = 10; the data correspond to six broods, two of size 20, three of size 19, and one of size 16.Bars in panel B show results taken from[9]; the data correspond to four broods of sizes 79, 77, 71, and 53.The width of the bins are chosen so that the expected number of counts in each bin is not smaller than 5.The probability assigned to each bar is proportional to the bar area.Symbols and error bars show the result of the best fit to the experimental data, together with their 95% confidence intervals: p = 1/15 in A, and p = 0, s = 20 in B. We simulated 10 3 independent realisations of the mating process.

FIG. 4 :
FIG. 4: Analytically computed heterozygosity during colonisation FIG.4: Analytically computed heterozygosity during colonisation (A), and in the equilibrium state (B).The lines shown from top to bottom correspond to: s = 10, s = 5, s = 3, s = 2, and s = 1.Remaining parameters: the mainland heterozygosity is set to unity, scaled female migration rate M = 0.05, number of females in each populated island N = 100, probability that two eggs are fertilised by the same sperm package p = 0.1, number of islands k = 10.
FIG. 5: (A) Heterozygosity as a function of the distance from the mainland and of time (single realisation of the model described).Mainland is not shown.The data correspond to 10 5 generations after the initial 7 • 10 6 generations.The number of available mates is s = 10.(B) Durations of low-and high-heterozygosity phases (blue, and red) relative to their corresponding values for s = 1.C Equilibrium-state heterozygosity (black) relative to its corresponding value for s = 1.Remaining parameters are as in Fig. 4B.

FIG. 1 :
FIG.1: Heterozygosity in the colonisation phase (A-C), and in the equilibrium state (D-F).The results corresponding to Eq. (14) are shown as solid lines, and the results of computer simulations are shown as symbols.The number of available mates is s = 1 (blue), s = 2 (red), s = 3 (green), s = 5 (magenta), and s = 10 (black).Averages are over τ = 10 4 independent realisations of the process of colonisation of empty islands in A, that is over τ = 2 • 10 4 realizations in B, and C. In panel D, averaging is done over τ = 1.5 • 10 7 generations (the initial τ = 10 7 generations being discarded), and over three independent realisations of our model.In panel E, averages are over τ = 4 • 10 7 generations (the initial τ = 5.5 • 10 7 generations being discarded) and over four independent realisations of our model.In panel F, averaging is done over τ = 5 • 10 7 generations (the initial τ = 5 • 10 6 generations being discarded) and over five independent realisations of our model.Remaining parameters used: mainland heterozygosity H (0) = 1, number of islands k = 10.
FIG. 2: (A) Illustration of the method used to determine the duration of low-and high-heterozygosity phases.The heterozygosity represented in terms of the low and high phases is shown by the magenta line.The black line depicts the result of computer simulation.The data shown correspond to those in Fig. 5A in the main text, for island 10. (B), (D), and (F) Durations of low-and high heterozygosity phases relative to their corresponding values for s = 1 (blue and red, respectively).(C), (E), and (G) Equilibrium-state heterozygosity at island 10 relative to its corresponding value for s = 1.The results of computer simulations are shown as symbols, and the analytical results are shown as solid lines.The parameters in B, and C correspond to those in Fig. 1D.The parameters in D, and E correspond to those in Fig. 1E.The parameters in F, and G correspond to those in Fig. 1F.
On average, M females migrate per generation from one island to [12]1:(A) Physical structure of the marine habitats of Littorina sax-FIG.1:(A)Physicalstructure of the marine habitats of Littorina saxatilis: mainland (red), islands (green), and skerries (blue).(B)Principalcomponents of allozyme population differentiation in L. saxatilis (data from[12], the presumably selected locus Aat − 1 is excluded).Populations are classified as mainlands (red), islands (green), or skerries (blue).(C) Stepping-stone model of a section