The effect of landscape structure on the evolution of two alternative dispersal strategies

Dispersal is an important event for most organisms at least once in their life cycle. The evolution of dispersal can be influenced by local adaptation, landscape structure, and perceived temporal and spatial variation. The interaction between local adaptation, landscape heterogeneity, temporal variability and rules of dispersal may be more complex than previously assumed. Therefore, we sought to understand the influence of emigration rules and landscape structure on emerging dispersal rates and traits. Here, we implemented an individual-based model (IBM) of trait evolution in scenarios characterized by different landscape structures and different degrees of spatial heterogeneity and global temporal variation. Individuals could evolve two traits coding for their environmental niche (position of niche optimum and niche width), and two traits determining nearest-neighbor dispersal: an individual emigrates with a probability defined by the first trait (random emigration), but emigrates with certainty if the fertility expected in the patch of residence falls below a threshold specified by the second trait (habitat-dependent emigration). We note an interaction effect between dispersal strategy and spatial variance—lower emigration under habitat-dependent than under random emigration if spatial heterogeneity is low, but eventually a reversal of this ranking if heterogeneity becomes large. Landscapes with sharp transition of habitat attributes result in a high degree of spatial sorting, while fractal landscapes do not. Emigration rates are overall lowest, when spatial variation is highest. We conclude that emergent emigration rates are influenced more by landscape structure and spatio-temporal heterogeneity than by the emigration strategy. With the ongoing land use change more research into this topic could help highlight the difficulties species might face under the change from landscapes characterized by gradual transition zones to landscapes dominated by abrupt ecotones, the latter typical for agricultural and urban settings.

) while also being influenced by the level of local adaptation itself (Kisdi et al. 2020).
The simplifying assumption of global dispersal is widespread in modeling approaches. It assumes that dispersal distances are large in relation to spatial structure of the landscape, but disregards the tripartite composition of the dispersal process. Individuals or propagules have to leave their natal patch (emigration), move between patches, and immigrate into a new patch (Bowler and Benton 2005;Crook et al. 2020). Emigration can be population dependent (e.g., density-dependent) or statedependent (e.g., fitness or life stage-dependent) (Bona et al. 2019;Cronin et al. 2020). Habitat choice is an important immigration strategy and was shown to evolve faster than local adaptation (Kisdi et al. 2020). Habitat choice can also accelerate local adaptation (Camacho et al. 2020;Jacob et al. 2018). Overall, dispersal can have opposing effects on local adaptation: It can either increase fitness by enhancing local adaptation via immigration of individuals that choose patches matching their niche attributes Kisdi et al. 2020). It can also reduce mean fitness, due to the immigration of maladapted individuals, which changes mean trait values and/or trait variation. Finally, dispersal can increase fitness by increasing the genetic variation, which increases evolutionary potential (Bridle et al. 2019). Which of those effects dominates in a population depends on the spatial structure of the landscape under consideration (Bridle et al. 2019;Jacob et al. 2019).
The most apparent way in which landscape structure can impede dispersal is landscape fragmentation (Barnes et al. 2015;Jacob et al. 2020). If suitable patch types are too far apart to be easily reached, selection for dispersal can weaken (Kubisch et al. 2013;Sinai et al. 2019). Even in non-fragmented landscapes, landscape (viz. patch type) heterogeneity can influence the evolution of dispersal because an organism's mobility determines the environmental variation it encounters during its lifetime (Bridle et al. 2019;Kaemingk et al. 2019;Schiffers et al. 2014). If individuals disperse far (search widely) and frequently in heterogeneous landscapes, they will encounter a wider variety of environmental conditions. This high variety of environmental conditions can lead to diminishing fitness, especially in specialized species. Landscape heterogeneity is thus usually considered to be detrimental to the evolution of frequent dispersal (Bridle et al. 2019;Kaemingk et al. 2019;Kubisch et al. 2013;Schiffers et al. 2014;Sinai et al. 2019). In principle, landscape heterogeneity could also be beneficial since it increases the chance to encounter new patches where the conditions fit the individual's niche optimum better. Without informed or habitat-orientated immigration this is, however, not likely to happen if the favored patch type is itself rare (Clobert et al. 2009;Kubisch et al. 2014). Consequently, landscape homogeneity will generally lead to increased dispersal probabilities since dispersal is less penalized by an increased possibility to encounter unsuitable patches. However, in this case, population density could be similar across the whole landscape and therefore decrease the expected differences in fitness expectation. This would eventually leave avoidance of kin-competition as the primary benefit of dispersal (Poethke et al. 2007).
Dispersal can also be a bet-hedging strategy against adverse conditions imposed by temporal variability. Particularly, under extreme conditions dispersers may find a patch that is, in this particular period (e.g., year, season), more similar to the optimal conditions of an individual than its natal patch, even though that habitat might not provide optimal conditions under average conditions. We show in a previous study that-especially in patches that deviate far from the landscape mean-individuals evolve high dispersal probability when spatial variation is high; this bet-hedging strategy leads to higher fitness expectations of dispersers than philopatrics in extreme years (Sieger and Hovestadt 2020). The benefit of such a strategy could further be advanced by selectively emigrating whenever fitness expectations become low. Obviously, for generalist individuals that do not experience a pronounced fitness-decline in a wider array of patch conditions, these arguments are generally of less relevance. Schiffers (2014) showed in a simulation study that not only the heterogeneity at large, but also the geometry of a landscape can influence the evolution of dispersal. In checkerboard landscapes, with distinct edges between patches, the grain size (i.e., the size of a cluster of patches with the same attribute) highly influenced the speed of range expansion, as well as other population characteristics. In gradient landscapes the evolved dispersal distance depended on the length of the gradient: steeper gradients favored shorter dispersal distances, while shallower gradients selected for longer distances. Furthermore, edge effects (Kaemingk et al. 2019), as well as patch connectivity (Cronin et al. 2020;Fobert et al. 2019;Masier and Bonte 2020;Schwarzmueller et al. 2019) and structure of the natal habitat (Ducros et al. 2020) have been found to shape dispersal syndromes. Overall, this leaves the impression that landscape structure can strongly influence dispersal.
The insights and findings mentioned above suggest that the interaction between local adaptation, landscape heterogeneity, temporal variability and rules of dispersal may be more complex than previously assumed. Here, we examine such interactions and especially investigate the influence of emigration rules on evolving dispersal traits and rates. In particular, we hypothesize that: (i) emigration decreases with increasing spatial variation, regardless of emigration strategy; (ii) existence of bigger habitat clusters, viz. stronger spatial autocorrelation leads to more emigration (into similar habitat); (iii) steep transition zones with (clear) habitat edges promote spatial sorting of dispersal phenotypes with lower emigration from ecotones, and (iv) habitat-dependent emigration results in lower overall emigration than random emigration. We expect the disadvantages of emigration to be alleviated when emigration only occurs when fitness expectations in the natal patch are low.

Overview
For this study, we expanded the metapopulation model of annual, haploid individuals already described in Sieger and Hovestadt (2020) to a set of simulations in differently structured landscapes and implemented different dispersal strategies, using the programming language Julia (Bezanson et al. 2012). Each patch of a landscape (gridcell) is characterized by a certain mean habitat value (e.g., a certain cover type) so that the landscape exhibits spatial variability in patch features. These values can be interpreted as reflecting spatial variance in habitat conditions like mean temperature or annual precipitation, but every other continuous environmental variable varying over space is just as likely (e.g., soil nitrogen content, water oxygen content, pH, or salinity), even though such parameters are less likely to rapidly fluctuate in time. As explained in more detail below, habitat attributes additionally underlie global temporal variability, so that conditions in any patch are variable in time. We assume, however, that such variability is completely synchronized across the landscape as might be true for the regional impact of actual climatic conditions. Any patch houses one population of haploid organisms as described in Sieger et al. (2019). Individuals can disperse to other patches, with dispersal determined by either a single (coding for random dispersal) or two heritable traits, with the second coding for fertility-dependent emigration; we use fertility expectation (habitat-dependent expected number of offspring) as a proxy for fitness here. The dispersal traits and both niche traits mutate during inheritance, thus allowing for adaptation to simulated conditions. Further details will be explained below.

Landscapes and scenarios
The metapopulation covers a spatially heterogeneous landscape of 64 by 64 habitat patches wrapped into a torus. We use either "tiled" landscapes or fractal landscapes. The tiled landscapes are characterized by habitat clusters (aggregations of several patches) that share the same habitat value and have sharp edges between them. Such landscapes were created with the NLMR-package (Sciani et al. 2018) in R 3.5.3 (R Core Team 2018), using the nlm_random_ rectangular_cluster() tool. The second were created with an algorithm for autocorrelated landscapes developed by Chipperfield et al. (2011), with the level of autocorrelation specified by the Hurst index. With both algorithms, the patches had habitat values between -1 and 1. Here, the neighboring patches have similar but not exactly the same habitat values, and show a gradual transition in patch conditions. Consequently, borders between different patch types are more pronounced in the "tiled" landscape type, since the edges between habitat clusters are clear cut. On the other hand, neighborhoods (a patch and its eight surrounding neighbors) in fractal landscapes usually have a non-zero standard deviation of habitat attributes. In tiled landscapes, all patches of a habitat cluster share the same environmental mean, which means that any neighborhood inside the habitat cluster has a standard deviation of zero; only patches directly at the edge between habitat clusters have a non-zero neighborhood standard deviation.
We standardized each of the generated landscapes to the mean value of 0, by calculating the mean of the landscape and subtracting it from each patch's value. The landscapes used as fractal landscapes with small patch cluster sizes (Hurst index of 0.3) were the same as the ones used in Sieger and Hovestadt (2020). To achieve similarly variable landscapes for the other three types of landscape structure, we created more landscapes than needed with the above-mentioned algorithms (fractal with a Hurst index of 0.9 and tiled with bigger [minimum habitat cluster length minlv: 21 patches, maximum habitat cluster length maxlv: 32 patches] or smaller habitat clusters [minlv: 2 patches, maxlv: 12 patches]) and selected five for each category that matched the standard deviation ( σ S ≈ 0.32 ) of the five landscapes used in Sieger and Hovestadt (2020). To reduce variance that was unaccounted for, between simulation experiments we stored five replicated realizations for each of the four landscape types, i.e., 'fractal' or 'tiled' and 'small' or 'big' habitat clusters (see Fig. 1 for examples) for use in all the simulation scenarios explained below. To examine the influence of the relationship between temporal and spatial heterogeneity, we further created instances of the 20 landscapes with different spatial heterogeneity. In the original set of landscapes, the average spatial standard deviation was σ S = 0.32 . To achieve different relationships, we kept the temporal heterogeneity the same, but multiplied each patch's habitat attribute with either 2, 4 or 8 ( σ S ∈ 0.32, 0.64, 1.28, 2.56 ). This resulted in 80 distinct landscapes.
Each landscape additionally experiences global temporal environmental heterogeneity: at the beginning of every time step t a single random value, drawn from a normal distribution with mean = 0 and σ T = 1, is added to any patch's mean environmental value to form the current environmental value of each patch. Consequently, local conditions vary in time, thus modulating the individuals' fitness expectations. For repeatability and to avoid unaccounted variance we created three different time series of environmental heterogeneity, which are utilized in all of the simulation experiments described below. All following landscape scenarios were thus based on 5 × 3 replicated simulation runs per type of landscape structure. Each landscape was used for modeling the metapopulation for three replicates, with the respective vector of temporal variation, each for the time span of 300 generations. This was done for both implemented dispersal strategies and a strong ( α = 2 ) or weak ( α = 4 ) trade-off (see below).

Life-cycle and population dynamics
After local patch conditions have been set to the new actual condition in each time step, newborn adult individuals first "decide" whether to disperse. In all simulation scenarios each individual carries a heritable and mutable trait d that codes for the probability to leave its natal patch. In some scenarios, however, a second dispersal trait f encodes an additional emigration fertility threshold: if the expected fitness (calculated according to Eq. 1; see below) falls below f the individual emigrates with certainty, otherwise it emigrates with its base probability p = d . In the latter case, an individual leaves the natal patch when a random number drawn from a uniform distribution U [0, 1] is lower than the individual's dispersal probability p . An emigrating individual either dies with a given dispersal mortality m = 0.1 or immigrates into a patch randomly selected from the eight neighboring patches. Depending on simulation scenarios, individuals therefore either exhibit two of the following emigration strategies: (1) random emigration, as implemented in many previous simulation or analytical studies of this type like (Comins et al. 1980;Kubisch et al. 2013Kubisch et al. , 2014Travis and Dytham 1999) with nearest-neighbor dispersal (referred to as NN), or (2) random nearest-neighbor dispersal combined with a patch (fitness)-dependent emigration criterion (referred to as HE). We did not include habitat choice of immigrants or scouting into our dispersal strategies.
After dispersal, density-independent but patchdependent reproduction of the N j adults i in each patch j takes place. The fit between environmental conditions and the individual i 's niche determines its reproductive success: the individual's expected fertility is defined by two heritable traits, the position of the niche optimum h i in the environmental space and its niche width (habitat tolerance) g i . Combined, the two traits define a normal distribution for the expected fertility around the niche optimum. The expected number of offspring for each adult i is thus calculated, with inclusion of a generalist-specialist trade-off term (Eq. 1), following Chaianunporn et al. (2015). The resulting, environment dependent expected number of offspring for adult i with traits h i and g i in patch j , L H j,t , g i at time t is calculated as: with R 0 the maximum possible offspring number. The specialist-generalist trade-off is calculated as: Example landscapes for each of the four landscape types generated (from scenarios with very high spatial heterogeneity, σ S ≈ 2.56 ). From left to right: fractal with big habitat clusters (Hurst index = 0.9), fractal with small habitat clusters (Hurst index = 0.3), tiled landscape with big habitat clusters (minimum habitat cluster length minlv: 21 patches, maximum habitat cluster length maxlv: 32 patches) and tiled with small habitat clusters (minlv: 2 patches, maxlv: 12 patches). Each landscape type group is presented by 5 replicates that were used in all simulation experiments. The legend shows the habitat value's deviation from the landscapes' mean (0) (see Initialization and scenarios) Note that larger values of α imply lower trade-off costs. When the difference between h i and H j,t is low or the value of g i is high, L i H j,t , g i is also high. The actual number of offspring (larvae) born by each adult i in patch j is generated by drawing from a Poisson-distribution with mean L i H j,t , h i , g i . After the birth of offspring, the adult population dies.
The total number of larvae L j, t produced in patch j at time t then undergo density-dependent mortality, with survival probability calculated according to the Beverton-Holt model: with a = R 0 −1 K ·R 0 and K the carrying capacity. This survival probability is used to allocate a random binomial factor to each individual offspring indicating whether it survives or not; the surviving larvae constitute the new adult population of the next generation. One time step t therefore equals one generation.
All four trait values of an individual are inherited from the parent and evolve by mutation and selection. Evolution of the niche optimum and dispersal traits are not penalized, but according to equation (Eq. 1) enlarging niche width underlies a trade-off of different strength (parameter α ), depending on the scenarios: It is either weak ( α = 4 , i.e., evolution of a higher tolerance diminishes maximum fitness at the optimum only slightly), or strong ( α = 2 ), leading to a bigger drop in maximum fitness when tolerance increases. The traits of each individual mutate separately and in each generation according to the following rules. A value drawn from a normal distribution with mean 0 and standard deviation 0.03 is added to the niche optimum inherited from the parent. The tolerance trait value must be restricted to a range g i > 0 ; for mutating the trait we thus multiply the parental trait with a value drawn from a uniform distribution between 0.97 and 1.03. The dispersal probability d and the fertility threshold f are changed additively by adding a value drawn from a normal distribution with mean 0 and standard deviation 0.001. Fertility traits of f < 0 will lead to the same results as f = 0 . Values for f ≤ 0 imply that the individual would never emigrate because of low fitness expectations. Values for d are allowed to take values outside the range [0, 1], but the dispersal routine implemented treats dispersal with d < 0 as d = 0 and values of d > 1 as d = 1 . Note that in the NN scenarios the fertility threshold f becomes a neutral trait as it does not affect the dispersal behavior. (3)

Initialization and scenarios
All 64 × 64 patches in a landscape were initialized with 100 individuals each. Each individual was initialized with a niche optimum drawn from a normal distribution with the environmental attribute of their patch of placement as mean and a standard deviation 0.2. The niche width (tolerance) was drawn from a log-normal distribution with σ g and µ g , which were calculated from the results for the last generation in the simulations of Sieger and Hovestadt (2020) for each trade-off strength, to speed up adaptive evolution. To calculate parameters µ g and σ g of the log-normal distribution from the evolved trait values g i , the means m and variance v of g i were inserted into the arithmetic moments of the log-normal distribution for mean (first moment, Eq. 4) and variance (second moment, Eq. 5): and The starting dispersal probability d for each individual is 0.2, while the fertility threshold f is drawn from a uniform distribution [0, R 0 ] , with R 0 = 10 . The carrying capacity K of each patch is 1000 individuals.

Analysis
Graphical presentations of results were created using R (R Core Team, 2018) with the 'tidyverse' package (Wickham 2019). For each patch in each landscape scenario, the population means of all trait values, as well as the mean population size and fertility (as a proxy for fitness) were calculated and stored every fifth generation. Local adaptation and expected fitness were calculated before dispersal so that fitness values relate to the situation in the individuals' natal patch. Values of the random dispersal trait d < 0 were set to d = 0 for analysis since values lower than 0 induce the same effect as d = 0.
Additionally, we fitted a linear mixed effects model (lmer from the lme4 package (Bates et al. 2015), R version 3.6.3 [R Core Team 2018]) to estimate effect of emigration strategy (HE, NN), landscape_type (fractal, tiled), habitat cluster type (big, small), spatial_heterogeneity ( σ S ∈ 0.32, 0.64, 1.28, 2.56 ) on arcsine transformed emigration rate, dispersal probability trait, and proportion of emigration caused by the fitness threshold trait. We included landscape number (1-5) time series number (1-3) as random effects. For the fitness threshold trait, we fitted a lmer model using only the results from the HE scenarios because the trait is not meaningful in the NN scenarios. The fixed effects here were the landscape type (fractal, tiled), habitat cluster type (big, small), spatial_ heterogeneity ( σ S ∈ 0.32, 0.64, 1.28, 2.56 ). Results of lmer models were not used for statistical significance tests as this is not meaningful in modeling approaches but to estimate variance components for each of the experimental variables. Variance components were extracted using the "insight" R library (Lüdecke et al. 2019) using the "get_variance()" function with fixed, random and residual options.

Results
Since the results from scenarios with a weak and strong trade-off were qualitatively similar, we only show the data from scenarios with a strong trade-off. A weaker tradeoff leads to evolution of higher tolerance trait values, subsequently to higher fitness expectations and therefore to higher fitness threshold values. This, however, did not lead to fundamentally different results and will therefore not be discussed in detail. In our simulations, average individual habitat optima match the average local conditions in their natal patch; the absolute deviation of the trait h from average patch conditions is near 0. Even in extreme patches no adaptation gap can be found after 300 generations. This is the case regardless of landscape type, habitat cluster size, or degree of spatial variation. This also contrasts with our previous findings for global dispersal, where certain adaptation gaps formed in extreme patches (Sieger and Hovestadt 2020). On average, the expected fitness of an individual (number of offspring) depends only on the current habitat value and neither on the habitat cluster size nor on the landscape type (Additional file 1: Fig. S1). The evolving niche width was also unaffected by the different scenarios and landscape structures and only depended on the specialist-generalist trade-off strength (Additional file 1: Fig. S2).
The realized emigration rate is generally higher in tiled landscapes than in fractal landscapes, especially when habitat clusters are big. With increasing spatial variation, emigration rate decreases. This decline is least in tiled landscapes with big habitat clusters, i.e., the effect of spatial variation on emigration rate is comparatively small. In fractal landscapes with big habitat clusters and both fractal and tiled landscapes with small habitat clusters, the decline is much steeper with increasing spatial heterogeneity. The lowest overall emigration rates emerge in fractal landscapes with high Fig. 2 Mean realized emigration rate across dispersal strategy in the last generation of each simulation (NN, HE), landscape type, and cluster size. Blue indicates random nearest neighbor emigration (NN), green results for habitat-dependent emigration (HE). Light colors show the result for the fractal, dark shades for the tiled landscapes. The left column corresponds to big habitat clusters, the right to small habitat clusters. Global temporal variation stays the same ( σ T = 1 ). Emigration rate decreases with increasing spatial variation (from left to right) spatial variation and small habitat clusters (Fig. 2, right  panel).
In scenarios with random nearest neighbor emigration (NN), evolution of the emigration trait d was highly influenced by the different landscape and spatio-temporal scenarios. First, emigration rate declined as spatial patch variation increased. Second, the magnitude of this decline depended on both the landscape type and the habitat cluster size. Generally, the decline was less pronounced in tiled compared to fractal landscapes, resulting in evolution of substantially higher emigration rates in the former landscapes if spatial variation is large. The lowest rates emerged in fractal landscapes with small habitat clusters and high spatial variation, where almost all patches had a mean emigration rate d ≈ 0 . Except for scenarios with very low spatial variation (0.32), emigration rate was also higher in tiled landscapes (Fig. 4, top). In NN scenarios, the emergent emigration rate depended only on the emigration probability trait values d , which were therefore almost identical in value.
The spatial sorting of the emigration rate ( Fig. 3) also was visible in the random emigration trait d, especially with higher spatial variation and in tiled landscapes with big habitat clusters (Additional file 1: Fig. S3).
In the scenarios with patch-dependent emigration (HE), the random emigration trait d generally evolved to lower values than in scenarios with random emigration only (NN), regardless of habitat cluster size and landscape type (Fig. 4, top). However, with increasing spatial variation, the difference between HE and NN scenarios decreased. In scenarios with small habitat clusters and high spatial variation, the trait d was actually lowest in the NN scenarios with fractal landscapes, rather than in the HE scenarios. The trait value of the fertility threshold f was selected towards mean values of around f ≈ 0.75 , whereas average fertility values were fertility ≈ 6.4 . With increasing spatial variation, the fertility threshold decreased, with overall lower values in fractal landscapes as compared to tiled landscapes. Small habitat clustering also led to lower trait values of f , in either landscape type. Only with very low (0.32) spatial variation did the difference between landscape type and habitat cluster sizes become negligible (Fig. 4, center). However, the effects of landscape structure and spatial heterogeneity on f were generally less pronounced compared to their effects on the dispersal trait d.
A systematic difference between the observed overall emigration rate and the emigration probability trait d indicated the surplus of emigration that was caused by emigration due to expected fertility being lower than f (with some stochastic variation caused by the stochasticity of the random baseline emigration). Considering this, we estimated the proportion of all emigration events that were caused by expected fitness below the fertility threshold. With increasing spatial variation, emigration was proportionally more often caused by low expected fertility, especially in fractal landscapes (Fig. 4, bottom). However, it is important to keep in mind that this does not equal high emigration rates: the total emigration rates in those scenarios were comparatively low. In the unlikely event of emigration, it was, Fig. 3 Emerging emigration rates in each patch in an example landscape (the same landscapes and panel configuration as in Fig. 1). Top row shows the results of scenarios with random emigration (NN), second row with habitat-dependent emigration (HE). The value is lowest at the ecotones, especially in tiled landscapes (panel three and four), indicating spatial sorting however, more likely that emigration was caused by low fitness expectations than by random emigration.
For the emigration rate, the dispersal probability trait and the proportion of emigration caused by the fitness threshold trait, the proportion of variance explained by the fixed effects (dispersal strategy, landscape type, spatial heterogeneity) was drastically higher than that explained by the random effects (landscape replicate, times series; Table 1). However, the largest variance component was the residuals variance with 62-78%. This was not surprising because the residual variance includes Fig. 4 Top: evolved mean emigration probability trait value d (coding for random emigration) of each patch in the last generation under random nearest-neighbor dispersal (NN) or patch-dependent emigration (HE) across levels of spatial heterogeneity (low-very high). Center: evolved mean fitness threshold values of each patch in the last generation for scenarios with patch-dependent emigration (HE) across levels of spatial heterogeneity (low-very high). If an individual's expected fertility in the natal patch falls below the threshold trait, it will disperse in those scenarios (HE). Typical fitness values are Fertility ≈ 6.4 . Bottom: proportion of emigration events that were caused by expected fertility lower than the individuals' fertility threshold (see text for details). Panel configuration and colors as in Fig. 2  the between-patch effects within landscapes. Given the spatial sorting of the dependent variables (Fig. 3), a high explanatory power of patch identity was expected. For the fitness threshold trait, residuals also explained a large proportion of variance, but fixed effects explained only 5% while random effects explained 42%. This means that the trait is highly influenced by inter-patch effects, landscape identity and time series identity.

Discussion
In line with previous studies, our simulations demonstrate emergence of decreasing emigration rates as spatial variance compared to temporal variance increases. With increasing spatial variation, the difference in patch attributes between single patches and even more so between habitat clusters becomes more pronounced. Therefore, it is increasingly unlikely for individuals to be adapted to both the patch conditions in the natal and the (potential) new patch. Consequently, emigration into another patch becomes increasingly penalized due to the decrease in fitness expectations for the more or less maladapted immigrants (and their offspring), thus selecting for reduced dispersal (Hastings 1983;McPeek and Holt 1992). We find that the evolution of the tolerance trait g is dominated by the global temporal variance, but hardly affected by the spatial heterogeneity (so that we do not see selection for larger values of g in more heterogeneous landscapes), similar to our previous study with global dispersal (Sieger and Hovestadt 2020). Instead, selection reduces dispersal to avoid the risk and implicit cost of arriving in non-suitable patches. The latter effect might be overturned if local carrying capacities would be (really) small. In that case kin-competition would promote dispersal more strongly (Poethke et al. 2007) and thus might induce selection for larger tolerance g. Our simulations demonstrate, however, that landscape structure plays a role in how strongly selection on dispersal is affected by an increase in spatial patch variance. This was also supported by the high variance in dispersal due to inter-patch effects. We also note an interaction effect between dispersal strategy and spatial variance: lower emigration rates evolve under patch-dependent than under random emigration if spatial heterogeneity is low, but eventually a reversal of this pattern occurs if spatial heterogeneity becomes large. Here we find higher emigration rates under random than under patch-dependent emigration. Following a similar logic, we found that in fractal landscapes, neighboring patches likely have a different environmental mean, especially if autocorrelation is weak. Indeed, almost no neighborhoods had zero deviation between patches in fractal landscapes, whereas in tiled landscapes all patches of a habitat cluster share an identical environmental mean. This means that selection against dispersal with increasing global spatial variance is stronger in fractal landscapes since the dispersing individuals typically have a lower fertility than philopatric individuals. In tiled landscapes, at least in the centers of habitat clusters, no such penalty for leaving the natal patch exists since an emigrating individual likely encounters similar patch conditions as in its natal patch. The discrepancy between fractal and tiled landscapes is less prominent, however, when spatial correlation is weak. Here, the proportion of edge patches increases in the tiled landscapes (habitat clusters are smaller), so that the probability of a neighboring patch to be different from the natal patch is larger. Our findings thus confirm our second hypothesis that landscapes with large habitat clusters promote evolution of higher emigration rates, in particular if patches are arranged in distinct clusters of similar patch type. Comparable results were also found by several previous studies (e.g., Gros et al. 2006;Shaw et al. 2014).
In addition to the inter-landscape differences of emigration rate, we also observed, in line with our expectation, edge effects and spatial sorting, in particular when landscapes have a tiled patch structure. Because of the aforementioned steep habitat gradient at cluster edges (ecotones), lower emigration at edges is selected for as emigration across ecotones is severely penalized for locally adapted individuals. The relevance of this argument is supported by the emergent dispersal rates in tiled landscapes, where neighboring habitat clusters have rather similar environmental mean. This is for example visible in Fig. 1, third column, where the environmental mean of several habitat clusters is depicted by yellow shades (bottom right of the landscape). The similarity in color shows the similarity in environmental mean. At the same position in Fig. 3, the decline in emigration rate is less pronounced than at borders between habitat clusters with a bigger between-habitat cluster difference. In such circumstances, the above-mentioned penalty is smaller and more emigration between habitat clusters takes place. Comparable empirical evidence was reported by Vespa et al. (2018), who found that edge effects on dispersal were less pronounced when the difference between native forest and monoculture plantations was lower. They show that when plantations are new and therefore the difference between the plantations and the native forest is very high, seed dispersal into the plantations decreases. With increasing plantation age, the contrast between the plantations and the forest diminishes, leading to increased seed dispersal from the forest into the plantations. Similarly, dispersal at treeline ecotones is lower in both the boreal and the alpine environment than in forests at lower altitude or latitude (Crofts and Brown 2020;Kambo and Danby 2018;Ribeiro et al. 2019). The dispersal across the boundary between forest and alpine meadows encompasses a more drastic change in patch conditions than dispersal inside the forest, where patch conditions are more homogeneous. For specialized forests, species dispersal from forest to wooded corridors is also lower (Paal et al. 2020), which is comparable to the situation in our simulations, since all individuals are locally adapted and can therefore be considered specialists of their native patch.
When comparing the random (NN) with the patchdependent (HE) emigration scenarios, we recognize that in most scenarios lower emigration rates emerge under the latter; this can-in principle-be traced to the better ability of "informed dispersal strategies" to equalize fitness expectations in space (Poethke et al. 2015). However, we also find that selection on HE dispersal is less sensitive to the increasing spatial heterogeneity, so that at very high spatial heterogeneity similar (in the tiled landscapes) or even larger emigration rates (in fractal landscapes; see Fig. 2) evolve in the HE compared to the NN scenarios. This overall effect is mirrored by the substantial sensitivity of the random dispersal trait d (see Fig. 4, top) to increasing spatial heterogeneity, but the nearly complete insensitivity of the fertility threshold trait f (see Fig. 4, center). As a consequence, the proportion of dispersal events that were induced by low fitness expectation increases with increasing spatial heterogeneity, in particular in the fractal landscapes (see Fig. 4, bottom).
Our findings could also have implications for conservation and restoration efforts. With ongoing land use change, sharp edged patches-that are typical for intense agriculture-become more prevalent (Liira et al. 2008). This could mean that the trait distribution inside the habitat clusters also changes, from no spatial sorting of traits, like in the fractal landscapes, to the spatial sorting found in the tiled landscapes. This change in spatial distribution of traits in habitat clusters could enhance the isolation of remaining natural patches due to selection for decreasing dispersal probabilities at the edges of habitat clusters. This would be especially damaging for species that rely on several types of patches to thrive, for example anurans. It was already shown that anurans particularly suffer from the impact of habitat fragmentation (Homola et al. 2019;Ribeiro et al. 2019). The same mechanisms that inhibit dispersal across sharp ecotones could consequently hinder the movement of anurans between their larval and adult habitats.
We also show that small-scale spatial heterogeneity (neighborhood relationships) has a strong influence on local selection resulting in substantial variance within scenarios. In other words, even at small spatial scales, evolution of dispersal strategies is already affected by spatial variance, often enough more than by the largescale differences between landscapes.
The huge influence of times series identity on the evolution of the fitness threshold trait is not surprising, when considering the high variability in realized fitness in the three different time series. The impact of time series on fitness subsequently influences the selection for the fitness threshold. Time series with an overall low fitness also lower the selection for a higher fitness threshold because the possibilities to reach this higher fitness are lower.
In our simulations, we only implemented short distance dispersal, i.e., dispersal to patches in direct vicinity of individuals' natal patches; this implies the assumption that the dispersal capabilities of species are limited in relation to the scale of spatial heterogeneity (possibly also because of selection for limited dispersal distance). Limited dispersal distance is prevalent in various ecosystems and across ecotones, e.g., in halophytes in salt marshes (Polo-Ávila et al. 2019), southern Atlantic forest trees (Vespa et al. 2018), boreal forest plants and trees (Paal et al. 2020;Trant et al. 2018) and anurans (Ribeiro et al. 2019). However, rare long-distance dispersal is important the colonization of new habitat, e.g., in eastern larch (Larix laricina, Trant et al. 2018), salt marsh plants (Polo-Ávila et al. 2019) and lesser prairie-chicken (Tympanuchus pallidicinctu, Earl et al. 2016). It was also consequential in the speciation of neotropical lizards (Sheu et al. 2020). In a previous study (Sieger and Hovestadt 2020), we showed that long-distance (global) dispersal can lead to bet-hedging benefits-in particular for individuals adapted to extreme and rare patches-against the effects of temporal variance, when the spatial variation of a landscape is high. In the current simulations, there is little selection of this sort on dispersal because short distance dispersal typically takes individuals to rather similar patches; bet-hedging benefits could possibly emerge by dispersing across distinct patch borders, but other than in our previous study this is not fundamentally more likely to happen for individuals adapted to extreme and rare patches. It may thus be a fruitful approach in future simulations to implement dispersal with an evolving dispersal kernel or dispersal switching between short-and long-distance dispersal. It could be particularly interesting to create simulations that combine random local with patch or fitness-dependent long-distance dispersal: Searching for a very different-and usually distantpatch could be favorable if the natal patch becomes very unsuitable for the individuals' phenotype under the conditions of an extreme year. Such emigration strategies could be coupled with habitat-matching immigration to enable individuals to assess their expected performance in potential future patches and choose the most suitable one (cf. Camacho et al. 2020).
The presented simulation scenarios could additionally be enhanced by including several other aspects that were not included. In our simulations, no sexual reproduction took place, possibly slowing down the emergence of optimal trait combinations due to the lack of recombination. It was shown that sexual reproduction can be an important aspect in withstanding climatic stress. Trant et al. (2018), for example, show that in black spruce sexual reproduction adapts faster to climatic conditions than clonal reproduction. This could mean that the influence of spatial or temporal heterogeneity on the evolution of the environmental niche and dispersal could change if sexual reproduction were to be included. We also did not consider the other two phases of dispersal, immigration and establishment, which also play an important role in colonization of new habitat. It was, e.g., shown that successful establishment was the key factor for successful dispersal at treeline ecotones (Crofts and Brown 2020;Kambo and Danby 2018;Paal et al. 2020). In our model individuals that emigrate either immediately die or successfully immigrate and establish, the risks of all stages of dispersal are implicitly included. The evolution of the environmental niche and dispersal could of course also be influenced by competition with or facilitation by other species, which was not accounted for in our study. Vespa et al. (2018) argue that the increased dispersal of forest tree species into plantations with increasing age of the plantations is most likely due to the return of bats and birds into the plantations and surrounding forest, which facilitate the dispersal of more diverse seeds. Nevertheless, our simulations suggest that there is also selection against dispersal at ecotones independent of possible mechanistic limitations.

Conclusions
Overall, we show that the landscape structure is paramount in the evolution of dispersal, especially when considering the differing degrees of spatial variance compared to temporal variance, but also when accounting for the patterns in which habitat is distributed across the landscape. With the ongoing land use change more research into this topic could also help shed light on the difficulties species might face under the change from landscapes characterized by gradual ecotones to landscapes with sharp edged ecotones, as they particularly occur in agricultural and urban settings.
Abbreviations IBM: Individual-based model; NN: Scenarios with random nearest neighbor emigration; HE: Scenarios with nearest neighbor emigration dependent on the patch.
Additional file 1: Fig. S1. Local adaptation in the niche optimum trait (blue line) for all degrees of spatial variation and all landscapes scenarios. In each row, spatial variation increases from left to right, while global temporal variation stays the same (σ T = 1). Individuals are adapted perfectly to local conditions (red line). Fig. S2. Local adaptation in the tolerance (niche width) trait for all degrees of spatial variation and all landscapes scenarios. Each iteration of the simulation is shown in its own color, however, iteration number 1 (red) is often not visible due to plotting order. In each row, spatial variation increases from left to right, while global temporal variation stays the same (σ T = 1). Tolerance is not influenced by either landscape type, patch size or degree of spatial variation and only slightly by the respective time runs of the iterations. Fig. S3. Emerging emigration probability trait values d in each patch in an exemplary landscape (the same landscapes and panel configuration as in Fig. 1 in the main text). Top row shows the results of scenarios with random emigration (NN), second row with habitat dependent emigration (HE). The value is lowest at the ecotones, especially in tiled landscapes (panel three and four), indicating spatial sorting.