Using mathematical modelling to investigate the adaptive divergence of whitefish in Fennoscandia.

Modern speciation theory has greatly benefited from a variety of simple mathematical models focusing on the conditions and patterns of speciation and diversification in the presence of gene flow. Unfortunately the application of general theoretical concepts and tools to specific ecological systems remains a challenge. Here we apply modeling tools to better understand adaptive divergence of whitefish during the postglacial period in lakes of northern Fennoscandia. These lakes harbor up to three different morphs associated with the three major lake habitats: littoral, pelagic, and profundal. Using large-scale individual-based simulations, we aim to identify factors required for in situ emergence of the pelagic and profundal morphs in lakes initially colonized by the littoral morph. The importance of some of the factors we identify and study - sufficiently large levels of initial genetic variation, size- and habitat-specific mating, sufficiently large carrying capacity of the new niche - is already well recognized. In addition, our model also points to two other factors that have been largely disregarded in theoretical studies: fitness-dependent dispersal and strong predation in the ancestral niche coupled with the lack of it in the new niche(s). We use our theoretical results to speculate about the process of diversification of whitefish in Fennoscandia and to identify potentially profitable directions for future empirical research.

Introduction

The diversity of species on Earth continues to provide inspiration for scientists studying speciation and the origins and maintenance of biodiversity. What makes these processes extremely complex and difficult to understand is that different evolutionary and ecological factors controlling their dynamics act simultaneously and often have opposing effects. The complexity of speciation processes implies that mathematical modeling can potentially play a very important role in uncovering its general rules and patterns. By now we have an impressive array of models and modeling techniques that shed light on the conditions, probability, waiting time to and duration of speciation, the degree of genetic and phenotypic divergence between the emerging species, and the way different resources (including space) are partitioned between the sister species. Models also explain the effects of different parameters and factors (such as the strength of selection, rates of mutation, recombination, migration, population size, number of loci, distribution of allelic effects, etc.) on the dynamics of speciation[1-9].

Most of this work has focused on models of speciation aims for both generality and mathematical simplicity. These models are very useful and insightful in uncovering general rules and patterns of speciation, adaptive radiation, and biological diversification. However, their generality almost necessarily implies that these models are very difficult to apply to any particular systems and species studied by empirical biologists. Therefore it is very important to supplement general models of speciation with those tailored for specific biological cases.

Studying models tailored for particular case studies can be very useful from several perspectives[9]. First, mathematical models emerging from these projects do lead to a better understanding of the evolutionary dynamics of the studied specific systems. Second, although the relevant models are case-specific, they contribute towards building the general theory of speciation, e.g. by supporting or undermining the generality of particular observations and patterns. Third, the process of building a mathematical model even for a particularly well-studied empirical system usually reveals the lack of biological understanding or crucial empirical data needed to make appropriate modeling assumptions or specify parameters. This can greatly stimulate further empirical work to remove these limitations.

By now a relatively small number of such models have been developed for some of the best studied systems. These include models aiming to capture the dynamics of non-allopatric speciation of cichlids in a lake[10-13] and palms on an oceanic island[14], hybrid speciation in butterflies in Central America[15], ecomorph formation in marine snails in Sweden[16], pulmonate snails[17], and parallel adaptation in threespine stickleback[18].

Here, we continue this work by focusing on a young and well-documented empirical system – the lacustrine European whitefish (Coregonus lavaretus (L.)) in Fennoscandian postglacial lakes. Ten to twenty thousand years ago Fennoscandia was covered by a thick ice sheet[19]. The ice sheet retracted 8000–10000 years ago[20,21], creating a landscape dominated by inter-connected lakes and rivers and thus providing the opportunity for postglacial immigration of cold-water adapted freshwater fish, such as whitefish. After the glacier melting and land uplifting, this region has a high density of lakes, with deep lakes having three major habitats (littoral, pelagic and profundal) for fishes. Almost all lakes in this region contain whitefish, but most lakes have only one littoral morph, some lakes two morphs (littoral and pelagic), and a few large and deep Fennoscandian lakes three morphs (littoral, pelagic and profundal)[22-24]. These different morphs of whitefish show habitat specific patterns in their resource use, body size, gill raker number, and life-history traits[22,23]. The littoral morph has the largest body size and intermediate number of gill rakers, the pelagic morph the smallest body size and highest gill raker number, and the profundal morph a small body size and the lowest gill raker number[25-27]. It should be noted that the denotations of the three morphotypes as the littoral, pelagic and profundal morphs do not imply that they solely reside in the respective habitats, but rather that their prime adaptation is related to the utilization of the specified habitat and its associated trophic niche. Genetic evidence suggests a rapid divergence of these morphs from an ancestral monophyletic lineage, driven by natural selection on gill rakers and body size[22,26,27]. Besides historical contingency, ecological opportunity mediated by interactions, such as resource competition and predation, may play an important role in diversification of whitefishes with regard to ecomorphology and life-history[24,28,29].

What is remarkable about this system, compared to other postglacial lakes with pelagic and littoral morphs, is the presence of up to three distinct morphs in some lakes. This raises the question: could they have diverged within the lake, or have different morphs evolved in different lakes and came into contact later? Our goal here is to use individual-based simulations to better understand conditions for and various factors (abiotic and biotic) controlling the divergence of Fennoscandian whitefish into three principal habitats of subarctic lakes. Our particular focus will be on the effects of selection for local adaptation, gene flow, carrying capacity of the habitat, and predation.

Methods

Model

Environment

We consider finite sexual diploid populations inhabiting isolated lakes. Each lake has three different ecological niches (habitats) with its own whitefish population: littoral, pelagic, and profundal. Each lake can have up to four predator species feeding on whitefish: pike (Esox lucius), perch (Perca fluviatilis), burbot (Lota lota), and brown trout (Salmo trutta) (Table 1)[28,30,31]. Time is discrete, the unit of time is one year, and the generations are overlapping.

Table 1

Habitats and the maximum prey size for different whitefish predators[28,30,31].

Predator	Max prey size	Habitats
Pike	25 cm	Littoral
Perch	15 cm	Littoral
Burbot	20 cm	Littoral, Profundal
Brown trout	15 cm	Littoral, Pelagic

Individual whitefish

Each individual is characterized by three phenotypic traits: size , the number of gill rakers , and whether or not it is sexually mature. We discretize size into four stages: 0, 1, 2 and 3, roughly corresponding to the size of 2 cm, 15 cm, 20 cm, and 30 cm in whitefish. Size stage 0 describes newly born individuals.

We assume that immature fish can grow while mature fish stop growing and invest all their available energy into reproduction. If an immature fish survives to the end of the year, it either matures, with a probability depending on its current size , or grows to the next size class , with probability . Newly born individuals can only grow but not mature (i.e., ), whereas individuals reaching the largest size always mature (i.e., ). Note that in our model maturation always happens at the end of the year reflecting the need to accumulate energy for gonad development e.g.[32]. Figure 1 illustrates these assumptions.

Figure 1

Life-history of a whitefish. The fish grows in size until it matures. The probabilities of maturation at stage 1 and 2 are and , respectively. If a whitefish reaches stage 3, it definitely matures the next year and stops growing.

The studies on Coregonus clupeaformis[33,34] and salmonids[35] show that maturation rate and the number of gill rakers are polygenic traits. Correspondingly, we assume that maturation rates and are additively controlled by diallelic loci each. (In numerical simulations ). The corresponding allelic effects are scaled so that and stay between zero and one. The number of gill rakers is also additively controlled by loci but each locus has multiple discrete alleles. These alleles are subject to step-wise mutation and the effects are scaled so that trait can take only integer values and mutation changes only by plus and minus one. We assume that all genes are physically unlinked. Mutations happens with probability per gene per reproduction.

Selection

We assume that the whitefish population is subject to density-dependent viability selection due to intraspecific competition, habitat-specific stabilizing selection on gill raker number, size-dependent mortality due to predation as well as fertility selection due to maturation rates and body size differences.

Condition

The gill raker number controls the food available to fish in a given environment. We define the condition of a whitefish with gill rakers in niche asHere is the optimum gill rakers number in niche and is a parameter measuring the strength of stabilizing selection on towards this optimum. For whitefish, we set at , and for the littoral, pelagic, and profundal habitats, respectively. These values are representative for well-established morphs and are close to those observed in whitefish morphs LSR, DR, and SSR discussed above[24,27].

Density-dependent selection

We assume that the population is subject to density-dependent mortality with the survival rate of individuals with gill rakers and of size in niche being

Here is the total number of individuals of size in niche , and constant parameter is the population size at which a population of perfectly adapted individuals (i.e., individuals with perfect condition ) would have a survival rate of (see the Appendix). We will interpret parameters as measures of carrying capacity of the corresponding sizes in the corresponding niches. The larger the number of competitors and the worse is the individual condition , the smaller is the survival rate . Our formulation follows the Beverton-Holt model[36] and implies that density-dependent competition occurs only between individuals of the same size in the same ecological niche. The latter assumption reflects that fact that resource partitioning among whitefish morphs is very strong e.g.[23,25] and that intra-morph resource competition is much higher and the most pronounced between individuals of nearest age cohorts and size classes using same resources at the same habitats[37,38].

Predation

Maximum gape of a predator restricts the sizes and shapes of prey that can be eaten. We assume that each predator species in our system is characterized by a maximum prey size . To describe predation, we posit that predator removes a random proportion of surviving whitefish with sizes smaller than or equal to its gape size . Values are parameters in our model. When whitefish are bigger than the maximum prey size of a predator, they are not predated upon. The maximum prey size of the four predators for whitefish are approximated in Table 1. Note that each predator is found in only a subset of habitats.

Dispersal

Individuals can change their ecological niches within the lake. We assume that, each year, after reproduction each fish chooses the ecological niche to go to with probabilities proportional to their overall survival in each niche (accounting for both density-dependent selection and predation). One interpretation of this assumption is that each fish samples different niches before deciding on the one to stay. Fitness-dependent migration has been used in a number of earlier ecological models[39-42] but has apparently been neglected in speciation modeling. Note that after each dispersal event, the niche population sizes change, so niche “attractiveness” can change as well. In numerical implementation of the model, to avoid a bias due to a sequence of events during the dispersal, individuals are chosen at random (without replacement, so that individuals move, at most, once per time step) to make a dispersal move.

Reproduction

We assume that mating is assortative by the size and the niche. [This modeling choice can also accommodate a scenario where fish do not mate within their physical niches but rather mate at a niche- and size-specific mating ground or mating time.] Each surviving female chooses a mate at random from an appropriate set of males. She then produces a random number of offspring sampled from a Poisson distribution with size-specific means and . All offspring are born in the same niche as their parents.

Life cycle

We assume there is a yearly sequence of events in the life of every fish. It starts with random death of eggs which is followed by density-dependent survival which is followed by predation mortality. The surviving fish grow or reproduce. Then all fish, except the newborn, have the opportunity to disperse among the niches within the lake, before this cycle starts again.

Numerical simulations

Scenarios

We simulated two scenarios of diversification differing in initial conditions. The Colonization-L scenario reflects the current view of the post-glacial colonization of lakes by large littoral fishes[24]. In this scenario, the initial population of size 3,000 eggs is introduced in the littoral niche to which it is adapted. Its modal trait values are: , , and , that is, fish grows to the “large” stage 2 only. The Colonization-G scenario is similar as Colonization-L scenario, but all initial colonizers are assumed to grow bigger, i.e. reach the “giant” stage three (, , and ).

In both scenarios, the initial population harbors genetic and phenotypic variation. The standard deviation of the gill raker number is five (which is close to values observed in natural populations). In the genes controlling the probabilities of maturation, one allele has frequency 95% and another has frequency of 5%. After introduction, the population then evolves for years. We then evaluate if and how many different morphs emerge and are maintained in the lake.

Parameter values

In numerical simulations, besides the initial conditions, we also varied the predation intensity . Specifically, for each of the four predators we set at , , and . The value of parameter measuring the strength of stabilizing selection was set to , which corresponds to moderately strong selection.

In our model, fertility parameters are set to respectively, following an exponential relationship with length[43,44]. Note that these numbers should be interpreted not as an actual number of eggs produced by a fish but rather as a number of offspring surviving to the moment when they are subject to selection.

To estimate carrying capacities (used in Eq. (2)) we used available data assuming that extant fish are adapted to their environment, i.e. that they have the optimal number of gill rakers and reach maturation at a proper size ( for the littoral niche, for the pelagic niche, and , for the profundal niche). We also assumed that extant fish in three niches are reproductively isolated and that fertility is a function of the female size[43,44].

We set the equilibrium densities of the largest morph in each niche to 1,000, 6,000, and 150 individuals for the littoral, pelagic, and profundal niche, respectively. These values are proportional to the observed densities in each niche[45,46].

From these values and empirical data, we estimated the corresponding equilibrium densities for other morphs (Table 2) and the values of corresponding carrying capacity parameters (Table S1 in the Appendix).

Table 2

Equilibrium population densities N for fish of different sizes in different habitats in the baseline model.

Stage	Littoral	Pelagic	Profundal
0 (2 cm)	15125	5517	656
1 (15 cm)	3025	6000	66
2 (20 cm)	756	0	150
3 (30 cm)	1000	0	0

The carrying capacities of different niches vary among different lakes. A particularly important factor affecting the likelihood of diversification is the carrying capacity of the profundal niche which can vary because of the size of the lake, its shape, and the nutrient deposition rate. In this study we used three different carrying capacities for the profundal niche resulting in equilibrium densities equal to one, two, and four times the values given in Table 2. We kept carrying capacities of the other two niches constant.

Overall, the results we present here explore different combinations of parameters, assuming assortative mating and fitness-dependent dispersal, for each of which we did 10 independent stochastic runs. We also performed a similar investigation of four related sets of models: 1) with random (rather than size-assortative) mating and random (rather than fitness-dependent) dispersal, 2) with size-assortative mating and random dispersal, and 3) with random mating and fitness-dependent dispersal. We also investigated our model with 4) smaller initial standard deviations in gill raker number and lower variation in the maturation loci. We never observed the emergence of the profundal morph in these models so we do not present the corresponding results.

Evolutionary outcomes

To interpret our numerical results, we say that a given niche has been successfully colonized if in the last year the population of newborn both 1) has gill rakers number adapted to their niche and 2) they mostly grow to a size which characterizes the niche in real lakes.

The former assumption was formalized as the requirement that the average gill raker number is close to the optimum value, specifically: . (Recall that is a parameter measuring the strength of stabilizing selection.) The latter assumption was formalized as requirements on the evolved average maturation rates and . Specifically,

if (i.e., fish mostly matures at size 1), then we say the fish is pelagic;

if , (i.e., fish mostly matures at size 2), then we say the fish is profundal;

if , (i.e., fish mostly reaches size 3), then we say the fish is littoral.

Results

Here we present our main results on the conditions under which each niche was colonized and by how many morphs. [Our results for the entire parameter space studied can be found here: 10.18710/PI8PJQ]

When we look across the entire parameter space, the majority of lakes fall into four different compositions of morphs present: littoral morph only, pelagic morph only, littoral and pelagic morphs, and all three morphs present (Tables 3 and 4). In most cases, the system reached a stochastic equilibrium state on the time scale of a few hundred generations (see the link above for some examples). Thus, the speed of diversification in numerical simulations is comparable with the results of ref. [47] whose estimates of the time of divergence in whitefish is on the order 1000–2000 generations in three subarctic watercourses. It is also comparable with that in earlier speciation models (e.g., refs. [7,11,48]). Below we will break down the conditions under which each outcome was more likely to be observed.

Table 3

Number of simulations with different compositions of morphs present for different relative carrying capacities (1 = baseline level, see Table 2), of the profundal niche in Colonization-L scenario.

Relative carrying capacity	1.0	2.0	4.0
None	2	2	1
Monomorphic	479	495	509
Littoral	439	477	501
Pelagic	39	17	8
Profundal	1	1	0
Dimorphic	329	313	295
Littoral Pelagic	329	313	279
Profundal Littoral	0	0	16
Profundal Pelagic	0	0	0
Three morphs	0	0	5

The data are summed up across all predation combinations.

Table 4

Number of simulations with different compositions of morphs present for different relative carrying capacities (1 = baseline level, see Table 2), of the profundal niche in Colonization-G scenario.

Relative carrying capacity	1.0	2.0	4.0
None	0	0	0
Monomorphic	425	484	520
Littoral	422	482	519
Pelagic	3	2	1
Profundal	0	0	0
Dimorphic	385	326	290
Littoral Pelagic	385	326	282
Profundal Littoral	0	0	8
Profundal Pelagic	0	0	0
Three morphs	0	0	0

The data are summed up across all predation combinations.

Only littoral morph (no diversification)

Survival of the littoral morph with no other morphs emerging was the most common outcome in our simulations (Tables 3 and 4 and Figs. 2 and 3). It typically required the presence of predation by trout. In this regime, the population density of the littoral morph is close to the carrying capacity in the littoral niche but is relatively low in two other niches in which we observe mostly migrants from the littoral niche. Lakes with only the littoral morph are common in Fennoscandia[24]. Their abundance in the geologically younger region is one of the reasons that we hypothesize that the littoral morph is ancestral to Fennoscandia lakes.

Figure 2

The frequencies of different outcomes in the Colonization-L scenario for different predation rates. 30 simulations for each parameter combination. The numbers are also reflected in the size of the font used. (a) Only littoral morphs. (b) Both littoral and pelagic morphs. (c) Both littoral and profundal morphs. (d) Only pelagic morph. (e) All three morphs. The data are summed up across all carrying capacities of the profundal niche.

Figure 3

The frequencies of different outcomes in the Colonization-G scenario for different predation rates. 30 simulations for each parameter combination. The numbers are also reflected in the size of the font used. (a) Only littoral morphs. (b) Both littoral and pelagic morphs. (c) Both littoral and profundal morphs. (d) Only pelagic morph. No cases with all three morphs were observed. The data are summed up across all carrying capacities of the profundal niche.

Two morphs

The emergence of the pelagic morph and the maintenance of the littoral morph was a frequent outcome in both colonization scenarios (about ). This outcome occurs more often when trout predation was absent or small (Tables 3 and 4 and Figs. 2 and 3). In this regime, littoral and pelagic niches are filled by the corresponding morphs close to carrying capacities. In the profundal niche, we observed migrants from the littoral niche. There were also a few cases () where the divergence occurred toward the profundal morph, instead of the pelagic (Tables 3 and 4 and Figs. 2c and 3c). In this regime, littoral and profundal niches are filled by the corresponding morphs close to carrying capacities. In the pelagic niche, we observed migrants from the littoral niche. All these cases were observed when the profundal environment was large and no burbot predator was present. This latter outcome is usually not observed in Fennoscandian lakes.

Single pelagic morph

Although lakes with only the pelagic morph present are not observed in Fennoscandia[24], they were observed in our simulations. This outcome implies the extinction of the littoral morph after it gives rise to the pelagic morph. This outcome was observed in both colonization scenarios under high predation by pike and perch but much more common in the Colonization-G scenario (Tables 3 and 4 and Figs. 2d and 3d). In this regime, the population density of the pelagic morph is close to the carrying capacity in the pelagic niche but is relatively low in two other niches in which we observe mostly migrants from the pelagic niche.

Three morphs

The cases with all three morphs present were very rare in our simulations - cases only. In this regime, the niches are filled by the corresponding morphs close to carrying capacity. They all happened in the Colonization-L scenario under large carrying capacity for the profundal niche (Table 3) when predation from pike and perch was strong but that from trout and burbot was weak (Fig. 2e). Lakes with all three morphs are rarely observed in Fennoscandia[24].

Failed adaptation

Four simulations resulted in populations with intermediate traits according to our metrics. These failed adaptations all involved high predation. Additionally, in one simulation with strong predation, the population went completely extinct.

Overall, both monomorphic and polymorphic lakes emerged from our simulations. The outcomes of Colonization-L and Colonization-G scenarios were similar, but Colonization-L scenario was more prone to diversification. The reason is that Colonization-L scenario started with the gill raker numbers typical for the littoral morph but with the size close to those for the pelagic and profundal morphs. This simplified the emergence of these two morphs. The lakes with all three niches colonized were observed but very rarely.

Discussion

We used individual-based simulations to study the likelihood of within-lake ecological and morphological diversification in Fennoscandian whitefish in the post-glacial time frame (about years). Our simulations show a common emergence of lakes where, in addition to the originally colonizing littoral morph, either the pelagic or profundal morphs have evolved and become established. The former outcome was much more frequent than the latter. We also observed the presence of lakes with all three morphs present, although very infrequently. However, this result fits with the observations from the wild which have revealed only a handful of trimorphic systems.

We started our simulations with a small founder population adapted to the littoral niche and no individuals adapted to the two other available niches: pelagic or profundal. Colonizing the pelagic niche required individuals to evolve increased gill raker number and smaller size. Colonizing the profundal niche required individuals to evolve decreased gill raker number and smaller size. The formation of new morphs typically happened on the time-scale of a few hundred generations. This is similar to other models with some initial genetic variation and strong selection[7,48,49].

Speciation theory tells us though that the mere existence of “empty” ecological niches does not guarantee that they will be colonized and “filled” by locally adapted organisms, especially if there is a possibility of gene flow from the ancestral niche. Some additional factors must usually be in place to simplify survival in the new environment and reduce the effects of deleterious gene flow preventing local adaptation. Several such factors turned out to be very important in our model as well.

First, populations must have sufficiently large initial genetic variation. In our simulations presented above, we set the initial standard deviation of the gill rakers at . With a lower standard deviation (and, correspondingly, lower initial genetic variation), the state with both littoral and pelagic morphs was possible, but the profundal morph never emerged. Second, diversification required mating to be assortative with respect to both size and ecological niche. Third, dispersal had to be fitness-dependent. Both these assumptions result in reduced effects of gene flow from the ancestral niche. With random mating and/or dispersal, empty niches remained largely empty. Colonization of the profundal niche was mainly observed when its carrying capacity was the largest. Biologically, this describes lakes that are large and deep. Large carrying capacity simplified survival in the new niche by effectively reducing within-niche competition. These four assumptions alone were still not enough to ensure diversification. The fourth important factor was the elevated predation pressure in the littoral niche. Colonization of the pelagic niche largely required the absence of brown trout which is the only pelagic predator in model assumptions. Correspondingly, colonization of the profundal niche largely required the absence of burbot which is the only predator that can be found there. By starting to utilize the niche with no predators, the fish gets an immediate additional survival advantage which does not require evolving any genetic adaptations. Moreover, the absence of predators in a niche implies that fish does not have to grow large and thus could evolve earlier maturation. The latter in turn would lead to faster population growth.

The appearance of the pelagic morph in addition to the littoral one was much more common than that of the profundal morph, which is in accordance with findings from empirical studies[24]. This happens because in our model, reflecting the situation in most lakes, the carrying capacity of the pelagic niche was much larger than that of the profundal niche making the pelagic niche much easier to colonize. Diversification was promoted if initial colonizers were littoral fish of “large” (Colonization-L scenario) rather than “giant” size (Colonization-G scenario), which happens because the size of the former is closer to those of the two other morphs (that are smaller than the littoral) and thus much smaller evolutionary change is required. However, body size is very plastic in whitefish and influenced by resource availability and strength of competition (e.g.[29,38]). The Colonization-G scenario may nevertheless be relevant as large size likely is beneficial for migration. Following this scenario, “giant” colonizers may initially have experienced untapped resources and fast growth rates, but subsequently turned into the “large” category due to reduced growth rates from increasing fish density and intraspecific resource competition. The modeling results corroborate with a recent empirical study from northern Fennoscandian whitefish indicating that body size and number of gill rakers are both targeted by natural selection[27].

Next, we discuss whether the factors identified in our model are present in natural populations of whitefish. Large genetic variation in colonizing individuals may have derived from the ancestral populations in the postglacial refugia, hybridization events during the loss of their pre-glacial habitats, and ongoing hybridization in contemporary times. Hybridization has an important role in historical and contemporary evolution of whitefish[50,51]. Whitefish is well-known for its high capacity for both fast forward and reverse speciation, which suggest a role of hybridization behind the high genetic variation harbored by the whitefish[22,26]. In the current modeling approach the minimum variation of five gill rakers of the ancestral morph was needed for the divergence to other morphs. Such variation of gill raker number is typical for generalist whitefish throughout much of the distribution suggesting ecological relevance of the model assumptions. However, the very few lakes with all three morphs in the simulation model overlooks the early individual specialization to all three habitats by the ancestral colonizing morph[52]. Such step is very likely a crucial initial step towards full-fledged ecological speciation described as five model simulations leading to occurrence of all three morphs. While the proportion of simulations leading to three morphs is around 0.6%, it could be indeed comparable to natural conditions where this morph is present only at very few large lakes having wide profundal habitat coverage and many predators in littoral and pelagic niches[28,53]. Furthermore, frequent occurrence of pelagic and littoral morphs (36–40% of simulations) obviously refers to frequent divergence to habitats providing main energy sources (pelagic phytoplankton and littoral algae/macrophytes) to the lake food webs. In contrast, the profundal habitat is an unproductive habitat and dependent on energy inputs (settling materials) from the pelagic and littoral zones. These pelagic and littoral inputs vary among systems due to lake morphology and productivity and this complexity makes whitefish divergence to this habitat rare[24,53].

Size assortative mating is common in fish[54-56]. Body size and size at sexual maturation of whitefish morphs are highly correlated with their resource use. Specifically, specialization to littoral benthic prey leads to large size and late maturation, pelagic prey to smallest size and fastest maturation, and finally profundal benthic prey to intermediate body size and late maturation[22,23,37,38]. Such differences in body size and maturation would be a logical step towards size assortative mating, if similar sized mature individuals among morphs are rarely present at the spawning sites. At least in the pelagic morph this seems likely as the very high mortality of small pelagic morph suggest that very few individuals are able to survive and grow to the size when other morphs mature[37,38]. Reproductive isolation may be further strengthened by temporal and spatial differences in time and place of spawning[57,58]. Such differences in spawning may arise via the habitat specific differences in temperature that regulates both prey resource availability and initiation of spawning[32,59]. In Fennoscandia, the littoral benthic resources are at the highest from early to mid-summer, pelagic zooplankton at mid to late summer and profundal benthic resources at late season[32,37,38]. The littoral reaches the highest water temperature in summer, but cools down at the earliest time followed by the pelagic and profundal[60]. Differences in realized water temperature among the whitefish morphs may be a major driver of divergence in Fennoscandian whitefish[60,61], as divergent temperature regimes alter the maturity status of all three morphs. There is also some evidence of spatial divergence of spawning sites among the morphs and in many salmonids such spawning site fidelity maintains population divergence[62-64]. In our model, we assumed that individuals exhibit certain habitat preferences. In fish, those can take different forms such as natal environment imprinting, condition dependent, density dependent, and predator avoidance[65,66].

In the model, the lack of pelagic (brown trout) and profundal (burbot) predation was a needed prerequisite for the rise of both pelagic and profundal whitefish morphs. Such conditions where strong predation occurs solely on littoral habitats of large lakes seems unlikely in contemporary conditions[28,30]. However, in the large lakes without polymorphic whitefish, brown trout and other predators (pike, burbot, perch) indeed used mostly littoral habitat[53]. It is likely that the rise of pelagic whitefish morph induces brown trout shift to pelagic habitat use and piscivory[31,53]. A similar process could be present with regard to burbot, where emergence of the profundal morph provides a new forage fish. However, burbot is a dark active predator that frequently use diel bank migration to feed on more abundant littoral prey resources[67]. Thus, the empirical data suggest that modeling requirements for predation could be met in nature, but such conditions are rare.

Due to computational considerations, our model has some obvious limitations. For example, we assumed that fish laid a relatively small number of eggs. In reality, fish can produce a very large number of eggs; much larger than what we could simulate using an individual-based approach. We expect that with larger number of eggs, selection will be more efficient and divergence may occur faster than what is observed here (cf.[68]). The initial conditions were centered only on scenarios of “giant” and “large” colonizers, but it is very difficult to determine the actual body size of the ancestor starting to diverge. However, increasing intraspecific competition for resources is likely after the colonization of a new lake, which tends to shrink the body size supporting shift from giant to smaller body size. Predation by multiple species had a strong effect on whitefish divergence, but in the wild there are lake systems with whitefish morphs without main predators or very low amount of predation[27]. However, the majority of lakes with three morphs have abundant predator populations and predation is likely to have strong influence on life-history divergence of prey. It is also likely that prey and predators co-evolve during the divergence process[69] but individual-based modeling of multi-species evolutionary processes is inherently difficult. The process of building the model has also revealed important gaps in the available data on whitefish populations and their environments. In particular, having more precise estimates of the population sizes, predation rates, and fecundity of different morphs would increase the power of our model.

In our model, the body size was subject to direct selection by predation and also due to fertility and maturation rate differences. We did not consider explicitly body size adaptation to the ecological niche. However since body size correlates with gill raker number[23,25,70], our model partially captures this effect.

There are a number of additional directions our model could be extended. For example one could study evolutionary dynamics on much larger temporal scales than used here and using a finer grid of numerical values to identify threshold values of parameters separating different dynamic outcomes. For simplicity we assumed constant predation pressure. In general, evolutionary changes in the the prey population can affect the predator density and the strength of predation experiences by the prey. We leave studying the effects of this feedback for future work.

The sizes and depth of lakes in Fennoscandia vary a lot as well as their productivity with evident implications to species divergence[24]. In our simulations, we considered only three major lake habitats and assumed pre-determined carrying capacities for each habitat. Lake morphometry and productivity in addition to deglaciation history likely fine-tunes the divergence processes further[24], but these are far too complex scenarios to model conclusively. Further modelling effort is needed to understand vertical pelagic divergence documented for example to Alpine whitefish[71] and North American ciscoes[72]. Also, it is well established that a small number of loci of large effect are more conducive to speciation than a large number of loci with small effects[7,49,73]. We have chosen a small number of loci to run our simulation to facilitate the process of adaptive radiation. With larger number of loci, we would not expect much diversification to be observed in our model. Given that diversification in whitefish does happen, we expect that the underlying traits are controlled by a small number of loci with large effects.

Overall, our modeling supported the possibility of divergence to three lake habitats during the postglacial time-frame, although such cases were rare and required a large profundal habitat without predators. Divergence to pelagic and littoral morphs was much more frequent and occurred with various levels of predation. The modeling effort indicated how little actually we know about reproductive isolation mechanisms as well as about habitat-specific carrying capacities and essential population dynamic parameters such as age-specific mortality in whitefish. Nevertheless, it would be interesting to apply our modeling approach to a true deep water species, such as Lake Baikal Cottus and various members of Salvelinus, in very large lakes using species-specific parameters[74-76]. Recent empirical studies have highlighted the important role of deep water habitats which have been previously neglected in monitoring efforts[77,78]. In these respects, mathematical modelling could provide an efficient predictive tool for finding new lakes with three morphs, as well as identification of key factors contributing to their disappearance from the postglacial lakes.

Supplementary information - appendix.