Models assessing the prospects of plant species at the landscape level often focus primarily on the relationship between species dynamics and landscape structure. However, the short-term prospects of species with slow responses to landscape changes depend on the factors affecting local population dynamics. In this study it is hypothesized that large herbivores may be a major factor affecting the short-term prospects of slow-responding species in the European landscape, because large herbivores have increased in number in this region in recent decades and can strongly influence local population dynamics.
The impact of browsing by large herbivores was simulated on the landscape-level dynamics of the dry grassland perennial polycarpic herb Scorzonera hispanica. A dynamic, spatially explicit model was used that incorporated information on the location of patches suitable for S. hispanica, local population dynamics (matrices including the impact of large herbivores), initial population sizes and dispersal rate of the species. Simulations were performed relating to the prospects of S. hispanica over the next 30 years under different rates of herbivory (browsing intensity) and varying frequencies of population destruction (e.g. by human activity).
Although a high rate of herbivory was detected in most populations of S. hispanica, current landscape-level dynamics of S. hispanica were approximately in equilibrium. A decline or increase of over 20 % in the herbivory rate promoted rapid expansion or decline of S. hispanica, respectively. This effect was much stronger in the presence of population destruction.
Browsing by large herbivores can have a dramatic effect on the landscape dynamics of plant species. Changes in the density of large herbivores and the probability of population destruction should be incorporated into models predicting species abundance and distribution.
Rapid changes in the landscape in recent years have resulted in increased levels of habitat fragmentation for many plant species. These changes have provoked discussion about the prospects of such species in the future agricultural landscape (e.g. Saunders et al., 1991; Collinge, 1996; Bastin and Thomas, 1999; Lindborg and Eriksson, 2004). Several authors have emphasized the importance of describing species dynamics at the landscape level to estimate the future prospects of species (e.g. Eriksson, 1996; Husband and Barrett, 1996; Bastin and Thomas, 1999; Hanski, 1999).
Modelling studies that simulate species dynamics at the landscape level focus primarily on the impact of changes in landscape structure (e.g. With et al., 1997; Hanski and Ovaskainen, 2000; Herben et al., 2006; Alados et al., 2009). However, for species that respond slowly to landscape change (e.g. long-lived species with limited dispersal; e.g. Mildén et al., 2006), changes in landscape structure may not be the most important factor influencing species dynamics. Instead, the factors that affect local population dynamics are probably crucial in influencing landscape-level species dynamics over the short term.
Browsing by large herbivores is one of the most important factors affecting local population dynamics of plant species (e.g. Bergelson and Crawley, 1992; Augustine and Frelich, 1998; Russell et al., 2001; Rooney and Waller, 2003). Browsing can strongly influence local population dynamics by affecting the components of the plant life cycle, such as seedling survival (Paige and Whitham, 1987; Knight et al., 2008), plant seed production (Knight et al., 2008; Ehrlén and Münzbergová, 2009; Lin and Galloway, 2009) and the probability of flowering in the next season (Knight et al., 2008; Ehrlén and Münzbergová, 2009). In addition to these negative effects, large herbivores can have positive effects on long-distance dispersal (reviewed in Nathan et al., 2008) and thus on species colonization. Large herbivores can also positively affect plant population growth rate by enhancing seedling recruitment (reviewed by Maron and Crone, 2006). The number of large herbivores, such as roe deer, has been increasing in the agricultural European landscape in recent decades (Meriggi et al., 2008). Increased herbivory and dispersal rates due to a higher number of large herbivores can have both negative and positive effects on the prospects of plant species in the landscape.
In the present study, we estimated the prospects of a grassland polycarpic perennial herb, Scorzonera hispanica, at the landscape level, incorporating the effect of browsing by large herbivores on the species dynamics. Specifically, the aims of the present study were to model the landscape dynamics of S. hispanica in northern Bohemia (Czech Republic) and to simulate the prospects of this species in the near future. Although S. hispanica is considered endangered in the Czech Republic, it is common in the study area. In the area, the species occurs on clearly delimited patches. Some of these patches are suitable but unoccupied, as identified by means of a sowing experiment indicating that S. hispanica is dispersal-limited (Münzbergová, 2004). Chýlová and Münzbergová (2008) demonstrated that this species prevails in grasslands established for at least 60 years, indicating that the dynamics of the species are quite slow. S. hispanica does not form a permanent seed bank (Münzbergová, 2004). Recolonization is thus only possible by means of long-distance dispersal. All the above properties indicate that S. hispanica fulfils the criteria for possessing metapopulation dynamics (Freckleton and Watkinson, 2002). Information on landscape-level dynamics of this species can thus be generalized to other species fulfilling the same criteria with slow response to landscape changes. The identified patterns could thus be generalized for many grassland and forest-understorey long-lived perennial herbs.
Browsing by large herbivores has been observed in most S. hispanica populations. We have also observed the destruction of habitats of S. hispanica within the study region, due primarily to ploughing, the construction of solar power stations or rooting by wild boars. Therefore, in the present study, we evaluated the effects of both browsing by large herbivores and population destruction on the landscape-level dynamics of S. hispanica, as they both may influence its future prospects.
To understand the future dynamics of S. hispanica in the landscape, we asked the following questions: (1) What are the future prospects of S. hispanica in the current landscape and under the current rate of herbivory (i.e. browsing by large herbivores)? (2) What is the effect of herbivory on the future prospects of S. hispanica? (3) What are the combined effects of herbivory and population destruction on the prospects of S. hispanica?
To answer these questions, we parameterized a model of landscape dynamics for S. hispanica based on available information on the distribution of suitable habitats, the local population dynamics (including the current rate of herbivory and the risk of population destruction), dispersal ability and current population sizes. We then simulated the prospects of the species after 30 years under a wide range of herbivory rates and with different levels of risk of population destruction. We assumed that the landscape would not change dramatically over such a short period and that the response to landscape changes would be slow. Under these assumptions, browsing by large herbivores is expected to be the primary factor influencing the prospects of S. hispanica. Model credibility was tested by performing sensitivity analyses of the model parameters.
MATERIALS AND METHODS
Study species and study area
Scorzonera hispanica L. (Asteraceae) is a rare, allogamous, polycarpic perennial herb inhabiting the dry grasslands of central and southern Europe. It has a single rosette and a single flowering stalk with one to seven yellow flowerheads. It is occasionally cultivated for its edible rootstock (Chater, 1976). The fruits of S. hispanica are achenes with a pappus. The presence of the pappus enables dispersal by wind and exozoochory. The species does not form a persistent seed bank; the seeds, which do not germinate, decompose within 2 years (Münzbergová, 2004).
In the Czech Republic, S. hispanica is a native species and is considered endangered. It occurs in central and northern Bohemia and in southern Moravia, occupying calcareous dry grasslands (alliance Bromion erecti of Ellenberg, 1988). To model the prospects of S. hispanica under different rates of herbivory and population destruction, we focused on a typical agricultural landscape with a common occurrence of both S. hispanica and large herbivores. All study populations are browsed by ungulates. No other type of herbivory has been observed. Roe deer, mouflon and wild boar are very common in the landscape, whereas fallow deer and red deer occur only rarely (Municipality Litoměřice, Department of Environment). Only browsing by roe deer, common herbivores of numerous plant species in both grassland and forest-understorey (e.g. Gill et al., 1996; Jepsen and Topping, 2004; Hewison et al., 2007), has been observed in S. hispanica populations. However, we consider that the other large herbivores in the landscape can also occasionally browse S. hispanica.
The study area (4·39 × 4·39 km) was situated in northern Bohemia in the Czech Republic (50°33′26″N, 14°12′45″E, to 50°31′21″N, 14°17′3″E). Calcareous dry grasslands are typical of the landscape. These grasslands form distinct patches surrounded by shrubs and large agricultural fields. These formerly maintained grasslands are now unmanaged and therefore experience very slow succession of shrubs and trees. Population sizes range from three to 2500 flowering individuals. Genetic variability in the field is high; Nei's genetic diversity values range from 0·04 to 0·32, indicating that all populations are genetically variable (Münzbergová and Plačková, 2010). Large herbivores favour the flowering stalks of S. hispanica. Our long-term field observations indicate that the flowering stalks of S. hispanica are browsed extensively without signs of leaf herbivory on the browsed individuals or on the surrounding vegetation (Z. Münzbergová, pers. obs.).
Field data collection
All dry grassland patches (73 in total, from 48·7 to 214 396·3 m2) in the study area were located within the region studied by Chýlová and Münzbergová (2008). In their study, a digital map of dry grassland patches was created and the presence of S. hispanica and 65 other species (Supplementary Data Table S1) were recorded at each patch. We added data from 12 populations (patches) of S. hispanica outside the study area to the present dataset to increase sample size. All external patches were ≤30 km from the study area and ranged in size from 882·9 to 62 365·9 m2. The external patches all hosted the same dry grassland vegetation (i.e. Bromion erecti, Ellenberg, 1988) as the patches within the study area. At each external patch, we recorded the presence of the 65 selected species of dry grassland vegetation. We counted the number of flowering S. hispanica individuals at all patches. We surveyed 85 patches of dry grassland, 35 of which hosted S. hispanica. The external patches were used to improve the predictive power of models of patch suitability and of the herbivory rate at each patch. External patches were not used to simulate the prospects of S. hispanica in the study area.
To model the impact of large herbivores on S. hispanica landscape-level dynamics, we incorporated the effect of herbivory on performance of S. hispanica into transition matrix models of the local population dynamics of the species. We used a set of eight transition matrices containing three size classes (seedling, large vegetative and flowering individuals) to simulate local population dynamics. These eight matrices were constructed for a previous study (Münzbergová, 2006) and included data collected between 2001 and 2004 in three populations over three transition intervals (population nos. 16, 18 and 20 in table 1 in Münzbergová, 2006). Population size ranged between 1632 and 2464 individuals, with at least 150 individuals marked in each population; see Münzbergová (2006) for additional details. Two populations (nos. 16 and 20) are found within the present study area; the third is nearby and occurs in the same type of habitat. This latter population is among the 12 external populations described above. We considered these populations to be representative as they contain a sufficient number of individuals for studying population dynamics and exhibit habitat conditions typical of other populations in the area. Two populations have been largely stable over the last 10 years. However, all marked plants in the third population were destroyed by wild boars during the last transition period; therefore, no transition matrix could be built from these data.
|Patch suitability||Herbivory rate|
|Local habitat characteristics|
|Ellenberg indicator values|
|Forest 1 km||x||n|
|Shrub 1 km||x||–||0·359|
|Forest 0·5 km||x||n|
|Shrub 0·5 km||x||n|
|Patch suitability||Herbivory rate|
|Local habitat characteristics|
|Ellenberg indicator values|
|Forest 1 km||x||n|
|Shrub 1 km||x||–||0·359|
|Forest 0·5 km||x||n|
|Shrub 0·5 km||x||n|
Key: +/– represents the positive/negative effect of characteristics included in the final model (significant values at P < 0·05 are in bold type), n indicates characteristics not included in the model and x indicates characteristics excluded from the test.
The plants used for matrix construction experienced browsing by large herbivores; however, browsing intensity was not quantified. It was thus necessary to identify those transitions within the matrices that were affected by herbivory and to replace these transitions by probabilities with a quantified rate of herbivory. Z. Červenková (unpubl. res.) found that only flowering stalks were browsed; there was very little herbivore damage to vegetative plants. Z. Červenková (unpubl. res.) also estimated the impact of large herbivores on performance of flowering S. hispanica in a field experiment (Supplementary Data Appendix S1). Specifically, she protected selected plants from browsing using cages and compared the performance of intact and browsed plants. She found that herbivory decreased the seed production and the production of clones by flowering plants. No other impact of herbivory was found. In addition, Münzbergová (2004) demonstrated that the recruitment and survival of seedlings and adult S. hispanica plants are not affected by the above-ground biomass at the localities. Seedling recruitment and survival are also unaffected by the presence of open spaces in the vegetation (Z. Münzbergová, unpubl. res.). These findings suggest that neither biomass removal nor an increase in canopy openness due to herbivory affect the reproductive success of S. hispanica. We therefore focused only on the impact of large herbivores on flowering individuals in the present study.
To estimate the rate of herbivory of flowering stalks, we collected data on intensity of browsing from 21 S. hispanica populations of varying size in 2009 and 2010. Ten populations were within the study area, and 11 were external. In all 21 populations, we recorded the total number of browsed and intact flowering plants. In populations comprising fewer than 150 flowering individuals, we recorded browsing data from all flowering plants present. In the larger populations, we collected data from approx. 150 flowering plants sampled along randomly selected transects. The study was conducted at the end of the flowering period (mid-July), when herbivory on S. hispanica ends but faded flowerheads are still present on the stalks (the herbivores browse flowering, not mature, flowerheads).
To estimate the suitability of individual grassland patches for S. hispanica and to identify the factors affecting the herbivory rate, we recorded data on 26 characteristics at each patch. These included both local habitat characteristics and characteristics describing landscape structure (Table 1), e.g. the location of individual patches in the landscape. To obtain data on the local habitat characteristics, we first constructed digital elevation models (DEMs) with a 5-m grid size. DEMs were derived from digital contours (1 : 10 000, 2-m vertical distance between contours) provided by the Czech Office for Surveying, Mapping and Cadastre. DEMs were constructed for the entire study area and the 12 external patches in ArcGIS 9·2 (Environmental Systems Research Institute Inc., Redlands, CA, USA). Based on these models, we created grids of slopes and potential direct solar irradiation (PDSI) for the 21st day of the month from December to June using ArcGIS 9·2, and created grids of topographic wetness index (TWI) using SAGA GIS 2·0·4. (SAGA User Group Association, Hamburg, Germany). For each patch, we then calculated the logarithm of the total area, the mean values of slopes, PDSI (from December to June) and TWI (Table 1).
Other local habitat characteristics were calculated using the presence of the 65 selected species from our species list (Supplementary Data Table S1). First, we calculated the Beals index, which expresses the probability of a species presence at a patch using the number of joint occurrences with other species (Beals, 1984; Münzbergová and Herben, 2004). We used the presence of all plant species from the species list in all patches for this calculation. Second, we calculated Ellenberg indicator values of light, temperature, moisture, nutrients, soil reaction and continentality (Ellenberg, 1988; Table 1) for each patch using all species recorded at the patch.
To obtain the parameters describing landscape structure around the patches, we calculated the nearest distance between each patch and shrubs, forests, roads and villages; we also recorded the amount of shrub and tree cover surrounding each patch (Table 1). We used digital maps of shrubs, forests, roads and villages for these calculations. The digital map of shrubs and forests based on NATURA 2000 mapping was provided by the Agency for Nature Conservation and Landscape Protection of the Czech Republic. The digital maps of roads and villages were created by combining information on the latest online cadastral and orthophotomaps provided by the Czech Office for Surveying, Mapping and Cadastre and by the Czech Environmental Information Agency, respectively. Using the digital maps, we first calculated the areas of both shrubs and forests within both a 500-m and 1-km radius of each patch. We then calculated the distance between each focal patch and (1) the nearest shrub, (2) the nearest forest, (3) the nearest road and (4) the nearest village (Table 1). Distances were calculated between centre points of each patch to the boundaries of these objects using ArcGIS.
To estimate patch suitability, we tested for the effects of habitat characteristics on the occurrence of S. hispanica. We used a generalized linear model (GLM) with a binomial distribution of the dependent variable (presence/absence of S. hispanica) in this test. We excluded data on landscape structure from our independent variables, as they related to patch availability not to patch suitability. We used data from all 85 patches, 35 of which hosted S. hispanica. To simplify the model (correlation matrix in Supplementary Data Table S2), we used step-wise both-direction regression starting with the maximal model. We identified those habitat characteristics that best explained S. hispanica presence using the Akaike information criterion (AIC; Crawley, 2002). Based on this model, we calculated the probability of S. hispanica presence at each patch (Crawley, 2002). These probabilities were used to identify suitable unoccupied patches for S. hispanica (see ‘Simulation plane’ below).
Our investigation of the factors determining S. hispanica herbivory rate involved a small number of observations (21 populations). Therefore, we primarily selected those landscape and habitat characteristics expected to influence herbivory rate (Table 1). Specifically, we used the area of forests and shrubs surrounding each patch within a radius of 0·5 and 1 km and the distances to the nearest forest, shrub, village and road as possible factors influencing the behaviour of large herbivores and the subsequent herbivory rate (e.g. Welch et al., 1990; Tufto et al., 1996; Hewison et al., 2001; Nilsen et al., 2004; Coulon et al., 2008). We also evaluated factors related to site vegetation, including S. hispanica population size, TWI, PDSI in June (i.e. in the growing season), the Ellenberg indicator value for nutrients and the Beals index (Table 1).
To identify the characteristics influencing herbivory rate, we used the mean rate of herbivory in 2009 and 2010 as the dependent variable (herbivory rates did not differ significantly between years, data not shown). The total numbers of browsed and intact plants over both years were used as the dependent variables with a binomial distribution in a GLM. As we had three similar measures of the impact of forest and shrubs, we made three partial tests, using: (1) the amount of cover of shrubs and forests within 1 km, (2) the amount of cover within 500 m and (3) the proximity of the nearest shrub and forest. In each test, we performed step-wise, both-direction logistic regression (using AIC; Crawley, 2002) starting with the maximal model. We then chose the best model (i.e. the model with the lowest number of independent variables and the highest explanatory power). Using the selected model, we calculated the predicted rate of herbivory (Crawley, 2002) at each patch in the area. All analyses were performed in S-Plus Professional Release 2 (MathSoft, Inc., Seattle, WA, USA).
Model description and estimation of model parameters
To simulate the dynamics of S. hispanica in the landscape, we used a dynamic, spatially explicit landscape-level model presented in previous studies by Münzbergová et al. (2005), Herben et al. (2006) and Mildén et al. (2006), following similar methods. This model does not assume equilibrium between species extinction and colonization; this assumption is important because disequilibrium species dynamics have often been observed in long-lived species following rapid changes in landscape structure (e.g. Matlack, 1994; Eriksson, 1996; Brunet et al., 2000; Lindborg and Eriksson, 2004; Herben et al., 2006). The model uploads (1) the information on location and size of suitable patches (habitats) for a species, (2) the initial habitat occupancy including local population sizes, (3) a set of matrices simulating local population dynamics and (4) the coefficients of dispersal curves (exponential and/or hyperbolic functions) and proportion of seeds dispersed independent of distance (for model details see Supplementary Data Appendix S2).
Suitable patches were identified on the grid (5-m cell resolution) by the probabilities of S. hispanica presence. Patches were classified as either suitable or unsuitable by finding the lowest calculated probability in the set of patches that, in the actual study area, host S. hispanica. We then considered all patches of the same or higher probability to be suitable for the species assuming that S. hispanica occurred on suitable patches only. We also found probability thresholds using methods recommended by Liu et al. (2005). As the threshold we took either the prevalence of the model-building data or the average predicted probability of the model-building data. Compared with the original threshold, these two thresholds identically identified three more unoccupied patches as unsuitable for S. hispanica. The lower number of suitable unoccupied patches was then used to estimate the sensitivity of the model to landscape structure. However, the changed model provided similar results to the original and is not discussed further.
Local population dynamics, herbivory and population destruction
Suitable patches were classified according to the predicted herbivory rate into 11 categories of habitat quality, corresponding to proportions of browsing of 0–100 % at 10 % step intervals. We then used the eight available transition matrices of the Münzbergová (2006) study to build 88 additional matrices (11 from each matrix). Each set of eight matrices included the rate of herbivory corresponding to specific herbivory rate (ranging from 0–100 %, 10 % step intervals). Specifically, in each of the eight matrices, we substituted those transitions significantly affected by browsing with the weighted mean of transition values in browsed and intact plants found by Z. Červenková (Supplementary Data Appendix S1). The set of the eight matrices with specific rate of herbivory was assigned to each patch according its predicted herbivory rate.
Ploughing, construction and rooting by wild boars occasionally occurring in the landscape can cause destruction to varying extents. Therefore, we included the probability of population destruction in the model. We had no reliable estimates of the probability and extent of population destruction. However, we assumed that some individuals could survive during the destruction. We thus used the 88 transition matrices described above and decreased all transition probabilities by 90 % to obtain a set of 88 modified matrices. In this way 90 % of all individuals that would have survived into the next year did not survive. Modified matrices were used at various frequencies to original ones. In all simulations, except those modelling the impact of population destruction on prospects of S. hispanica, we used a frequency of one disturbance matrix per 29 original matrices, i.e. one population destruction per 30 years per population. This proportion was chosen based on the observation of Münzbergová (2006) and our subsequent monitoring of the populations in the area.
Initial population size
The numbers of seedlings and vegetative individuals at each patch were calculated from the numbers of flowering individuals (counted in the field) according to the mean stable stage distribution occurring under a specific rate of herbivory. However, stable stage distributions are reached only in populations with stable local population dynamics. To estimate the sensitivity of the model to this assumption, we used half the numbers of seedlings and vegetative individuals calculated from the stable stage distribution. The results of this alternative model were, however, very similar to the original and are therefore not reported further.
To simulate density-dependence, we estimated the maximum population density at any patch, based on the number of S. hispanica individuals at each patch and patch size. The calculated maximum was 0·97 individuals m−2. Based on our field experience, we assumed that the maximum density of seedlings was four times higher than the maximum density of vegetative or flowering individuals. Thus, seedlings had one-quarter the competitive effect of flowering and vegetative individuals. We also performed a sensitivity analysis of this parameter, using 0·97/3 or 0·97 × 3 individuals m−2 during the simulations.
We assumed no incoming diaspores to the simulation plane, as the study area was somewhat isolated from other S. hispanica populations (the nearest flowering population was 4·3 km from the area border). Therefore, low numbers of incoming diaspores could be expected. Outgoing diaspores during the simulations were considered to be lost. S. hispanica was expected to disperse by wind and exozoochory. Wind dispersal was modelled as distance-dependent using a negative normalized exponential function (Münzbergová et al., 2010):Bullock and Clarke, 2000; Nathan et al., 2002). We thus modelled S. hispanica landscape dynamics under several scenarios of dispersal ability (including different exponential curves). We then checked whether higher/lower dispersal ability influenced the results of the model. Specifically, dispersal coefficient α was calculated as 1/D, where D is the mean dispersal distance of the seeds calculated from the formula (e.g. Augspurger, 1986; Soons and Heil, 2002; Tremlová and Münzbergová, 2007) Münzbergová (2004). Due to the lack of data on the range of terminal velocities, we chose a range of 1·78 m s−1 ± 33 %. Based on the ranges of all parameters (w, h and t), we calculated a dispersal distance range of 0·24–22·80 m and a mean dispersal distance of 2·51 m. In the simulations, we used the mean dispersal distance. Minimum and maximum dispersal distances were used to perform sensitivity analyses of the dispersal parameter.
We did not have an estimate of the proportion of seed dispersed via animal fur in the field as obtaining realistic estimates of such a value is difficult (Nathan et al., 2008). Therefore, we assumed that only 0·1 % of all seeds were dispersed by exozoochory and by rare events (see also Münzbergová et al., 2005). As herbivores attack individuals during flowering and not when the seeds are mature, the proportion of damaged flowerheads with mature seeds is very low. Therefore, endozoochory was not considered in our simulations. Exozoochory was assumed to affect primarily the long-distance dispersal of S. hispanica. It was modelled as independent of distance. Although this assumption seems to be unrealistic, it was used in previous studies (e.g. Münzbergová et al., 2005; Mildén et al., 2006). In our study, dispersal was modelled within small study area. This does not suggest that dispersal is independent of distance at any scale, but rather that the animals can easily cross the whole model landscape within a short period of time. The sensitivity analysis of distance-independent dispersal was performed using 1 % and 0·01 % of the dispersed seeds.
All forecasts were run for 30 steps (30 simulation years), with each forecast replicated 100 times. We ran the simulations for up to 30 years, assuming that the landscape would not change dramatically over such a period. However, running the model for 100 years provided qualitatively very similar results, with the time frame having no impact on our conclusions (data not shown).
We estimated the impacts of the range of the model parameters on the prospects of S. hispanica under different herbivory rates. First, we simulated the prospects of the species under the current rate of herbivory predicted for individual patches in the area. We then simulated a gradual (at 10 % intervals) decrease (or increase) in the predicted rate of herbivory over the entire study area, until the herbivory rate (i.e. the number of browsed flowering individuals) decreased to 0 % (or increased to 100 %) in all patches. In this way, we obtained 17 different simulations of different herbivory rates.
To assess the effect of the frequency of population destruction on the prospects of S. hispanica, we simulated different rates of population destruction under different herbivory rates. The frequencies of population destruction ranged between 0·5 and 10 disturbances per population per 30 years. We obtained 11 different simulations of the impact of disturbance regime on the prospects of S. hispanica under 17 different rates of herbivory.
Determinants of patch suitability and herbivory rate
The step-wise regression identified 11 of 17 local habitat characteristics as significant predictors of the presence of S. hispanica in a patch (Table 1). These include PDSI in various months, the Beals index and the Ellenberg indicator values. Based on the model, we identified 31 patches as suitable for S. hispanica, eight of which were unoccupied. Five of 12 landscape and local habitat characteristics were selected in the step-wise regression as significant predictors of the rate of herbivory in a patch (Table 1): patch area, TWI, shrub cover within 1 km of the patch, proximity of the nearest village and proximity of the nearest road. The predicted rate of herbivory ranged from 40 to 100 % (mean 77 %, median 80 %) among single patches.
Impact of herbivores at the landscape scale
Our simulations revealed a strong effect of large herbivores on the long-term prospects of S. hispanica. Simulation using the predicted herbivory rate showed an equilibrium in the number of S. hispanica individuals in the area (mean population size after 30 years = 86 371, s.d. = 23 632; initial population size = 78 462, Fig. 1A). Under high rates of herbivory, S. hispanica tended to go extinct; under low rates, population size increased substantially. Similarly, herbivory rates had a negative effect on patch occupancy (Fig. 1B). Using the current rate of herbivory in our simulations, the patches hosting small populations experienced higher turnover of patch occupancy than the patches hosting large populations (Fig. 2).
Sensitivity analyses of dispersal parameters (Supplementary Data Figs S1a, b and S2a, b) revealed an effect of both wind dispersal and exozoochory on the number of occupied patches under low rates of herbivory at the end of the simulations (Supplementary Data Fig. S2a, b). Sensitivity analysis of maximum population density (Supplementary Data Figs S1c and S2c) showed a strong effect on the total number of individuals under low rates of herbivory (Supplementary Data Fig. S1c). In both cases, parameter effects disappeared under high rates of herbivory, which indicates that increased herbivory reduces the positive effects of longer dispersal rates and higher carrying capacity (Supplementary Data Figs S1 and S2).
Our simulations also revealed that not only herbivory rates but also frequencies of population destruction had a strong effect on the landscape-level dynamics of S. hispanica. Rates of population destruction higher than those observed in the field (i.e. one per 30 years in a patch) led to a considerable decrease in the number of individuals and the number of occupied patches in the area (Fig. 3). The pattern was observed for all rates of herbivory. A rapid decline of the number of S. hispanica individuals was observed when high frequencies of population destruction were combined with high rates of herbivory.
In this study, we demonstrated the importance of browsing by large herbivores on the landscape-level dynamics of S. hispanica. Despite the negative effects of browsing on the performance of S. hispanica, the landscape-level dynamics of S. hispanica are currently at equilibrium. The future prospects of this species, however, depend on the prospects of large populations. We also found that the potential effects of large herbivores on landscape-level dynamics may be considerable. Simulated declines or increases in the rate of herbivory throughout the landscape by more than 20 % often led to the rapid expansion or decline of S. hispanica. A similar effect of herbivory rate was observed in the sensitivity analyses of various model parameters. These findings indicate that herbivores can be among the major drivers of landscape dynamics of long-lived perennial herbs.
The results of the simulations demonstrated relatively high turnover of the local populations, especially those that are small. This indicates that despite being a long-lived perennial, S. hispanica exhibits features of metapopulations, as we expected. The metapopulation framework (Hanski, 1999) is thus a suitable approach for modelling the dynamics of S. hispanica and other similar species. More specifically, the high turnover of the small populations and high survival of the large ones suggests that S. hispanica is a likely representative of species with mainland–island metapopulation dynamics (Harrison, 1991). From a conservation point of view, survival of large populations is crucial for survival of the whole metapopulation. These large populations may, however, be threatened in the landscapes by factors such as human-induced population destruction as well as an increased rate of herbivory.
Effect of herbivores on landscape-level dynamics
We expected S. hispanica to decline by the end of simulations using the current rate of herbivory because of (1) the negative effects of large herbivores on performance of S. hispanica, (2) its slow dynamics and dispersal-limitation (Münzbergová, 2004) and (3) the increased fragmentation of the landscape over the last 60 years (Chýlová and Münzbergová, 2008). Unexpectedly, our simulations of the prospects of S. hispanica following 30 years of current landscape conditions (herbivory and disturbance) suggest that the total number of individuals in the study area is largely stable. We can, however, expect slight declines in future habitat occupancy. S. hispanica was unable to establish new large populations during the simulations due to the high risk of extinction of small populations (as demonstrated also in the field by Münzbergová, 2006).
The maintenance of landscape dynamics of S. hispanica near equilibrium (except the extinction of small populations) can be explained by the type of its local population dynamics. Population dynamics are very stable over time (even with the current high rate of herbivory) due to the high survival probability of individuals and occasional clonal reproduction. Nevertheless, a simulated 20 % decline in the current herbivory rate markedly increased the number of seeds produced, resulting in higher seed dispersal and more successful colonization. In contrast, a 20 % increase in the current rate of herbivory led to a serious population decline. This decline was due primarily to the large changes occurring in the most abundant populations. Sensitivity analyses of the dispersal parameters revealed that habitat occupancy depended partly on the estimation of the dispersal parameters. However, the dispersal parameters did not influence overall decline/increase of S. hispanica under higher/lower than current rates of herbivory. In addition, changes in dispersal rates had only negligible effects on landscape-level prospects of S. hispanica under high rates of herbivory. This indicates that the negative effects of herbivory on performance of species, specifically on generative reproduction (e.g. Knight et al., 2008; Ehrlén and Münzbergová, 2009; Lin and Galloway, 2009), can be much stronger than the possible positive effects of herbivores as dispersal agents (Fischer et al., 1996; Nathan et al., 2008). The effect of herbivory was enhanced significantly when combined with population destruction (resulting from large disturbances). The importance of population destruction to species landscape dynamics has been previously demonstrated, for example by Münzbergová et al. (2005). Such a clear negative effect of herbivory and population destruction on species dynamics is caused by the absence of any positive effects of these activities on plant performance. If overcompensation (Paige and Whitham, 1987) or enhanced seedling recruitment (Gomez, 2005) was found in the case of S. hispanica, the effect of herbivory rate and population destruction on its landscape-level dynamics would be much less clear.
The current equilibrium state of the S. hispanica metapopulation may reflect several factors. One possibility is that the expansion of S. hispanica has been constrained in the past (e.g. by cattle grazing and mowing). This hypothesis is supported by the fact that many S. hispanica populations occur in former pastures (Münzbergová, 2004). In such a scenario, current patch occupancy should reflect high past landscape connectivity. Alternatively, landscape connectivity may still be the same, but with an ongoing increase in the rate of herbivory and the frequency of large disturbances in the area. This would result in reduced growth rates of local populations and in higher probability of extinction of small populations.
Patterns of herbivory
First, it is important to note that the proportion of browsed individuals in a population does not necessarily relate to the frequency of visits by large herbivores, as the herbivores can readily pass over patches without browsing. We found less herbivore damage in larger patches and in patches surrounded by high shrub cover. This pattern corresponds to the results of previous studies of the habitat preferences of large herbivores. For example, roe deer prefer small patches (Aulak and Babinska-Werka, 1990; Welch et al., 1990) with rich ground vegetation (Welch et al., 1990), and the density of deer increases with increasing habitat heterogeneity (Kie et al., 2002) and the density of habitat edges (e.g. Tufto et al., 1996; Saïd et al., 2005; Miyashita, 2008). According to Lamberti et al. (2006), roe deer prefer open habitats (e.g. orchards and fields) to woodlands and scrublands. This latter observation suggests that the presence of shrubs or forest may decrease the attractiveness of patches to herbivores.
Another important factor affecting behaviour of large herbivores is human activity in the landscape, especially developed areas and roads. Several studies have found that these factors negatively impact deer density (e.g. Hewison et al., 2001; Coulon et al., 2008). In our study, populations near villages were browsed less heavily whereas the proximity of roads had the opposite effect. The increased rate of herbivory near roads may result from the use of roads as corridors by large herbivores. The affinity of large herbivores to individual patches may have been also affected by landscape topography and surface (Coulon et al., 2008). We found a higher rate of herbivory in the drier patches above valleys than in the wetter patches closer to valley bottoms.
Estimated parameters and model credibility
Several things should be kept in mind in interpreting our simulation results. The sensitivity analyses indicated that dispersal ability and carrying capacity had an effect on the total number of S. hispanica individuals under low, but not high, rates of herbivory. Similarly, these parameters had a greater effect on habitat occupancy under low than under high herbivory rates (Supplementary Data Figs S1 and S2). Similar results were found when using smaller initial population sizes and running the simulations for 100 years (data not shown). These results indicate that our conclusions regarding the effects of large herbivores on the species landscape dynamics are independent of the parameter estimates.
To simulate the prospects of S. hispanica in the landscape, we set the initial patch occupancy at 74 %, corresponding to the observed occupancy of this species in the study area. This occupancy level is significantly higher than the 32 % found by Münzbergová (2004), who used sowing experiments to identify patches suitable for, but unoccupied by, S. hispanica. However, in the study by Münzbergová (2004) patches were distributed over a larger area (approx. 400 km2) than in our study (approx. 20 km2), but within the same landscape. In addition, in the study of Münzbergová (2004), the predictions of patch suitability were based on seedling establishment, which does not necessarily reflect patch suitability for adult plants (Ehrlén et al., 2006). In our study, suitable patches were identified using a combination of abiotic conditions and species composition. These types of factors have been shown previously to explain species distribution (e.g. Dupré and Ehrlén, 2002; Münzbergová, 2004; Ehrlén et al., 2006; Chýlová and Münzbergová, 2008). As suggested, for example, by Tájek et al. (2011) and observed in the present study, the combination of these two types of factors provides the best predictions of habitat suitability for a species. Specifically, we identified drier, shaded, basic, nutrient-rich patches as more suitable for S. hispanica in dry grasslands. Patch preference may reflect both specific abiotic conditions and past land use (particularly as pastures). Nevertheless, the sensitivity analyses revealed that a reduction in the number of suitable patches had little effect on the model results.
In our simulations, local population dynamics was assumed to be the same among patches (except for the impact of herbivory) and largely stable over time. We simulated demographic stochasticity representing random changes in local population dynamics over time. There was no indication for a gradual change in environmental conditions. In our simulations we thus assumed that among-year variation in the local population dynamics was caused only by environmental stochasticity. Environmental stochasticity was simulated by drawing a random transition matrix (for each population, in each step) from a set of matrices. The matrices were very similar and thus their random sample had little effect on the local population dynamics. Although differences between patches and changes in local population dynamics over time could occur, we argue that these factors are unlikely to strongly affect our conclusions. First, all our populations occur within a small area under very similar habitat conditions, minimizing potential differences among populations. Second, our simulations extended only 30 years into the future, making it unlikely that habitat conditions will change dramatically. Third, it is likely that the responses of S. hispanica to any landscape changes would be slow, and thus minimal over this period. Therefore, browsing by large herbivores is probably the primary factor influencing the prospects of S. hispanica.
Finally, we assumed a stable rate of herbivory at each patch during the individual simulations. However, the incidence of browsing by large herbivores fluctuates between years due to changes in their abundance. For example, 22 % more roe deer were recorded in the landscape in 2008 than in 2007 (Municipality Litoměřice, Department of Environment). However, we suggest that fluctuations in herbivory rate are unlikely to alter the main conclusions of the model.
Our field observations indicate that over 60 % of flowering S. hispanica individuals are damaged by large herbivores in most populations each year. Our simulations, however, suggest that current dynamics of S. hispanica are approximately in equilibrium under the current rate of herbivory and frequency of large disturbances (one per 30 years per population). The simulation results also revealed a higher survival probability of large populations than that of small ones. Therefore, under current landscape conditions, the prospect of S. hispanica in the landscape depends heavily on the prospects of large populations.
Simulations of the effect of herbivory rate on the dynamics of S. hispanica indicated that a decline or increase in the herbivory rate of more than 20 % over the entire landscape could lead to a rapid expansion or decline of the species. This effect is predicted to be much stronger under the additional occurrence of disturbance. These results confirm our hypothesis that browsing by large herbivores can have dramatic effects on the landscape dynamics of species if important components of the life cycle are strongly affected by these herbivores.
Finally, as concluded in other studies, our study suggests that the probability of population destruction should be incorporated into models predicting changes in species distributions. Incorporating the effect of large herbivores and population destruction into models of species landscape dynamics should be a major endeavour of future metapopulation studies.
Supplementary data are available online at www.aob.oxfordjournals.org and consist of the following. Appendix S1: details of collection and testing of the data on the impact of large herbivores on performance of S. hispanica. Appendix S2: detailed description of the landscape-level model. Table S1: list of selected dry grassland species. Table S2: correlation matrix of habitat characteristics recorded at patches of dry grasslands. Fig. S1: sensitivity analysis of the model: the effect of herbivory rates on the total number of S. hispanica individuals. Fig. S2: sensitivity analysis of the model: the effect of herbivory rates on patch occupancy by S. hispanica.
We thank A. Lampei Bucharova, R. Leimu, J. Knappova, H. Skalova and the anonymous referees for their comments on a previous version of this paper. We thank the Czech Hydrometeorological Institute for providing climate data. This study was supported by GAUK 64709, GACR P504/10/0456 and also in part by long-term research development project no. RVO 67985939.