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.


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.
Individual whitefish. Each individual is characterized by three phenotypic traits: size s, the number of gill rakers x, and whether or not it is sexually mature. We discretize size s 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 m s depending on its current size s, or grows to the next size class + s 1, with probability − m 1 s . Newly born individuals can only grow but not mature (i.e., = m 0 0 ), whereas individuals reaching the largest size 3 always mature (i.e., = m 1 3 ). 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.
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 m 1 and m 2 are additively controlled by L diallelic loci each. (In numerical simulations = L 4). The corresponding allelic effects are scaled so that m 1 and m 2 stay between zero and one. The number of gill rakers x is also additively controlled by L loci but each locus has multiple discrete alleles. These alleles are subject to step-wise mutation and the effects are scaled so that trait x can www.nature.com/scientificreports www.nature.com/scientificreports/ take only integer values and mutation changes x only by plus and minus one. We assume that all genes are physically unlinked. Mutations happens with probability µ = − 10 5 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 x controls the food available to fish in a given environment. We define the condition ω x ( ) j of a whitefish with x gill rakers in niche j as Here θ j is the optimum gill rakers number in niche j and σ 2 is a parameter measuring the strength of stabilizing selection on x towards this optimum. For whitefish, we set θ at 26,36, and 18 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 x gill rakers and of size s in niche j being Here N s j , is the total number of individuals of size s in niche j, and constant parameter K s j , is the population size at which a population of perfectly adapted individuals (i.e., individuals with perfect condition ω = 1 j ) would have a survival rate of . 0 5 (see the Appendix). We will interpret parameters K s j , as measures of carrying capacity of the corresponding sizes in the corresponding niches. The larger the number of competitors N s j , and the worse is the individual condition ω j , the smaller is the survival rate ν j . 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 g k . To describe predation, we posit that predator k removes a random proportion π k of surviving whitefish with sizes smaller than or equal to its gape size g k . Values π k 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][40][41][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 b b , 1 2 and b 3 . 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: = x 26, = m 0 1 , and = m 1 2 , 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 ( = ). 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 10,000 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 0%, 30%, and 50%. The value of parameter σ measuring the strength of stabilizing selection was set to 4, which corresponds to moderately strong selection.
In our model, fertility 3 are set to {4, 16, 64} 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 K s j , (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 ( = = m m 0 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 K s j , (Table S1 in the Appendix). 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 × × = www.nature.com/scientificreports www.nature.com/scientificreports/ 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 x 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 m 1 and m 2 . Specifically, (i.e., fish mostly matures at size 1), then we say the fish is pelagic; (i.e., fish mostly matures at size 2), then we say the fish is profundal; (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: https://doi.org/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.
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. two morphs. The emergence of the pelagic morph and the maintenance of the littoral morph was a frequent outcome in both colonization scenarios (about 40%). This outcome occurs more often when trout predation was   www.nature.com/scientificreports www.nature.com/scientificreports/ 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 (24) 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 -5 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 . 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 10,000 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 5. With a lower standard deviation (and, correspondingly,  www.nature.com/scientificreports www.nature.com/scientificreports/ 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 www.nature.com/scientificreports www.nature.com/scientificreports/ 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 www.nature.com/scientificreports www.nature.com/scientificreports/ 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][55][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][63][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 Scientific RepoRtS | (2020) 10:7394 | https://doi.org/10.1038/s41598-020-63684-3 www.nature.com/scientificreports www.nature.com/scientificreports/ 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][75][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.