Response of forage plants to alteration of temperature and spring thaw date: implications for geese in a warming Arctic

Abstract

Changes in summer temperatures in Arctic Alaska have led to longer and warmer growing seasons over the last three decades. Corresponding with these changes in climate, the abundance and distributions of geese have increased and expanded over the same period. We used an experimental approach to assess the response of goose forage plants to simulated environmental change. We subjected Carex subspathacea, a preferred goose forage growing on the Arctic Coastal Plain (ACP) of Alaska, to manipulations of temperature and timing of spring thaw to measure potential effects in terms of plant nitrogen concentration, aboveground biomass, and total nitrogen availability. Carex subspathacea responded to warming in a dynamic fashion. Increases in temperature led to decreases in leaf nitrogen concentration but increases in aboveground biomass. The increase in biomass was stronger than the decline in nitrogen concentration such that total nitrogen availability was increased with temperature for the first 35–40 d of the season. Grazing removal accounted for only minimal offtake of biomass, and we found no indication that grazing maintained elevated levels of nitrogen concentration longer in the season as reported in other studies. Based on demonstrated relationships in the literature between forage nitrogen concentrations and gosling growth rates, we conclude that there is currently abundant high‐quality forage available across the ACP. This finding fits with recent evidence of high gosling growth rates and increasing trends in goose abundance on the ACP. Our results suggest that with climate warming of a few degrees, nitrogen concentration of forage may decrease, but forage biomass and total nitrogen availability will increase. Our data suggest that nitrogen concentration will not fall below the minimum threshold required by geese in the near future. As such, we suggest that there is currently no bottom‐up limitation to goose numbers on the ACP.

Introduction

The effects of climate change are expected to be magnified in polar regions compared with temperate latitudes (Serreze and Francis 2006). Average summer temperatures in Arctic Alaska have increased approximately 2°C over the last 30 yr and are predicted to warm by an additional 2°C over the next 50 yr, leading to longer and warmer growing seasons (Hinzman et al. 2005). These changes have already resulted in a complex series of changes in biological processes (Hinzman et al. 2005, Van Hemert et al. 2015). Rapid permafrost degradation is leading to changes in foraging habitats, and, simultaneously, goose populations are changing in abundance and distribution (Jorgensen et al. 2006, Flint et al. 2008, Tape et al. 2013, Amundson et al. 2019).

Arctic-nesting geese, and their young in particular, require high-quality forage to meet nutritional demands for rapid growth in the short breeding season. On the Arctic Coastal Plain (ACP) of northern Alaska, geese appear to have benefited from recent climate warming (Hupp et al. 2017) and all species have shown increases in abundance corresponding with the period of temperature increase (Amundson et al. 2019). Further, this pattern appears to be supported at the broader scale across the circumpolar Arctic (Jensen et al. 2008, Madsen et al. 2011). Previous studies have concluded that inadequate forage during the brood-rearing period is a constraint on population dynamics in geese (Sedinger et al. 2016). However, the observations of high rates of gosling growth in the Arctic (Sedinger et al. 2001, Hupp et al. 2017) imply that forage is not limiting and support the idea that these increases in population size are the result of bottom-up drivers of population processes (Van Hemert et al. 2015).

The specific response of plants to warming appears to be species and location dependent (Doiron et al. 2014). However, multiple studies support a broad general pattern where increases in temperature result in an increase in plant biomass and a reduction in plant nutrient concentration (Chapin et al. 1995, Arft et al. 1999, Doiron et al. 2014, Lameris et al. 2017). This matters because, with their simple digestive system, geese require forage with relatively high nitrogen concentration (Sedinger 1997, Sedinger et al. 2016). Whereas recent warming in the Arctic is associated with an increase in plant primary productivity (Bhatt et al. 2010, Epstein et al. 2012), there likely exists a temperature-induced threshold beyond which the foliar ratio of carbon to nitrogen exceeds what is suitable for goose digestibility. The relationship between warming and overall plant response as it relates to goose populations is complicated because timing of herbivory has been shown to have a larger influence on plant nutrient concentration than timing of spring thaw (Post and Pedersen 2008, Olofsson et al. 2009, Beard et al. 2019).

The primary goal of this study was to assess the response of coastal salt marsh habitats to changes in the timing of spring thaw and temperature to better understand expected future conditions for Arctic-nesting geese. We used an in situ field experimental approach to manipulate the temperature and timing of spring thaw on monotypic stands of Carex subspathacea (hereafter subspathacea), an important goose forage plant in Arctic Alaska. We increased temperature via open-top-chamber greenhouses and decreased temperature with shade fabric that reduced solar radiation. We also modified spring thaw timing through snow manipulations by removing snow from plots to advance the growing season and adding snow to plots to delay soil thaw and associated green up. We quantified seasonal patterns in the relationships between temperature and spring thaw treatments relative to aboveground biomass, percent nitrogen, and available nitrogen (biomass × %N). We also assessed cumulative impacts of consistent warming and cooling treatments on the same plots over three years.

Methods Study site

We conducted our study from mid-May to mid-August, 2011–2013, in the Smith River estuary near the coast of the Beaufort Sea in the Teshekpuk Lake Special Area (TLSA) on the Arctic Coastal Plain (ACP) of Alaska (70.9° N, 153.2° W; Fig. 1). The TLSA annually supports 75,000–100,000 geese of four species (black brant Branta bernicla nigricans, greater white-fronted geese Anser albifrons, lesser snow geese Anser caerulescens carulescens, and cackling geese Branta hutchinsii) during the flightless wing molt, making it the most significant goose molting area in the circumpolar Arctic (Flint et al. 2008). The TLSA also supports increasing numbers of breeding geese (Amundson et al. 2019). Distributions of molting black brant have shifted from primarily inland, freshwater lakes toward coastal, brackish wetlands in recent decades (Lewis et al. 2011, Flint et al. 2014).

The Smith River estuary, which we consider to be representative of the coastal portion of the TLSA, is delimited by low-lying areas that are <1 m above sea level and dominated primarily by subspathacea and Puccinellia phryganodes; both species typically have high nitrogen content and are preferred forage for geese (Gadallah and Jefferies 1995, Person et al. 1998, 2003). These halophytic meadows abut Lowland Moist or Wet Sedge Tundra (Jorgenson and Heiner 2003) at least a meter higher in elevation, which is dominated by sedges and underlain by characteristic ice-wedge polygonal ground and ice-rich permafrost (Kanevskiy et al. 2013). Grazing lawn habitats along the Smith River and other estuarine environments on the ACP are expanding due to inundation, sedimentation, and subsidence (Tape et al. 2013). Virtually unused by molting black brant in the 1970–1980s, the Smith River estuary now supports more molting and post-molt black brant than most other wetland complexes in the TLSA, and increasing numbers of molting and brood-rearing white-fronted geese and snow geese (Flint et al. 2008, Lewis et al. 2010, Lewis et al. 2011). The climate in this area is characterized by a short growing season (June–August), high variation in daily temperature during spring and summer (range = −3.2° to 23.4°C, data from this study), and low average annual precipitation (˜140 mm; Gallant et al. 1995).

Experimental manipulations

We established two experimental study designs to examine the effects of in situ manipulation of environmental conditions on forage plant nutrient composition and production. The first design (temperature-thaw) involved temperature and timing of spring thaw manipulations paired with weekly vegetation sampling in which new plots were selected annually. The basic goal of the temperature-thaw manipulations was to simulate historic conditions (i.e., later thaw, colder summers) and future conditions (i.e., earlier thaw, warmer summers) to compare with current unmanipulated conditions. In the second design (multi-year), we manipulated temperature and maintained the established plots for three years, but restricted vegetation sampling to three occasions per year. The primary design allowed for detailed analysis of temperature and spring thaw effects on forage characteristics within a given year whereas the multi-year designed allowed us to assess the potential cumulative effects of temperature manipulations and restricted grazing.

The temperature-thaw study was a blocked design consisting of three temperature treatments (exclosed, greenhouse, and shadehouse) crossed with three spring thaw treatments (advanced, natural, and delayed) for a total of nine plots at each site (Fig. 1). Grazing effects were added to this design by establishing plots that were identified but otherwise unmanipulated and were available to be grazed by geese. Exclosed plots were fenced with 2-cm plastic mesh deer netting attached to a structure (120 × 120 × 60 cm) constructed of 1.25-cm pvc pipe; these structures were designed to exclude grazing without altering temperature or precipitation regimes. Greenhouse plots were surrounded with semi-clear fiberglass open-top-chamber greenhouses (OTCs) that have been used previously to increase ambient temperatures without impeding normal precipitation input (Henry and Molau 1997, Hollister and Webber 2000, Klein et al. 2005, Sullivan and Welker 2005). Our OTCs were conical in shape, measuring 68 cm in height with a diameter of 200 cm at the base and 120 cm at the top. Shadehouse plots were shaded with 10% red shade cloth (ChromatiNet) suspended 30 cm above ground by a structure constructed from 2-cm steel pipe with a footprint of 365 × 380 cm. Red shade cloth reduces the spectrum of blue, green, and yellow light and increases the red and far-red light spectrum, simulating solar radiation typical of an overcast day (Stamps 2009). Thus, this approach was intended to decrease solar input and heating without restricting precipitation input. Spring thaw treatments were conducted by manipulating snow depth associated with plots during late spring, approximately two weeks prior to breakup and before considerable snowmelt had occurred on the landscape. For advanced spring thaw plots, we removed snow from an area of 5 × 5 m area to promote soil thaw and expose the vegetated surface to sunlight. For delayed spring thaw plots, we piled snow to a height of 1 m over an area of 5 × 5 m area in order to slow the rate of soil thaw and limit sunlight penetration. Snow depths were unmanipulated in natural spring thaw plots.

We used data loggers (HOBO U-12, Onset) to record soil and air temperature in all plots on 15-min intervals. Soil temperature sensors were positioned at a depth of 6 cm below the soil surface; air temperature sensors were located inside vented solar radiation shields and positioned at a height of 10 cm above the soil surface. Temperature sensors and temperature treatment structures were put in place when snowmelt had exposed surface vegetation, but prior to substantial soil thaw. Large-scale flooding associated with spring breakup delayed the placement of treatment structures and temperature sensors for some sites.

Plots measured 120 × 120 cm and were sub-divided into sub-plots measuring 20 × 20 cm; these sections were separated from each other by a 6 cm buffer (Fig. 1). On approximately 7-d intervals, beginning at the onset of the growing season, we randomly selected a section from each plot for vegetation sampling. We removed the 20 × 20 cm section to a depth of 15 cm. Turves from control plots were examined for evidence of grazing, and grazing intensity was indicated as none, light, moderate, or heavy. All aboveground biomass was then clipped from the turve using scissors, rinsed in freshwater, and transported to camp where samples were dried at 50°C for 48 h and stored for later processing. Turves were returned to the substrate within their respective plots after vegetation sampling. Following completion of the field season, we sorted vegetation samples to separate live from dead material, dried samples to a constant mass at 50°C, and then weighed samples on a digital scale to the nearest 0.001 g. A sub-sample (0.5–1.0 g) of live material from each sample was analyzed for percent nitrogen using a C-N analyzer at the University of Alaska Anchorage stable isotope laboratory. We replicated the experiment at two sites in 2011 and at three sites in 2012–2013 with new sites and plots being selected each year.

For the multi-year design, we established six sites, each consisting of three plots (one plot for each temperature treatment) that were maintained over the 3-yr study period. Temperature sensors and treatment structures were put in place when snowmelt had exposed surface vegetation, but prior to substantial soil thaw. Sampling methods for multi-year plots were the same as described above, except that multi-year plots were sampled only three times per year. All data from this study are available in Meixell and Flint (2021).

Statistical analysis Overall approach

Our overarching objective was to quantify the effects of our experimental manipulations on aspects of forage plants that are important to geese. As such, we measured percent nitrogen, aboveground live biomass, and total nitrogen availability. All three of these metrics are known to have non-linear seasonal patterns of variation under normal environmental conditions. While biomass increases in an asymptotic fashion throughout the growing season, percent nitrogen peaks early in the season and then declines. Total nitrogen availability is the net result of these processes and typically rises early in the season and may then stabilize or decline. These inherent non-linear patterns add significant complexity to assessing the effects of our manipulations. Our primary interest was in assessing deviation from these basic seasonal patterns as caused by our manipulations. Therefore, we used polynomial regression as fit to the entire data set for a given variable relative to date to describe the seasonal non-linear process and we calculated the residuals from these regressions. We then assessed patterns of variation within these residuals relative to year, spring thaw, and treatment effects. The fact that we used residuals as our dependent variable functionally linearized the data. The interpretation of the effect size from these analyses then becomes a deviation from the inherent non-linear seasonal pattern.

In our study design, we used alterations of the timing of spring thaw as well as temperature treatments (i.e., greenhouse and shadehouse) to simulate potential changes in climate conditions. While our treatments are clearly defined and applied consistently, the realized effect of a given treatment is uncontrolled. For example, while the greenhouse treatment generally raised the temperature within the structure, the actual amount of change was variable across sites, years, and days. We used a two-stage modeling design to analyze these data. First, we directly considered the effects of our spring thaw and temperature treatments on the residuals of percent nitrogen, live biomass, and total nitrogen availability. If a model with temperature treatment was supported, we then conducted a second stage of modeling that replaced the temperature treatment effect with an index of soil and/or air temperature. Because the realized temperature treatment effect was variable, there was potential for the continuous temperature variable to explain more variation than the treatment class variable. That is, the temperature variable could explain variation within and among treatment classes. As our indices to temperature (both soil and air), we calculated cumulative thaw degree days (TDD; Van Wijk et al. 2012) for each plot’s air and soil temperature with a threshold of 0°C for plant growth along the Arctic coast (Dennis et al. 1978, Miller et al. 1980). Because cumulative TDD is a continuous summed quantity, there was a non-linear exponential increase with date. Therefore, we used polynomial regression with date across all plots to estimate the overall function for soil and air TDD (separately) and calculated the residual TDD for each date and plot. We used these residual air TDD and soil TDD as our continuous indices to temperature.

Temperature

To assess the effectiveness of our experimental treatments, we calculated the difference between the TDD for manipulation treatments (i.e., greenhouse and shadehouse) and exclosure plots for each sampling day. If we assume that exclosures represent natural conditions at each location and time, the difference between greenhouses/shadehouses and exclosures represents the realized effect size of each manipulation treatment. We calculated the average of these metrics for each of the spring thaw treatments (advanced, natural, delayed) for both air and soil temperatures.

A small proportion of temperature sensors (11 of 204) logged sporadic erroneous readings, presumably caused by a poor connection between the sensor cable and the data logger. In most cases, erroneous readings were highly unrealistic (e.g., <−10°C or >40°C), spanned an interval of <3 h, and accounted for <10% of all temperature measurements over the season for a given plot. To rectify data containing erroneous readings, we used temperature data from a plot within the same site while matching spring thaw and treatment as closely as possible. We filtered temperature records from the plot in question if a given reading differed from the reference sensor by >2°C, then filled in the datapoints between existing readings linearly. In rare cases where erroneous temperature readings spanned >12 h and for sensors that failed entirely (n = 4), we substituted missing temperature data with those collected from a nearby plot with comparable spring thaw and treatment. All sensors with complete failure occurred in 2011; one was from a plot in the primary study, and three were from multi-year plots.

Nutrient dynamics

We examined variation in vegetation metrics (i.e., nitrogen concentration, biomass, total nitrogen) using a linear mixed models procedure in SAS with a maximum-likelihood approach to fit models and assessed relative support for candidate models using Akaike’s information criterion corrected for sample size (AIC_c; Burnham and Anderson 2002). Plots within a given site were likely not independent, so to control for intra-plot correlation we included year-site as a random effect in all models.

Temperature and spring thaw effects

To quantify the effects of temperature and spring thaw on plant nutrient dynamics, we limited the data set to plots that were ungrazed (i.e., exclosure, greenhouse, and shadehouse plots). The dependent variable was the deviation from the overall non-linear regression. We conducted separate analyses for each of three dependent variables: percent nitrogen, live biomass (g), and total available nitrogen. For live biomass and available nitrogen, we scaled data to g/m². Available nitrogen is the product of percent nitrogen and live biomass, representing the amount of nitrogen in a square meter. To account for likely sources of variation in nutrient dynamics and to limit the number of models considered, we conducted model selection in a two-stage approach. In stage 1, we considered a suite of models assessing support for variation in vegetation metrics relative to treatment and spring thaw manipulations. In stage 2, we replaced the class-level treatment variables and directly assessed the effect of soil and air temperature on vegetation metrics.

In stage 1, we sought to identify the top approximating model structure for spring thaw manipulations. We started with a base model that included year, day of year, treatment (treatment represented the classes of greenhouse, shadehouse, and exclosed from grazing), and interactions between date, year, and treatment. We then added spring thaw treatment as a three-class factor (advanced, natural, and delayed) with interactions. Because the goal of the spring thaw manipulations was to simulate variation in breakup timing, we added a model that substituted the continuous variable thaw date (thaw) for the spring thaw factor. Thaw was the day of year for which cumulative soil TDD exceeded 0.10. The goal of this phase was to assess support for a spring thaw effect and assess if that effect was functionally equivalent to alteration of spring thaw date. Stage 1 considered three models, and the best supported structure for spring thaw effect was used in stage 2.

Our objective in stage 2 was to quantify the effects of temperature treatments on vegetation characteristics (i.e., biomass and nitrogen concentration). We considered three alterations of the best model from stage 1. We first assessed the best model from stage 1 without any treatment effects. We then replaced the treatment effects in the stage-1 best model with our index to airTDD and soilTDD. We then considered models in which the effect of temperature varied with day of year.

We then considered one final model designed to demonstrate the degree of annual variation in our data. All models run previously in stages 1 and 2 included year effects and interactions. We took the best model from stage 2 and removed all year effects. The goal was simply to demonstrate the level of variation that was being explained by the categorical year variable, even after controlling for spring thaw and temperature effects.

Grazing effects

To quantify the effects of grazing on plant nitrogen dynamics and to estimate offtake of biomass, we conducted analyses using only exclosure and control plots. Modeling was conducted using the same approach and suite of models as described above for the temperature-effect analysis, the difference being that treatment represented a grazing effect with two classes: grazed and exclosed.

Multi-year effects

Multi-year plots were sampled on only three occasions per year, so we assessed seasonal variation using a three-level categorical variable (i.e., early, middle, and late). Because date was categorical in this analysis, we used the direct measurements as opposed to residuals as the dependent variables. The goal of the multi-year plots was to assess whether the effect size of the treatments (greenhouse/shadehouse) changed through time. As such, we were primarily interested in the interaction term of treatment and year, or the treatment proxy of TDD and year. There was no spring thaw manipulation associated with the multi-year plots. We modeled variation in forage nutrient dynamics on multi-year plots using a set of six models (three pairs, one with and one without interaction term) to assess the support for the interaction term.

Predicting effects of increased temperature

To demonstrate the effect size of temperature, we calculated the difference between the model predicted values of percent nitrogen, biomass, and total nitrogen using the best model with a 1°C increase in temperature on each day from day 1 to day 50 as compared to a model with a zero temperature deviation. Recall that our temperature index was calculated as a residual from the overall temperature relationship. Therefore, the estimate with no temperature effect would be the predicted value in an average or normal year. We then used a temperature index that deviated by one degree for each day. This would be the equivalent of a year where the temperature was one degree above average every day of the year. The difference between the two lines can be interpreted as the net effect of that temperature change. The realized relationships are non-linear because temperature was included in a 3-way interaction with year and date.

Results Temperature variation

The patterns of temperature varied substantially among years, spring thaw, and temperature treatments (greenhouse, shadehouse). The initial year of study (i.e., 2011) was relatively cold, particularly in the latter half of the season (Fig. 2). The opposite pattern was seen in 2012, as temperatures started relatively normal and then were above average in the latter part of the season ending with the highest TDD total. The pattern in 2013 showed a relatively consistent increase in temperature that was higher early in the season but ended between the previous years. There were stronger manipulation effects on relative temperature in 2013. The effect size of the greenhouses relative to the exclosures was greater than the shadehouses (Fig. 3). The spring thaw manipulation was ineffective in 2012. While the manipulations were applied in a consistent fashion each year, the realized effects were highly variable (Fig. 3).

Temperature effect

Percent nitrogen, biomass, and total nitrogen all showed non-linear patterns with date (Fig. 4). Percent nitrogen peaked about 18 d after initiation of sampling, then declined gradually for the remainder of the season. Biomass increased rapidly, then growth appeared to slow during the latter half of the season. Because total nitrogen availability is the product of percent nitrogen and biomass, total nitrogen increased to a peak at about 35 d followed by a gradual decline thereafter.

In the primary spring thaw and temperature analyses, the best models to describe the residual variation in all three metrics (i.e., percent nitrogen, biomass, and total nitrogen) showed similar structure (Table 1). All three metrics varied by year, spring thaw treatment, day of year, and temperature index. Only percent nitrogen supported substituting date of thaw for our spring thaw treatment. In all cases, the net effect of increased temperature was a decline in percent nitrogen (Fig. 5). In some cases, the decline was nearly linear throughout the season; in other cases, the nitrogen concentration declined initially, then stabilized. For biomass, the net effect of increased temperature was an increase in biomass (Fig. 6) that primarily occurred early in the season and peaked between 20 and 30 d after initiation of sampling. The net effect of increased temperature on total nitrogen availability was similar to that for biomass (Fig. 7) with an initial increase in total nitrogen availability that peaked about 20 d after initiation of sampling. In most years, this increase in nitrogen availability was diminished by day 40, after which increased temperature had a negative effect on nitrogen availability.

Table 1 Model selection results for percent nitrogen, aboveground biomass, and total nitrogen availability for the primary temperature and spring thaw manipulation experiment.

Model structure†	−2 LL‡	No. parameters§	AICc¶
Percent nitrogen
yr, treat, day, day × yr, day × treat, yr × treat	390.4	16	422.4
yr, spring, treat, day, day × yr, day × spring, day × treat, yr × spring × treat	340.4	36	412.4
yr, thaw, treat, day, day × yr, day × thaw, day × treat, thaw × yr × treat	365.2	22	409.2
yr, thaw, day, day × yr, day × thaw, thaw × yr	457.7	12	481.7
yr, thaw, treat, day, day × yr, day × thaw, day × treat, yr × thaw × treat	365.2	22	409.2
yr, thaw, air, day, day × yr, day × thaw, day × air, thaw × air × yr	376.2	15	406.2
yr, thaw, soil, day, day × yr, day × thaw, day × soil, thaw × soil × yr	387.6	15	417.6
yr, thaw, air, day, day × yr, day × thaw, day × air × yr, thaw × air × yr	375.6	17	409.6
yr, thaw, soil, day, day × yr, day × thaw, day × soil × yr, thaw × soil × yr	384.6	17	418.6
thaw, air, day, day × thaw, day × air, thaw × air	474.4	12	498.4
Biomass
yr, treat, day, day × yr, day × treat, yr × treat	4706.7	16	4738.7
yr, spring, treat, day, day × yr, day × spring, day × treat, yr × spring × treat	4601.0	36	4673.0
yr, thaw, treat, day, day × yr, day × thaw, day × treat, thaw × yr × treat	4671.2	22	4715.2
yr, spring, day, day × yr, day × spring, spring × yr	4685.6	16	4717.6
yr, spring, treat, day, day × yr, spring × day, day × treat, spring × yr × treat	4601.0	36	4673.0
yr, spring, air, day, day × yr, day × spring, day × air, spring × air × yr	4563.2	22	4607.2
yr, spring, soil, day, day × yr, day × spring, day × soil, spring × soil × yr	4589.0	22	4633.0
yr, spring, air, day, day × yr, day × spring, day × air × yr, spring × air × yr	4557.0	24	4605.0
yr, spring, soil, day, day × yr, day × spring, day × soil × yr, spring × soil × yr	4588.5	24	4636.5
spring, air, day, day × spring, day × air, spring × air	4617.3	16	4649.3
Total N
yr, treat, day, day × yr, day × treat, yr × treat	877.6	16	909.6
yr, spring, treat, day, day × yr, day × spring, day × treat, yr × spring × treat	777.8	36	849.8
yr, thaw, treat, day, day × yr, day × thaw, day × treat, thaw × yr × treat	839.3	22	883.3
yr, spring, day, day × yr, day × spring, spring × yr	853.8	16	885.8
yr, spring, treat, day, day × yr, day × spring, day × treat, yr × spring × treat	777.8	36	849.8
yr, spring, air, day, day × yr, day × spring, day × air, spring × air × yr	767.5	22	811.5
yr, spring, soil, day, day × yr, day × spring, day × soil, spring × soil × yr	781.2	22	825.2
yr, spring, air, day, day × yr, day × spring, day × air × yr, spring × air × yr	751.7	24	799.7
yr, spring, soil, day, day × yr, day × spring, day × soil × yr, spring × soil × yr	775.5	24	823.5
spring, air, day, day × spring, day × air spring × air	825.7	12	849.7

Notes:

All models contained the random effect of site-year. The dependent variables were the residual from a polynomial regression of each sample measurement against day of year. These dependent variables can be interpreted as the deviance of a given sample from the overall average, controlled for date. The model with the lowest AIC_c value for each dependent variable is noted in bold.

^†Yr = year of study with three classes: 2011, 2012, and 2013. Spring = spring thaw treatment where we attempted to alter timing of initiation of plant growth; this variable has three classes: advanced, natural, and delayed. Day = day of year of sampling. Treat = temperature treatment; this variable has three classes: greenhouse, shadehouse, and exclosed from grazing. Thaw = date when cumulative soil thaw degree days (TDD) exceeded 0.1. Air = residual from a polynomial regression of air TDD measured against day of year for all samples; this metric can be interpreted as the deviance of TDD for a given sample from the overall average, controlled for date. Soil = residual from a polynomial regression of soil TDD measured against day of year for all samples; this metric can be interpreted as the deviance of TDD for a given sample from the overall average, controlled for date.

^‡−2 log-likelihood of the model.

^§Number of parameters in the model including two parameters for the variance and the random effect.

^¶Small sample corrected Akaike information criterion.

Grazing effect

We detected no effect of grazing on percent nitrogen after controlling for year and date effects (Table 2). Thus, percent nitrogen followed the same pattern regardless of grazing treatment. For biomass, the best supported model included a treatment effect which was realized as a gradual offtake of biomass with date. Biomass increased as the season progressed, but the effect size of grazing was the removal of 0.46 (95% CI 0.65, 0.28) grams of dry mass per square meter per day in the grazed plots. Thus, the grazed plots increased in biomass more slowly than the exclosed plots. The total nitrogen availability results mirrored those for biomass.

Table 2 Model selection results for percent nitrogen, aboveground biomass, and total nitrogen availability for the multi-year temperature manipulation experiment.

Model structure†	−2 LL‡	No. parameters§	AICc¶
Percent nitrogen
yr, date, treat, yr × date, yr × treat	48.3	17	86.6
yr, date, treat, yr × date	49.4	13	77.8
yr, date, air, yr × date, yr × air	47.9	14	78.8
yr, date, air, yr × date	52	12	78.1
yr, date, soil, yr × date, yr × soil	52.2	14	83
yr, date, soil, yr × date	56.4	12	82.5
Biomass
yr, date, treat, yr × date, yr × treat	1431.3	17	1469.6
yr, date, treat, yr × date	1432.5	13	1460.9
yr, date, air, yr × date, yr × air	1431.3	14	1462.1
yr, date, air, yr × date	1435.7	12	1461.8
yr, date, soil, yr × date, yr × soil	1434.8	14	1465.7
yr, date, soil, yr × date	1439.8	12	1465.9
Total N
yr, date, treat, yr × date, yr × treat	290.9	17	329.2
yr, date, treat, yr × date	291.8	13	320.3
yr, date, air, yr × date, yr × air	286.9	14	317.8
yr, date, air, yr × date	293.2	12	319.3
yr, date, soil, yr × date, yr × soil	29.7	14	322.6
yr, date, soil, yr × date	294.5	12	320.6

Note:

The model with the lowest AIC_c value for each dependent variable is noted in bold.

^†Yr = year of study with three classes: 2011, 2012, and 2013. Treat = temperature treatment; this variable has three classes: greenhouse, shadehouse, and exclosed from grazing. Date = plots were sampled three time per year, this variable describes timing of sampling as a class with three categories: early, middle, and late. Air = residual from a polynomial regression of air thaw degree days (TDD) measured against day of year for all samples; this metric can be interpreted as the deviance of TDD for a given sample from the overall average, controlled for date. Soil = residual from a polynomial regression of soil TDD measured against day of year for all samples; this metric can be interpreted as the deviance of TDD for a given sample from the overall average, controlled for date.

^‡−2 log-likelihood of the model.

^§Number of parameters in the model including two parameters for the variance and the random effect.

^¶Small sample corrected Akaike information criterion.

Multi-year treatments

There was little change in percent nitrogen, biomass, and total nitrogen as measured across three years with consistent temperature manipulations (Fig. 8). There was only marginal evidence that the effect of our treatments varied across years for the multi-year plots. For all three vegetation metrics, models containing an interaction between year and temperature index received minor support (Table 3). These models indicated that the net effect of temperature on biomass and total nitrogen increased between the first and second years and remained elevated in the third year (Fig. 9).

Table 3 Model selection results for percent nitrogen, aboveground biomass, and total nitrogen availability for the grazing manipulation experiment.

Model structure†	−2 LL‡	No. parameters§	AICc¶
Percent nitrogen
yr, day, day × yr	222.5	8	238.5
yr, graze, day, day × yr, day × graze, yr × graze	220.7	12	244.7
yr, spring, graze, day, day × yr, day × spring, day × graze, spring × yr × graze	194.2	26	246.2
yr, thaw, graze, day, day × yr, day × thaw, day × graze, thaw × yr × graze	213.1	17	247.1
yr, spring, graze, air, day, day × yr, spring × day, day × graze, air × day, air × graze, spring × yr × graze	179.3	29	237.3
yr, spring, graze, soil, day, day × yr, day × spring, day × graze, soil × day, soil × graze, spring × yr × graze	183.7	29	241.7
spring, graze, air, day, day × spring, day × graze, air × day, air × graze, spring × graze	255.9	15	285.9
Biomass
yr, day, day × yr	787.7	8	803.7
yr, graze, day, day × yr, day × graze, yr × graze	766.4	12	790.4
yr, spring, graze, day, day × yr, day × spring, day × graze, spring × yr × graze	724.3	26	776.3
yr, thaw, graze, day, day × yr, day × thaw, day × graze, thaw × yr, × graze	760.4	17	794.4
yr, spring, graze, air, day, day × yr, day × spring, day × graze, air × day, air × graze, spring × yr × graze	678	29	736
yr, spring, graze, soil, day, day × yr, day × spring, day × graze, day × soil, soil × graze, spring × yr × graze	691.8	29	749.8
spring, graze, air, day, day × spring, day × graze, day × air, air × graze, spring × graze	742.6	15	772.6
Total N
yr, day, day × yr	1586.1	8	1602.1
yr, graze, day, day × yr, day × graze, yr × graze	1565.2	12	1589.2
yr, spring, graze, day, day × yr, day × spring, day × graze, spring × yr × graze	1508.0	26	1560.0
yr, thaw, graze, day, day × yr, day × thaw, day × graze, thaw × yr × graze	1553.6	17	1587.6
yr, spring, graze, air, day, day × yr, day × spring, day × graze, air × day, air × graze, spring × yr × graze	1470.6	29	1528.6
yr, spring, graze, soil, day, day × yr, day × spring, day × graze, soil × day, soil × graze, spring × yr × graze	1487.5	29	1545.5
spring, graze, air, day, day × spring, day × graze, air × day, air × graze, spring × graze	1530.4	15	1560.4

Notes:

The dependent variables were the residual from a polynomial regression of each sample measurement against day of year. These dependent variables can be interpreted as the deviance of a given sample from the overall average, controlled for date. The model with the lowest AIC_c value for each dependent variable is noted in bold.

^†Yr = year of study with three classes: 2011, 2012, and 2013. Spring = spring thaw treatment where we attempted to alter timing of initiation of plant growth; this variable has three classes: advanced, natural, and delayed. Day = day of year of sampling. Graze = grazing treatment; this variable has two classes: available for grazing and exclosed from grazing. Thaw = date when soil thaw degree days exceeded 0.1. Air = residual from a polynomial regression of air thaw degree days (TDD) measured against day of year for all samples; this metric can be interpreted as the deviance of TDD for a given sample from the overall average, controlled for date. Soil = residual from a polynomial regression of soil TDD measured against day of year for all samples; this metric can be interpreted as the deviance of TDD for a given sample from the overall average, controlled for date.

^‡−2 log-likelihood of the model.

^§Number of parameters in the model including two parameters for the variance and the random effect.

^¶Small sample corrected Akaike information criterion.

Discussion

Our study showed that predicted increases in temperature and growing season length in the Arctic are likely to result in increased biomass and a reduction in nitrogen concentration in goose forage plants. Our results showing increases in biomass and reductions in nitrogen concentration with increases in temperature match the results from previous work on other plant species in the Arctic (Chapin et al. 1995, Arft et al. 1999, Doiron et al. 2014, Lameris et al. 2017); however, we found that the increase in biomass was stronger than the reduction in nitrogen concentration such that total nitrogen availability on an area basis increased. The patterns were not consistent with date such that temperature-derived increases in biomass occurred early in the growing season, but then peaked and biomass at the end of the season was only marginally increased due to temperature. This pattern matches the results of Livensperger et al. (2016) who showed that advanced spring thaw and warming advanced the timing of primary productivity but did not necessarily increase total annual production. Total nitrogen availability was initially increased by warmer temperature, but then peaked and declined quickly. The initial increase in total nitrogen was driven by the increase in biomass, whereas the subsequent decline in total nitrogen was driven by the late season decline in nitrogen concentration. The net result was that total nitrogen availability was increased by warming for the first 35–40 d of the season, but warming had a negative effect on nitrogen availability for the remainder of the season. Thus, while warming increased nitrogen availability early in the season, it appeared to increase the rate of seasonal senescence such that aboveground nitrogen concentration declined below levels from unmanipulated plots late in the season. The patterns we observed in total nitrogen availability fit with results of experimental warming of other plant species in the Arctic (Chapin et al. 1995, Doiron et al. 2014, Lameris et al. 2017).

Multi-year effects

Our multi-year plots maintained consistent treatments across three growing seasons. It appears that the patterns we observed in single-year treatments were maintained in the multi-year plots. There was some evidence that the effects of temperature may be stronger after the first year but the support for a treatment by year interaction was equivocal. If there is a cumulative multi-year effect of warming, it appears to be primarily manifested as an increase in the effect of temperature on biomass. Other long-term warming studies found that positive warming effects on biomass only lasted for the first four years of treatment, possibly due to soil nutrient depletion (Arft et al. 1999). The relatively short duration of our experiment would not have the ability to detect such effects. Nitrogen concentration was relatively consistent across the three years, whereas biomass and total nitrogen tended to increase through time. There was an increase in the effect size of 1°C temperature increase after the first year. It is unclear whether this effect size increase represents a true multi-year effect or simply annual variation in effect size. That is, the observed increase may have been a function of the relatively lower temperature in 2011 compared with subsequent years. Ruess et al. (1997) concluded that nitrogen availability was ultimately regulated by belowground processes and the multi-year decline observed in other studies is likely a realization of belowground depletion. Our plots occurred on highly water-saturated soils, and as such, our manipulation of relatively small plots may not have influenced subsurface processes independent of the surrounding habitats as a result of lateral water and nutrient flow. As such, our plots may not reflect the true multi-year effect as would be applied to the entire area.

Grazing effects

We did not detect a positive feedback grazing effect as was demonstrated by Person et al. (1998) where grazing resulted in a net increase in nitrogen concentration of subspathacea compared with ungrazed plants on the same date. This was likely a consequence of the fact that we detected very little grazing pressure. In the study by Person et al. (1998), a high proportion of aboveground biomass was removed by grazing and this may have facilitated the response they observed. In our case, we did not detect substantial grazing pressure and most grazing appeared to occur late in the season after nitrogen content had already peaked. At the Colville River Delta, Alaska, USA (150 km east of our study area), Sedinger et al. (2001) and Hupp et al. (2017) found very little grazing offtake, similar to what we measured at the Smith River. Person et al. (1998) and Beard et al. (2019) argued that grazing response delayed senescence and maintained peak nitrogen content longer into the season than for ungrazed plots. It may be that a similar response does not occur after the peak and during the seasonal decline in nitrogen content (Beard et al. 2019). Alternatively, we may not have had sufficient grazing pressure to stimulate a detectable response in the plants. In our sample, we showed similar aboveground biomass in our grazed plots to what Person et al. (1998) observed in their ungrazed exclosures, further supporting that grazing pressure in our sample was minimal.

Annual variation

Considerable annual variation in all three plant metrics remained in our data even after controlling for temperature and spring thaw effects (Table 1). It is not clear what the source of this variation might be but there are a number of factors that likely play some role. Ruess et al. (2019) demonstrated that belowground processes are directly related to aboveground nutrient metrics for subspathacea. These belowground effects may be related to immediate conditions as well as carryover from previous years (Ruess et al. 2019). The soils beneath our plots were highly water saturated, which likely influences soil respiration and mineralization rates, but it is unknown if soil water content varied among years. Periodic storm surges likely altered soil salinity among years. Precipitation may have played some role, but our plots were within 0.5 m of sea level such that immediate precipitation appeared to have little direct effect on soil water content observed during sampling. Overwinter precipitation may have played a role through its influence on spring flooding which inundated our plots each spring. The high water saturation we observed also suggests that hydraulic connectivity with the surrounding area potentially reduced the effect of our relatively small manipulations on belowground processes associated with a specific treatment.

The use of greenhouse and shadehouse to alter temperatures was highly dependent on weather conditions. That is, clear days tended to increase the effect size of the treatments, particularly the greenhouses. There was substantial annual variation in temperature and weather patterns among years which likely influenced our results. The first two years of the study showed opposing patterns in terms of relative temperature trends within seasons with temperature decreasing in 2011 relative to the three-year average and increasing in 2012 relative to average (Fig. 2). Given that tissue nitrogen concentration peaks early in the season, the actual timing of temperature variation likely affects realized outcomes such that not all temperature increases may be equivalent in terms of their effects on plant nutrient dynamics. We conclude that factors beyond temperature and timing of spring thaw have a substantial influence on annual variation in plant biomass and nutrient content, but it is unknown what these processes may be.

Effect of vegetation changes relative to goose foraging

Under a climate warming scenario, our data suggest that goose forage plants will increase in biomass and decrease in nitrogen content. This largely agrees with previous temperature manipulation studies on other plant species in the Arctic (Chapin et al. 1995, Arft et al. 1999, Doiron et al. 2014, Lameris et al. 2017). In our study, the increase in biomass was stronger than the decrease in nitrogen concentration and this resulted in a net increase in total available nitrogen. Doiron et al. (2014) argued that an increase in available nitrogen would not compensate for reduction in percent nitrogen because goslings were already consuming as much as possible, and therefore, any reductions in nitrogen concentration should result in decreased rates of gosling growth. However, Sedinger et al. (2001) showed that black brant goslings on the ACP with access to forage of higher biomass, but lower nitrogen concentration, grew more quickly than goslings on the Yukon–Kuskokwim Delta. Our data, when compared to data from Person et al. (1998) collected on the same putative species on Yukon–Kuskokwim Delta, suggest little difference in nitrogen concentration, but a substantially higher biomass and total nitrogen availability in the Arctic. This pattern of higher forage availability in the Arctic agrees with a recent study that has continued to validate high gosling growth rates on the ACP (Hupp et al. 2017). The projected decline in nitrogen concentration resulting from a 1°C increase in temperature would still have our values for nitrogen concentration above the 2.7% that Sedinger et al. (2001) reported from areas further east on the ACP. As such, we suspect that predicted decreases in nitrogen concentration resulting from small increases in temperature are not likely to reduce gosling growth. Compared to Person et al. (1998), we found minimal offtake from foraging, implying that forage availability on the ACP is not limiting. Similarly, Sedinger et al. (2001) found high biomass availability and minimal offtake on the Colville River Delta in the mid-1990s and Hupp et al. (2017) showed that while goose abundance in the same area increased considerably since the Sedinger et al. (2001) study, forage availability and rates of gosling growth remained unchanged.

It is unknown what the lower threshold for nitrogen concentration is before a potential forage is functionally unavailable to geese, particularly goslings. On the Yukon–Kuskokwim Delta where black brant gosling growth is apparently limited by availability of subspathacea, there is an abundance of Carex ramenskii (Ruess et al. 1997, Person et al. 2003). Ruess et al. (1997) measured biomass of Carex ramenskii at >500 g/m² and nitrogen concentration of about 1.25%. Given that black brant on the YKD show low rates of gosling growth as a result of forage limitation, we can assume that the abundant Carex ramenskii is below the nitrogen concentration threshold required by goslings. Sedinger et al. (2001) showed high gosling growth rates in the Arctic where subspathacea was present at 80–100 g/m² and nitrogen content of 2–3% was available. Therefore, we conclude that the minimum threshold for nitrogen content occurs somewhere in the 1.25–2.00% range. Our overall average peak nitrogen concentration was 3.2%, and our manipulations suggest a decline of about 0.3% resulting from a 1°C increase in temperature. Thus, we predict that under a climate warming scenario of <4°C, forage quality will remain adequate for geese while biomass will increase. As such, forage availability would seem to able to support a larger goose population on the ACP than currently exists.

Multi-year projections

Our study area was located in the TLSA, which was given special area designation, in part, because of the high concentrations of molting geese (Flint et al. 2008). This area also has some of the highest breeding densities for geese across the ACP (Amundson et al. 2019). Thus, our study area was located in a region with relatively high numbers of geese, yet we measured minimal offtake due to grazing. Further, Tape et al. (2013) concluded that permafrost degradation and subsidence were likely increasing the areal extent of subspathacea grazing lawns along the coast of the ACP. Whereas plots of subspathacea that were exclosed from grazing on the Yukon–Kuskokwim Delta showed a 25-fold increase in biomass and 50% reduction in percent nitrogen after three years (Person et al. 2003), we found little change in biomass or nitrogen concentration following three years of grazing exclusion, even when combined with continuous greenhouse treatments. Thus, we did not find any evidence of conversion of subspathacea to a growth form similar to Carex ramenskii, as occurs on the Yukon–Kuskokwim Delta in western Alaska (Person et al. 2003). As such, the grazing system along the ACP appears to be far more stable than on the Yukon–Kuskokwim Delta and frequent grazing is not required to maintain the areal extent of forage. Our data indicate that under warming conditions, the initial response is likely that available forage biomass will increase without a substantial loss of nitrogen content. Therefore, we suggest that goose populations are not constrained by bottom-up process in Arctic Alaska and we expect further increases in available forage under projected levels of warming for the near future.

Acknowledgments

This work is part of the U.S. Geological Survey Changing Arctic Ecosystem Initiative and was supported by funding from the Wildlife Program of the USGS Ecosystem Mission Area. Many field technicians logged many hours with cold wet hands to collect these samples or spent tedious hours sorting dried samples, including the following: K. Carlson, K. Crawford, A. Driskill, D. Gerik, C. Hansen, B. Linschoten, B. Mason, T. McIntyre, A. Meyers, M. Johnson, P. Stacey, S. Stark, L. Tennant, C. Twellmann, and T. Spivey. Field logistics were supported by the pilots B. Gill and M. McCrary of 70 North Air Service in Deadhorse, Alaska, as well as T. Shoemaker from Arctic Air Alaska in Fairbanks. M. Rogers and A. Brownlee at the University of Alaska Anchorage Stable Isotope Lab conducted the nitrogen content analyses. We thank J. Pearce for his support in conducting the work and for comments on the manuscript. R. Ruess as well as several anonymous reviewers provided valuable comments on the manuscript. We thank G. Brower and Commissioners of the North Slope Borough Planning and Community Services Department and the North Slope Borough, Department of Wildlife Management for their assistance and feedback on this research. We also thank the members of the NPR-A Subsistence Advisory Panel for their comments on how to best conduct our work in the Arctic to reduce disturbance to wildlife and subsistence activities. We thank the Arctic Field office of the Bureau of Land Management for logistical and technical support. Any use of trade, firm, or product names is for descriptive purposes only and does not imply endorsement by the U.S. Government.

Data Availability

Data are available from U.S. Geological Survey: https://doi.org/10.5066/P9WMFBVN.