Extinction dynamics in experimental metapopulations
See allHide authors and affiliations

Edited by Simon A. Levin, Princeton University, Princeton, NJ (received for review July 8, 2004)
Abstract
Metapopulation theory provides a framework for understanding population persistence in fragmented landscapes and as such has been widely used in conservation biology to inform management of fragmented populations. However, classical metapopulation theory [Levins, R. (1970) Lect. Notes Math. 2, 75107] ignores metapopulation structure and local population dynamics, both of which may affect extinction dynamics. Here, we investigate metapopulation dynamics in populations that are subject to different migration rates by using experimental metapopulations of the annual plant Cardamine pensylvanica. As predicted by classical metapopulation theory, connected populations persisted longer than did isolated populations, but the relationship between migration and persistence time was nonlinear. Extinction risk sharply increased as the distance between local populations increased above a threshold value that was consistent with stochastic simulations and calculation of metapopulation capacity [Hanski, I. & Ovaskainen, O. (2000) Nature 404, 755758]. In addition, the most connected metapopulations did not have the highest persistence levels. Stochastic simulations indicated an increase in extinction risk with the highest migration rates. Moreover, calculation of population coherence [Earn, D. J. D., Levin, S. A. & Rohani, P. (2000) Science 290, 13601364], a metric that predicts synchronous cycles, indicated that continuous populations should cycle in phase, resulting in an increased extinction risk. Determining empirically the optimal migration level to improve survival chances will be challenging for any natural population. Migration rates that would not increase migration above the threshold value would be ineffectual, but migration rates that would homogenize local densities could increase the risk of coherent oscillations and enhance extinction risk.
Habitat fragmentation is implicated in the decline of many endangered wild species (14). To offset the problems associated with habitat fragmentation, conservation biologists have called for the construction of corridors to facilitate migration between fragmented habitats (5), but the construction of corridors and their attendant increase in migration rates is not without risks (6). Predicting how migration affects population persistence is important for species conservation. Metapopulation theory provides a framework for understanding population persistence in fragmented landscapes and as such has been used in conservation biology to inform management of fragmented populations. Originally conceived by Levins, metapopulation dynamics are determined by the balance between local colonization and extinction events (7). When colonization events exceed extinctions, metapopulations are predicted to persist, because migration allows an extinct local population to become recolonized (79). Absent from this theory is any consideration of how differences in patch suitability (area, distance to other local patches) and local dynamics (such as reproductive rate and intraspecific competition) influence metapopulation persistence. However, recent theoretical work has included more realistic aspects of metapopulations. Hanski and Ovaskainen (10) examine how differences in patch area and in metapopulation structure affect predicted metapopulation persistence (refs. 1114; for a review, see ref. 15). These formulations of metapopulation theory still ignore local population dynamics. Local dynamics may greatly affect metapopulation persistence, because local extinction risk may be a function of local density. For example, intense local competition may induce local population cycling, which may enhance local extinction. Moreover, if the network of local populations becomes highly synchronized because of high migration rates, global extinction risk may increase (1619). Incorporating realistic estimates of both patch structure and local dynamics should therefore enhance our ability to predict metapopulation dynamics.
Here, we report on population persistence in experimental plant metapopulations that were subject to different migration rates. To simulate natural conditions, we allowed the main demographic processes, competition and dispersal, to occur naturally, but we varied the distance between adjacent subpopulations, which indirectly altered migration rate. As migration altered local density, differences in migration rate had the potential to alter local dynamics. This approach allowed us to determine the interplay between migration and local competition in a naturalized setting.
Materials and Methods
Study Species. Experimental populations were initiated with the annual plant Cardamine pensylvanica (Brassicacea), a selffertile weed that has no seed dormancy, produces explosively dispersed seeds, and can go from seed to seed in ≈2.5 months in controlled conditions (20). Because reproduction and dispersal occurred naturally after populations were initiated, the experimental plants completed their life cycle without further manipulation; thus, the role of local competition and dispersal in population persistence could be directly assessed.
Experimental Design. The experimental arrangement consisted of a set of linked local populations [i.e., metapopulations sensu Levins (7)] where the degree of connection between the local populations varied from completely connected to completely isolated (Fig. 1). In our experiment, local populations were adjacent (continuous treatment), separated by 23.2 cm (medium treatment), separated by 49.5 cm (long treatment), and separated by partitions so that there was no seed dispersal between populations (isolated treatment) (Fig. 1). Local populations consisted of a pair of pots (each 18 cm deep × 4 cm in diameter, filled with vermiculite), one of which contained bare soil. Pots were set up in an alternating zigzag pattern, a bare soil pot (open circle in Fig. 1) being placed next to each local population (hatched circle in Fig. 1). The bare soil pot provided a receptacle for dispersing seeds, thus constituting the next generation and providing a way to maintain populations continuously over multiple generations. Each metapopulation was composed of four local populations and therefore consisted of a narrow array of eight pots (two pots × four local populations for the entire metapopulation). Once 80% of the plants in all of the metapopulations had dehisced their seeds, all pots harboring plants were simultaneously removed (to impose consistent discrete generations) and were replaced with bare soil pots. We then counted all individual plants within each pot to determine the number in each local population. We repeated this procedure for 15 generations or until all four local populations within a metapopulation went extinct. The generation time was ≈10 weeks, so the experiment was run for ≈3 years. At generation 0, each local population was initiated with one seed taken from >100 individuals from each of three longestablished greenhouse populations of C. pensylvanica. The experimental metapopulations were maintained in growth chambers (1,800 microeinsteins·m^{2}·s^{1}, 24°C, 18 h daylight, watered four times daily, nutrients applied once daily) at the University of Vermont.
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 firstdegree 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 nonnormally 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., densityindependent) 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_{t}e ^{bnt }, 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.6462.499], b = 0.098 (95% C.I., 0.0700.122), a = 0.292 (95% C.I., 0.2150.406), and c = 2.471 (95% C.I., 2.3032.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 _{0} e ^{bnt }, where w _{0} 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 _{0} e^{dt}e ^{bnt }, 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.3822.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 (KruskalWallis 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 (1619). 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.
Acknowledgments
We thank Drs. Brian Beckage, Neil Buckley, Sebastian Lavergne, Ophelie Ronce, Peter Thrall, and two anonymous referees for reviewing the manuscript. This work was supported by National Science Foundation Grant DEB9527986 (to J.M.).
Footnotes

↵ * To whom correspondence should be addressed. Email: jane.molofsky{at}uvm.edu.

↵ † Present address: Laboratoire Génome, Populations, Interactions, Centre National de la Recherche Scientifique Unité; Mixte de Recherche 5171, Bâtiment 24, Université; de Montpellier II, 34095 Montpellier Cedex 5, France.

Author contributions: J.M. designed research; J.M. performed research; J.B.F. analyzed data; and J.M. and J.B.F. wrote the paper.

This paper was submitted directly (Track II) to the PNAS office.

Abbreviations: AIC, Akaike Information Criterion; AR, autoregressive; C.I., confidence interval.
 Copyright © 2005, The National Academy of Sciences
References

↵
Bierregaard, R., Lovejoy, T., Kapos, V., Dossantos, A. & Hutchings, R. W. (1992) Bioscience 42 , 859866.

Quinn, J. F. & Hastings, A. (1987) Conserv. Biol. 1 , 198208.

Robinson, S. K., Thompson, F. R., Donavon, T. M., Whitehead, D. R. & Faaborg, J. (1995) Science 267 , 19871990.
 ↵

↵
Simberloff, D. & Cox, J. (1987) Conserv. Biol. 1 , 6371.

↵
Simberloff, D., Farr, J. A., Cox, J. & Mehlman, D. W. (1992) Conserv. Biol. 6 , 493504.

↵
Levins, R. (1970) Lect. Notes Math. 2 , 75107.
 ↵
 ↵
 ↵
 ↵
 ↵

↵
Earn, D. J. D., Levin, S. A. & Rohani, P. (2000) Science 290 , 13601364. pmid:11082064

Holt, R. D. (1993) in Species Diversity in Ecological Communities, eds. Ricklefs, R. & Schluter, D. (Univ. Chicago Press, Chicago), pp. 7788.
 ↵

↵
AlShebaz, I. A. (1988) J. Arnold Arbor. 69 , 85166.

↵
Hilborn, R. & Mangel, M. (1997) The Ecological Detective: Confronting Models with Data (Princeton Univ. Press, Princeton).