Extinction dynamics in experimental metapopulations

Data Analysis

Synthetic Measures. Metapopulation dynamics were determined by summing the number of individuals for each local population at each generation. From the resulting time series, we calculated metapopulation persistence time (set to a maximum value of 16 for metapopulations that persisted for the duration of the experiment) and quantified metapopulation cycling by using a first-degree autoregressive (AR) analysis. For this latter analysis, we truncated all dynamics at generation 9 to remove a possible bias in AR coefficient estimation that might result because of fewer data points in metapopulations undergoing earlier extinctions. We then ran an autoregressive integrated moving average (ARIMA) statistical model (one degree of autoregression, one degree of differencing) on these truncated dynamics to separate oscillations in metapopulation density from a possible trend in the evolution of density over time. The previous calculations yielded 16 independent estimates (4 for each treatment) for both persistence times and AR coefficients. Because these statistics were non-normally distributed, we calculated the difference between all possible treatments and then performed permutation tests to determine significance. For example, when comparing persistence times between the continuous and the isolated treatments, we computed the difference in mean persistence time. We then performed all possible permutations of persistence times between treatments, estimated the difference in mean persistence time in the permuted samples, and compared the difference with the actual value. The P value in this unilateral test is the proportion of permutations that yielded a higher (lower) difference than the observed value, when this observed difference is positive (negative).

Two other synthetic measures of population performance were calculated: metapopulation capacity (10) and population coherence (16). Metapopulation capacity incorporates both the size and location of each local population and has been shown to correlate with persistence time (10). In our system, local populations had the same area, so the only landscape feature that varied between treatments was the distance, and hence migration rate, between populations. Technically, we estimated metapopulation capacity as λ, the leading eigenvalue of the matrix of migration between pairs of populations (elements on the diagonal of this matrix being set to 0).

Population coherence occurs when local dynamics become synchronous. Synchronous oscillations can result in an enhanced extinction rate because local populations that cycle in phase will experience a simultaneous decline. This synchronous cycling will reduce the possibility of a rescue effect and make populations more susceptible to extinction. Earn and colleagues (16) developed two distinct criteria that can be used to predict coherent oscillations. The first criterion can determine whether populations can enter coherent oscillations; it incorporates both local dynamics and migration. Technically, the first criterion is calculated as e ^μ|λ|, where μ is the Lyapunov exponent of the local recruitment function and λ is the second leading eigenvalue of the dispersal matrix. Coherent oscillations are possible when this quantity is <1. When coherent oscillations are possible, the second criterion can be used to determine whether they are inevitable. Technically, this criterion is calculated as the product of r, the maximum (i.e., density-independent) local recruitment rate, and λ, the second leading eigenvalue of the migration matrix. If the quantity r|λ| is <1, coherent oscillations are inevitable.

These calculations require a knowledge of migration rates and an estimate of the recruitment function within local populations. The methods by which we obtained these estimates are described in the following section.

Demographic Models. We constructed a demographic model that incorporated the effects of intraspecific competition and migration on metapopulation dynamics. First, we fit a set of nested competition models on the observed dynamics by using maximum likelihood techniques (21). The parameters of the model were fit on data from the isolated treatment because local dynamics in this treatment were not influenced by migration. Second, we determined migration from an independent experiment that directly measured seed dispersal.

Intraspecific competition. The competition function predicted the average local density of a generation as a function f of the local density of the previous generation: f(n_t ) = rn_te ^-bn_t , where r is the growth rate and b is the parameter that determines the strength of competition. We assumed that local densities followed a negative binomial distribution. This distribution is highly flexible, with the Poisson and geometric distributions being two of its limiting cases. We defined the variance of the distribution of density as a function of the mean: var(n_{t +1}) = a(1 + f(n_t ) ^c ). The parameter c allows us to model situations where the variance increases more than linearly with the mean. Simpler models where the variance is constant or varies linearly with the mean were tested against this more complex model. These models all yielded higher values for the Akaike Information Criterion (AIC) and were therefore discarded.

The seeds used to initiate the first generation of our experiment were harvested from plants growing in environmental conditions different from those in our experiment. Thus, we dropped the first generation and fit our demographic models on the resulting truncated dynamics; best fit values obtained from this procedure are r = 2.009 [95% confidence interval (C.I.), 1.646-2.499], b = 0.098 (95% C.I., 0.070-0.122), a = 0.292 (95% C.I., 0.215-0.406), and c = 2.471 (95% C.I., 2.303-2.650). The equilibrium density predicted from these parameter estimates is 7.12. With twice this number, average individual fertility is reduced to <25% of its maximum value, which indicates strong competition.

To use data from all four treatments, we also fit a competition model in which fertility (the number of fruits produced by an individual) is a function of local density; data consisted of the number of plants and number of fruits produced in each generation. Again, average fertility was predicted by an exponential competition model w_{t +} = w ₀ e ^-bn_t , where w ₀ is the fertility in the absence of competition. Fertility was assumed to follow a negative binomial distribution. Additionally, to account for possible changes over time of demographic parameters, we allowed fertility, in the absence of competition, to vary exponentially with time. The model then becomes w_{t +1} = w ₀ e^dte ^-bn_t , where d quantifies the rate of change in fertility over time.

Migration. The migration model was fit on data from an independent experiment where nine individual plants were allowed to disperse their seeds, and the distance between each dispersed seed and the mother plant was measured. We fit an exponential distribution to these data, which allowed us to predict the rate at which dispersal probability decreases with distance. We obtained an average dispersal distance of 22.58 cm (95% C.I., 22.38-22.78 cm). From this analysis, we estimated migration rates between adjacent populations. This analysis revealed large differences among treatments (3.91%, 1.40%, and 0.44% for the continuous, medium, and long treatments, respectively).

Simulations. The competition and migration models were then combined, and simulations were run in which competition occurred as in the isolated treatment but where migration varied. Rates of dispersal between populations were estimated from the fitted dispersal function and standardized to sum to 1. For each generation, we predicted the number of adults that were recruited from a local population by using our probabilistic competition model and then randomly assigned the recruited adults to the four populations of the simulated metapopulation according to the predicted migration rates. Extinction probabilities were estimated from 2,000 simulations that were run for 15 generations.

Results

Synthetic Measures. In our replicate metapopulations, extinction dynamics depended on the amount of fragmentation (Figs. 2 and 3). All four metapopulations of the isolated and long treatments went extinct before the experiment ended. In contrast, two medium metapopulations and one continuous metapopulation persisted until the end of the experiment at 15 generations (Figs. 2 and 3). To assess the statistical significance of these results, we compared the extinction dynamics of the three connected metapopulation treatments to that of the isolated treatment. We used the isolated treatment as a reference because local populations in this treatment have independent dynamics. Using maximum likelihood methods, we fit a negative binomial distribution to the actual distribution of local population persistence times in the isolated treatment. The probability that a local population survives for at least 15 generations is 0.003. We then estimated the distribution of survival times for metapopulations that comprise four local populations that do not exchange migrants. The probability that one or more metapopulations (out of four) survives for 15 generations (as in the continuous treatment) is 0.012, and the probability that two or more metapopulations survive for 15 generations (as in the medium treatment) is 0.00005. Thus, the isolated and long treatments had higher extinction rates than did the medium and the continuous treatments. Globally, persistence times differed among the treatments, with connected metapopulations persisting longer (Kendall rank correlation between persistence and distance between populations, r = -0.49, P = 0.008). Pairwise permutation tests indicated that medium metapopulations persisted longer than long and isolated metapopulations (P = 0.042 and P = 0.037, respectively), but the most connected metapopulations (continuous treatment) did not persist the longest; there was no significant difference between the continuous and long metapopulations (Fig. 3).

Empirically derived estimates of metapopulation capacity indicated that metapopulation capacity increased with decreasing distance between populations (metapopulation capacity = 0.0, 0.1, 0.31, and 0.70 for isolated, long, medium, and continuous treatments, respectively). These values reflect differences in recolonization probabilities: the average number of generations that elapsed between local extinction and subsequent recolonization was 1.8, 1.3, and 1.1 for the long, medium, and continuous treatments, respectively.

It is interesting that metapopulation capacity predicts that the continuous treatment will persist longer than the medium treatment, which is not consistent with our experimental results. This discrepancy is most likely due to how increased migration influences local dynamics. An autoregressive integrated moving average (ARIMA) model indicates a change in local dynamics that correlates with migration rate. First, only 2 of the 16 estimated AR coefficients were significantly <0 (P < 0.05), indicating cycling, and these coefficients occurred in the continuous treatment. In addition, we found that overall AR coefficients were significantly more negative in the continuous treatment than in the other treatments (permutation test, P = 0.04). Additional confirmation comes from the criteria developed by Earn and colleagues (16). We found that coherent oscillations are possible in all treatments but inevitable (i.e., |λ|r drops below 1) only in the continuous treatment. Coherent oscillations can increase extinction risk and may explain why continuous metapopulations had observed and predicted lower persistence levels than did medium metapopulations (16).

Demographic Models. Demographic models fit on the isolated treatment indicated a significant nonlinear density dependence (Fig. 4). The experimentally measured dispersal kernel indicated clear differences in the proportion of seeds dispersing into the three connected metapopulation treatments (Fig. 5). Applying the demographic data and dispersal data in our simulation model, we find a strong and nonlinear effect of migration on extinction risk. Extinction risk suddenly decreases when distance between adjacent populations drops below 50 cm (Fig. 6). We found a qualitatively similar behavior for all of the scenarios we tested, but the magnitude of the drop and the distance at which it occurs varies considerably with the parameter values used. Comparing predicted and observed extinction dynamics then provides a way for us to evaluate a posteriori the validity of our model.

Predicted extinction rates were 74% for the isolated treatment (probability that all four metapopulations go extinct, as observed, is 0.3000), 22% for the long treatment (probability that all four metapopulations go extinct as observed, 0.0024), 10% for the medium treatment (probability that two or more metapopulations go extinct, 0.0523), and 9% for the continuous treatment (probability that three or more metapopulations go extinct, 0.0027). Although we achieved a good fit to observed extinction risk with the isolated treatment, the overall risk of extinction was greatly underestimated in the three remaining treatments. This low prediction of extinction risk could be due to a bias in our estimate of demographic parameters. To test this possibility, we determined the joint 95% C.I. around our estimates of r and b and ran simulations by using values of r and b that lie within this interval. We obtained better predictions of extinction risk with lower values of r. For example, with r = 1.35 and b = 0.066, we obtained extinction risks of 0.76, 0.39, 0.31, and 0.35 for the isolated, long, medium, and continuous treatments, respectively. Note that the correspondence between our predictions and our observations is as good for the isolated treatment and much better for the medium and continuous treatments. This improvement in our predictions indicates that fitting our competition model on the isolated treatment probably overestimated fertility (as quantified by r). The fact that fertility varies across treatments is confirmed by the analysis of a probabilistic model where the number of fruits produced by a plant is predicted as a function of local density. This model can be fit on all four treatments and revealed marginal differences in fertility (Kruskal-Wallis test, P = 0.047), with fertility being lower in the continuous treatment than in the other three treatments.

Surprisingly, taking into account a possible bias in parameter estimation does not resolve the discrepancy between predicted and observed values in the long treatment: whatever estimate of r and b we choose within the 95% C.I. still underestimates extinction risk. Observations on the long treatment suggested that plant fertility (as measured by fruit production) declined over the course of the experiment. To investigate this possibility, we let fertility in the absence of competition vary with time. We then compared the AIC of this more complex model with that of the previous model, where fertility is constant. We found a significant decrease in fertility in the long treatment only (d = -0.06, AIC = 2853.6 for the complex model, AIC = 2858.1 for a model with d set to 0). We then incorporated this decrease into our simulation model and found an extinction probability of 0.72 for the long treatment, which is compatible with our observations. The cause of this decrease remains unclear, particularly because we did not find a similar decline in the isolated treatment.

Discussion

As predicted by classical metapopulation theory (7), migration strongly affected persistence times: connected populations persisted longer than isolated ones. However, the relationship between migration and persistence time was nonlinear. Using a simulation model, we find that there is a sharp increase in extinction risk as distance between local populations increases above a threshold value (≈50 cm) (see Fig. 6). Metapopulation capacity, which incorporates migration but ignores local dynamics, reproduces beautifully this nonlinear relationship (see Fig. 6). Yet our experimental data suggest that extinction risk does not always decrease when migration increases. Rather, in our experiment, the medium metapopulations had the highest persistence levels. This discrepancy between the predicted persistence times and the observed persistence times can be resolved if we examine population coherence (16). Calculation of coherence values for our four treatments indicates that although all fall into the region of parameter space in which coherent oscillations are possible, only the continuous treatment falls into the region where coherent oscillations are inevitable. This result is consistent with our observation that in the continuous treatment, two of the four AR coefficients are significantly negative, and that in general, AR coefficients are more negative in the continuous treatment than in the other three treatments. Thus, synchronous cycles may well occur in this treatment. When metapopulations cycle in phase, extinction risk can increase because all local populations experience bad years at the same time (16-19). Stochastic simulations, where demographic parameters were obtained from the isolated populations and where migration varied, predicted for some parameter values an increase in extinction risk when migration becomes extremely high (see the solid thin curve in Fig. 6). Taken together, these results suggest that the lower persistence times in the continuous treatment result from the interplay of local dynamics and migration. However, caution is needed in attributing different persistence times to the medium and the continuous metapopulations, because the difference in persistence times was not statistically significant, probably because of small sample sizes. Therefore, we have insufficient statistical power to show that coherent oscillations in the continuous metapopulations lowered persistence time. More generally, our simulations indicate that, for our study system, the magnitude of this effect on the short time scale is low. Detecting it would therefore require much more experimental work. From our simulations, we found that such experimental work would not succeed unless >10 replicate metapopulations were followed for at least 50 generations.

In any case, we found convincing experimental evidence that the relationship between migration and persistence is nonlinear and some evidence that high migration rates modify local dynamics. These findings make accurate measurement of migration rates imperative in conservation studies: corridors that would not increase migration above the threshold value would be ineffectual; corridors that would homogenize local densities could increase the risk of coherent oscillations and enhance extinction risk. Determining empirically what migration level is required to improve survival chances will be challenging for any natural population. However, we established that synthetic measures such as metapopulation capacity (10) and population coherence (16) can adequately predict various aspects of the relationship between migration and persistence. Because these methods only require data collected on short time scales (e.g., migration between patches that can be estimated by using population genetics methods, patch capacity that can be assessed by using demographic censuses, and local recruitment), they should provide useful tools for efficiently managing endangered natural populations.