Treatment-based Markov chain models clarify mechanisms of invasion in an invaded grassland community

Lisa Castillo Nelis, J. Timothy Wootton


What are the relative roles of mechanisms underlying plant responses in grassland communities invaded by both plants and mammals? What type of community can we expect in the future given current or novel conditions? We address these questions by comparing Markov chain community models among treatments from a field experiment on invasive species on Robinson Crusoe Island, Chile. Because of seed dispersal, grazing and disturbance, we predicted that the exotic European rabbit (Oryctolagus cuniculus) facilitates epizoochorous exotic plants (plants with seeds that stick to the skin an animal) at the expense of native plants. To test our hypothesis, we crossed rabbit exclosure treatments with disturbance treatments, and sampled the plant community in permanent plots over 3 years. We then estimated Markov chain model transition probabilities and found significant differences among treatments. As hypothesized, this modelling revealed that exotic plants survive better in disturbed areas, while natives prefer no rabbits or disturbance. Surprisingly, rabbits negatively affect epizoochorous plants. Markov chain dynamics indicate that an overall replacement of native plants by exotic plants is underway. Using a treatment-based approach to multi-species Markov chain models allowed us to examine the changes in the importance of mechanisms in response to experimental impacts on communities.

1. Introduction

Most ecological systems are dynamic; and yet we usually study them with static snapshots of population abundance. While these methods are good for examining overall short-term trends, to make long-term predictions and understand how important processes respond to environmental impacts, it may be helpful to explore changes in underlying dynamics. Multi-species Markov chain models (hereafter referred to as Markov models) provide one way to characterize ecological dynamics and can be estimated relatively easily from empirical data. Here, we explore ways in which applying such approaches to studying dynamics may improve insight into how communities and species respond to environmental impacts such as species invasion.

Markov models characterize dynamic relationships among system components by modelling transitions in ecological states over time. The probability of each transition to the next state given the current state is independent of past states. For example, given that the current state is bare space (i.e. no plants), the Markov chain gives the probability of transitioning into all other states in the next time step. Markov models can use data at fixed points in space over time, are relatively easy to parametrize and some studies have shown that they accurately predict composition in both the systems they were parametrized in and in novel situations (Wootton 2001b, 2004).

Waggoner & Stephens (1970) and Horn (1975) began using multi-species Markov models to describe the successional dynamics of forest communities in the 1970s. Since then, others have expanded the use of Markov models to other systems, including coral reefs, grasslands, intertidal communities and behavioural studies (Austin 1980, Augustin et al. 2001, Caswell 2001; Wootton 2001a,b, 2004; Hill et al. 2002; Jang et al. 2008). Recently, the use of Markov models has been expanded by the rigorous experimental testing of their prediction abilities in intertidal systems (Wootton 2001b, 2004) and by the creation of indices that allow comparison between or within systems (Wootton 2001a; Hill et al. 2004). Additionally, Markov models take into account many underlying interactions, including indirect interactions, such as apparent competition, that are frequently left unexplored in community analysis (White et al. 2006). Here, we further apply Markov chain models to explore the effects of treatments on the underlying dynamics of a multiply-invaded grassland community to predict the long-term effects of treatment on community composition, and to use treatments to predict long-term outcomes in novel situations. Hence, our analysis extends prior approaches comparing different empirically derived transition matrices of life stages in a single species (Caswell & Werner 1978; Levin et al. 1987; Harvell et al. 1990; Walls et al. 1991), examining effects of perturbations to empirically derived transition matrices of life stages of single species (reviewed in Caswell 2001; Morris & Doak 2002) and contrasting transition matrices among different communities (Tanner et al. 1994; Wootton 2001a; Hill et al. 2004; Wootton et al. 2008).

Historical records of invasions often reveal negative impacts of invaders on native species (Schofield 1989; Clarke et al. 2005; Telfer et al. 2005). When multiple simultaneous invasions occur, the possibility of facilitation or even ‘invasional meltdown’ is present, especially when an exotic organism and its exotic pollinator or dispersal agent co-occur (Simberloff & Holle 1999; Parker & Haubensak 2002; Ness & Bronstein 2004; Wonham et al. 2005; Gurevitch 2006; Simberloff 2006).

We understand relatively little about the mechanisms affecting multiple interacting invasions. What is certain is that the interactions resulting from these mechanisms are likely to be extremely complex. For example, in the situation of an exotic herbivore and exotic plants co-occurring in a native grassland, the herbivore-mediated mechanisms involved could include disturbance, selective herbivory and epizoochorous seed dispersal (seeds dispersing by sticking to the skin of an animal vector). These mechanisms can in turn affect interactions between invasive and native plants (Edwards et al. 2000; Parker et al. 2006; Bailey et al. 2007; Peters 2007).

To explore the mechanisms among the exotic herbivore and native and exotic plants further, Markov chain transition matrices can be calculated. Even though transition matrices do not explicitly define mechanisms, they reflect mechanisms (Wootton 2001b). For example, transitioning into bare space can reflect disturbance and/or herbivory, transitioning from bare space can reflect colonization, dispersal and recruitment limitation and transitioning between species reflects direct and apparent competition. When these transition matrices are compared among treatments in conjunction with knowledge of the natural history of the system, we have a powerful tool for analysing underlying mechanisms and understanding any given community in more depth than is possible from simply examining overall responses. This depth of understanding can aid in our general understanding of community dynamics, and when studying at-risk systems like those with multiple invasions, can even be used to aid conservation and restoration.

Here, we explore how the use of Markov chain models can aid understanding, long-term community prediction and information for conservation. The study system is a mixed native and exotic grassland with no native mammals, but an introduced mammalian herbivore.

2. Study system

This study took place in Vaquería valley, Robinson Crusoe Island, Chile. The island (33°37′ S, 78°53′ W) is 667 km west of mainland Chile (Tulke 1979). Robinson Crusoe is part of the Juan Fernández Archipelago and has an area of 93 km2 (González-Ferrán 1987). The Juan Fernández Archipelago has high plant endemism (Greimler et al. 2002a), but has no native mammals, reptiles or amphibians (Skottsberg 1920).

Since human discovery in 1574, the plant composition of the islands has changed dramatically (Haberle 2003). Its native grassland areas have been invaded by exotic species, including rabbits (Zunino 1989; Jaksic 1998) and many plants, some of which are epizoochorous (Greimler et al. 2002b; Danton 2004). By introducing exotic plants and mammals, humans have created a new biological community with potentially interacting invasive organisms. Associated with these invasions, native plant species richness has declined (P. Danton 2007, personal communication).

The Vaquería valley, our study site, is composed of gully-separated patches of grassland, predominately composed of Nassella laevissima (Greimler et al. 2002a), a native bunchgrass. In 1987, the valley was fenced to keep out cattle, goats and horses as part of a larger programme to restrict grazing of large exotic mammals (I. Leiva 2007, personal communication). Among the small mammals present in the valley (rats, mice, cats and rabbits), rabbits appear to cause the most impact to vegetation at our study site (L. C. Nelis 2003, personal observation) and were therefore selected as the focal herbivore in our study.

Based on characteristics of our study system, we generated the following a priori hypotheses for the pattern of effects we expected with the manipulation of rabbits.

  1. Exotic plants with epizoochorous seeds will decrease where rabbits were excluded because of decreased seed dispersal.

  2. Native plants will increase in the absence of rabbit grazing because the community has developed in the absence of mammalian herbivores.

  3. Exotic plants will decrease in the absence of grazing because there will be no selective herbivory by rabbits on natives to reduce competition.

  4. Soil disturbance will favour exotic plants because of selective grazing on native seedlings, susceptibility of natives to disturbed conditions owing to their lack of shared history with digging mammals and higher dispersal ability typical of exotic species.

  5. Rabbits will have stronger effects on disturbed sites because herbivory causes higher mortality rates on seedlings and because of the effects of differential seed dispersal where competition from residents is weak.

We examined these hypotheses using both Markov chain models and traditional population-level analyses. The experimental system was composed of rabbit presence and absence treatments crossed with disturbed and not disturbed treatments. Each treatment was modelled independently using Markov models to compare the underlying dynamics among treatments.

3. Material and Methods

(a) Data collection and model parametrization

In October 2004, we constructed eight replicated experimental rabbit-exclusion blocks. Each block was designed with treatments at two levels: half-block and subplot. The half-block treatment was rabbit presence or absence. Rabbit absence half-blocks were fenced with a metre high chicken wire buried 20 cm deep. Rabbit presence half-blocks had no fence but had posts and were disturbed to imitate the burying of a fence. Each rabbit presence or absence half-block was then subdivided into four 10 × 10 m plots. Each plot was subsampled by 12 0.5 × 0.5 m permanent subplots in a stratified arrangement (Bergelson et al. 1993), half of which were disturbed by having the soil turned over in 2004, and half of which were left in the state in which we found them. This layering of half-block and subplot treatments allowed us to examine an experimental system of rabbit presence crossed with disturbance with four treatments: rabbit disturbed, rabbit not disturbed, no rabbit disturbed and no rabbit not disturbed. There were a total of 768 subplots, or 192 per treatment.

We created the disturbed treatments by turning over the soil to a depth of 10 cm within the 0.5 × 0.5 m subplots using shovels. While this type of disturbance most closely resembles human agricultural disturbances, it does resemble rabbit digging in that a gap in established plants is created. The gap size experimentally created, 0.5 × 0.5 m, is intermediate between the two types of digging disturbances rabbits create: tunnel construction and small excavations. Tunnel construction results in a hole with varying depth surrounded by an area cleared of plants and covered with loose soil from the excavation; a plant gap up to 1 m2 in area results. Small excavations tend to be triangular in shape, 2–3 cm deep at the point of the triangle and result in 10–20 cm2 free of plants. While not mirroring rabbit disturbances exactly, the disturbances we created allowed us to monitor results of competition among plant species for colonization of gaps, and subsequently allowed us to follow the plant community recovery after an episodic soil disturbance.

We collected data from each permanent subplot using a 0.5 × 0.5 m quadrat divided into 25 10 × 10 cm squares. Censuses was taken during 2004–2006 in December, when most plants were producing seeds but had not yet begun to senesce. These permanent subplots were marked with two numbered flags placed into the ground to orient our portable quadrats. Within each quadrat square, we noted each plant species present. We chose this method instead of obtaining count data for three reasons: (i) growth forms of some plants (e.g. bunch grass) make counting stems less valuable, (ii) we can use the number of squares in which a species was present within a subplot as an estimated count, (iii) it is time efficient and allows us to sample a higher percentage of the grassland for rare plants without losing spatial definition. To check that we had sampled enough subplots to find the majority of plant species, we used data from our 2004 field season to create species accumulation curves by plots (Colwell 2005). As we found the curves flattened with increased sampling showing that we had sampled sufficiently, we continued these methods in subsequent years.

(b) Natural history

To test our hypothesis that rabbits preferentially carry exotic seeds, we examined six hunted rabbits borrowed from hunters in 2004 in the Vaquería valley. The stomachs were cut open and the contents were visually examined for seeds. The fur was also carefully checked for seeds by using a fine-toothed comb.

To further examine the possibility of epizoochorous seed transport, we stuffed two commercially purchased rabbit pelts with rocks and newspaper and dragged them on roughly oval transects around each experimental block in 2005. Every 20 m we stopped, counted, identified and collected seeds combed out of the rabbit pelts, and then continued the transect. Overall, 893 m of transects were conducted across the valley.

To determine the input from the seed bank, we also germinated soil samples. In 2005, soil samples were collected outside of experimental blocks. The samples were taken to a depth of 50 cm. The cores were then divided into 10 cm depth samples, and these samples were placed into Petri dishes in a germinator (18–20°C, 12 h of light, medium humidity) to grow so plants could be identified.

(c) Model analysis

To use the Markov chain approach, quadrat squares were assigned to mutually exclusive ecological states. We selected states with both biological interest and sufficient observations of changes to parametrize the model. Most of our 57 identified plants had too few occurrences to be parametrized as individual states. Conversely, some individual plants had so many transitions that any grouped ecological state they were in reflected only their own transitions. Therefore, after exploring the general structure of the dataset, the most ecologically informative grouping with sufficient sample sizes that we found is shown in table 1. These definitions are the rules by which plants were divided into ecological states.

View this table:
Table 1.

Ecological states for Markov chain analysis.

Once we had determined ecological states, we analysed the final census (2006) data using split-plot MANOVA (JMP 5.1), because group proportions are not independent and because rabbit treatments were split at the half-block level necessitating a split-plot analysis. Each quadrat square was identified as only one state, following the same rules established for Markov states (table 1). The dependent variables were estimated count data for each ecological state transformed by (ln(data) + 1). Estimated counts were the total number of times one or more plants from a given ecological state were present within a subplot. Fixed factors estimated were rabbit, disturbance and disturbance * rabbit. Because multiple disturbance treatments were located within single-rabbit treatments, we used a split-plot design in which we estimated the random effect half-block[rabbit] and used this as the error term for testing rabbit effects, and estimated disturbance * half-block[rabbit] and used this error term to test for disturbance effects.

For the Markov Chain analysis, the analysis level used is a quadrat square; 10 × 10 cm squares within a subplot. There were 4800 quadrat squares sampled each year for each treatment. Once each quadrat square was assigned to a mutually exclusive ecological state for each year, we separated our data by treatment. We then calculated the transition frequencies between each ecological state within each treatment over our three censuses, giving us four transition matrices of our rabbit crossed by disturbance treatments pooled over two time steps. To probe the possible effects of non-independence of quadrat squares from localized interactions or the occasional individuals that were larger than a single quadrat square, we also analysed only non-adjacent quadrat squares (nine of 25 quadrat squares analysed per subplot). We calculated transition frequencies for each treatment as described above and ran a regression between the full and reduced data transition frequencies.

To determine whether transition patterns of our four treatments were significantly different, we used a log-linear model analysis (Sokal & Rohlf 1995; Caswell 2001; Hill et al. 2002). We fit hierarchical models accounting for ecological state at the prior census, rabbit treatment, disturbance treatment and all interactions among these terms on ecological fate of a site, using a fully saturated model and models lacking combinations of interaction terms among rabbit treatment, disturbance treatment and ecological state (Hill et al. 2004). To assess statistical significance of these effects, we calculated maximum-likelihood G-statistics for each model and tested the difference between the statistics for models with and without the terms of interest, using a χ2 distribution and degrees of freedom equal to the difference in degrees of freedom of the nested pairs of models (Microsoft Excel 2004 for Mac; Hill et al. 2004).

We analysed the data for time-varying disturbance and rabbit effects by separating transitions by the starting year of their census interval and using log-linear analysis with terms for rabbit treatment, disturbance treatment, start year and their interactions, compared with a model without interactions between time and treatments. We also examined correlations among transitions from different census periods. Because all disturbed plots initially started as bare space, we could not probe time-varying transitions from other starting states in a complete factorial design, so we targeted subsets of transitions. These included comparison of the full matrix of transitions in undisturbed plots to explore time-varying rabbit effects, and comparison of transitions starting from bare space for each rabbit/disturbance treatment combination.

We examined large-scale patterns of change in the magnitudes of transitions among treatments, using coloured versions of our matrices to aid in visualization. We also calculated most of the indices proposed in Hill et al. (2004), though here we focus only on two of the most relevant to our questions: replacement by a species, and recurrence times. Replacement by a species, or replacement ability, is the probability pij that ecological state i can replace other ecological states j in one time step. To calculate replacement ability, we usedEmbedded Image where b is bare space, and s is the number of plant states (s = 8). Note that species replacement can happen directly through competitive displacement, or indirectly if the first plant dies and a second plant colonizes the same space (Tanner et al. 1994).

Recurrence times are also informative about our system. Recurrence time is the average time elapsed between a point leaving state i and returning to it (Iosifescu 1980; Hill et al. 2004), and is calculated asEmbedded Image where wi is the ith element of the eigenvector associated with the dominant eigenvalue standardized to sum one. For the mean biological recurrence times in a stationary community,Embedded Image

For the eigenanalysis, we calculated the eigenvector associated with the dominant eigenvalue for each treatment matrix and standardized the eigenvector of each treatment to sum to one. Besides their use in the above recurrence time analysis, these standardized eigenvectors are the prediction for the long-term plant composition of the ecological community, known as the stationary distribution (Caswell 2001).

Additionally, it is often of interest for restoration to know what the original community prior to invasion was like, but we usually lack this information. Here, we made predictions about the composition of the original community by modifying our models to examine (i) the effect of rabbits on native plants in the absence of exotic plants and (ii) the composition of the system in the absence of all exotics (Wootton 2001b, 2004). To do this, we eliminated all squares involving exotic plants and/or rabbits and recalculated our transition probability matrices, which is mathematically equivalent to removing columns and rows of the transition probability matrix and proportionally reassigning the probabilities of removed elements to the remaining elements in each column. Using the new probability matrix without exotics, we calculated the stationary distribution as above, giving us the long-term composition predictions without exotic species. For this analysis we combined disturbed and undisturbed states into one pooled matrix, so we could examine only the effect of rabbits.

Comparing the results from the unmanipulated matrix with the matrix with exotic plants removed provides a picture of the effects of rabbits on native plants, whereas examining the rabbit-exclusion matrix with exotic plant transitions removed provides a picture of the pre-invasion community. We are assuming that removal of some species from the model does not disproportionately change the transitions among the remaining species. This assumption has been upheld in intertidal and lake systems (Wootton 2001b, 2004; 2009, unpublished data). We therefore extrapolate, but cannot prove, that it is valid in our grassland system.

4. Results

Plants identified over the 3 years of this study (electronic supplementary material, appendix S1) were grouped into ecological states that were mutually exclusive for the Markov study (table 1). The MANOVA (electronic supplementary material, appendix S2) was significant for rabbit (p < 0.0001), disturbance (p < 0.0001) and disturbance * rabbit (p < 0.0001). Follow-up analyses with ANOVA (electronic supplementary material, appendix S2) revealed significant effects of disturbance (figure 1); native plants (Nassella, Native, Juncus) did worse with disturbance, and exotic plants (sticky, non-sticky and Briza) did better. Bare space, sticky, Briza, Juncus and Nassella were all significant for disturbance * rabbit (figure 1). Curiously, no significant effect of rabbit treatment was detected on a species by species case, perhaps because compensatory responses in plants are not accounted for in ANOVA but are in MANOVA (Sokal & Rohlf 1995), and because the experimental design had less power to detect rabbit effects.

Figure 1.

Mean 2006 plant abundance by treatment. Nassella, Native and Juncus are native plants, the other groups contain exotics and mix is a mixture of native and exotic plants. The statistical analyses were conducted on ln(data) + 1, therefore the means shown are exp((ln(data) + 1)) − 1. **p < 0.05 (disturbance); *p < 0.10 (disturbance); ••p < 0.05 (disturbance * rabbit).

(a) Patterns of transitions

Patterns of transitions in ecological states were significantly affected by rabbit treatment, disturbance treatment and all their interactions (table 2; electronic supplementary material, appendix S3). The overall species effects confirmed that different species exhibit different life histories and ecologies. Terms involving disturbance and rabbit treatments demonstrated that treatments affected the patterns of ecological state replacements, that effects of rabbits depended on whether disturbance was present and that many treatment responses were species-specific. The regression of the transition frequencies calculated with all quadrat squares compared with transition frequencies calculated from reduced data (non-adjacent quadrat squares) showed a strong correlation (slope = 1.01 (±0.01 s.e.), intercept = −0.001 (±0.002 s.e.), adjusted r2 = 0.994); therefore, any non-independence arising from local interactions did not appear to strongly affect the analysis.

View this table:
Table 2.

Results of log-linear model analysis testing different treatment effects on transitions in ecological state from one census to the next. Model codes represent the highest levels of interactions among variables (S, initial state; F, subsequent fate; R, rabbit treatment, D, disturbance treatment) present in the model.

Rabbit effects varied with time in the no-disturbance plots (G = 535.14, d.f. = 112, p < 0.0001), but these differences had little impact on predicted long-term community composition (compositional similarity between years = 0.90 without rabbits, 0.92 with rabbits). Disturbance effects were smaller in the second time step than the first: for plots starting as bare space, the correlation between transitions in disturbed and non-disturbed plots increased from 0.745 in 2004–2005 to 0.892 in 2005–2006. However, disturbance as a factor was still significant when only 2005–2006 data were analysed (G = 64.08, d.f. = 14, p < 0.0001), indicating effects of disturbance lasted more than 1 year.

We examined the treatment-specific matrices for patterns of differences among transitions using colour-coded matrices to aid in visually interpreting patterns of matrix structure and changes with treatment (figure 2). The largest values across all matrices were generally associated with self–self transitions (matrix diagonal), transitions to exotics (lower half of the matrix) and transitions to the mixed native and exotic state (bottom row), but other specific transitions had high values in the different treatments. For example, transitions to native states (Nassella, Native and Juncus) were highest in the absence of rabbits and disturbance, as we expected. In contrast to our expectations and to assumptions underlying current management plans, overall transition probabilities to natives tended to decline when rabbits were removed in disturbed plots, suggesting that simply removing rabbits may not benefit native species.

Figure 2.

Coloured-coded transition matrices where column transitions to row. More intense colour indicates stronger transitions. Nassella, native, and Juncus are native plants, sticky, non-sticky, and Briza are exotics, and mix is a mixture. Treatments shown are (a) rabbit disturbed, (b) rabbit not disturbed, (c) no rabbit disturbed and (d) no rabbit not disturbed.

Summary matrix indices (electronic supplementary material, appendix S4) showed that the mixed native–exotic ecological state tended to maintain itself over time (figure 3a) and replaced other states (figure 3b). Exotic states replaced others more frequently than native states (figure 3b), with highest probability in disturbed sites. In undisturbed plots, the native Nassella also has a high replacement ability (figure 3b). Within rabbit treatments, biotic mean recurrence times show that all native plant groups recurred faster (i.e. they did better) when undisturbed, while all exotic species recurred faster with disturbance (figure 3c).

Figure 3.

Indices to compare transition matrices. (a) A measure of competition, the winner from a native–exotic mix. (b) Average replacement ability of a functional group. (c) Recurrence time at a site following disappearance for each functional group.

Transitions from and to bare space, which reflect the underlying mechanisms of colonization and disturbance, respectively, responded to herbivory and disturbance (figure 4). In undisturbed plots, rabbits increased transitions from native species to bare space as expected, but tended to reduce transitions from exotic species to bare space. In the presence of disturbance, however, rabbits strongly increased transitions to bare space from exotics and the native species Nassella, but reduced transitions to bare space by natives. Contrary to our hypothesis, rabbits decreased colonization of bare space by epizoochorous, or sticky-seeded, plants in both the presence and absence of disturbance (figure 4). Surprisingly, rabbit presence was extremely beneficial for the colonization of two native species: Juncus colonization of bare space benefited from rabbits in both disturbed and undisturbed plots, and Nassella colonization increased in disturbed sites when rabbits were present.

Figure 4.

Rabbit effects shown as proportional transition change (the difference between transition probabilities over the mean) on (a) transition to bare space and (b) colonization.

The predictions of community composition, or stationary distributions, were constructed for both current conditions and novel conditions (figure 5). When the rabbit impact under current conditions is used to make the prediction, we find that rabbits reduce the abundance of rare natives, increase Briza, an exotic species, and increase Juncus, a native species (figure 5a). Projected effects of rabbits on bare space and Nassella vary with disturbance, increasing them in disturbed plots but having minor effects in undisturbed plots. These predictions, made using data from 2004 to 2006, strongly correlate with the final measured composition in 2006 (r = 0.99).

Figure 5.

Prediction by Markov chain models of community composition in (a) the current parametrized situation with both rabbits and exotic plants, (b) the current parametrized situation and novel situations without exotic plants and (c) the current parametrized situation and in the novel situation with no exotic plants or rabbits.

Removing exotic plants from the Markov chain models with and without rabbits and projecting long-term composition indicated that the presence of exotic plants weakens the effect of rabbits on native plants (figure 5b). Bare space, however, showed an opposite trend for the effect of rabbits in the presence versus the absence of exotic plants. Note that disturbance treatments were not separated for this analysis.

Most useful for conservation, when all invasive organisms were removed from the prediction, we are able to model the pre-invasion grassland composition and compare it with the current situation (figure 5c). As expected, all native groups probably exhibited higher abundance in the absence of exotics, with Nassella dominating the community.

(b) Natural history

We found that European rabbits do carry exotic seeds epizoochorously. Of the six rabbits we were able to examine, none had seeds in their guts, only grass clippings. On their pelts, five rabbits carried seeds from exotic species (four, Acanea argentea, one, Avena barbata), and no rabbits carried native seeds. The number of seeds ranged from two to 164. We also found that exotic seeds attached more frequently than native seeds to stuffed rabbit pelts even though the native N. laevissima had the highest biomass in the valley. The rabbit pelt transects picked up 693 seeds from nine exotic species (dominant were 383, A. argentea and 253, Galium aparine) and 100 seeds from three native species (dominant was 93, N. laevissima).

There was no evidence that the seed bank affects plant recruitment in this system. Of the close to 100 soil samples we attempted to germinate, only a few seeds were found, only three samples germinated plants and they were all from 0 to 10 cm depth. Two of the plants were N. laevissima, the dominant native bunch grass, and one was a dicot that died before it produced a real leaf.

5. Discussion

Our results support only two of our original five hypotheses: that disturbance favours exotic plants, and that rabbit effects are the strongest in disturbed areas. Exotic plants with epizoochorous seeds did not increase in the presence of rabbits, apparently because the effect of rabbit grazing was a more important factor than seed dispersal limitation. Other mechanisms besides rabbit dispersal may be adequate dispersal agents for these epizoochorous plants. Native plants did not uniformly increase and exotic plants did not uniformly decrease when rabbits were removed because rabbits did not specialize on evolutionarily naive natives. Perhaps the shared evolutionary history between rabbits and some of the exotic species resulted in strong adaptations by rabbits to handle plant defences (Nelis 2008).

Examining only abundance patterns to study this system, we learned that there is an effect of rabbit presence across plant groups. However, individual plant groups were more strongly affected by disturbance and the interaction of disturbance with rabbit presence than by only rabbit presence. Although we know rabbits do carry exotic seeds on their coats, sticky exotic plants did worse in the presence of rabbits. While this is useful information, by adding the use of treatment-based Markov chain models, our understanding of the community and the multiple ongoing interactions becomes much deeper.

For example, our use of treatment-based Markov chain models shows a pattern of local short-term replacement of natives by exotics in all treatments. Transition probabilities in all treatments tend to be the highest entering a mixed state, and when a mixed state transitions to another group, it tends to transition to exotic species. Furthermore, exotic species compete much better than natives in disturbed areas, and rabbits continually propagate disturbance with their digging. Even though it seems that rabbits help Nassella and Juncus colonize bare space (figure 4b), these native species do not seem to persist well in the long term, especially in disturbed areas.

In addition to examining the mechanisms in which we were interested a priori, we also discovered some unanticipated interactions. For example, we had not expected that the removal of rabbits in the disturbed habitat would decrease native recruitment to such an extent, nor is this pattern clear from looking at mean ecological state abundances. We had also not suspected that Briza was so dependent on rabbit presence that the long-term prediction would show a strong positive affect of rabbits across disturbance treatments. This ability of Markov models to point out overlooked interactions and community patterns is one of this method's greatest strengths, and in this case gave us important conservation information.

Our prediction of community composition shows that to restore the study system to a more native state, management targeted at eliminating rabbits will be insufficient because even without rabbits exotic species out-compete native species in previously disturbed habitats. The effects of disturbance persist for an extended period of time; there was little recovery of native species in disturbed habitats over 2 years, the effects of disturbance persisted through the second transition time step and calculated recurrence times for native plants in disturbed areas are very high. In addition to eliminating rabbits, the park must focus on reseeding and exotic plant removal in previously disturbed habitats where natives do not recruit strongly on their own. While Markov models may not exactly predict future community composition in this system, they have been shown to be accurate in intertidal systems (Wootton 1994, 2001b) and can be used to educate management decisions.

Also, the modified model in which we eliminated exotic plants and rabbits from the long-term composition predicts that invasive species as a group have strongly reduced all native species, and that this reduction has been proportionately higher for currently rare species than for Nassella. While removing both rabbits and exotic plants from the model cannot exactly reverse-engineer the original conditions, the composition predictions give restoration work a reasonable starting goal. Further analyses to assist in conservation could be done, such as finding the sensitivities of factors to minimize the difference between the estimate of the original conditions and the current stationary distribution.

With the insight into mechanisms that we have gained by using Markov chain models to examine this multiply-invaded system, we can now design more focused experiments or analyses to tease apart the importance of specific mechanisms. For example, we know that rabbits positively impact the colonization of bare space by the native species Juncus in both disturbed and undisturbed conditions. Whether this pattern is caused because Juncus is not a preferred food item and is indirectly positively impacted by rabbits, because rabbits somehow improve its dispersal or because Juncus is able to strongly benefit from nutrient deposition by rabbits could be easily determined now that the pattern has been discovered.

In conclusion, Markov chain models were extremely useful for examining the underlying dynamics in our experimental treatment-based system using spatially explicit data over multiple years. The models were relatively easy to parametrize, they accurately delivered a wealth of information about community dynamics in a multiply-invaded system, there were readily available tools such as indices for further analysis and Markov chain models allowed us to make predictions in current and novel situations that helped inform conservation at our field site. Further studies using a similar treatment-based approach of Markov chain models are called for to explore differences among the underlying dynamics of different communities.


We would like to thank J. M. Fariña, without whose help this project would not have been possible. We are also grateful to the Corporación Nacional Forestal de Chile for their support from both Región V and the Parque Nacional Archipiélago Juan Fernández. Thanks to G. Dwyer, P. Amarasekare, J. Bergelson and C. Pfister for their scientific and editorial help throughout this project, to C. Johnson, M. Polk, J. Meza, I. Leiva and V. Reyna for logistical support, to our army of field and laboratory assistants who are too numerous to name, to an anonymous reviewer, M. Spencer, and H. Caswell who immeasurably improved this manuscript, and to all of those friends and colleagues who contributed their support to the authors. This project was funded by the National Science Foundation GRF and NSF grants 011780, 10456110, 0452687 and 0708462, by the ARCS Foundation Chicago, by the University of Chicago Hinds Fund, by Sigma Xi, by Rotary One, by the Southeast Chicago Rotary Club and by personal funds.


  • Present address: Department of Biology, Stanford University, Gilbert Building, Room 109, 371 Serra Mall, Stanford, CA 94305-5020, USA.

    • Received August 28, 2009.
    • Accepted October 6, 2009.


View Abstract