Climate change has the potential to influence the persistence of ecological communities by altering their stability properties. One of the major drivers of community stability is species diversity, which is itself expected to be altered by climate change in many systems. The extent to which climatic effects on community stability may be buffered by the influence of species interactions on diversity is, however, poorly understood because of a paucity of studies incorporating interactions between abiotic and biotic factors. Here, I report results of a 10-year field experiment, the past 7 years of which have focused on effects of ongoing warming and herbivore removal on diversity and stability within the plant community, where competitive species interactions are mediated by exploitation through herbivory. Across the entire plant community, stability increased with diversity, but both stability and diversity were reduced by herbivore removal, warming and their interaction. Within the most species-rich functional group in the community, forbs, warming reduced species diversity, and both warming and herbivore removal reduced the strength of the relationship between diversity and stability. Species interactions, such as exploitation, may thus buffer communities against destabilizing influences of climate change, and intact populations of large herbivores, in particular, may prove important in maintaining and promoting plant community diversity and stability in a changing climate.
Late Pleistocene climate change resulted in compositional turnover in many palaeo-communities [1–3], and ongoing and future climate change is expected, similarly, to drive the development of novel ecological communities . Community stability, an important precursor to the persistence of communities, has, however, long been studied with an emphasis on biotic drivers, such as competition and species diversity [5,6], with abiotic factors receiving consideration mainly in the context of perturbations, such as drought . Increasing recognition of the importance of the interaction between abiotic and biotic factors in determining species distributional-  and community compositional responses to climate change [9,10] suggests the stability response of ecological communities to climate change may likewise be determined by the interaction between abiotic and biotic factors .
The earliest treatments of community stability were theoretical and focused on biotic interactions such as exploitation and interference among members of food webs across and within trophic levels [12–15]. Ensuing experimental investigations of community stability have, by contrast, focused nearly exclusively on influences of interference interactions on species diversity within a single trophic level [5–7,16,17]. Where interactions across trophic levels have been taken into consideration, however, they have revealed important effects on consumer-level stability of forage species diversity [18,19], effects of soil invertebrate fauna on plant species diversity , and of invertebrate herbivore abundance on relative plant species abundances . Top-down, exploitative effects of consumers on stability of resource-level species remain, however, unexamined. Such effects should be prominent because diversity confers stability  and consumers such as herbivores promote plant species diversity by mediating competitive interactions [23–25]. Such effects should likewise be responsive to abiotic variation such as climate change because herbivores also influence the response of plant species diversity to rainfall variation , and plant productivity and community compositional responses to warming , and variation in snowpack .
In 2002, I initiated a long-term exclosure and warming experiment in a simple, low-Arctic plant community near Kangerlussuaq, Greenland, to investigate the influence of exploitation by large herbivores (caribou, Rangifer tarandus; and muskoxen, Ovibos moschatus) on plant community response to climate change . Previous results from this experiment revealed a mediating influence of herbivory on plant community compositional response to warming: in the absence of large herbivores, the plant community responded to warming by shifting away from dominance by graminoids towards dominance by competitively superior dwarf shrubs, whereas the continued presence of large herbivores maintained the original composition of the plant community under warming . The interaction between warming and herbivore exclusion has subsequently been demonstrated to alter ecosystem function in this system by significantly increasing carbon uptake by the plant community . Here, I focus on the implications of the interaction between warming and herbivore exclusion for diversity and stability within the plant community.
2. Material and methods
(a) Experimental design
Six permanent, circular, large-herbivore exclosures, measuring 800 m2 and separated by up to 1 km, were erected early in the plant growing season in 2002, following baseline sampling of plant community composition on exclosure- and adjacent control sites. Three of these exclosures and their adjacent control sites were selected to receive a warming treatment, whereas the remaining three were set aside for future manipulations. The warming treatment began on two of these exclosure- and adjacent control sites prior to the onset of the plant growing season in 2003, and was expanded to the third pair of exclosure and control sites in advance of the plant growing season in 2004. The warming treatment involved use of open-topped, passive warming chambers constructed and deployed according to the protocols of the International Tundra Experiment , measuring a basal diameter of 1.5 m and height of 0.4 m; ambient, control plots inside and outside exclosures received no warming treatment . The warming chambers elevated near-surface temperatures by 1.4–2.0°C . The number of plots increased from six exclosed/ambient, six exclosed/warmed, seven grazed/ambient and seven grazed/warmed in 2003 to 12, 12, 13 and 13, respectively, in 2004 when the experiment was expanded to the third pair of exclosure- and control sites.
(b) Community composition, diversity and stability
Baseline community composition was quantified prior to the establishment of the exclosure treatment in June, 2002, using non-destructive point-frame sampling, and, thus, was quantified prior to establishment of sampling plots subsequently receiving warming treatments that began at the onset of the plant growing season in 2003. The point-frame used for baseline measurements in 2002 contained 10 sampling pins, whereas that used from 2003 and onward annually contained 20 sampling pins. Beginning in 2003, point-frame sampling was performed on all species on treatment and control plots annually at the peak of the growing season . Point-frame sampling was used because it allows for the assessment of cumulative responses to experimental manipulation without, presumably, altering the structure or growth of the vegetation being sampled. Although suitable for estimating proportional community composition, point-frame data can only, at best, be considered an index of biomass or productivity. Vegetation sampled by this method on non-experimental plots at the study site in 2003 and subsequently clipped, dried at 60°C for 24 h and weighed, indicated closer agreement between pin hits and biomass for graminoids (r = 0.81, p < 0.001) and deciduous shrubs (r = 0.80, p < 0.001) than for forbs (r = 0.71, p = 0.001) . Moreover, the best relationship between pin hits and biomass for graminoids and forbs was log-linear (y = axb; slopes = 0.13 and 0.48, respectively), whereas that for deciduous shrubs was linear (y = bx + a; slope = 1.09) . Rather than attempting to transform pin hit data into estimates of biomass (g m−2), as in an earlier study , I elected to work strictly with raw pin hit data in this study to avoid introducing additional uncertainty into my analyses. As a consequence of the log-linear nature of the relationship between pin hits and biomass for graminoids and forbs, it is possible that the pin hit data may underestimate actual biomass of graminoids and forbs, in comparison with that of deciduous shrubs, at higher levels of biomass. This suggests that the point-frame method may potentially underestimate any biomass increase in response to the warming and exclosure treatments , and, if so, result in an underestimate of variability—the inverse of stability—associated with biomass estimates. If this is the case, then the reduction in stability by the warming and exclosure treatments reported below may actually be more pronounced than detected by this method. Because the index of diversity used in this study is a simple count of functional groups and species encountered on each plot during point-frame sampling, it is unlikely that the point-frame method has biased the estimation of diversity.
For point-frame sampling, each plot was marked permanently at its corners to facilitate consistent orientation of the square point-frame during each sampling event. Plant community diversity was quantified during the first 2 years of annual point-frame sampling as richness of functional groups, including dwarf shrubs, comprising only two species (dwarf birch, Betula nana nana; and grey willow, Salix glauca), graminoids, forbs and moss; non-specific deciduous shrub litter was also recorded but not included in functional group richness. In the second and third years of the study (2004 and 2005), a major outbreak of a noctuid moth, Eurois occulta, severely reduced above-ground biomass across all plant functional groups in the community, without bias, and regardless of warming or exclosure treatments . The caterpillar outbreak, though unforeseen, offered an opportunity to monitor the recovery of the plant community to this major disturbance, an important measure of stability in long-term studies , in response to the ongoing experimental manipulations. Hence, following the major caterpillar outbreak of 2005, increased taxonomic resolution was applied by recording data at the species level for forbs, which had previously been recorded only at the functional group level. Analyses of diversity in this study were performed on data collected from 2005 through to 2011. Diversity in this study is, therefore, quantified for analyses conducted at the community level as the total number of species of deciduous shrubs and forbs, plus the number of additional functional groups (graminoids, mosses, lichens) recorded on each plot; and for analyses conducted at the forb component of the community as the total number of forb species recorded on each plot. Thus, this study does not mirror the more common design of diversity–stability experiments that have relied on manipulation of plant species richness on sown plots , but does address criticisms of that approach [33,34] by representing a test of the diversity–stability hypothesis in a natural, field setting [35,36].
(c) Analyses of total community dynamics
To test the response of each plant functional group, in a community context, to the exclosure and warming treatments over the course of this ongoing experiment, I used a multivariate general linear model (GLM; MANOVA) of proportional contributions of all plant functional groups to community composition, rather than multiple individual univariate GLMs (ANOVAs) performed on each functional group independently. This MANOVA approach was used to allow for possible interactions among dependent variables (in this case, plant functional groups) as well as the responses of multiple dependent variables to common factors. Moreover, MANOVAs are more robust to violations of the assumption of homogeneity of variances than are ANOVAs when the experimental design is balanced with regard to sample size across treatments. MANOVAs of proportional contributions of plant functional groups to community composition were performed on un-transformed, raw proportions of all community components, including dwarf birch, grey willow, graminoids, forbs, moss and litter, because the arcsine transformation was not warranted in this case . The structure of the MANOVAs included year (2003–2011), the exclosure treatment (exclosed or grazed) as well as the warming treatment (warmed or ambient) as fixed factors and the length of the annual warming treatment as a covariate to account for variation among years in the number of days of the warming treatment. Because the assumption of equal variances was not met for any of the plant functional groups, I compared Pillai's trace significance values for each of the independent variables with Wilk's Lambda significance values to determine whether unequal error variances were due to unequal sample sizes, and these significance values were identical in all cases. However, in the interest of reducing the probability of spurious detection of significance, I applied a Huynh–Feldt reduction in degrees of freedom by 25 per cent for the F-tests deriving from this MANOVA, as was carried out for the univariate GLMs performed on forb-only data described in §2c.
Community stability was defined in this study as the inverse of the coefficient of variation (CV) in total live community biomass, indexed by point-frame pin hits per plot, through time [5,38]. I analysed the relationship between community diversity and stability using a univariate GLM of the plot-scale CV (calculated across 2005–2011) of peak annual above-ground live biomass estimated using point-frame sampling as the dependent variable, and the exclosure and warming treatments as fixed factors, a site blocking variable as a random factor, and plot-scale mean community diversity across years (2005–2011) as a covariate, with post hoc pairwise comparisons of the treatment means to test their significance. A Levene's test revealed that the assumption of homogeneity of error variances was not in this case violated (F3,46 = 0.45, p = 0.72).
(d) Analyses of forb dynamics
Because the forb component of the plant community is the most species rich at the study site, additional analyses focused on dynamics and stability within this functional group only. In addition to examining the relationship between species richness and stability within the forb subcomponent of the plant community, as well as differences in the strength of this relationship among experimental treatments, three separate univariate GLMs were performed on the data on forb species richness collected from 2005 through to 2011. The relationship between stability of the forb component of the plant community and forb species richness was derived using estimates of the mean number of forb species counted on each of the study plots during the experiment (2005–2011), and the CV (the inverse of stability) of peak forb biomass (estimated based on point-frame sampling) calculated across the same years. One plot contained zero forbs and so was excluded from the analysis.
The three GLMs of forb species richness were performed in hierarchical fashion, and were aimed at: (i) investigating the dynamics of forb species richness through time after recovery from the caterpillar outbreak in response to the warming and herbivore exclosure treatments; (ii) examining whether the potential influence of the warming or exclosure treatments on forb richness was related to increases in deciduous shrub biomass (indexed by point-frame sampling) in response to those treatments ; and (iii) examining whether the potential influence of increases in deciduous shrubs on forb species richness reflected increases in shrub litter, which may suppress soil nutrient turnover and availability, thereby adversely affecting forb competitors , in response to the warming or exclosure treatments. In the first case, annual means of forb species number were derived from a univariate GLM of the number of forb species counted on each plot during annual point-frame sampling of above-ground biomass; this model included year, the exclosure treatment and the warming treatment as fixed factors, and the annual length of the warming treatment as a covariate. Note that a forb species count of, for example, ‘1’ in any given year on any given plot does not imply that a single individual was present on that plot in that year; in fact, observations of single species most often included several individuals. In the second case, the model included plot-specific forb species number as the dependent variable, the exclosure and warming treatments as fixed factors, year as a random factor and total deciduous shrub biomass (based on point-frame sampling) and the annual length of the warming treatment as covariates. In the third case, the model included plot-specific forb species number as the dependent variable and the exclosure and warming treatments as fixed factors, year as a random factor, and both the annual length of the warming treatment and litter biomass as covariates. This analysis was preceded by a GLM of litter biomass that included the exclosure and warming treatments as fixed factors and year as a random factor to determine whether litter biomass differed according to these treatments; annual length of the warming treatment was omitted as a covariate because it was not significant (p = 0.31), as should be expected because litter production is cumulative on the plots. The means of litter biomass by treatment in figure 3 derive from the GLM of litter biomass described earlier.
In performing the three GLMs of forb species richness, Levene's tests revealed that the assumption of homogeneity of error variances was violated in each instance (p ≤ 0.001 in all cases). This is because the data on counts of forb species followed a non-normal distribution. Because transformations typically applied to remedy departures from normality do not work consistently well for count data , I opted instead to apply a Huynh–Feldt reduction in the degrees of freedom by 25 per cent for the F-tests upon which the significance of the predictor variables was based. Only reduced degrees of freedom and associated p-values are reported. Additionally, I repeated all analyses of the forb species count data using generalized linear models with the same factors and covariates, specifying a Poisson distribution with log-link functions to account for non-normality in the forb data . In every case, results obtained using univariate GLMs with Huynh–Feldt reduced degrees of freedom were confirmed by Poisson-specified generalized linear models. Results of both approaches are reported below as F-values from the univariate GLMs followed by Wald's χ2-values from the Poisson generalized linear models for comparison.
Because multiple (in this case, three), hierarchical GLMs were applied to the data on forb species richness, it could be argued that the probability of detection of significant relationships might have been spuriously increased. However, rather than applying a sequential Bonferroni adjustment of significance tests, I favour a logical and judicious interpretation of the results deriving from this approach  because it is aimed at improving a mechanistic understanding of the response of a dependent variable, forb richness, to experimental manipulation. For instance, if forb richness is shown in the first analysis to decline in response to warming, the second analysis is designed to test whether this decline relates to an increase in deciduous shrubs, and, if so, the third analysis is designed to test whether the effect of deciduous shrubs on forb diversity reflects an increase in shrub litter in response to warming. The sequential nature of these analyses would in and of itself, therefore, preclude advancing to a subsequent test if the preceding test did not identify a significant relationship.
Data are archived online at the US National Center for Atmospheric Research Earth Observing Laboratory (http://data.eol.ucar.edu/codiac/dss/id=106.415).
3. Results and discussion
Across the 7 years of the study during which diversity was monitored (2005–2011), community stability displayed a positive association with community diversity (F1,27 = 10.1, p = 0.003; figure 1). A post hoc comparison revealed significantly greater mean stability on ambient than on warmed plots (F1,48 = 7.06, p = 0.01; figure 1, ‘A’ versus ‘W’), indicating that, overall, warming reduced the stability of biomass production across years at the community level. Further post hoc comparisons indicated a substantial reduction in community stability by the warming treatment on exclosed plots (mean difference =−0.066, p = 0.005; figure 1, ‘EA’ versus ‘EW’), and significantly lower stability on double-treatment plots compared with double-control plots (mean difference =−0.058, p = 0.01; figure 1, grazed-ambient (GA) versus exclosed-warmed, EW). These results indicate that the reduction in plant community stability by warming was further exacerbated by removal of the competition-mediating effects of exploitation of the plant community by large herbivores. In accordance with this pattern, the greatest community diversity and stability in biomass production occurred on GA plots, whereas the lowest community diversity and stability occurred on EW plots (figure 1).
Although the deciduous shrub component of the plant community on study plots at this site comprises only two species, dwarf birch (Betula nana nana) and grey willow (Salix glauca), the forb component includes up to seven to 10 species , although many fewer species co-occur on average within each plot. Following an initial rebound from the peak of the caterpillar outbreak in 2005, forb diversity subsequently declined and remained low on warmed plots, whether grazed or exclosed (figure 2a), reflecting a significant, overall reduction in forb diversity by the warming treatment (F1,240 = 3.77, p = 0.053; Wald's χ2 = 6.39, p = 0.01). The MANOVA of variation in functional group proportional contributions to total community biomass over the entire length of the warming/exclosure experiment (2003–2011) revealed a significant interaction between the warming and exclosure treatments on the proportions of dwarf birch (F1,290 = 4.78, p = 0.03) and grey willow (F1,290 = 4.59, p = 0.033), both of which increased in the community in response to warming and herbivore removal . Hence, the reduction in forb diversity in response to the warming treatment may, in part, relate to increasing dominance by dwarf shrubs on warmed and exclosed plots, a point revisited in greater detail below.
Within the forb component of the community, temporal stability of biomass production from 2005 through to 2011 displayed a positive, nonlinear relationship to diversity, quantified, as described earlier, as the number of species per plot (Spearman r = 0.80, p < 0.001; figure 2b). The nature of this relationship indicates, first, that stability of the forb component of the community varied by approximately 80 per cent across the range of diversity observed; and, second, that in species-poor systems such as those typified by Arctic plant communities, species losses may have disproportionate effects on stability compared with losses in more species-rich systems such as temperate meadows . Furthermore, the strength of the relationship between forb diversity and stability, quantified either as the slope of the regression between ln(species number) and (μ/σ)peak biomass, or as the Spearman correlation between species number and (μ/σ)peak biomass, was greater on grazed plots than on plots exclosed from large herbivores (figure 3a), suggesting that removal of large herbivores weakened the relationship between diversity and stability. As well, there is some evidence, albeit weak, that warming may have enhanced stability within the forb component of the community where large herbivores interacted with the warming experiment: the strength of the stability–diversity relationship was highest on grazed, warmed plots (figure 3a), and variability across years in biomass production by forbs (the inverse of stability) was lower on grazed, warmed plots than on exclosed, ambient plots (mean difference =−0.62, p = 0.02). These latter results suggest that large herbivores may promote plant community stability in a warmer climate.
As suggested earlier, two possibilities might explain the reduction in diversity and stability in response to warming and large herbivore removal. First, in Arctic plant communities such as the one in this study, removal of large herbivores may reduce species diversity within the plant community by promoting the expansion of deciduous shrubs , which can outcompete forbs for limiting nutrients such as nitrogen, and, through expansion and associated inputs of recalcitrant litter, inhibit their growth through immobilization of soil nitrogen or suppression of nutrient cycling rates [43,44]. Second, expansion of shrubs may result in suppression of soil microbial activity owing to shading and cooling of the soil surface , potentially further reducing soil nutrient availability, with adverse consequences for inferior competitors within the plant community . Several results are consistent with the former hypothesis, but data are lacking in this study to address the latter. For instance, the GLM of forb species richness across the study period revealed a significant warming effect (F1,244 = 3.86, p = 0.051; Wald's χ2 = 7.73, p = 0.005), and a significant relation to shrub biomass (F1,240 = 29.6, p < 0.001; Wald's χ2 = 3.78, p = 0.05). Furthermore, the GLM of litter biomass over the same period identified a significant effect of the warming treatment (F1,4 = 10.3, p = 0.033), but not of the exclosure treatment or its interaction with warming (p > 0.05 in both cases). Finally, forb species richness bore a significant, negative relationship to shrub litter across treatments (F1,240 = 17.1, p < 0.001; Wald's χ2 = 35.3, p < 0.001), and was highest on grazed/ambient plots, where litter biomass was lowest, and lowest on exclosed/warmed plots, where litter biomass was highest (figure 3b).
Just as understanding, the interplay between abiotic and biotic factors in community composition and dynamics is a central focus in ecology, understanding the interaction between environment and species interactions in community dynamics should be central to the study of climate change because of the implications both pose for species diversity and its influence on community stability. Most generally, and classically, community composition is regarded either as a consequence of species–environment interactions  or of species–species interactions . As this and other studies have demonstrated, species interactions are of fundamental importance to community response to climate change because they may, in the case of exploitation, stabilize interactions among competitors that are exploited, or, in the case of interference in the absence of exploitation, lead to competitive exclusion . This study highlights the added importance of the trophic dimension to understanding the contributions of environment and species interactions to community dynamics, as well as emphasizing the importance to plant community stability under climate change of maintaining intact populations of large herbivores.
Assistance with establishment and maintenance of the experiment and annual field sampling was provided by Jesper Bahrenscheer, Pernille Bøving, Sean Cahoon, Todd Costello, Nell Herrmann, Toke Høye, Syrena Johnson, Megan MacArthur, Christian Pedersen, Ieva Perkons, Mason Post, Taylor Rees, Chris Wilmers, Tyler Yenter and CH2MHill Polar Services. This study is a continuation of a project initially funded by the National Geographic Society's Committee for Research and Exploration, with additional funding from the Office of Polar Programs at the US National Science Foundation. Comments by Jedediah Brodie, Dan Doak, Jeff Kerby, Daniel Schindler, Os Schmitz, Mark Vellend, David Watts and two anonymous referees improved the manuscript.
- Received November 16, 2012.
- Accepted January 23, 2013.
- © 2013 The Author(s) Published by the Royal Society. All rights reserved.