Ecological succession and climate change are pushing tundra as well as Arctic and subarctic lowland plant communities toward increased woody vegetation cover. However, areas along the Hudson Bay Lowlands that are over-grazed by hyper-abundant lesser snow geese are experiencing drastic losses of grass, sedge, and woody cover. We assessed longterm changes in proportional ground cover and habitat patch characteristic at a subarctic supratidal marsh that was largely vacated by breeding snow geese over a decade ago. We found no evidence of habitat recovery. Rather, snow geese leave a legacy on the land that propagates degradation of habitat long after their direct removal of vegetation through foraging. Over a 35 year period, we documented a 46% reduction in graminoid cover and an 84% reduction in shrub cover, which led to smaller and more isolated patches of shrubs that many avian species depend upon for foraging and nesting. Recent experimental goose exclosures suggest that recovery of degraded habitat is possible, but habitat management at a large scale will require drastic reductions in lesser snow goose numbers.
Introduction
Changes in vegetation assemblages can occur via long-term ecological succession (Clements, 1916; Gleason, 1926), climatic forcing (Adler at al., 2006), competition for limiting resources (Tilman, 1982), and through herbivory and over-grazing (Bazely and Jefferies, 1986). Following retreat of the Laurentide ice sheet ∼8000–9000 years ago, the lowlands surrounding Hudson Bay were inundated by the historical Tyrrell Sea (McAndrews et al., 1982). Over geological time, isostatic uplift has exposed new shoreline areas. As soils accumulated on the exposed saline substrate, the general succession of vegetation at a given geographic location transitioned from saltwater marsh to freshwater marsh, shrubby marsh, treed fen (Larix and Picea spp.), and eventually treed peatland (Tarnocai, 1982; Klinger and Short, 1996; Glaser et al., 2004). Elsewhere in the Arctic and subarctic, recent evidence suggests that climate warming is leading to a rapid increase in plant growth, decreased bare ground, and increased shrubification of tundra plant communities (Elmendorf et al., 2012). Two different mechanisms are thus pushing tundra as well as Arctic and subarctic lowland plant communities toward increased woody vegetation cover.
Along areas of the Hudson Bay Lowlands used by lesser snow geese (Chen caerulescens caerulescens; hereafter snow geese), a very different pattern in vegetation change has been observed over the past 30–40 years. Driven by agricultural modifications and associated nutrient subsidies on their wintering grounds and migratory stopover areas (Boyd et al., 1982; Abraham et al., 2005b), the mid-continent population of snow geese has increased at an annual rate of 5–14% (Alisauskas et al., 2011). The late Robert L. Jefferies and colleagues clearly documented that overgrazing and grubbing (of roots and rhizomes) of vegetation by hyper-abundant snow geese has led to increased barren ground and extensive loss of graminoids and shrubs (e.g., Jefferies and Rockwell, 2002; Abraham et al., 2005a). Between 1973 and 1993, snow goose foraging led to the loss of 2454 ha of coastal marsh habitat at the long-term La Pérouse Bay study area (Jano et al., 1998) and >35,000 ha at other areas along the Hudson Bay Lowlands (Jefferies et al., 2006).
Severe snow goose grubbing removes the insulating layer of graminoid swards and exposes dark soils. The low albedo of exposed dark soils results in high soil temperatures and increased evapotranspiration, which in turn allows salts from deeper soils of the historical Tyrrell Sea to be drawn and concentrated at the surface. Increased soil salinity in turn leads to premature leaf drop in the dwarf shrubs (i.e. Salix and Betula spp.), resulting in higher plant mortality (Iacobelli and Jefferies, 1991). In following years these areas hold less snow and are the first to be exposed during spring snowmelt, with further grubbing by pre-laying snow geese exacerbating the problem. Past a threshold, hyper-salinity of the soil seemingly pushes the ecosystem to an alternative state (Hik et al., 1992; Srivastava and Jefferies, 1995) and affects its ability to recover (Handa et al., 2002). As the snow geese degrade their habitat, they disperse to healthier vegetated sites within the greater Cape Churchill region (Cooch et al., 2001), leaving behind increased areas of barren ground (Abraham et al., 2005b). Unfortunately, the widespread loss of vegetation has resulted in marginal habitat for other ground-nesting avian species that rely on intact and contiguous patches of grasses and shrubs for nesting cover and foraging (Rockwell et al., 2003, 2009).
Our objective here is to assess the transition of vegetation states at the La Pérouse Bay study area over 35 years, and determine whether vegetation conditions have improved or continued to deteriorate since the departure of the historically abundant snow goose colony that eventually sought out better quality habitats (although some geese still use the area). Given the experimental findings of Jefferies and colleagues, we might expect that the legacy snow geese have on soil dynamics will have inhibited the ability of graminoids and shrubs to recover from historical degradation (Hik et al., 1992; Srivastava and Jefferies, 1995; Handa et al., 2002). Within goose exclosures placed on degraded habitat, however, some recovery of graminoid vegetation has recently been observed (Abraham etal., 2012). Facilitated by climate change (Elmendorf etal., 2012) and deposition of soils during spring runoff, we alternatively hypothesize that graminoids (but not shrubs) have had enough time to begin initial recovery since last assessed at a large scale in 1999 when most snow geese were dispersing out of the La Pérouse Bay study area.
Methods
STUDY AREA
Study plots were located on coastal supratidal marsh habitat near La Pérouse Bay, approximately 30 km east of Churchill, Manitoba, Canada (58°52.3′N, 93°41.0′W), which is part of the western Hudson Bay Lowlands and now within the northern boundary of Canada's Wapusk National Park (Fig. 1). Vegetation of the study area is characterized by dwarf shrub species (i.e. Salix, Betula spp.), and the saltmarsh grass (e.g. Puccinellia phryganodes) and sedge species (e.g. Carex subspathacea) that snow geese prefer to forage on by grubbing and shoot pulling (Jefferies et al., 2003). With the increased foraging by snow geese on these graminoids, larger extents of hypersaline soils are now more common throughout the general area. For more details on soil and vegetation interactions of the study area, see Iacobelli and Jefferies (1991), Srivastava and Jefferies (1995), and Jefferies and Rockwell (2002).
VEGETATION SAMPLING
Classification of vegetation and habitat condition was conducted on 5 study plots, which were set up in a grid system of 50 m2 cells and included portions of a 1976 study site (7 ha) established to investigate the relationship between mating systems of Savannah sparrows (Passerculus sandwichensis) and habitat quality (Weatherhead, 1979). In addition to this long-term plot, four nearby study plots were established in 1999, representing heavily degraded habitats (2 plots, one of 3 ha and the other 5 ha) and marginally intact habitats (2 plots, one of 3 ha and the other 5 ha) that were geographically adjacent and representative of the habitats surrounding La Pérouse Bay.
Vegetation sampling was conducted on the long-term plot in 1976, and on all 5 study plots in 1999 and 2010. Duplicating the same sampling procedure used by Weatherhead (1979) and Rockwell et al. (2003), a total of 28 grid cells in 1976 and 92 grid cells in 1999 and 2010 were surveyed during the summer (see Peterson, 2012). A modified step-point method was used to classify the type of vegetation and habitat condition underfoot at each ∼ 1 m pace, along two diagonal transects within each grid cell (Evans and Love, 1957; Owensby, 1973; Rockwell et al., 2003; Abraham et al., 2005a). In 1999 and 2010, all shrubs, grasses, sedges, and forbs were identified to genus and some to species level. Shrub height was also recorded and condition was classified as either dead (100% dead branches), in poor condition (> 1/3 dead branches), or healthy and alive (<1/3 dead branches). Habitat status in 1976 fell into one of only six categories: barren ground, pond, sedge-short grass, mixed grass-short willow, Elymus (Leymus arenarius), and other Salix spp. (Weatherhead, 1979, Rockwell et al., 2003).
The objective of vegetation surveys in 1999 and 2010 was to assess whether habitat conditions had remained the same, improved, or degraded further on all plots. Given the differential detail of vegetation recording over time, habitat status was reduced into three classes for analysis: barren ground, graminoid cover, and shrub cover–a level of detail suitable for defining habitat quality for ground-nesting avian species on the study plots that could be affected by snow goose habitat degradation. The barren ground class consisted of all bare soils, ephemeral ponds and streambeds, algal mats, mosses, and forbs (i.e. Salicornia borealis, Senecio congestus, Atriplex glabriuscula), indicative of disturbance and hypersaline soils, and completely dead shrubs with no growth or ground cover. All grass and sedge species were included in the “graminoid cover” class, and the shrub cover class included all living or partially living shrubs (i.e. Salix spp., Betula glandulosa, and Myrica gale) (Rockwell et al., 2003).
DATA ANALYSIS
The total number of paces from both diagonal transects within each grid cell were summed, and the proportion of each habitat class was then calculated. These proportional data were then transformed using the logit transform: log(y/[1 -y]), with the addition of a small value (+0.0001) to the denominator and numerator of the equation for graminoid and shrub cover proportion (due to some values close to or equal to zero), and subtraction of a small value (-0.0002) from the proportion of barren ground (due to some values close to 1.0). This ensures that the transformed data does not include undefined values (-∞ and ∞) and that proportional estimates fall between 0 and 1. Given the monotonic nature of the logit transform, estimated coefficients are naturally interpreted within a linear model (as opposed to an arcsine transformation; Warton and Hui, 2011).
To determine if there were significant changes in proportion of the three habitat classes (barren ground, graminoid cover, and shrub cover) over time, and to determine the proper fixed effect (i.e. year) model structure for later analyses, an initial multivariate analysis of variance (MANOVA) was run in the program R (using the vector of proportional cover amongst the 3 habitat classes as the multivariate response). A comparison of habitat changes among 1976, 1999, and 2010 was conducted first on the long-term study plot (n = 28), since the paired study plots had not been established prior to 1999. All plots were then analyzed together as one single study area (n = 92) between 1999 and 2010.
Three separate model structures were compared to one another: a null model, a time trend model, and a model with year treated as a factor. To determine the best model structure for temporal change in habitat cover, we used Akaike's Information Criterion adjusted for sample size (AICc; Akaike, 1973; Burnham and Anderson, 2002). Based on the best-performing MANOVA model structure for temporal change in habitat cover, a generalized linear mixed effects model (GLMM) was developed using the lme4 package in R (with a Gaussian distribution and identity link; Bates et al., 2011). In each GLMM, year was treated as a fixed effect and grid cell nested within study plot was treated as a random effect. This accounted for grid cell and study plot heterogeneity in the model structure but allowed us to focus our attention on temporal change in vegetation cover and habitat conditions, which was of primary interest (Bolker, 2008).
To adequately investigate avian ground nesting habitat quality on all study plots, we also measured an index of mean patch size of all three habitat classes, estimated from the average number of sequential paces within a habitat class along both diagonal transects in a grid cell. Similarly, we also estimated the mean distance between shrub patches, measured from the average summed number of both graminoid and barren ground paces between points where live shrubs were recorded. This parameter was essential in addressing the connectivity of quality nesting habitat “islands,” where the shrub component is critical for many ground-nesting birds. Analysis of habitat patch size and distance between shrub patches was similar to the analyses of habitat cover described above; however, we log-transformed the mean step length data to meet assumptions of normality (Kutner et al., 2004). The Delta method (Seber, 1982) was used to calculate the variance and 95% confidence interval for each habitat parameter (e.g., proportion cover) derived from the statistical analyses of transformed data.
Results
LONG-TERM STUDY PLOT (1976, 1999, 2010)
Proportional Cover
The initial MANOVA analysis indicated that treating year as a factor was the best parameterization of temporal change in proportional cover of the three habitat classes (barren ground, shrub cover, and graminoid cover) on the long-term study plot (ΔAICc = 23.69 for year treated continuously and 294.68 for the null model). After implementing these findings into a GLMM to account for random variation amongst the grid cells (i.e., a model with a fixed-year effect and a random effect for grid cells), we found that the proportion of barren ground increased over the course of the study (by a multiple of 2.6 between 1976 and 1999, and by 20% between 1999 and 2010), while the proportions of graminoid and shrub cover both decreased (by 29% and 65% between 1976 and 1999, respectively, and by 24% and 53% between 1999 and 2010, respectively; Fig. 2 and Appendix Table A1).
Patch Characteristics
The proportional increase of barren ground and related decrease in graminoid and shrub cover over the long-term study occurred in tandem with the fragmentation of habitat. Initial MANOVA analysis of the mean patch size of each habitat class, and distance between “islands” of shrub cover yielded a top model structure with year treated as a factor (ΔAICc = 128.54 for year treated continuously and 440.07 for the null model). These results were then used to develop a GLMM, such as the one described above for the analysis of habitat class proportions. The index to mean patch size of barren ground increased by ∼1.5 m between 1976 and 1999, but increased fourfold (∼6.6 m) between 1999 and 2010. In parallel, the mean patch size of shrub cover decreased by approximately 6.5 m between 1976 and 1999, but only by 0.3 m between 1999 to 2010. The mean patch size of graminoids decreased by approximately 2.5 m between 1976 and 1999, with little to no change between 1999 and 2010 (Fig. 3, parta). As this habitat fragmentation occurred, the distance between islands of shrub patches actually decreased by ∼1.3 m between 1976 and 1999. Between 1999 and 2010, however, the fragmentation had become so extreme that the mean distance between islands of shrub patches increased by 5.1 m (Fig. 3, part b, and Appendix Table A2).
ALL STUDY PLOTS (1999, 2010)
Proportional Cover
Analyses across all study plots in 1999 and 2010 generally confirmed the long-term analyses at a larger scale. Given the two points in time, initial MANOVA analysis of temporal changes in proportional habitat cover gave equal weight to treating year as a trend or as a factor (ΔAICc = 71.87 for the null model). For consistency, we estimated changes in proportional cover using a GLMM with year treated as a factor and parameterized a random effect for grid cells nested within study plots. The proportion of barren ground across the entire study area significantly increased between 1999 and 2010 (by 31%) in parallel with a sharp decrease in both the graminoid and shrub cover proportions (37% and 61%, respectively; Appendix Table A3).
Patch Characteristics
Initial MANOVA analysis of the habitat patch parameters also showed that year treated as a trend or as a factor received equal support (ΔAICc = 53.10 for the null model). Using the same GLMM structure as that used for the proportional cover analysis, mean patch size of barren ground increased by 7.5 m across the 11-yr period. In parallel, both graminoid and shrub patch size decreased (by 0.39 m and 0.55 m, respectively), albeit in smaller increments than patches of barren ground (Appendix Table A4). The average distance between islands of shrub patches increased from 6.8 to 13.6 m, doubling the isolation distance (Table A4).
Discussion
We found that snow geese are driving change in vegetation assemblages along the Hudson Bay Lowlands that opposes the mechanisms of climate change and ecological succession (Tarnocai, 1982; Elmendorf et al., 2012). Although experimental goose exclosures on degraded habitat indicate that some recovery of graminoids is possible (Abraham et al., 2012), we found that largescale recovery has yet to occur. In fact, snow geese have left a legacy on the soil and vegetation that lasts at least a decade after taking their hungry bellies to greener pastures. Since 1999 when most breeding snow geese had departed the long-term La Pérouse Bay study area (Cooch et al., 2001), we found additional loss of graminoid and shrub cover and associated increases in barren ground that offer little habitat to endemic wildlife (Rockwell et al., 2003; Milakovic et al., 2003; Milakovic and Jefferies 2003).
FIGURE 3.
(a) Indices of mean patch size (m) for the ground cover categories (barren ground, graminoid cover, and shrub cover) on the long-term study plot in 1976 (dark gray), 1999 (light gray), and 2010 (medium gray). (b) Mean estimates of distances between shrub patches in 1976, 1999, and 2010.
The shrub assemblage that is important for breeding passerine and shorebirds declined most abruptly over time. On the long-term study plot alone, over 84% of shrub cover was lost between 1976 and 2010 (Fig. 2), and 61% was lost between 1999 and 2010 across all study plots (Table A3). Although graminoid cover declined on all plots over the study, it did not decline as rapidly as shrub cover. This may be because many of the graminoids are more tolerant of saline conditions than the shrubs (Jefferies et al., 1979; Handa et al., 2002; Rockwell et al., 2003). Persistence of graminoid patches suggests that some revegetation may be possible in locations with suitable freshwater conditions, unconsolidated sediment, and vegetative fragments (i.e., cover of residual graminoids and/or dead or dying shrubs), which can facilitate nitrogen fixing such that clonal propagation may take place (Handa and Jefferies, 2000). The potential for regrowth of graminoid assemblages has been observed within experimental exclosures on the study area, but to manage habitat at a large scale, snow goose numbers will need to be drastically reduced or deterred from using important habitats (Abraham et al., 2012).
The simultaneous loss of graminoids and shrubs led to an obvious increase in barren ground (see Fig. 2). Unfortunately, as patch size of barren ground increases, soil water salinity increases, creating a feedback that subsequently limits recolonization by shrub and graminoid assemblages (Srivastava and Jefferies, 1996). Thus, continued growth of barren ground patches will likely exacerbate in the near term even when snow geese are completely absent from an area due to the salinity and erosion mechanisms explained above. In the long term, however, invasion of these barren grounds by ruderal forbs adapted to saline soils (i.e. Salicornia borealis, Senecio congestus, Atriplex glabriuscula; Handa et al., 2002) could eventually leach some of the salts, and their decomposition could eventually create soils amenable to graminoid and shrub growth. In the meantime, we found that increased patch size of barren ground reduced the connectivity of shrub cover (Fig. 3, part b), which could negatively affect the breeding success of many groundnesting birds in the study area, especially the Savannah sparrow that uses low-lying dwarf shrubs in place of tall grasses as nesting cover in the northern limit of its range (Wheelwright and Rising, 2008). Rockwell et al. (2003) showed a 77% decline in average nesting densities of Savannah sparrows on the long-term study plot between 1976 and 1999, which was attributed to the extensive loss of shrubs during the same period. Given observed changes in the states of shrub cover and distance between shrub patches, we predict that these habitat metrics will continue to change in the near future at the rates most recently observed.
Continued monitoring of vegetation states and transitions (Bestelmeyer et al., 2009; Rumpff et al., 2011) will nevertheless be crucial to learning about the ability of the Hudson Bay Lowlands to recover if snow geese can be effectively managed (Abraham et al., 2012). If not, the expansion of habitat degradation to a larger scale across the coast of western Hudson Bay is likely and would have deleterious impacts on not only the plant community, but also on invertebrate and vertebrate biodiversity (Jefferies et al., 2006). Given these potential threats to lowland ecosystems in the north, we agree with the recommendations set forth by Abraham et al. (2012): (1) ground evaluations of vegetation changes should be conducted along the entire western coast and southern portion of Hudson Bay; (2) assessment of the potential for vegetation assemblages to recover from snow goose degradation should continue; (3) areas of highest potential for expansion of breeding snow goose populations should be identified; and (4) efforts should be invested in estimating the spatial carrying capacity for sustaining snow geese in North America, such that managers can judge the full extent to which Arctic and subarctic ecosystems might be vulnerable to snow goose degradation.
Acknowledgments
We would like to thank the Central and Mississippi Flyways, the Arctic Goose Joint Venture, the American Museum of Natural History Chapman Fund, Wapusk National Park, and Utah State University for funding. Koons was supported by NSF (DEB 1019613). We thank the late R. L. Jefferies for his contributions to our understanding of goose-plant interactions, P. J. Weatherhead for beginning this study, and D. Iles, who provided critical assistance in the field and with Figure 1.
References Cited
Appendices
APPENDIX
TABLE A1
Linear mixed model results for changes in proportional cover of barren ground, graminoids, and shrubs on the long-term study plot across 1976,1999, and 2010 on the logit-transformed scale. In all cases, 1976 represents the intercept, and coefficients for other years are estimated relative to the intercept. Proportional estimates and 95% confidence intervals are also shown. The grid-cell random effect shows the variance (Var) amongst the cells in the study plot on the logit scale (all rounded to the second decimal). The formula used in the lme4 package was: habitat class ∼ as.factor(year) + (1|cell_id).
TABLE A2
Liner mixed model results for changes in habitat patch size of barren ground, graminoid, and shrub classes, and distance between shrub patches on the long-term study plot across 1976, 1999, and 2010 on the log-transformed scale. In all cases, 1976 represents the intercept and coefficients for other years are estimated relative to the intercept. Habitat patch size estimates and 95% confidence intervals are also shown. The grid-cell random effect shows the variance (Var) amongst the cells in the study plot on the log scale (all rounded to the second decimal). The formula used in the lme4 package was: habitat patch size (or distance) ∼ as.factor(year) + (1|cell_id).
TABLE A3
Linear mixed model results for changes in proportional cover of barren ground, graminoids, and shrubs across all study plots between 1999 and 2010 on the logit-transformed scale. In all cases 1999 represents the intercept and the coefficient for 2010 is estimated relative to the intercept. Proportional estimates and 95% confidence intervals are also shown. The grid-cell random effect shows the estimated variance (Var) amongst grid cells in a given study plot on the logit scale (all rounded to the second decimal). The formula used in the lme4 package was: habitat class ∼ as.factor(year) + (plot|cell_id).
TABLE A4
Linear mixed model results for changes in habitat patch size of barren ground, graminoid, and shrub cover, and distance between shrub patches across all study plots between 1999 and 2010 on the log-transformed scale. In all cases, 1999 represents the intercept, and the coefficient for 2010 is estimated relative to the intercept. Habitat patch size estimates and 95% confidence intervals are also shown. The grid-cell random effect shows the estimated variance (Var) amongst grid cells in a given study plot on the log scale (all rounded to the second decimal). The formula used in the lme4 package was: habitat patch size (or distance) ∼ as.factor(year) + (plot|cell_id).
