Mutualisms (bidirectional, beneficial, interspecific interactions) occur across a large number of organisms and encompass a diversity of benefits (Bronstein, 2015). Despite arguments over the exact fraction of interactions in a community being mutualistic, there is a general consensus that mutualisms are fundamental interspecific interactions at the level of the community and ecosystem (Bronstein, 1994;Dodds, 1997). Mutualisms were traditionally classified as "obligate" or "facultative" although it is accepted now that mutualisms lie along a continuum (Bronstein, 2015;Ollerton, 2006). This continuum consists of a focal species with varying degrees of dependence on its mutualist partner from complete dependence (obligate mutualisms) to intermediate dependence (facultative mutualisms) to no dependence (independent growth, i.e., not a mutualism). In other words, dependence is how much of the focal species' growth is contingent on its own intrinsic growth rate versus on its partner mutualist (Janos, 2007;Moyano, Rodriguez-Cabal, & Nuñez, 2020).
Mutualisms, however, are context-dependent interactions. A focal species often increases the population densities of their partner through exchange of resources (e.g., carbon, nitrogen, pollen) or services (defense, pollen transfer) but the provision of these resources or services can be energetically expensive (Miller et al., 2002;Pyke, 1991;Southwick, 1984;Stanton & Palmer, 2011). These costs of mutualism can reduce focal species growth potentially turning the interaction exploitative. This shift in interaction outcome is shown in empirical work and meta-analyses where mutualisms turn parasitic (exploitative) based on several factors including species densities or climactic gradients (Aizen et al., 2014;Chamberlain et al., 2014;Drew & King, 2022;Koide & Schreiner, 1992). Density-dependent shifts in interaction outcome have also been predicted by a previous theory that models mutualisms as a bidirectional consumer-resource interaction (Holland & DeAngelis, 2009, 2010). Variation in space of mutualist interaction frequency and outcomes due to density effects can thus influence the spatial range dynamics of the partners.
Several examples exist of a species' range expansion being hindered or facilitated by mutualistic interactions based on partner density. Von Holle and Simberloff proposed the term "invasional meltdown" where "…a group of nonindigenous species facilitate one another's invasion in various ways, increasing the likelihood of survival and/or of ecological impact, and possibly the magnitude of impact" (Simberloff & Von Holle, 1999). For instance, in plant-fungal mutualisms, partner presence is important for establishment of invasive species in novel environments with fungi having been hypothesized as the driving factor for European plant invasions in the Holocene era either through formation of novel mutualisms or through co-invasion of both mutualist species (Nuñez & Dickie, 2014;Simonsen et al., 2017;Wilkinson, 1998). Based on partner presence (or absence) and their densities, range expansion of species into new habitats can be slowed down or even stopped both in plant-microbe and plant-pollinator mutualisms (Barthell et al., 2001;Chen et al., 2018;Nuñez et al., 2009;Parker, 1997;Stanton-Geddes & Anderson, 2011). However, much less is understood about how partner dependence along with the population dynamic effects of mutualisms can shape range expansions and their limits.
Specifically, there is no consensus on whether increased partner dependence speeds up or slows down range expansions. Two opposing ideas may explain these patterns. First, species that do not depend on a mutualist partner will spread fastest (Ideal Weed Hypothesis) as its own invasion is not contingent on its partner's ability to disperse (Baker, 1965). Second, conversely, greater dependence on a partner can facilitate the establishment of the focal species by increasing their growth and thereby persistence in newly colonized regions in space (Kubisch et al., 2014;Nara & Hogetsu, 2004). Indeed, pine tree species with a greater dependence on mutualist fungal partners have increased spread capacity (Correia et al., 2018;Moyano et al., 2021;Moyano, Rodriguez-Cabal, & Nuñez, 2020). Current theoretical models of mutualisms explain only a subset of these empirical results. They predict that mutualist species will evolve lower dispersal distances, reducing their ability to colonize new habitats and slow down their range expansion (Chaianunporn & Hovestadt, 2012;Kubisch et al., 2014;Mack, 2012). This phenomenon occurs as an expansion into new territory in the absence of a partner that can lead to reduced growth or repeated extinction events. These models, however, do not explicitly incorporate degree of dependence and hence cannot explain when more dependent species should spread faster (Moyano et al., 2021;Moyano, Rodriguez-Cabal, & Nuñez, 2020). Accomplishing this outcome requires a new theory under which greater or lesser dependence on a partner can lead to the fastest range expansion of a species.
Here, we address this gap by constructing a spatially explicit mechanistic model of a two species mutualism. We consider mutualisms along a spectrum ranging from no dependence (independent growth) to complete dependence (obligacy). We determine when partner dependence could enhance or slow down range expansions relative to independently spreading species. We explored how interspecific differences in the magnitude of benefits exchanged shaped the degree of dependence that resulted in fastest species range expansion. We next observed how differential partner dependence generated distinct spatial population distributions for a focal species. Finally, we identified when mutualisms could turn exploitative due to density effects, thereby leading to range limit formation of a focal species.
To study the role of mutualistic interactions on species growth and subsequent dispersal, we set up and analyzed a pair of integro-difference equations (IDEs). IDEs assume that the species' annual cycle has two distinct phases (growth and dispersal) that occur sequentially (Lutscher, 2019). The growth and interaction phase occurs from some time period t to t + T followed by dispersal at the end of the year that is (t + 1). Growth is modeled by coupled ordinary differential equations (ODEs) for the two species while the dispersal of each species is governed by their dispersal kernel. We tracked the speed at which species expand into new regions in space along with how their population size is spatially distributed. The general functional form of our IDEs are as follows:
where k P and k F are the dispersal kernels of the two species P and F and M 1 and M 2 are nonlinear growth functions that describe how the abundances of species P and F change over time.
The growth of species P and F at each point in space was modeled using ODEs. Our model used a bidirectional consumer-resource framework of interactions between P and F similar to Holland and DeAngelis (2010). Our model best modeled nutritional mutualisms (e.g., plant-fungal or legume-rhizome mutualisms) or plant-pollinator mutualisms where resources/benefits exchanged (carbon [C] and phosphorus [P], C and nitrogen [N], nectar and pollen transfer respectively) contribute directly to species growth rates. In these mutualisms, the costs incurred would be resources a species uses to make the C, N, P, or other products. We include dependence to study different mutualisms along a facultative to obligate continuum. In addition to their intrinsic growth and death rates, the species' population sizes are governed by the benefits of mutualism obtained and costs incurred to and from their partner. The equations are given by:
where r P and r F are intrinsic density-independent growth rates and d P and d F are death rates. The dependence of each species on its partner is denoted by D P and D F for P and F respectively with values between 0 (i.e., species growing independently) and 1 (complete dependence on partner for growth). When there is no partner dependence, P (the focal species here) logistically grows to its carrying capacity r P d P . Similarly, when
The second term of Equations ( 2a) and (2b) describe the benefit obtained by the focal species from its partner where α PF is the maximum benefit offered to species P by species F and α FP is the maximum benefit to F by P.
Benefit from the mutualism is mediated by the population size of the partner (F) and the dependence of the focal species (D P ). The total growth of a species (say P) is the sum of its own intrinsic growth and partner benefits weighted by its degree of dependence (D P ) resulting in some trade-off between intrinsic growth versus benefit uptake from partner. Such trade-offs between reproduction and root architecture have been observed in invasive forbs with different degrees of dependence on their mycorrhizal fungal partner (Seifert et al., 2009). For P, when F is small, the benefit obtained (increase in per capita growth rate) scales linearly with F and saturates with increasing F where h P and h F are half-saturation constants mediating the rate at which benefits are provided to partner species. Such saturating benefits capture physiological limits to uptake or handling of resources (Soberon & Del Rio, 1981;Wright, 1989) and have been observed in empirical systems such as the fig-fig wasp, ant-treehopper, ant-aphid systems and have been used theoretically for plant-pollinator dynamics (Breton & Addicott, 1992;Bronstein, 2001;Morales, 2000).
Although, mutualisms are associated with mutual benefit exchange, costs are incurred as well (Bronstein, 2001;Frederickson & Gordon, 2007;Morris et al., 2010;Stadler & Dixon, 1998). Costs can include resources used to produce rewards for partners, energetic costs to form structures to attract partners, or cheating behaviors of the partner resulting its exploitation. In some plant-fungi and plant-pollinator interactions, increased partner densities increases costs and reduces growth of the plant species thereby turning the interaction exploitative (Aizen et al., 2014;Koide & Schreiner, 1992). Hence, costs of mutualisms also depend on partner population densities yet, to our knowledge, there are no empirical data on specific functional forms. In our equations, the third term describes the cost incurred (reduction in per capita growth rate) by the focal species while providing a benefit. This cost increases linearly with F's (partner) density ( β PF e P ) when P's density is low but tends to zero when P >> F. For a fixed F, as P's density increases, the cost paid by each individual of P is increasingly diluted. e P , and e F are half-saturation constants mediating the rate costs are incurred. β PF is the coefficient of cost inflicted on species P due to F and β FP is the coefficient of the cost on species F by P. Combining the costs and benefits of the mutualism (second and third terms of the equation), the per capita growth rate of a species, say P (given by dP Pdt ), increases with increasing F until a maxima following which it starts decreasing with increasing F. This effect of increasing partner density turning P's per capita growth rate from positive to negative has been observed empirically (Gange & Ayres, 1999).
The growth functions that are denoted by the equations M 1 and M 2 in Equation 1 and are solutions to the differential equations (2a and 2b) that describe the growth of the two species for the time period T. Equations (2a) and (2b) are both analytically intractable and hence almost all analysis performed on this model is numerical. We used MATLAB (v. R2021a) for our simulations and solved our ODEs with the ode45 package.
Dispersal follows the growth and interaction phase. Species dispersal is governed by their dispersal kernel, a probability density function describing the probability of an individual dispersing to and establishing at a location "x" given it started at another location "y." Here, we assume that the dispersal kernels of both species are equal and follow a Gaussian kernel, one of the most commonly used functional forms for modeling dispersal (Gilbert et al., 2014). The dispersal kernels are given by:
where σ 2 P and σ 2 F are the variances of the dispersal kernels of species P and F. We assume that the variances of both species' dispersal kernels were the same (σ 2 P = σ 2 F = 0.25).
For the simplest version of our model where there exists no cost to providing the partner species with mutualistic benefits (i.e., β PF ¼ β FP ¼ 0) and the dependence of each species on the other is the same (D P ¼ D F ), we can analytically obtain expressions for range expansion speeds for the two species. We identified the degree of dependence on the mutualism that led to greatest range expansion.
We performed all analysis of our model with costs (β > 0) numerically. Incorporating mutualism costs makes the equations analytically intractable and hence all further results arise from simulations. We initialized the simulations with extremely low densities of individuals (P = F = 0.01) at the zero spatial coordinate of a 1-D line of length 1200 units (-600 to +600 units). The total length is subdivided into a total of 65,537 equally spaced nodes in order to have a high spatial resolution (see Appendix S1: Figure S1; for single species example of spatial spread). We allowed growth for T = 10 time steps following which we allowed for a dispersal event to occur for both species. The choice of 10 time steps of growth between dispersal events is arbitrary but altering the number of time steps does not qualitatively affect our results. Following dispersal, the growth and interaction phase occurs again. We allowed for 60 such iterations (each iteration took the system from time "t" to "t + 1"). This time duration was sufficient for the simulations to reach steady state that is, species have either attained a constant spread speed or remained stationary (the distance of the edge from the center is constant over consecutive time points). We located the range edge of our populations by identifying the farthest point in space where the density of the species was above a threshold (density > 0.05). To quantify the range expansion speed we took the difference of the location of the range edge between the final and penultimate time steps (between t = 59 and t = 60) in the simulation. We define range limits as those edges formed when the speed of range expansion is zero. This outcome occurred when a species' per capita growth rate was less than or equal to zero. To generate our results, we only considered parameter regions that exhibit mutualism at low density (i.e., α > β). In some cases, interactions shifted from mutualism to exploitation as densities changed through the simulation, consistent with empirical understanding of some mutualistic interactions. Across all simulations, we held constant the intrinsic death rates of the species (d P = d F = 0.1) and the half-saturation constant of benefits provided and costs associated with the mutualisms (h P , h F , e P , e F = 0.3). We varied several parameters including the magnitudes of benefit (0 ≤ α ≤ 2), cost (0 ≤ β ≤ 1), and intrinsic growth rates (0 ≤ r ≤ 1). Further, we considered a range of dependence values where 0, 0
Þ with a step size of 0.1 (full list of parameter values in Table 1). We did not explicitly define units for our model parameters as we aimed to build a system agnostic model. We could apply this model to large trees with wind dispersed seeds where dispersal every year could be several hundred meters, or it could be smaller bacterial or fungal species whose dispersal distances are a few centimeters/meters.
We considered four sets of simulations: two to address mutualism driven change in range expansion speed (Sets 1 and 2), and two to explore how spatial distribution and range limits can arise from mutualistic interactions (Sets 3 and 4). We observed Set 1: Symmetric species interactions (α PF ¼ α FP ,β PF ¼ β FP , D P ¼ D F ) to understand when increased dependence speeds up versus slows down range expansion speed. With no costs, we can derive results analytically as a function of all model parameters. With costs, we simulated for two different magnitudes of benefit (α PF = α FP = 0.6 (low benefits) and α PF = α FP = 1.2 (high benefits). We alter, in a pairwise fashion, growth parameters (α and β) to identify the degree of dependence that leads to fastest range expansion. We varied α (0
. We then measured the speeds of range expansion for the two species. We present results for cases where our focal species (P) invades into regions where its partner (F) has already invaded and grown. The results are qualitatively similar to when both species co-invade into new territory (Figure 1).
Set 2: We next considered asymmetric magnitudes of benefits exchanged (i.e., α PF ≠ α FP ) and analyzed how these differences could alter the degree of dependence that allows for greatest speeds of expansion. We varied α PF and α FP (0 ≤ α PF , α FP ≤ 2). We also explored how interspecific differences in growth rates and costs of
Table of parameters and their ranges. Parameter Meaning (units) Parameter range P Population density of species P (total numbers) … F Population density of species F (total numbers) … α ij Mutualist benefit offered by species j to species i (i,j = P or F
Cost of providing benefit to species i by species j (i,j = P or F; i ≠ j) mutualisms affect species range expansion speeds (Appendix S1: Figures S7 and S8). Here, we swept pairwise through two parameters while holding the others constant. We varied r P and r F (0 ≤ r P , r F ≤ 1) and swept through all possible pairs of dependence ( 0, 0
We repeated this while sweeping through α PF and α FP (0 ≤ α PF ,α FP ≤ 2) and β PF β FP (0 ≤ β PF , β FP ≤ 1). Here, our simulations were initiated with co-invasion of both species into new ranges (Figure 2).
Set 3: We simulated range expansion of mutualists with different degrees of dependence on each other while holding other model parameters constant. We looked at differences in the spatial distribution of both species at the range front and core for the different cases and analyze the population dynamics of both species in a nonspatial analog of the model (i.e., without spatial spread). Here too, we initialized our simulations with co-invasion of species P and F in the spatial model (Figure 3).
Set 4: We obtained analytical conditions when range limits for P were generated through its interactions with F. We used this expression to show how partner density could influence whether P continues to expand or whether its range limits form (Figure 4).
We found that increased dependence could speed up or slow down range expansions, and that intermediate dependence could lead to the fastest spread if species interactions were asymmetric. We also found that mutualisms that affected the degree of dependence could affect the range expansion speeds and spatial distributions of the interacting species. Furthermore, we identified the conditions when range limits of species could form as a result of its interaction with the partner species. We described each result in a subsection.
We see that mutualisms can speed up or slow down range expansion speed based on the magnitude of benefits exchanged. The speed of a species' range expansion wave is contingent on magnitude and symmetry of benefits exchanged as elaborated below (see Appendix S1:
Table S1 for summary). First, we consider two identical species exchanging symmetric benefits (in terms of magnitude). In the simplest model with no costs (β ¼ 0), we obtain an analytic expression for the spread speed of the focal species P (c P ; Appendix S5.1). We find that c P ¼ ffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi 2σ P 2 1 -D P ð Þr P T p . Because there is symmetry in benefits exchanged, these expressions hold for species F (c F ) as well. This speed for focal species (say P) arises when it co-invades into new territory with its partner F from low densities. Hence, it receives negligible benefits from F. In this co-invasion case, P spreads fastest whenever D P ¼ 0 (no dependence on F). With increasing D P , P's speed reduces and goes to zero if it is obligately dependent on F for growth. The same holds for F.
Alternatively, we can derive an expression for expansion speed when P spreads into a region where F has already spread to and is present at high density, that is, its carrying capacity. When P invades this region and forms a mutualism with F, we can calculate the spread speed of P (Appendix S5.1) to be:
Here, when the total benefits received from the partner are greater than its own intrinsic growth rate
, P spreads fastest when it is completely dependent on F (D P ¼ 1). Conversely, when
, P spreads fastest when D P ¼ 0. Thus, greater dependence of a species on their partner (species F) slows down P's expansion speed when the benefits (α) provided by each species to the other is low. Here, being dependent on a partner slows P down as it receives lesser benefits from F resulting in slower growth. Furthermore, as P provides fewer benefits, F's own population size reduces, thus limiting the cumulative amount of benefits it can provide to P. Conversely, greater dependence increases the speed of range expansion when benefits exchanged (α) are large. These results arose because the per capita growth rate of a species (say P) was mediated both by its own intrinsic growth rate and partner benefits (which is a function of partner population size, maximal benefits received (α), and dependence of P on F).
The dependence on F that maximizes P's growth will ultimately lead to its fastest spread. This pattern holds even when a cost of providing the mutualism is introduced to the model where P invades into regions with F (Figure 1).
We next introduced a cost to engaging in a mutualism (e.g., provisioning carbon, nectar rewards, oviposition sites) and consider two cases of spread: co-invasion of partners and invasion of focal species into region where partner is present. In both cases, we see that increasing dependence on a partner can slow down spread when the maximal benefits (α) exchanged are low. This is slightly different to our previous result (Equation 4; where β ¼ 0) in that regardless of co-invasion or invasion into partner's range, the focal species spreads slower at higher dependence if costs of mutualism are high. However, when α is high, range expansion speeds increase with increasing dependence on the mutualism following an initial dip. When dependence (say D P ) is low, the focal species acquires minimal benefits from its partner, whereas the costs incurred increase almost linearly with F's density. Thus, at low dependence, the effects of a cost of mutualism slows down a species relative to nonmutualists. But with increasing D P , the benefit from the partner outweighs the cost and results in increased speed of P. In summary, when magnitude of benefits exchanged between partners is symmetric (α PF ¼ α FP ), independence or complete dependence on partner leads to fastest speed of spread of a focal species regardless of co-invasion or invasion into partner species range. Our results are robust across a wide range of α, β ð Þ values (Appendix S1: Figures S2 and S3 for P and F, respectively).
We also varied (r, α) and (r, β) in a pairwise fashion. We found, in the former, that when r is low to medium, with increasing α, species spread fastest when they completely depend on their partner (Appendix S1: Figure S4). In the latter, complete dependence leads to fastest spread when species growth rates are low and the cost of mutualism is also low to medium in magnitude (Appendix S1: Figure S5).
We next consider asymmetric benefits, by relaxing the assumption that the benefits exchanged by the partners are equal (i.e., now α PF ≠ α FP ). As the benefits exchanged are often in different currencies (pollination, defense, nutrition provisioning), their impacts on the growth of species that receive the benefits can be very different. When benefits exchanged are unequal or asymmetric, the degree of dependence leading to fastest spread is no longer complete (D P ¼ 1) or no dependence (D P ¼ 0) on partner. Furthermore, we also assume that both species can have a different degree of dependence on their partner. In some cases, P spreads fastest when it depends to an intermediate degree on F while F completely depends on P (seen in top left region of Figure 2). This occurs when α PF has a large value (1:8 ≤ α FP ≤ 2:0) and is greater than some intermediate α FP (1:0 ≤ α PF ≤ 1:2). Similarly, when α PF has an intermediate value and α FP is large (1:8 ≤ α FP ≤ 2:0), P spreads fastest when it completely depends on F and F has some intermediate degree of dependence on P (see bottom right region of Figure 2c). A focal species (P) spreads faster when it depends less on its partner in cases where it receives greater benefits (α PF > α FP ).
Let us consider when α FP > α PF (here D P > D F ). D P is at its maximal value (=1) so P can make use of all the benefits from its partner. If D F is large, then F can do the same and grow even faster than P (as α FP > α PF ). However, if D F is too large, then there is an extremely large growth in population size of F, which also places a greater cost on P. This results in reduced growth and spread of P. Reduction in P's numbers will lead to lower total benefits offered to F thereby slowing down F's growth as well. Hence, P will spread fastest when F has an intermediate D F where it gains enough benefits to grow and spread with P while also not inflicting too high a cost on P (by growing too fast) and to prevent its extinction.
When α FP and α PF are symmetric and low ( ≈ 1), P spreads fastest when it spreads independently. This means that the dependence of each species on the other would be close or equal to 0 (D P , D F ≈ 0). Engaging in any type of mutualism in such scenarios slows down their speed of range expansion. When the benefits exchanged are symmetric and high (α PF , α FP ≥ 1:4), complete dependence of both species on the other leads to fastest speeds of range expansion (Figure 2). These results are in line with our expectations from the previous section where we study symmetric benefit exchanges of different magnitudes. We show our results for spread of species P but the results are symmetric for species F (Appendix S1: Figure S6).
We further varied pairs of values of r (r P and r F ) to understand when our focal species (P) spreads fastest. When r P ≈ r F ≈ 0, P and F spread fastest when they are completely dependent on each other (Appendix S1: Figures S7, right panel, bottom left; S8, left panel, bottom left). As partner growth rates increase relative to it (say r F >> r P ), focal species P spreads fastest when its dependence (D P ) increases while its partner continues to spread independently (i.e., D F ¼ 0; Appendix S1: Figure S8). In such cases, we find that P spreads fastest because it is dependent on F as it trails F through space. Presence of a partner in the new region enables the focal species to obtain benefits as it invades and hence grows and spreads faster.
Beyond affecting expansion speed, mutualisms alter the spatial distributions of the interacting species (relative to nonmutualist spread) across their ranges. We define this distribution of population size as the shape of the range expansion. We show that different degrees of dependence can lead to very different range expansion dynamics despite the two species having the same growth parameters. For a single species, the range expansion wave spreads at a constant velocity given by c ¼ ffiffiffiffiffiffiffiffiffiffiffiffi 2σ 2 rT p (where r is intrinsic growth rate of the species; σ is the variance of the species' dispersal kernel, and T is the length of time between two dispersal events). The population size at the range center (given by P ¼ r P d P ) is much higher than at the edge of the wave (P ! 0) (see Appendix S1: Figure S1).
When two species spread independently, the densities of the species at the core is far higher than at the edge (Figure 3a,d). The two waves spread with different speeds of range expansion proportional to the square root of their growth rates (Speed of P = 0.9878 units, F = 1.3973 units; r P ¼ 0:2, r F ¼ 0:4). P trails F with both species reaching their individual carrying capacities behind the front of the wave (Figure 3a). In Figure 3d, we analyze the nonspatial analog of this model studying their population dynamics over time where the equations defining the mutualism were solved at a single point in space. By analyzing the linear stability of the system, we see that there exists a stable fixed point reached by both species from a low density over time which is equal to its carrying capacity behind the range edge that is, ( r P d P , r F d F ) (Appendix S1: Section S4.2.1). The densities of the species at range core remain higher than the edge even in a completely dependent mutualism. When species engage in a completely dependent mutualism 3c,f), both species spread at the same speed (Speed of P = Speed of F = 0.7183 units). Behind the front, both species grow to the same population size. This is because the difference in their intrinsic growth rates does not affect how they grow due their complete dependence on their partner's benefits for growth. The population size equilibrium that is reached in the nonspatial model (Figure 3f) is a saddle point with both eigenvalues being real with opposite sign (Appendix S1: Section S4.2.2).
When there is an asymmetry in dependence between the species 3b,e), the shape of the spreading wave is no longer monotonic (strictly decreasing from core to edge) and there are oscillations behind the front of the wave for both species (Figure 3b). Species F travels faster than P and expands into new range by itself, whereas P trails and invades into regions already occupied by F (Speed of P = 0.6579 units, F = 1.0823 units). The trajectory of the populations in the nonspatial model in this case shows oscillations before reaching an equilibrium (Figure 3e). This equilibrium is a stable focus as its eigenvalues are complex conjugates with a real negative part (Appendix S1: Section S4.2.2). The oscillations arise due to overshoot in species growth followed by suppression of partner growth due to costs placed on them. This results in their reduced growth as well. Through these corrections, the population sizes of both species goes to a lower value and the oscillations dampen over time. At the edge of the population (in Figure 3b), F spreads faster than P but behind the front, P's density is much larger and drives F to lower densities than itself over time. Dynamics at the range edge of the spatial model mimic that of the early phases of the nonspatial model and dynamics at the range core mimic the latter phases of the nonspatial model.
Under certain conditions mutualistic interactions can also lead to formation of range limits for a focal spreading species. These range limits are formed the invading focal species (say P) cannot establish after dispersal, that is, its per capita growth rate at extremely low densities is nonpositive. Due to its dependence on its partner (F), the criteria for establishment and spreading through new regions is contingent on partner density. At the range edge of P, we can set P to 0 (as P is at extremely low densities in Equation 2a) to obtain an expression for per capita growth rate:
Furthermore, for invasion to be possible, dP Pdt > 0, so
Rearranging terms we find that,
is necessary for invasion to occur. This expression is a quadratic in F, and we can use it to identify whether expansion occurs or a limit forms for P across a range of F densities. Let us call this expression g(F). When g(F) is negative, invasion and spread is possible but spread will stop if g F ð Þ > 0 and there will be a range limit formed. In the case of co-invasion of P and F into new territory (F = 0), range expansion will always occur for P as long as D P ≠ 1. However, if F is of some finite nonnegative population size, whether P forms a range limit is contingent on the degree of dependence of the species on each other. If D F ¼ 0, a range limit does not occur for any D P value. If D F ≠ 0, then there can exist a range of densities for F for which invasion and spread of P is possible but beyond which there is a formation of a range limit (Figure 4). In Figure 4, although α PF >> β PF , the population size of F impacts whether there is expansion or a limit is formed highlighting the context-dependency of mutualistic interactions of species along their range.
Mutualisms have been shown to constrain and expand a species range (Fowler et al., 2023;Stephan et al., 2021). Here, we develop a spatially explicit population dynamic model of mutualism that explains how partner dependence alters the speed and shape of species range expansion. We find that increasing partner dependence can increase or decrease a species' range expansion speed based on the initial partner density and the magnitude and symmetry of benefits exchanged. Differential dependence on partners also influences population dynamics of species thereby altering the population distribution behind the expansion front. Finally, we see that mutualisms can turn exploitative due to density-dependent effects which generates species range limits. Our results show that accounting for mutualisms can provide insights into species' range expansion speeds and their persistence and distribution across space.
Resolving contradictions of the role of mutualisms on the speeds of range expansion described in previous literature
Our findings that increased dependence on a partner can either slow down or speed up range expansion bridges a conflict between previous theoretical and empirical studies. First, when we ignore mutualism costs (β ¼ 0), we derive an analytical expression and find that species co-invading into new regions will spread fastest when they are not dependent on their partner. This result agrees with previous theories predicting that mutualisms should lead to reduced dispersal and range expansion of species relative to species spreading independently (Chaianunporn & Hovestadt, 2012;Kubisch et al., 2014;Mack, 2012). Second, we find that when a species spreads into regions where its partner is already present, greater dependence leads to faster spread. For example, this happens when fungal species are dormant underground until they detect the presence of an invading plant to form a mutualism with (Callaway et al., 2004;Dickie et al., 2017;Nara, 2009). To our knowledge, we provide the first theoretical basis for empirical scenarios where greater dependence leads to faster expansion (Moyano et al., 2021;Moyano, Rodriguez-Cabal, & Nuñez, 2020). When costs are incorporated (β ≠ 0), we use numerical simulations to show that increased dependence speeds up or slows down range expansion based on magnitude of benefits exchanged (α) regardless of whether the species co-invade or spread into regions where their partner is present.
In contrast with the extremes above, we also find that intermediate dependence (facultative mutualisms) can lead to fastest spread, but only when benefits exchanged are asymmetric where a species receives greater benefits than it provides (α PF > α FP ). This is because, on the one hand, high dependence on a partner yields high growth benefits from the partner, which speeds up spread. On the other hand, higher dependence incurs higher costs, which leads to lower partner population sizes, slowing spread. Intermediate dependence balances maximizing benefits received from the partner while preventing density-mediated shifts from mutualism to exploitation. Existing data suggest that asymmetry in benefits between partners is true for most mutualisms (Bascompte et al., 2006). In this case, our model posits that intermediate dependence should most often drive fastest speeds of range expansion for a focal species. Indeed, recent work across >800 genera showed that plants with facultative mycorrhizal associations showed greatest invasion success into new ranges (Moyano, Dickie, et al., 2020). Thus, our findings highlight the importance of considering mutualisms along a continuum of independence to complete dependence as opposed to a binary.
Overall, consideration of partner dependence allows us to move past the verbal hypotheses that mutualisms should impact species range expansion rates to understand conditions under which they affect said rates (Svenning et al., 2014). Furthermore, it adds to previous literature that establishes how antagonistic interspecific interactions (such as predation, parasitism, and competition) alter the spatial distribution and abundances of the interacting species. While competition is generally thought to reduce expansion rates, both predation and host-enemy (parasite) interactions are shown to speed up or slow down spread due to a variety of mechanisms including predation rates, enemy release and spatial variation in prey abundances (Fagan et al., 2005;Holt et al., 2011;Keane & Crawley, 2002;Neubert et al., 2000;Okubo et al., 1989;Prakash & de Roos, 2002). Here we add to existing literature on how positive interspecific interactions speed up or slow down species spread based on degree of partner dependence.
Our results demonstrate that mutualisms can lead to several different spatial population distributions by impacting species population dynamics. In our model, a species spreading by itself has a higher population density in the range core (at carrying capacity) and declines monotonically toward the range edge (Figure 3a). We see the same pattern when both species completely depend on the other (Figure 3c). However, in other cases we find the opposite: that species densities are higher at the range edge than core (Figure 3b). In this scenario, F has a higher intrinsic growth rate than P and thus grows much faster from the initial low density (at the population edge). Furthermore, the intermediate dependence on P (D F ¼ 0:4) accelerates its own growth and suppresses the growth of P. Because benefit exchange between species in our model is density dependent, low population size leads to P's density being depressed at the edges due to exploitation by F (which is at higher densities). This result is analogous to "disease free halo" patterns observed in host host-parasite interactions when disease spread is density dependent (Antonovics, 2009;Bruns et al., 2019). Over time the dynamics reverse: the greater dependence of P on F (D P ¼ 1:0) means P slowly increases to a greater size than F (at the population core). Thus P manages to obtain benefits that allow it to spread behind F's front. Overall, this creates a pattern where the species' densities spike up near their respective range edges.
Our results also agree with empirical results across several taxa that species need not follow distance-abundance relationships (i.e., species have maximum abundance at the center of their range) (Dallas et al., 2017). This suggests that mutualisms are another biotic mechanism that confound distance-abundance relationships (for more hypotheses see Sagarin & Gaines, 2002), paralleling similar results in predator-prey and host-parasitoid systems (Frick et al., 2010;Hastings et al., 1997;Robinson et al., 2010). Finally, these spatial patterns in density in turn mean that partner benefits can vary between range core and edge. Such density-dependent effects are supported in meta-analyses of empirical host-symbiont mutualisms, which have found to change in the interaction sign between species (Chamberlain et al., 2014;Cunning & Baker, 2014;Drew & King, 2022). Our work demonstrates the consequences of this relationship in a spatially explicit setting, showing how mutualisms alter the spatial distribution of species.
In order to invade into new territory and expand their range, a species needs to be able to grow from low density in the presence of its partner species. The density of the mutualist partner species plays a very important role in the focal species' invasion ability as it confers benefits and inflicts costs. We find a nonlinear relationship between the partner species' density and the focal species' ability to invade. However, choosing a different functional form (e.g., fixed or saturating costs with partner density) can sometimes lead to qualitatively different outcomes in invasion success of a focal species with changing partner density (see Appendix S1: Section S3 for different benefit/cost functional forms). In the functional form considered here, we find that increasing density of the partner initially increases the per capita growth rate of the focal species in agreement with empirical evidence. However, we find that, above a threshold density, the partner imposes a large per capita cost on the focal species. The nonlinearity in benefits acquired from partners in our model results in range limit formation of a dispersing species, a pattern seen in plant-fungi mutualisms (Gange & Ayres, 1999). Past literature has shown repeated evidence that population decline shapes range limits (Moeller et al., 2012;Nuñez et al., 2009;Stanton-Geddes & Anderson, 2011). Our study shows that interaction sign (mutualist to exploitative), mediated by the effect of the species' relative population densities can also act as a factor causing these declines (see Table 2 in Sexton et al., 2009).
While our model shows how density-mediated shifts from mutualism to exploitation can impose species range limits, limits are often delineated through a combination of biotic and abiotic factors (which we do not consider in our model) (Lee-Yaw et al., 2016;Louthan et al., 2015). With global change, a species' range is often in a disequilibrium with suitable abiotic environments (Svenning & Sandel, 2013). Our model suggests that high dependence on partners at a species' leading edge should facilitate invasion while also buffering it against extinction in unsuitable environments (trailing edge) thus slowing down range shifts. Some evidence of this phenomenon exists in tree-fungi mutualisms, yet it is still not well understood (Lankau & Keymer, 2016). Conversely, low dependence could result in reduced partner benefits leading to faster extinction at the trailing edge resulting in accelerated range shift and contraction. Future evaluation of species range shifts must consider population dynamics and partner dependence to better predict their distributions and spatial overlap with partners (HilleRisLambers et al., 2013).
Our framework could be developed in several new directions. First, we assume that both species have identical dispersal kernels. One could explore the effect of dispersal differences both in distance and in shape of kernel. For example, seed dispersal of pine trees extends up to several hundreds of meters, while fungal spores tend to disperse smaller distances (Bullock et al., 2017;Galante et al., 2011). Absence of its mutualist could lead to the generation of range limits for pine trees. Second, in reality, mutualisms are of many varieties (e.g., dispersive, nutritional, pollination, defensive, etc.) resulting in different mechanisms and currencies of benefit (Hale et al., 2020). This can lead to different spatial dynamics for the species. Third, we assume our mutualism is highly specific. Several species that engage in mutualisms interact promiscuously with multiple species forming a mutualistic network (Bascompte & Jordano, 2013). Incorporating multiple species interactions along with explicit space would alter the community dynamics of the system that could consequently alter each species' spatial distribution. Finally in our model we do not assume any relationship between maximal benefits exchanged (α) and degree of dependence (D), but we compare all possible combinations α and D to identify which degree of dependence leads to greatest range expansion. In nature however, it is unlikely that dependence is uncoupled from the benefits of mutualism; real mutualisms are unlikely to have high dependence and low benefits.
Here, we show that mutualisms impact the speed and shape of a species' range expansion wave. We identify conditions in which increased partner dependence can speed up or slow down a species' range expansion speed. We further show that mutualisms influence the population dynamics of the species, thereby altering the shape of their range expansion wave. Finally we find that density-dependent shifts from mutualism to exploitation between partners can lead to formation of range limits of a focal species. We hope that our work, along with suggested future directions, will serve as interesting questions for other researchers to pursue and help develop this understudied field.
of 15 NAVEN NARAYANAN and SHAW 19399170, 2024, 1, Downloaded from https://esajournals.onlinelibrary.wiley.com/doi/10.1002/ecy.4171, Wiley Online Library on [09/01/2024]. See the Terms and Conditions (https://onlinelibrary.wiley.com/terms-and-conditions) on Wiley Online Library for rules of use; OA articles are governed by the applicable Creative Commons License
19399170, 2024, 1, Downloaded from https://esajournals.onlinelibrary.wiley.com/doi/10.1002/ecy.4171, Wiley Online Library on [09/01/2024]. See the Terms and Conditions (https://onlinelibrary.wiley.com/terms-and-conditions) on Wiley Online Library for rules of use; OA articles are governed by the applicable Creative Commons License