We constructed an individual-based model to simulate a metapopulation of two multi-locus, haploid species that occupy discrete patches located on a 2-dimensional grid of linear length L and toroidal shape (a doughnut) to avoid edge effects. Each patch on the grid is characterized by an environmental value Ei,i=1,…,L2, which are random variables distributed by the multivariate normal distribution:

where ET=(E1,…,EL2) and μT=(μ1,…,μL2) is a vector whose elements are the expected values of the environmental values, i.e., E(E_i) = μ_i. Here, Σ is the covariance matrix whose elements are Σij=σi2σj2ρij where σi2 is the variance of the environmental value at patch i and ρ_ij is the correlation between the environmental values at patches i and j, which we choose to depend on the Euclidian distance d_ij between those patches. Explicitly, we set ρij=ρdij where ρ∈[0, 1] is the correlation between the environmental values of patches for which d_ij = 1. Of course, d_ij = 1 is the smallest distance between any two patches in the grid. We note that the correlation decreases exponentially with the distance between patches, i.e., ρ_ij = exp(−d_ij/ξ) where ξ = 1/|ln ρ| is the correlation length of the environment.

A word is in order about the calculation of the Euclidean distance d between two points (i₁, i₂) and (j₁, j₂) in a rectangular grid with cyclic boundary conditions (toroid). Let us assume that the open grid is L₁×L₂, i.e., that there are L₁ patches in the horizontal direction and L₂ in the vertical direction (L₁ = L₂ = L for the grid considered in this paper), so that i₁, j₁ = 1, …, L₁ and i₂, j₂ = 1, …, L₂. The horizontal and vertical distances between these points are given by the equations d_h = min (|j₁−i₁|, L₁−|j₁−i₁|) and d_v = min (|j₂−i₂|, L₂−|j₂−i₂|), from where we can readily calculate the Euclidean distance, viz., d=dh2+dv2 .

To avoid a profusion of parameters we assume that the patches are statistically identical, i.e., μ_i = μ and σi2=σe2 for i = 1, …, L². With this assumption we can set μ = 0 without loss of generality, since a different choice of μ would amount to a uniform shift on the environmental values and so it would be inconsequential. Although we set σe2=2 in most of our simulations, we have also analyzed the effect of σe2 on species' coexistence. We note that for ρ = 0 the environmental values E_i are statistically independent normal random variables with mean zero and variance σe2. Most of our analysis will focus on the uncorrelated environment, but we have also analyzed the possibility of coexistence in environmentally autocorrelated landscapes (i.e., ρ>0). As a brief technical note, we mention that in the case the matrix Σ is symmetric and positive definite we can readily produce samples of the random vector E by setting E = μ+Σ^½X where XT=(X1,…,XL2) is a random vector whose components are statistically independent standard normal random variables (Wasserman, 2004). The main difficulty here is the calculation of the square root of the matrix Σ, which can be done using its spectral decomposition.

Following Scheiner et al. (2020), the phenotype Z_i of an individual located at patch i at the time of development was determined by 40 haploid loci as

where R_k are the allelic values at the m_r = 20 non-plastic or rigid loci (i.e., loci whose phenotypic expression does not ontogenetically react to the environmental value), P_k are the allelic values at the m_p = 20 plastic loci (their phenotypic expression depends on external environmental cues that influence development) and ϵ is a normally distributed environmental effect with mean 0 and variance σ_ϵ = 1/10. Here, b is the plasticity parameter that takes on the value b = 0 for the non-plastic species (species 1) and b = 1 for the plastic species (species 2). There is no lack of generality in this choice because the allelic values P_k (and R_k as well) are real-valued variables and so any other choice of the plasticity parameter can be reset to b = 1 by a proper rescaling of P_k.

The initial allelic values for all loci were also independently drawn from a normal distribution with mean 0 and variance 1/10. Hence, the sum of allelic effects for each set of loci is a normal random variable of mean 0 and variance 20/10 = 2, which matches our typical choice for the variance of the environmental values, viz., σe2=2. For a given genotype, the phenotype Z_i at patch i is a linear function of the environment value E_i and so ∑k=1mrRk is the intercept and b∑k=1mpPk is the slope (Scheiner, 2013; Scheiner et al., 2020). In the initial setup, the expected values of these quantities are zero.

Selection is only for viability, and the survival probability of an individual at patch i depends on its phenotype and the cost of plasticity. Here, we assume a Gaussian fitness model (Scheiner et al., 2020):

where E_i is the optimum phenotype that coincides with the patch's environmental value (i.e., stabilizing selection with a moving optimum), w² is inversely proportional to the strength of stabilizing selection and c≥0 determines the cost of plasticity. Here, we set w² = 1 without loss of generality. In fact, we can easily eliminate the parameter w² from the model by rescaling the adaptive non-plastic allele values Rk′=Rk/w and the patches environmental values Ei′=Ei/w. The adaptive plastic allele values P_k do not change. Since Ei~N(0,σe2), we have Ei′~N(0,σe2/w2). Hence, the effect of w² is simply a rescaling of the variance of the patches environmental values σe2. In other words, increasing the strength of selection (i.e., decreasing w²) is equivalent to increasing σe2 and hence to increasing the roughness of the landscape. Equation (3) includes maintenance costs of plasticity (DeWitt et al., 1998) because there is a proportional reduction in survival when c>0 even if plasticity is not expressed as it is the case for species 1. However, in this paper we assume that individuals belonging to species 1 do not carry plastic alleles, so effectively b = 0 and c = 0 for them.

We note that there are no intra and interspecific interactions during the viability selection process, whose net effect is to decrease the population of both species. After passing the viability selection process, the surviving individuals compete among themselves to repopulate their patch, as described next.

We assume that when the two species occupy the same patch they interact and there is interference or scramble competition. Given the species abundances N_1i and N_2i in patch i = 1, …, L² after viability selection, the total number of offspring N1i′ and N2i′ produced by the individuals of each species is determined by Ricker's equations (Ricker, 1954):

where e^r is the maximum growth rate in a low-density population, K_max is the carrying capacity of each species when alone (assuming a₁₁ = a₂₂ = 1), and a_ij is the per capita effect of species j on species i (Godfray et al., 1991). For simplicity, here we assume that both species have the same maximum growth rate and equilibrium population size when alone in a patch.

In the absence of migration, so the individuals are confined to their birth patches, we consider the scenario where the non-plastic species (i.e., species 1) outcompetes the plastic species (i.e., species 2). This scenario is achieved by setting r∈(0, 2), a₂₁>a₁₁ and a₁₂<a₂₂. Furthermore, since we do not want to distinguish between the two species when only one species is present in the metapopulation, we set a₁₁ = a₂₂ = 1. To investigate the possibility of species coexistence (or lack thereof) when the ecologically inferior species 2 displays plasticity we set a₂₁ = 3/2, a₁₂ = 1/2 and r = 0.6 in our simulations. Assuming a relatively large r is reasonable in the case of expanding species as, e.g., insects (Frazier et al., 2006) and plants (Appendix in Franco and Silvertown, 2004). We emphasize that our choice for the competition matrix a excludes the possibility of mutual invasion, which is a standard requisite for species coexistence (Pásztor et al., 2006). In fact, since a₂₁ = 3/2>1 = a₁₁ species 1 can invade a resident population of individuals of species 2 at equilibrium, but since a₁₂ = 1/2 < 1 = a₂₂ species 2 cannot invade a resident population of individuals of species 1 at equilibrium. It is instructive to note that det(a) = 1/4>0, so our competition matrix offers a counterexample to the fallacious statement that det(a)>0 implies negative frequency dependence (i.e. rare advantage) and hence ensures mutual invasibility (TBox 9.1 in Pásztor et al., 2006).

We note that Ricker equations (4a) and (4b) yield real values for the number of offspring of each species in patch i = 1, …, L² and, in fact, the conditions that guarantee the superiority of species 1 over species 2 presented above for the single-patch situation are valid only when species numbers or abundances N_1i and N_2i are real variables. It turns out that transforming those real variables into integer variables that are necessary for our individual-based simulations may produce spurious results, such as permanence of a few individuals of species 2 in patches dominated by individuals of species 1 or the stability of patches dominated by individuals of species 2 against the invasion of a few individuals of species 1. Here, we circumvent this difficulty by taking the ceiling function in Equation (4a) (i.e., the least integer greater than or equal to N1i′) and the floor function in Equation (4b) (i.e., the greatest integer less than or equal to N2i′), which biases the competition in favor of species 1. In fact, this procedure impairs considerably the ability of species 2 to colonize a vacant patch i even when there is no competition (i.e., N_1i = 0) because a single founder of species 2 cannot produce more than one offspring for our choice of growth rate (r = 0.6). For instance, if there is a single adult individual of species 2 in an otherwise empty patch i (i.e., N_2i = 1 and N_1i = 0), then Equation (4b) yields N2i′<er≈1.8. Since floor(1.8) = 1, a single founder of species 2 cannot populate a vacant patch. This effect is mitigated when the migration rate is large since in this case there is a good chance that several individuals of species 2 migrate together to the same patch. However, this is actually a convenient scenario for our purposes since the more ecologically impaired species 2 is, the more remarkable the finding that plasticity can guarantee its permanence in the metapopulation.

In addition, we set K_max as a hard upper bound to the number of offspring N1i′ and N2i′ in patch i. In other words, whenever Nli′>Kmax we set Nli′=Kmax for l = 1, 2. This procedure is actually inconsequential because the populations at each patch approach the support capacity K_max from below because of our choice of the growth rate (viz., r = 0.6). For Ricker's growth equation, overshooting and the possibility of limit cycles happens for r>2 only (see, e.g., Franco and Fontanari, 2017). In the Supplementary material, we present several instances of the time evolution of the species abundances that support our claim that N_li<K_max for l = 1, 2.

The ecological competition procedure described above determines the number of offspring of each species N1i′ and N2i′ in patch i = 1, …, L². We re-emphasize that although Equations (4a) and (4b) produce real values for the species abundances, we take the integer values of those abundances using the floor and ceiling functions with the care to bias the competition in favor of species 1 (see Section 2.3). Now we need to specify the phenotypes of the N1i′ offspring of species 1 and of the N2i′ offspring of species 2 in patch i. We assume that the individuals that passed the viability selection sieve reproduce asexually (see Supplementary Section 7 for a brief discussion of the effect of recombination) and that the mother of each offspring is chosen randomly, with replacement, among the survivors. We recall that the numbers of surviving individuals of species 1 and species 2 in patch i are N_1i and N_2i, respectively, and that all survivors have the same probability of being chosen as mothers regardless of their fitness. Hence, selection works at the level of survival (viability selection) only and not at the level of the (genetic) differences of reproduction of survivors. In this way the reproductive output reflects the ecological dynamics and not the population composition at each generation.

The differences between mother and offspring are due solely to mutations in the m_r non-plastic loci and in the m_p plastic loci, which were implemented as follows. Each allele of the offspring can mutate with probability u_r or u_p depending on whether it is a non-plastic or a plastic allele. (Here, we assume u_r = u_p = 5/1, 000, which gives a genome-wide mutation rate U = 0.2.) Once a mutation occurs, say at the plastic locus k, we add a normal random variable ξ of mean zero and variance 1/100 to the existing allelic value which then becomes P_k+ξ. This is Kimura's continuum-of-alleles model (Kimura, 1965). As usual, generations were discrete and non-overlapping. During the development of an individual in a particular patch, we ignored any potential influence of parental phenotypes as, e.g., transgenerational plasticity (Uller, 2008).

In sum, the offspring generation of species l = 1, 2 in patch i is obtained by selecting with replacement Nli′ individuals from the N_li survivors of the selection sieve. The selected survivors are referred to as mothers. The phenotype differences between offspring and mothers are the mutations in the non-plastic and plastic loci.

Each individual within each patch can migrate to one of the eight surrounding patches (Moore neighborhood) with probability p_mig, and its destination is equally likely to be any of the eight patches. We were careful in keeping track migrant and non-migrant individuals that remained in their natal patch (Hassell et al., 1995). The flow of migrants between patches takes place simultaneously and it may result in some patches becoming empty or exceeding the carrying capacity. After arriving at their destination patches, the migrants as well as the residents of those patches pass the viability selection sieve as described in Section 2.2. Again, some patches may become empty at this stage.

As originally defined by Levins (1969), metapopulation dynamics consists of the extinction and colonization of local populations. Early models analogous to Levins' showed that two competitors could coexist globally even if coexistence was impossible in a single patch (Levin, 1974; Slatkin, 1974; Nee and May, 1992). However, here we do not impose a random extinction probability; only local adaptation and gene flow may interact to promote genetic variation and coexistence in the metapopulation. Actually, real metapopulations may contain local populations that never go extinct (Schoener and Spiller, 1987). A cautionary note: Kawecki and Ebert (2004, p. 1232) rightly pointed out that “local adaptation is about genetic differentiation”, but warned to minimize non-genetic effects such as plasticity (thus considering it a “nuisance parameter”) when studying local adaptation. However, at the metapopulation level studied here the total phenotypic variation for plastic species 2 is the result of the variation in the reaction norm intercepts (first term on the right side of Equation 2) and slopes (second term on the right side of Equation 2), both of which have a genetic basis (Scheiner, 1993; Sommer, 2020). Therefore, local adaptation is better understood as how close the mean phenotype matches the patch's environmental optimum value.

In this paper we consider only an expanding population scenario. More pointedly, the initial population was located on a randomly selected patch of the 2-dimensional grid at carrying capacity K_max for the two species, each at equal frequency (i.e., K_max/2 individuals from each species), and all other patches were empty. We recall that the carrying capacity of the metapopulation is KmaxL2 individuals, so there is plenty of room for expansion from this initial setup. For each species independently, there was an equilibration period of 2,000 generations at the seed patch before the colonization of the empty patches started. We initialize the allelic values R_k, k = 1, …, m_r and P_k, k = 1, …, m_p in Equation (2) for each individual as described before. In the equilibration period the two species evolve independently in the seed patch, i.e., there is no interspecific (as well as intraspecific) competition since after viability selection one guarantees that there will be exactly K_max/2 offspring of each species. Thus, Equations (4a) and (4b) are not used in the equilibration period. In sum, during equilibration viability selection decreases the population of each species by eliminating the less fit individuals and reproduction resets the population of the seed patch to its original size.

After the equilibration period, the colonization of empty patches starts. The order of events is migration, phenotype determination, viability selection, ecological competition and reproduction. The sequence of these five events comprises one generation. Note that this sequence of events guarantees that a given individual undergoes the processes of selection, competition and reproduction within the same patch and that only their offspring have the possibility to migrate to neighboring patches.

A word is in order about the ecology that our model describes. Consider a particular patch, say patch i, at a moment just after migration, so its population consists of the offspring that stayed in patch i and those that migrated to patch i. We recall that the model assumes that only the offspring migrate. To reach the reproductive age, these offspring must pass the selection sieve in patch i. Those who passed this sieve become adults: they are the survivors, which amount to N_1i individuals of species 1 and N_2i individuals of species 2. The survivors compete among themselves in patch i to secure the resources to support their potential offspring. This competition is described in a coarse-grained manner by Equations (4a) and (4b), which output the number of offspring that each species can give rise to and sustain in patch i, viz., N1i′ and N2i′. At this point, we could argue that the survivors produce an infinite number of offspring but only N1i′+N2i′ of them survive because of resources limitation. Alternatively, we could argue that the survivors produce exactly the number of offspring determined by Equations (4a) and (4b). This last interpretation is the one adopted in population dynamics (Godfray et al., 1991; Pásztor et al., 2006), from where we have borrowed those Ricker-like equations. In any case, assuming one or the other scenario would not affect the outcomes. Next, the mothers of the N1i′+N2i′ offspring are chosen randomly with replacement among the survivors of each species. Behind the coarse-grained approach is the assumption that adults of the same species are indistinguishable with respect to their competitive and reproductive abilities. Finally, each offspring decides if it will stay in patch i or move to one of the neighboring patches.

Individual-based simulations were independently implemented in Fortran and in MATLAB (2020) algebra environment using tools supplied by the Statistics Toolbox. Simulation results were double-checked by different authors to avoid any potential error. The results presented here are based in the Fortran code because it has speed advantages over MATLAB. The variable parameters were: L (landscape dimensionality), σe2 (variance of the environmental values E_i), ρ (environmental correlation), K_max (patch's carrying capacity), c (plasticity cost), and p_mig (migration probability). For each set of conditions, we run 1,000 independent simulations (a random landscape for each simulation). The metapopulation dynamics was run for at most 2,100 generations and we used the last 100 generations to average over the quantities of interest (e.g., the abundance of each species) in the equilibrium regime. If one of the two species fixed before that upper limit, we halted the dynamics. Otherwise, we considered that coexistence was achieved. However, in the study of the single-species metapopulation dynamics all runs reached the upper limit of 2,100 generations. In the Supplementary material, we present many instances of the time evolution of both species (e.g., Supplementary Figure 12), which show that the running time of 2, 000 generations is sufficient to guarantee that the metapopulation dynamics reaches the equilibrium regime.

Since the quantities used to characterize coexistence at equilibrium are averages over patches (typically L² = 400), last generations of the colonization stage (100) and runs (typically 500 runs result in coexistence), the number of samples used to estimate their mean values is very large, resulting in error bars smaller than the sizes of the symbols used in the figures. However, in order to assess the variability of the equilibrium variables described next, in Supplementary Section 8, we offer a variety of scatter plots for selected values of the model parameters.

In this paper we aim at the characterization of the metapopulation in the equilibrium regime, defined as the regime between generations t = 2, 000 and t = 2, 100. In the Supplementary material, we present results for the time evolution of both species in a variety of scenarios, but here we consider the equilibrium regime only. We focus on the following four variables.

To facilitate the interpretation of these variables, in Supplementary Section 3, we offer snapshots of the grid where the abundances the two species in each patch is shown in a color scale.

As the uncertain heterogeneous environment poses an adaptive challenge to both species through the viability selection sieve, it is instructive to study the metapopulation dynamics separately for each species before considering the competition between them. In addition, for the runs that do not result in coexistence, the equilibrium of the metapopulation is described by the single-species dynamics. As before, the initial single-species population was located on a randomly selected patch of the 2-dimensional grid at carrying capacity K_max and there was an equilibration period of 2,000 generations before the individuals were allowed to migrate to the neighboring patches.

We consider now the general setup where the two species are first let to reach equilibrium independently of each other in the seed patch and then are allowed to compete and migrate to the neighboring patches. Of course, the focus here is on the runs that led to coexistence since the runs that do not lead to coexistence were already fully characterized in the previous subsection.

Figure 4 summarizes the effects of the environment on the probability that a run results in coexistence, which is measured by Γ (Figure 4A), and on the fraction of patches that harbor the two species, which is measured by 〈〈Π〉〉 (Figure 4B). To a good approximation the effect of the environment is represented by the single variable σe2(1-ρ), which means that ρ can be absorbed in σe2 and we can study the uncorrelated landscape only without loss of generality. In other words, increasing the correlation between patches is equivalent to decreasing the variance of environmental values in an uncorrelated landscape. The important message from Figure 4 is that the plastic species 2 is extinct in a quasi-homogeneous or smooth environment (i.e., for σe2(1-ρ)≈0). We note that in this region there are no data for 〈〈Π〉〉 because no run resulted in coexistence.

Interestingly, increase of the environment roughness has only a limited effect on the probability of coexistence Γ, which quickly levels out and remains unaffected by further changes on σe2 (Figure 4A). The probability that a patch exhibits coexistence 〈〈Π〉〉 displays a more interesting behavior (Figure 4B). For smooth environments, most patches are occupied by species 1 only, but as the environment roughness increases, those patches begin to harbor both species. The slow decrease of 〈〈Π〉〉 we observe for large σe2 is due to the appearance of patches occupied by species 2 only (see Supplementary Figures 9–11).

Figure 5 shows the environmental effect on the relative abundances of both species. For smooth environments, species 2 is present in a few patches only (Figure 4B) and so its relative abundance 〈〈n₂〉〉 must necessarily be small, even if its density is high in the patches where it is present. In fact, the relative abundances are informative only when 〈〈Π〉〉≈1, in which case they represent the proportions of each species within a patch. The low density of species 1 for rugged environments is an indication that there may be patches occupied by species 2 only, which supports our explanation for the decreasing of 〈〈Π〉〉 for increasing σe2. We recall that within-patch coexistence can happen only if the density of species 1 is below the threshold 1/a₂₁ = 2/3, which is indicated by the dashed horizontal line in Figure 5A. Otherwise, we can observe among-patch coexistence only, in which species 2 occupies patches characterized by extreme environment values that are not suitable to species 1.

Figure 6 shows the effect of migration on species coexistence for an uncorrelated landscape (ρ = 0). Increasing the migration probability p_mig has an effect similar to increasing the environment ruggedness. As pointed out in our study of the single-species dynamics, migration affects the adaptation of species 1 but has little to none influence on the adaptation of species 2 when phenotypic plasticity is non-costly. Hence, the increase of the abundance of species 2 with increasing p_mig shown in the figure is a result of the effect of migration on the abundance of species 1 which in turn affects species 2 in the ecological competition stage.

The parameters σe2, ρ and p_mig influence mainly the adaptation of the non-plastic species 1. The plasticity cost c, however, affects the plastic species 2 only and Figure 7 shows its effect on species coexistence. For the migration probability considered (p_mig = 0.3), species 1 cannot prevent invasion (and, consequently, coexistence) but for large c species 2 cannot take advantage of the maladaptation of species 1. We note that it is the presence of species 1 that drives species 2 to extinction, since species 2 alone can thrive for large c by turning off the plastic alleles (Figure 3). The data missing for c = 0.4 is because none of the runs resulted in coexistence.