Evolution of dispersal under spatio-temporal heterogeneity

This is a PDF file of an article that has undergone enhancements after acceptance, such as the addition of a cover page and metadata, and formatting for readability, but it is not yet the definitive version of record. This version will undergo additional copyediting, typesetting and review before it is published in its final form, but we are providing this version to give early visibility of the article. Please note that, during the production process, errors may be discovered which could affect the content, and all legal disclaimers that apply to the journal pertain.

whether they can successfully invade the population of residents. Details of this approach are described in the Appendix.

Spatial and temporal heterogeneity
Let us now especially consider the case of two qualities, q ∈ {H, L}, representing patches with High 158 and Low productivity, respectively. In each patch, the transition of quality occurs independently of the other patches. The transition is governed by a time-homogeneous Markov chain, in which the transition 160 probabilities P (j ← i) (i, j ∈ {H, L}) describing the probability that the quality of a patch changes from where 0 < α, β < 1. Equation (1) can be re-written as J o u r n a l P r e -p r o o f Journal Pre-proof (see eq. (2) of Massol and Débarre (2015)) with the following transformation from (α, β)-space to (π, τ )-164 space: which maps the region of feasible parameters, {(α, β) 0 < α, β < 1}, to a pentagon-like region (Fig. S1), 166 (π, τ ) 0 < π < 1, max − π 1 − π , − 1 − π π < τ < 1 .
The new parameter π represents the equilibrium proportion of high-productivity patches, and τ represents the temporal autocorrelation of patch quality between two consecutive censuses (see Massol and Débarre 168 (2015); see also our Appendix B for details).

Model parameters 170
In Appendix A we show that the metapopulation fitness, and thus the evolution of dispersal, depends on the fecundities F k and population densities n k only through their product, F k n k . In the case of two 172 qualities, this dependence occurs through the productivity ratio f = F H n H (F L n L ). Without loss of generality we can assume f > 1, so that high productivity patch type has a larger product, F H n H > F L n L .

174
Concerning arrival-bias parameters λ k , only their relative magnitude matters (if multiplied with the same constant, the arrival probabilities ϕ k remain unchanged). In the case of two patch qualities, metapopu-176 lation fitness depends on the parameters λ H and λ L only through λ = λ H λ L , which measures arrival bias to high-productivity patches relative to the low-productivity ones. 178 An explicit expression of the metapopulation fitness of mutants can be derived (details in Appendix C and the electronic appendix). That expression uses the parameters f and λ together with other parameters, 180 listed in Table 1.
Model parameters 0 < p ≤ 1 dispersal survival probability f > 1 productivity ratio 0 < π < 1 proportion of high-productivity patches −1 < τ < 1 temporal autocorrelation λ > 0 arrival bias (λ > 1 represents bias towards high-productivity patches) Evolving strategies 0 ≤ m ≤ 1 dispersal probability of unconditional strategy 0 ≤ m H ≤ 1 dispersal probability from high-productivity patch of conditional strategy 0 ≤ m L ≤ 1 dispersal probability from low-productivity patch of conditional strategy  182 We will first study the evolution of unconditional dispersal, i.e., assume that the patch quality has no effect on the dispersal strategy of an individual, so that m 1 = m 2 = m.
which is evolutionarily attracting from all initial strategies (Fig. 1bc). For λ = 1 our eq.(7) agrees with 200 eq.(18) of Cohen and Levin (1991) and eq.(13) of Massol and Débarre (2015). Fig. 2a-d illustrates the dependence of the singular strategy on parameters depicting spatiotemporal heterogeneity, discussed in 202 more detail in Section 3.2. For comparison purposes, Fig. 2e-h show analogous results for conditional dispersal, which will be discussed in Section 4. 204 Whether the singular strategy is uninvadable (evolutionarily stable strategy, ESS) or not, can be determined based on the second derivative of metapopulation fitness, ∂ 2 ∂m 2 mut R((m mut , m mut ); (m, m)) mmut=m .

206
According to Theorem D.3, the sign of this second derivative is given by the sign of −τ . Therefore, if τ > 0 (positive autocorrelation), the evolutionarily attracting singular strategy is uninvadable, and dis-208 persal evolution will converge to the singular strategy (Fig. 1b). If, τ < 0 (negative autocorrelation), the evolutionarily attracting singular strategy is not uninvadable, and it is thus a branching point (Fig. 1c). 210 In Fig. 2a, the black thick line τ = 0 separates the parameter regions of these two cases. The branching threshold (τ = 0) corresponds to a case when the quality of a patch changes fully randomly (i.e., inde-212 pendent of the previous state). This threshold has already been found by Massol and Débarre (2015) for the case λ = 1.

214
When m * is a branching point (evolutionarily attracting and τ < 0), dispersal evolution is first expected to approach m * , but then the population becomes dimorphic, and disruptive selection will cause the 216 two strategies present in the population to evolve away from each other. In other words, evolutionary branching will happen. The expression for metapopulation fitness derived in Appendix C gives 218 the metapopulation fitness of a rare mutant in the environment set by a monomorphic resident only. Therefore, in Appendix D.7, we present results from numerical simulations illustrating that evolutionary 220 branching in this model typically leads to the coexistence of highly-dispersing individuals and almost sessile individuals. . Singular dispersal strategies of unconditional dispersal m * (left panels) and singular dispersal strategy components m * H of conditional dispersal (right panels) (a-c,e-g) with respect to autocorrelation τ and proportion of high-productivity patches π, when f = 2.5, (d,h) with respect to τ and productivity ratio f , when π = 0.5. Other parameters: p = 0.95, λ = 1. Purple curves in (b-d) correspond to branching points, black curves to evolutionarily stable strategies. The conditional dispersal strategy component m * L = 0.
J o u r n a l P r e -p r o o f Journal Pre-proof

Spatiotemporal heterogeneity promotes the evolution of dispersal
According to (5) and (6) (See also Appendix D.2), the zero-dispersal strategy m = 0 is evolutionarily The first part of this expression is non-positive and the second one is non-negative. Positive dispersal can thus evolve for intermediate π, 232 when p and f are large enough.
When patches have high productivity most of the time (π max < π ⩽ 1) or when they have low productivity 234 most of the time (0 ⩽ π < π min ), there is not much spatial heterogeneity in the model, and the zerodispersal strategy is evolutionarily attracting (Fig. 2a, light gray shading). The thresholds π min and π max 236 given by (D.6) in the Appendix do not depend on τ , and thus the boundaries separating the light-gray area of zero-dispersal and positive dispersal in Fig. 2a are vertical lines. As illustrated in Fig. 2b, the 238 singular dispersal strategy typically reaches its maximal value for intermediate π. Spatial heterogeneity, in the form of substantial proportions of both patch types, thus promotes the evolution of dispersal in 240 this model.
The productivity ratio f is another measure of spatial heterogeneity. When f ≈ 1, there is very little 242 spatial heterogeneity, as patch productivities are then similar. In such a situation we haveD 0 < 0, and positive dispersal does not evolve. Theorem D.4 in Appendix D.6 shows thatD 0 is an increasing function 244 of the productivity ratio f . Therefore, the zero-dispersal strategy is evolutionarily attracting only for 1 < f <f , and positive dispersal can evolve for f >f , in which the expression forf is given in (D.7). This 246 effect is illustrated in Figs. 2d and 3. Increasing the productivity ratio f increasesD 0 and thus promotes the emergence of dispersal.

248
The effect of f on the singular dispersal strategy m * , however, is more complicated. The following properties are proved analytically in Theorem D.4: When λ ⩽ 1, the singular dispersal strategy m * increases 250 when f increases (Fig. 3ab). The same is true, when λ > 1 and p is small enough (Fig. 3c). However, when λ > 1 and p ≈ 1, m * can decrease when f increases. We can thus conclude that spatiotemporal 252 heterogeneity, in the form of substantial differences in the patch qualities, mostly promotes the evolution of dispersal.
The autocorrelation parameter τ can be interpreted to measure temporal heterogeneity (Appendix B). If τ is close to 1, patch qualities change very seldom, so that there is very little temporal heterogeneity.

256
On the other hand, if τ is close to −1, patch qualities change almost every time, which means strong temporal heterogeneity. The quantityD 0 does not depend on τ . IfD 0 > 0, positive dispersal evolves 258 under our assumption τ < 1. Autocorrelation τ thus does not affect the emergence of positive dispersal, but has an effect on the singular strategy m * . According to Theorem D.4, m * is a decreasing function 260 of τ , and m * → 0 as τ → 1. Furthermore, there exists such a threshold value τ thresh , that when τ is small enough, τ < τ thresh , full dispersal (m * = 1) may evolve (Fig. 2c). These effects of τ on the singular 262 strategy m * were already found by Massol and Débarre (2015) for λ = 1. Theorem D.4 thus provides J o u r n a l P r e -p r o o f Journal Pre-proof a mathematical proof on this dependence for general λ. The analytical expression for τ thresh is given in 264 (D.10). In Fig. 2a, this is the thick black curve separating the areas of positive singular strategy (color shading ranging from blue to red) and full dispersal m * = 1 (dark gray shading). We can thus conclude 266 that temporal heterogeneity promotes the evolution of dispersal.
The parameter plot of the singular dispersal strategy m * with respect to the proportion of high-productivity 268 patches π and autocorrelation τ in Fig. 2a illustrates that spatiotemporal heterogeneity promotes dispersal. In Appendix D.8 we show that depending on other parameters, such parameter plots may look 270 qualitatively different. Fig. S3 shows all qualitatively different types of such parameter plots for p < 1, and Fig. S5 for p = 1. Furthermore, Fig. S4 illustrates when each of these cases occurs depending on 272 the choice of the other remaining parameters, productivity ratio f , dispersal survival probability p, and arrival bias λ.

274
In this section we have observed that spatiotemporal heterogeneity promotes the evolution of dispersal.
In particular, the evolutionarily attracting singular dispersal strategy m * reaches a maximum for inter-

Larger dispersal survival promotes dispersal
It is relatively easy to see thatD 0 increases with p when f > max{1, λ}. Therefore, increasing p promotes 288 the emergence of positive dispersal. Furthermore, according to Theorem D.4, the singular strategy m * increases with dispersal survival p (Fig. 4a), which is again a generalization of a result by Massol and  According to Theorem D.4, when λ is large, λ > f , bothD 0 < 0 and D(1) < 0 hold, which means that dispersal evolves to zero, m * = 0 from all initial conditions. Using continuity arguments, ifD 0 > 0 for 294 some values of λ < f , there exists λ * < f at whichD 0 = 0. Furthermore, the singular strategy is decreasing with respect to λ for λ < λ * at least when λ ≈ λ * . Numerical explorations illustrate (Fig 4b) that the 296 singular strategy can either be a decreasing function of λ for all λ < λ * , or a non-monotonic function of λ (Fig. 4b). 298

Evolution of conditional dispersal
Next, we will study the evolution of conditional dispersal, i.e., assume that the dispersal strategy of a

304
According to Theorem E.1, the fitness gradient for conditional dispersal is given by (see Appendix E.1) The coefficients are given by

308
The evolutionary dynamics of conditional dispersal can be understood by analyzing phase-plane plots (Fig. 5) showing the isoclines (nullclines) of the components of the fitness gradient, • If p = 1, the two isoclines (11) are identical and dispersal evolution first converges to the isocline 314 but is neutral along it (Fig. 5bc). The isoclines may reside outside the feasible domain, in which case m H evolves to 1 and m L evolves to 0 (Fig. 5d).
According to Appendix E.3, in case 0 < m * H < 1, we have  Fecundity F k measures reproductive success, and population densities n k alter the probability of surviving competition. Patches with the larger value of F k n k have high productivity. Naively one could have 324 expected that it would be better to avoid patches of low productivity, so that m L would evolve to some positive value and m H would evolve to zero. However, the local conditions individuals observe consist not 326 only of F k and n k , but also of the amount of competitors therein. We will consider the heuristic reasons behind our analytical result in Discussion.

Spatiotemporal heterogeneity and evolution of conditional dispersal
According to (9) and (12), positive dispersal from high-productivity patches, m * H > 0 evolves, if B 0 > 0, 330 or equivalently τ < 1 and Substantial spatial heterogeneity is present when π is intermediate and f is large. We observed in 336 Section 3 that such conditions favor unconditional dispersal. In contrast, the evolution of the conditional dispersal component m * H is promoted by decreasing π: Assuming f p > 1, which is necessary for positive 338 dispersal to emerge, the quantity B * 0 is a decreasing function of π. In particular, lim π→0 B * 0 = f p−1 holds. Positive dispersal for small π thus evolves, when f > 1 p.
which means that m * H = 0 for π ≈ 1 (Figs. 2ef and 6). More precisely, we have m * The motivation for an individual to emigrate from a high-productivity patch is the hope to immigrate into a low-productivity patch to avoid competition, which has a very low probability when π ≈ 1. Here 344 we observe that the existence of low-productivity patches together with their possibility to become highproductivity patches promotes the evolution of conditional dispersal.

346
In contrast with unconditional dispersal m * , the singular strategy component m * H of conditional dispersal can either be a monotonically decreasing or non-monotonic function of π. Detailed conditions for the 348 two cases are given by Theorem E.3: For τ < 0, the condition (16) is only a necessary condition for the singular strategy m * H to be a non-352 monotonic function of π (Fig 6b).
According to Theorem E.3, zero dispersal m * H = 0 evolves for low productivity ratios, 1 ⩽ f ⩽f , in which 354f is given in (E.18). When the productivity ratio is large enough, f >f , the strategy m * H is positive, and increases with f (Fig. 2h). The productivity ratio f can thus have qualitatively different effects on 356 unconditional and conditional dispersal, since the unconditional dispersal strategy m * was observed to be either increasing or non-monotonic with respect to f (Figs. 2d and 3). Spatiotemporal heterogeneity, in 358 the form of substantial differences in the patch qualities, promotes the evolution of conditional dispersal. Proportion of high-productivity patches, π Singular dispersal, m * H τ = 0.5 Figure 6: The singular conditional dispersal strategy component m * H with respect to the proportion of high-productivity patches π for (a) different values of productivity ratio f (b) different values of autocorrelation τ . When τ ⩾ 0, the singular strategy m * H is a non-monotonic function of π, when λ > 1 p and f >f (panel a, f = 3, 5, 30). Panel b) illustrates, that when τ < 0, the condition λ > 1 p and f <f is only a necessary condition for m * H to be a non-monotonic function of π. When τ is small, the value of π at which m * H would reach its maximum (marked in red) may lie outside of the parameter domain. Parameters (a) Temporal heterogeneity (small τ ) promotes evolution of unconditional dispersal as well as conditional

Larger dispersal survival promotes dispersal
The expression B 0 given by (14) is an increasing function of 370 p. Furthermore, the condition B * 0 > 0 can be written as .
According to Theorem E.3, the singular dispersal strategy component m * H is an increasing function of p 372 (Fig. 7a). (Note, however, the neutrality of dispersal evolution for p = 1, explained in Section 4.2, proved in Theorem E.2.) We can thus conclude that increasing dispersal survival p promotes also the evolution 374 of conditional dispersal.

The effect of arrival bias λ may be non-monotonic 376
In Section 4.3 we noticed that positive dispersal may evolve, m * H > 0, if λ is not too large. The expression B 0 given by (14) is a decreasing function of λ. Therefore, we can write the condition B * 0 > 0 as (1 − p)π and f p > 1.
Increasing the arrival bias λ thus hinders the emergence of dispersal. Its effect on the singular strategy m * H may, however, be non-monotonic. According to Theorem E.3, if 1

Summary
We extended previous models (Cohen and Levin, 1991; Massol and Débarre, 2015) of the evolution of 384 unconditional dispersal in which patch qualities are spatially heterogeneous and temporally fluctuating, to include a case in which arrivals can be biased towards (λ > 1) or against (λ < 1) the high-productivity 386 patches. Increasing λ thus makes competition in the high-productivity patches stronger. The extended model shows richer dynamics including evolutionary bistability and non-monotonic parameter depen-388 dence with respect to the productivity ratio f . We also studied the evolution of conditional dispersal. Surprisingly, we found that conditional dispersal from low-productivity patches never evolves. Further-390 more, evolutionary branching occurs for τ < 0 in the unconditional case, while it never occurs in the conditional case. Phenotypic plasticity can be expected to hinder evolutionary branching. Evolutionary 392 bistability can occur only for unconditional dispersal.
The effects of parameters on the evolutionarily singular dispersal strategies are summarized in Table 2. For 394 some of the parameters, their qualitative effect is similar for both unconditional and conditional dispersal. First, increasing the dispersal survival probability p increases both m * (unconditional dispersal) and m * H 396 (conditional dispersal). This result is rather intuitive, but note that increasing dispersal survival has been shown to decrease dispersal in some models (Comins et al., 1980;Gandon and Michalakis, 1999;398 Heino and Hanski, 2001). Second, increasing the temporal autocorrelation decreases both m * and m * H . Temporal heterogeneity (decreasing τ ) thus promotes dispersal. Third, the singular dispersal strategy 400 either decreases when the arrival bias λ is increased, or is non-monotonic with respect to λ. For other parameters (f and π) there are differences in their qualitative effect.

402
The productivity ratio f is one of the parameters with qualitatively different effects for unconditional and conditional dispersal: Although in both cases f needs to be large enough for dispersal to emerge, 404 the conditional dispersal strategy m * H is always increasing with respect to f , whereas the unconditional dispersal strategy can also be non-monotonic with respect to f . This can happen only when the cost of 406 J o u r n a l P r e -p r o o f Journal Pre-proof dispersal is small, p ≈ 1, and there is arrival bias towards high-productivity patches, λ > 1.
The proportion of high-productivity patches π has two qualitatively different effects on unconditional and 408 conditional dispersal. First, the unconditional dispersal strategy m * is always non-monotonic with respect to π, whereas the conditional dispersal strategy m * H is either decreasing or non-monotonic with respect 410 to π. Second, the unconditional dispersal strategy m * = 0 for π ≈ 0 or π ≈ 1, whereas if conditional dispersal m * H evolves for some π, we have m * H > 0 for π ≈ 0. The reason for the latter difference is 412 as follows. Positive unconditional dispersal can evolve, if the average benefit of dispersal from highproductivity patches and from low-productivity patches is positive. When one patch type dominates 414 the population, the selection pressure in dominant patches almost exclusively determines m * (and the rare patches matter very little). Dispersal from the dominant patch type is not beneficial, because one 416 may die during dispersal, and survivors arrive with high probability into a similar patch they left from. This is the reason why no unconditional dispersal evolves when π is close to either 0 or 1 (see Fig. 2b).

418
Such evolution to no dispersal in a homogeneous environment has been observed before (Hastings, 1983;Holt, 1985;Cohen and Levin, 1991;Parvinen, 1999;Gyllenberg et al., 2002;Parvinen, 2006). In case of 420 conditional dispersal, the mechanism described above explains also why m * H = 0 for π ≈ 1 and m * L = 0 for π ≈ 0. For conditional dispersal, each strategy component affects dispersal behavior in the corresponding 422 patch type only. Dispersal from the rare patch type may evolve, if the living conditions in dominant patches are better. When high-productivity patches are rare, π ≈ 0, individuals dispersing from high-424 productivity patches to low-productivity patches experience less competition, for which reason positive dispersal m * H > 0 can evolve. In our model, however, we always have m * L = 0, reasons for which are 426 discussed in Section 5.2.
increasing OR non-monotonic, when λ > 1 and p ≈ 1 increasing 0 Larger autocorrelation, τ decreasing decreasing 0 Larger proportion of highproductivity patches, π non-monotonic decreasing OR non-monotonic, when f > 1 p 2 and λ >λ > 1 0 Larger arrival bias λ decreasing OR non-monotonic decreasing OR non-monotonic 0 Evolutionary branching when τ thresh < τ < 0 and 0 < m * < 1 never Evolutionary bistability when λ > 1, τ < τ thresh , and f is intermediate never respectively, and our results tell us m * L = 0. For unconditional dispersal m H = m L , so both expressions simply become m * . In the conditional dispersal case, the average dispersal rates (19) depend not only 432 on the evolved dispersal rate from high-productivity patches m * H , but also on the proportion of highproductivity patches, π. Although m * H in Fig. 8a decreases with respect to π, both patch-average and 434 juvenile-average dispersal rates are non-monotonic with respect to π, reaching a maximum at an intermediate proportion of high-productivity patches (Fig. 8). In both unconditional and conditional cases, 436 average dispersal rate is large when both high-and low-productivity patches are present in substantial magnitudes and when the productivity ratio f is large. In this sense, we can say that spatiotemporal (thin curve) can be decreasing (panel a) with respect to π, but both the patch-average dispersal rate (thick blue curve) and the juvenile-average dispersal rate (thick orange curve), given by (19), are non-monotonic with respect to π. Parameters: a) λ = 1 b) λ = 2.5. Common parameters: τ = 0.1, As expected, the arrival bias λ has an effect on the evolutionarily singular strategies m * (unconditional) 440 and m * H (conditional), provided that spatiotemporal heterogeneity is present. If we would have f = 1 or τ = 1, positive dispersal would not evolve, so that arrival bias (λ ≠ 1) alone does not promote dispersal.

442
As explained above, m * and m * H can be either decreasing or non-monotonic with respect to λ. However, including the arrival bias can also change the qualitative effect of other parameters on m * and m * H . First, 444 for λ = 1, the unconditional dispersal strategy is increasing with respect to f , but for λ > 1 the effect of f can be non-monotonic. Second, for λ = 1, the conditional dispersal strategy m * H is decreasing with 446 respect to π, but for λ > 1 the effect of π can be non-monotonic (Fig. S8). Furthermore, evolutionary bistability of unconditional dispersal is possible only for λ > 1. This is illustrated by Fig S4, from which 448 we also observe that the boundaries of the parameter regions (red and purple) allowing for evolutionary bistability cross at λ = 1. Analogously, the parameter region for which full dispersal never evolves (light 450 blue) has a corner at λ = 1. Consequently, unbiased dispersal λ = 1 is a special case, as an infinitesimally small change of λ allows richer evolutionary scenarios of unconditional dispersal evolution.

Why conditional dispersal from low-productivity patches does not evolve
One of our major findings is m * L = 0, irrespective of parameter values. One may naively expect that 454 dispersing offspring from a low-productivity patch could be adaptive especially when τ ∼ 1, but this is not the case. Here we intuitively explain why.

456
Let us consider one complete life-cycle starting at the moment immediately after emigration before dispersal survival selection, immigration, and patch-quality transition. At that moment, there are three 458 different types of juveniles in the population; H-juveniles are those juveniles who decided to stay in a currently high-productivity patch, L-juveniles are those juveniles who decided to stay in a currently Note that this result holds for any α > 0 and β > 0 (or equivalently, for any τ < 1). Even if patch 488 quality rarely changes (i.e., τ ∼ 1), dispersal from high-productivity patches can evolve but dispersal from low-productity patches never evolves from (m H , m L ) ∼ (0, 0).

Comparison with previous literature & future directions
Our model and the result for the unconditional dispersal without arrival bias (λ = 1) are identical to the 492 model of "juvenile dispersal with local density regulation" by Massol and Débarre (2015). We extended the model to include arrival bias (λ ≠ 1) and found richer dynamics such as evolutionary bistability, as 494 summarized in Section 5.1.
We also studied the evolution of conditional dispersal and found a novel and simple result that dispersal 496 from low-productivity patches never evolves. As Massol and Débarre (2015) have shown for the case of unconditional dispersal, our results for the conditional dispersal may also depend on the life cycle 498 assumptions (see also Johst and Brandl (1997)).
For a case of evolution of conditional dispersal with spatial heterogeneity but without temporal fluctua-500 tion, McPeek and Holt (1992) studied a two-patch model and found the selectively neutral combination of (m H , m L ). An analogous result was found by Parvinen (1999), Fig. 3b therein. These results corre-502 spond to our neutral isocline with τ = 1 and p = 1 (Fig.5b).
McPeek and Holt (1992) also studied cases with spatiotemporal fluctuation, but the strategy they investigated was assumed to be conditional on the 504 J o u r n a l P r e -p r o o f Journal Pre-proof patch index, not on the patch quality. Thus our results with spatiotemporal heterogeneity have no direct correspondence to their results. 506 We introduced arrival bias λ in immigration as a model parameter. The differences in arrival bias could be caused by active decisions by dispersers when they encounter potential patches to immigrate into.

508
Alternatively, patches with different qualities could possess different vegetation, so that the attachment probabilities of seeds dispersing by wind would be different. Although one might expect arrival bias 510 toward high-productivity patches (λ > 1) to be advantageous, there is typically less competition in the low-productivity patches. When λ < 1, dispersers may thus experience less competition. Alternatively, types of patches for settlement. In our current model, we have considered only patch-quality-dependent emigration strategy, but considering patch-quality-dependent immigration strategy and its coevolution 516 with emigration strategy (i.e., Weigang (2017)) will further enrich our understanding of evolution of dispersal traits.

518
If we assume the local population sizes to be finite, kin selection will be present, which is expected to promote dispersal (e.g. Hamilton and May, 1977). Concerning conditional dispersal, one would thus 520 expect positive dispersal to evolve from both patch types. That is, our prediction m L = 0 might not hold for models with finite local population sizes. However, depending on parameter values, dispersal from 522 one patch type may evolve to zero, while dispersal from the other patch type evolves to a positive value, analogous to our result in Section 4. 524 We have studied a particular life-cycle assumption in our model, but we could consider many of its variants. As Massol and Débarre (2015) stressed, the order of events in life-cycle assumptions can 526 dramatically change the evolutionary outcomes, and our prediction here may not be robust against changing those details.

528
For mathematical simplicity, we have assumed in our model that there are only two different patch qualities and that each patch follows the same transition rule between those quality states. However, for a

536
In summary, we have studied evolution of dispersal in a heterogeneous environment with arrival bias. We have found that arrival bias introduces richer evolutionary dynamics including bistability. We have 538 also found that evolutionary outcomes are strikingly different between unconditional and conditional dispersal. We believe that those results shed light on ecological factors that impact dispersal evolution.  In the initial phase of potential invasion by mutants, practically all patches are occupied by residents only. In patches of type k, the adult population density is n k , and adults reproduce with fecundity γF k .

664
Juveniles disperse with probability m k . Since the proportion of patches of type k is π k , the average amount of resident emigrants is

A.2 Resident immigrants
Emigrants will survive dispersal with probability p. The parameters λ k describe the relative attractivity 668 of patches of quality k, so that emigrants arrive into a patch of type k with probability ( If all patches are equally attractive (λ k = 1), we obtain from (A.2) thatλ k = 1 for all k, and the probability ϕ k to arrive into a patch of quality k is equal to the proportion π k of such patches. In such case each 676 patch will receive the same amount of immigrants, γI res , independent of its quality.
According to (A.3), each patch of type k receives γλ k I res immigrant residents, and there are (1−m k )γF k n k 678 natal residents, so the total population density in patches of type k after immigration is (1 − m k )γF k n k + γλ k I res (A.4)

680
After immigration, each patch experiences the possibility of transition of patch type. Consider a patch that was of type k at time t, and is of type j at time t + 1. Let 1 γ S j←k denote the settlement probability 682 in such patches. Since n j is the population density after competition, we have S j←k = n j (1 − m k )F k n k +λ k I res (A.5)

Journal Pre-proof
The settlement probability thus depends both on the current and previous quality of the patch.

A.4 Metapopulation fitness
In the moment after immigration, but before transition, consider a mutant juvenile present in a patch 686 of type k. This mutant and all its descendants remaining in this patch will be called the mutant colony. Over the years, this mutant colony will send emigrants. Let R k denote the total amount of emigrants 688 sent by the mutant colony founded by the focal mutant juvenile. The patch type after transition will be j with probability P (j ← k). In such case, the settlement probability for the mutant juvenile is 1 γ S j←k .

690
It will reproduce with fecundity γF j . The proportion m mut,j of the offspring will disperse and survive dispersal with probability p. On the other hand, the proportion 1 − m mut,j of the offspring will remain in 692 this patch, and will eventually produce R j emigrants. We thus obtain Note that the scaling factors γ in this equation cancel each other. These relations form a system of linear 694 equations, from which R 1 , . . . , R N can be solved, at least numerically. A dispersing juvenile will arrive at a patch of type j with probability ϕ j , and the metapopulation fitness is From the last expression of (A.6) we observe that the relative fecundity and local population density always appear as a product, F j n j or F k n k . Therefore, these parameters affect dispersal evolution only 698 through their product.

A.5 Metapopulation fitness from adult perspective 700
From a viewpoint of an adult, a focal adult mutant in a patch of type k produces F k juveniles, among which the fraction 1 − m mut,k stays in the natal patch whose type might change from k to j. Thus, the 702 number of emigrants produced by the focal adult and its descendants is given bỹ To combine them to obtain metapopulation fitness, we again calculate the number of emigrants produced 704 by a mutant disperser that has just landed in a patch of type j. The patch type might change from j to k, and the juvenile settles there to become an adult in a patch of type k. Such an adult will have a 706 metapopulation fitnessR k . Thus we have Actually the two methods yield the identical expression of R, irrespective of the explicit forms of P (j ← k) 708 and S j←k .

J o u r n a l P r e -p r o o f
Journal Pre-proof B Spatial and temporal heterogeneity for two patch qualities 710 Let us now especially consider the case of two qualities, q ∈ {1, 2}. In the main text, label 1 corresponds to High productivity, and label 2 corresponds to Low productivity, but the following calculations hold 712 for any two patch-type situations. The transition is governed by a time-homogeneous Markov chain, in which the transition probabilities P (j ← i) (i, j ∈ {1, 2}) describing the probability that the quality of a 714 patch changes from i to j are where 0 < α, β < 1. Let the random variable X(t) denote the quality of a focal patch at time t, and 716 p i (t) = P (X(t) = i) the probability that the quality of the focal patch is i at time t. These probabilities satisfy Since the eigenvalues of the transition matrix are 1 and 1 − α − β, we obtain Therefore, after some transient, the probability distribution of quality of each patch becomes independent of its initial condition and of that of other patches as well. Consequently, 722 the probabilities that a given patch at a given census time is in quality 1 and 2 respectively converges to π 1 and π 2 , where they satisfy 724 π 1 π 2 = P (1 ← 1) P (1 ← 2) P (2 ← 1) P (2 ← 2) with π 1 + π 2 = 1. The solution is which we call π in the main text.

726
The meaning of parameter τ in eq.(3) in the main text can be explained by following the argument by Rodrigues and Gardner (2012), as follows. Let us introduce a dummy variable S t , which is equal to 1 if the 728 patch quality at t-th census time is 1 and is equal to 0 if the quality is 2. Then, the correlation coefficient between S t and S t+1 of the same patch is deemed as a quantity that measures how similar/different the 730 patch quality is between two adjacent census times, which takes a value between −1 (perfect negative correlation) and 1 (perfect positive correlation). Let τ denote this correlation coefficient. By definition Each term is evaluated as J o u r n a l P r e -p r o o f Journal Pre-proof where we used (B.5). This expression makes intuitive sense because larger α and β mean that the patch quality frequently changes and hence that the stability of patch quality τ is low. By solving Eqs.(B.5, 738 B.9) with respect to α and β, we obtain Therefore, our model is parameterized either by (α, β) or by (π 1 , τ ). In the latter case, its domain is a 740 curved pentagon-like region represented as as illustrated in Fig. S1.  Figure S1: Dependence of τ and π on α and β The limit of no temporal heterogeneity, τ → 1, can be obtained by letting the transition probabilities α 744 and β tend to zero, while keeping their relative magnitudes constant: α = ϵπ and β = ϵ(1 − π) and ϵ → 0. Alternatively, we can assume that the patch qualities do not change, α = β = 0, and assume that the 746 patch proportions are π and 1 − π. Both approaches result in the same expression of metapopulation fitness, which means that our results for τ → 1 are the same as for τ = 1.

748
The limit τ → −1, can be obtained by letting α → 1 and β → 1, so that π → 1 2. That limit, however, is not consistent with starting with α = β = 1. A Markov chain with such a transition matrix has the 750 equilibrium (1 2, 1 2), but it is not asymptotically stable. From any other initial condition (x, 1 − x), the Markov chain fluctuates between the states (x, 1 − x) and (1 − x, x). In other words, it is not meaningful

C.1 Arrival probabilities
Concerning the relative attractivity parameters λ k , only their relative magnitude matters. According to 756 (A.2), we haveλ .
In the case of two patch types, metapopulation fitness thus depends on the parameters λ 1 and λ 2 760 only through λ = λ1 λ2 , which measures the attractivity of high-productivity patches relative to the lowproductivity patches.

C.3 Metapopulation fitness
The pair of linear equations (C.5) can be solved explicitly. Substituting the solution into (A.7), and

D Evolution of unconditional dispersal
We will first study the evolution of unconditional dispersal, i.e, assume that the patch quality has no 776 effect on the dispersal strategy of an individual, so that m 1 = m 2 = m and m mut,1 = m mut,2 = m mut .

D.1 Fitness gradient 778
The fitness gradient is the first derivative of the metapopulation fitness (C.6) with respect to the mutant dispersal strategy, evaluated when the mutant dispersal strategy is equal to that of the resident.
Since the coefficientÃ is positive, the direction of selection is completely determined by Z + mY , which 784 is linear with respect to m. Therefore, there exists at most one singular strategy m * = −Z Y , provided that 0 < m * < 1. Before investigating the properties of the singular strategy, consider the attractivity of 786 the boundaries of the strategy space.

788
The fitness gradient D(0) is not defined, but we can determine the direction of selection at zero dispersal by calculating the limit of (D.2) when m → 0, The denominator of (D.4) is clearly positive. The sign is thus determined by the numeratorD 0 given by (6): The first part ofD 0 is either zero, when p = 1, and negative otherwise. The second part ofD 0 has the multiplier π(1 − π), which is zero when there is no spatial heterogeneity, π = 0 or π = 1. This multiplier 794 reaches its maximum for π = 1 2 . The sign of the second term is determined by f − λ. If f > λ, the second term is positive for intermediate π, and zero dispersal may be repelling. If 1 < f < λ, zero dispersal is  productivity patches π and autocorrelation τ . As expected, dispersal does not evolve, when π ≈ 0 or π ≈ 1 (areas with light-gray shading). Dispersal evolves to a positive singular strategy for intermediate 800 values of π (color shading ranging from blue to red). The curves separating these two areas are straight vertical lines, because the expressionD 0 does not depend on autocorrelation τ . By solvingD 0 = 0 for π 802 we obtain π min,max = 1 2 By solvingD 0 = 0 for f we obtain
Consequently, potential evolutionary endpoints are Proof. The second derivative at singularity is where the components C, B 1 and B 2 are listed below.
The componentD 0 is the numerator of the scaled selection gradient for zero dispersal, given in (D.4).

838
It is positive, when the zero-dispersal strategy is evolutionarily repelling, which is a necessary condition for the singular strategy to be evolutionarily attracting. Furthermore, f > λ is a necessary condition for 840 the zero-dispersal strategy is evolutionarily repelling. Therefore, assuming that the singular strategy is evolutionarily attracting, we have (f − λ)D 0 > 0.

842
The component J o u r n a l P r e -p r o o f Journal Pre-proof does not affect the sign of (D.12).

844
The component is clearly non-negative. It is zero, when p = 0, π = 0, π = 1, or τ = 1, which all are such special cases that 846 dispersal does not evolve. Finally, Since f > λ when the singular strategy is evolutionarily attracting, the expression B 2 is certainly positive, when τ ⩾ 0. The expression is increasing with respect to τ . Therefore it is enough to prove that B 2 > 0 850 when τ = τ thresh defined in (D.10). We obtain which is positive, when the zero-dispersal strategy is repelling. • The singular strategy decreases with autocorrelation τ . (Fig. 2c). 856 • The singular strategy increases, when dispersal survival p is increased. (Fig. 4a) • Increasing the relative fecundity f increasesD 0 and thus promotes the emergence of dispersal. The 858 singular strategy m * increases with f when λ ⩽ 1 (Fig. 3bc). For λ > 1 the singular strategy m * mostly increases when f increases, but m * decreases when f increases when p is sufficiently close 860 to 1 (Fig. 3a).

862
Numerical explorations illustrate that the singular strategy can be decreasing or non-monotonic with respect to λ (Fig. 4b). 864 Proof. In order for dispersal to evolve, we must haveD 0 > 0, for which a necessary condition is f > λ. AssumingD 0 > 0 and Y < 0 we have the following effects of parameters on the expression − Z Y :

866
• Autocorrelation τ : When τ = 1, we have Z = 0 and Y < 0, so that dispersal evolves to zero. It is easy to see that J o u r n a l P r e -p r o o f Journal Pre-proof which means that when 0 < m * < 1, the singular strategy decreases with autocorrelation τ .
• Dispersal survival p: which means that increasing p promotes the emergence of dispersal. Furthermore, The component determining the sign of (D.21) thus decreases when τ is decreased. Substituting the lower bound τ = τ thresh from (D.10) results in which is positive, when the singular strategy is evolutionarily attracting. Therefore, the singular strategy m * increases with p.

878
• Productivity ratio f : With the help ofD 0 > 0 we obtain which means that increasing f promotes the emergence of dispersal.
The component determining the sign of (D.25) thus decreases when τ is decreased. Substituting 882 the lower bound τ = τ thresh from (D.10) results iñ in whichW 0 does not depend on p and ThereforeW 0 − pW 1 decreases when p is increased. Substituting p = 1 results iñ J o u r n a l P r e -p r o o f

Journal Pre-proof
We thus reach the following conclusion: If λ ⩽ 1 thenW 0 −W 1 ⩾ 0, so that the singular dispersal 886 strategy increases when f increases. If λ > 1 then the singular dispersal strategy usually increases when f increases, but can decrease, when p ≈ 1. In particular, for p = 1 we have which is positive for λ < 1 and negative for λ > 1. 890 In the special caseW Therefore, the singular strategy m * increases with f at least when λ = 1.  As discussed in subsection 3.2, the parameter plot of the singular dispersal strategy m * with respect to the 900 proportion of high-productivity patches π and autocorrelation τ in Fig. 2a illustrates that spatiotemporal heterogeneity promotes dispersal. However, depending on other parameters, such parameter plots may 902 look qualitatively different. Fig. S3 shows all qualitatively different types of such parameter plots for p < 1. Below we will explain the details of all such cases ranging from a-g. Furthermore, Fig. S4 illustrates when 904 each of these cases occurs depending on the choice of the other remaining parameters, productivity ratio f , dispersal survival probability p, and arrival bias λ.  Figure S3: Qualitatively different parameter plots with respect to autocorrelation τ and proportion of high-productivity patches π, when p < 1. a) Dispersal evolves to zero (m * = 0) for all combinations of (τ, π), b) For large τ , dispersal evolves to zero (m * = 0), for τ close to -1 there is bistability: depending on the initial strategy, dispersal evolves either to m * = 0 or m * = 1, f-g) The zero-dispersal strategy is attracting only when π ≈ 0 or π ≈ 1. For intermediate π, a singular strategy 0 < m * < 1 exists, when τ is large enough. For τ close to -1 full dispersal evolves (m * = 1), e) As in f), but full dispersal does not evolve for τ close to -1 c-d) As in f-g), but with a region of bistability as in b).
In c and f the boundary between m * = 1 is below τ = 0, while in d and g, a part of it is above τ = 0.
a) Dispersal evolves to zero (m * = 0) for all feasible combinations of π and τ (Fig. S3a). This occurs at least when the productivity ratio is smaller than the arrival bias, f < λ (grey area in Fig. S4), 908 and the parameter range becomes wider when the dispersal survival probability is decreased. b) Bistability: dispersal evolves either to zero or to full dispersal (m * = 0 or m * = 1) in Fig. S3b. The 910 zero-dispersal strategy m * = 0 is at least locally evolutionarily attracting for all feasible combinations of π and τ , and full-dispersal strategy m * = 1 is locally evolutionarily attracting when τ is small 912 enough, τ < τ thresh , in which the analytical expression for τ thresh is given in (D.10). This occurs (purple area in Fig. S4), when λ > 1 and f is increased enough from the case a. More precisely, by 914 substituting π = 1 2 into D(1) given by (5) and (6) and taking the limit τ → −1, we obtain that D(1) > 0 in the bottom corner of Fig. S3, when On both sides of this region, a region of bistability exists. For π ≈ 0 and π ≈ 1 only zero-dispersal evolves (Fig. S3cd), red area in Fig. S4. Necessarily λ > 1. Figures S3c and d differ in that respect 922 that in Fig. S3c τ thresh < 0 for all π, whereas in Fig. S3d τ thresh > 0 for some π. The parameter regions in which cases c and d occur are plotted in red in Fig. S4, with a dashed curve separating 924 the cases c and d. By substituting τ = 0 into D(1) and finding its maximal value with respect to π, we observe that case d occurs when e) Positive dispersal evolves for intermediate π, but full dispersal never evolves (m * < 1), Fig. S3e. This case occurs, when λ < 1 and f is increased from the case a, shown as a blue region in Fig. S4.

928
By solving when the discriminant in (D.6) is positive, we obtain that positive singular strategies exist for intermediate π, providing analytically the curve between the grey (a) and blue (e) areas (λ < 1) and purple (b) and red (c) areas. in which cases f and g occur are plotted in green in Fig. S4, with a dashed curve separating the cases f and g. The curve separating the red and green areas is obtained by solving numerically 938 when τ thresh meets the boundary of the feasible parameter region (B.11) at π min or at π max given by (D.6). According to our numerical explorations, these two conditions occur at the same time. λ > 1, λ = 1 and λ < 1, the threshold τ thresh satisfies τ thresh > 0 (Fig. S5b), τ thresh = 0 (Fig. S5c), and 944 τ thresh < 0 (Fig. S5d), respectively.  Figure S5: Qualitatively different parameter plots with respect to autocorrelation τ and proportion of high-productivity patches π when p = 1. a) Dispersal evolves to zero (m * = 0) for all combinations of (τ, π), b-d) A singular strategy 0 < m * < 1 exists, when τ is large enough. For τ close to -1 full dispersal evolves (m * = 1), e) As in d), but full dispersal does not evolve. Panels b-d differ in respect of whether evolutionary branching is possible or not. Panel f is the same as Fig. S4f, but here the dots correspond to parameter combinations in panels a-e.
Earlier we obtained the fitness gradient for unconditional dispersal (D.2) by differentiating (C.6) when 954 m 1 = m 2 = m and m mut,1 = m mut,2 = m mut . It can also be obtained from the fitness gradient for conditional dispersal (E.3), by adding its components together and setting m 1 = m 2 = m. Therefore, the connection 956 between (6) and (E.3) is the following: (E.4)

958
According to Theorem E.1, the fitness gradient isoclines are straight lines (eq. (11) in the main text).
Next we show that all qualitatively different types of phase-plane plots of conditional dispersal evolution 960 are as shown in Fig. 5. Consequently, there will be no dispersal from less productive patches, i.e., m L evolves to zero.  where Proof. All qualitatively different types of phase-plane plots are shown in Fig. 5, illustrating that condi-966 tional dispersal converges to {m * H , 0}, where m * H is given in (E.6). Next we go through all cases.
• When p = 1, we obtain from (E.3) In this case, the coefficients of the isoclines are which means that the isoclines (E.5) are equal. According to (E.2) and because the coefficients 970 (E.7) are non-negative, it is clear that strategies converge towards the isocline. At the isocline, both components of the fitness gradient are zero, so dispersal evolution is neutral along it (Fig. 5bc).

972
• Consider the remaining case p < 1. To complete the proof concerning convergence, we show that the isocline for D H = 0 lies above the isocline for D L = 0 which means that singular strategies in the 974 interior of the strategy space do not exist. Together with the fact that the isocline for D L = 0 lies below the origin, the result then follows (Fig. 5defg). (E.10) The expression (E.10) is clearly positive, if 0 ⩽ τ < 1 and p < 1. However, the expression f pτ (marked 982 with *) is negative, when τ < 0. Consider therefore the case τ < 0 with p < 1 in more detail. In such case, we use −τ 1−τ ⩽ π ⩽ 1 1−τ (recall the stationary fraction of high-productivity patches satisfies 984 π = α α+β = α 1−τ = 1−τ −β 1−τ ). From π ⩽ 1 1−τ we obtain (1 − τ )(1 − π) ⩾ −τ . Furthermore, π(1 − π) reaches J o u r n a l P r e -p r o o f Journal Pre-proof its maximum at π = 1 2 , and its minimum at π = −τ 1−τ or π = 1 1−τ . Therefore so that (1 − π)π(1 − τ ) 2 ⩾ −τ . As a result, the expression in square brackets satisfies We have shown that the isocline for D Equality in (E.13) holds only in the special case p = 1 investigated above (Fig. 5b). When p < 1, (E.13) holds with strict inequality. The isoclines cross only at the origin, and conditional dispersal evolves to 996 zero (Fig. 5a), so that also (E.6) holds in this case. The neutral contour line of a pairwise invasibility plot (PIP) is thus a vertical line. This is the boundary case between strict ESS and branching, meaning that evolutionary branching does not happen. Analogous to the evolution of unconditional dispersal, the parameter plots of the singular dispersal strategy component m * H with respect to the proportion of high-productivity patches π and autocorrelation 1060 τ (as in Fig. 2e) may have qualitatively different forms depending on the other parameters. Fig. S7 shows all qualitatively different types of such parameter plots. Below we will explain the details of all such 1062 cases ranging from a-f. Furthermore, Fig. S8 illustrates when each of these cases occurs depending on the choice of the other remaining parameters, productivity ratio f , dispersal survival probability p, and 1064 arrival bias λ.  Figure S7: Qualitatively different parameter plots with respect to autocorrelation τ and proportion of high-productivity patches π, when p < 1. a) Dispersal evolves to zero (m * H = 0) for all combinations of (τ, π) b-f) Positive dispersal evolves for 0 < π <π. In b, d and e, m * H is non-monotonic with respect to π (at least for τ > 0), whereas in c and f, m * H is decreasing with respect to π. In b and c full dispersal does not evolve, whereas it evolves in d-f. Panel d differs from e and f in the positioning of the dark grey region corresponding to the parameter combinations of (π, τ ) for which full dispersal evolves. Parameters: p = 0.8. a) Dispersal evolves to zero (m * H = 0) for all feasible combinations of π and τ (Fig. S7a). Based on 1066 results presented in Section 4.3, this scenario occurs when f < 1 p (grey area in Fig. S8). Again, this parameter range becomes wider when the dispersal survival probability is decreased, but it is 1068 considerably smaller than for unconditional dispersal (Fig. S4).

E.3 No evolutionary branching for conditional dispersal
b) Positive dispersal can evolve, 0 ⩽ m * H < 1 and m * H is a non-monotonic function of π (Fig. S7b): 1070 Dispersal evolves to zero, m * H = 0 forπ ⩽ π ⩽ 1, in whichπ is given by (15). Positive dispersal evolves for 0 < π <π, but full dispersal does not evolve. This scenario occurs, when f > 1 p, but f will disperse (which we call "L-dispersing-juveniles"). However, an H-and L-dispersing-juvenile leave the 1136 same number of progeny in the future because survival probability p is independent of juvenile's place of origin, and hence we call both of them simply "dispersing juveniles". Let v H , v L , v D be reproductive values 1138 of a single H-juvenile, of a single L-juvenile, and of a single dispersing juvenile, respectively. Since the fraction 1 − m H of H-born-juveniles become H-juveniles and the fraction m H of them become dispersing 1140 juveniles, and since the fraction 1 − m L of L-born-juveniles become L-juveniles and the fraction m L of them become dispersing juveniles, from the preservation law of reproductive values it follows that (F.1) Next, consider the moment immediately after step "3. Emigration" but before step "4. Immigration" in Section 2.1. There are those juvenile who survived mortality in dispersal, and we call them "surviving 1144 dispersing juveniles", the reproductive value of each of which is denoted by v S . From the preservation law of reproductive values, it follows that Next, consider the moment immediately after step "4. Immigration" but before step "5. Transition" in Section 2.1. There are only two types of those juvenile. They are those juveniles who are in a currently 1148 high-productivity patch (which we again call "H-juveniles") and those juveniles who are in a currently low-productivity patch (which we again call "L-juveniles"). When we consider a "surviving dispersing 1150 juvenile", the chance that he arrives at a currently high-productivity patch to become an H-juvenile is λ H π {λ H π + λ L (1 − π)} = λπ {λπ + (1 − π)}, and the chance that he arrives at a currently low-productivity 1152 patch to become an L-juvenile is (1 − π) {λπ + (1 − π)}, so from the preservation law of reproductive values we have We next consider the moment immediately after step "5. Transition" but before step "6. Competition" in Section 2.1. After a transition of patch productivity, there are four types of juveniles in the population.

1156
They are, those juveniles who are in a patch whose productivity was high before transition and is high after transition (which we call "HH-juveniles"), those juveniles who are in a patch whose productivity 1158 was high before transition and is low after transition (which we call "HL-juveniles"), those juveniles who are in a patch whose productivity was low before transition and is high after transition (which we call 1160 "LH-juveniles"), and those juveniles who are in a patch whose productivity was low before transition and is low after transition (which we call "LL-juveniles"). Reproductive values of each of those juveniles, 1162 denoted by v HH , v HL , v LH and v LL , satisfy Finally, we consider the moment immediately after step "2. Reproduction" but before step "3. Emigration" 1164 once again and consider the relationship between v HH , v HL , v LH , v LL and v HB , v LB to complete a series of recursions on reproductive values. For that purpose we need to count how many H-born-juveniles a 1166 single HH/LH-juvenile bears on average, and count how many L-born-juveniles a single HL/LL-juvenile bears on average. Those quantities have already essentially been calculated in Section C, and by using 1168