Introduction

Populations of muskoxen (Ovibos muschatus) in high arctic Greenland are locally migratory. The animals move inland in autumn and spend winter on high ground where there is generally less snow than on the coastal area where they spend summer (Forchhammer and Boertmann, 1993). They calve inland in March/April and only return to the coast, where there is substantially greater biomass of for-age, after the melt in May/June (Thing et al., 1987). Females, therefore, have to meet the energy demand of early lactation by drawing on energy reserves acquired and stored since the previous summer (White et al., 1989; Kerby and Post, 2013). Additionally, bulls increase energy uptake before the rut in August/September, to meet the additional energy requirements for mating (Forchhammer, 1995). It follows that weather variables such as date of snow melt and summer temperature may influence the nutritional status of the animals in summer and also supply cues for migration to lower altitudes. For instance, Forchhammer et al. (2005) demonstrated that the movement of muskoxen into and out from the Zackenberg Valley, Northeast Greenland, displayed a negative density dependence and suggested that the strength of density dependence might be related to variation in forage biomass. However, as climate changes are gradually advancing spring in the High Arctic (Høye et al., 2007; Mortensen et al., 2014), it is likely that the high arctic muskox populations are experiencing relaxed competition for forage and, hence, earlier immigration cues, as vegetation can be accessed progressively earlier. Finally, females may be regarded as a resource formales during the rut, where, in particular, the operational sex ratio (number of females per male) significantly influences the spatial movement of males (Forchhammer and Boomsma, 1998). Thus, to expand on our current understanding on the muskox population dynamics during spring migrations, we examine the combined effects of variation in climate (temperature, snow cover), forage availability (length of plant growth season), and male access to females (females per male) on the average seasonal values of immigration rate, abundance, and density dependence, respectively, in a high arctic muskox population. Specifically, we predict that earlier snowmelt could induce higher immigration rates, as earlier snowmelt will allow earlier access to food resources, prompting an incentive to be first on-site. Additionally, we expect that the density dependence will be associated with the sex ratio in the valley, as the ratio of females to males decreases, competition among males in the rut will increase (Forchhammer and Boomsma, 1998).

Methods

Data and Analysis

Data on the total number, age/sex class, and location of muskoxen were obtained from weekly censuses of the entire area from June to August each year (Schmidt et al., 2012) The census were conducted by a transect walk, encircling the valley, ensuring vantage points to see all of the valley All encountered muskoxen were classified and counted, using binoculars. A weekly census could be cancelled or postponed, due to bad weather conditions or low visibility

For short-term changes in population size, such as those considered here on an intraannual, weekly timescale, we may assume that birth and death rates equal zero. Furthermore, given that the animals were not individually recognizable, the model considered only net migration (i.e., I — E).To embrace both the apparent nonstationary trend in muskoxen observed (Fig. 1) and potential density dependence in the net migration of muskoxen into the Zackenberg Valley, we applied first-order autoregressive integrated moving average model (arima(1,0,0)) with an external regressor (xreg) (R Core Team, 2014). With xreg defined as week number (t)—that is, representing a vector with the length(X_t )—we get

where X_t is the ln-transformed population size week t. The slope a is a within-year average estimate of direction and size of net migration. Whereas positive values of a indicate overall immigration, negative values indicate an emigration over weeks. The intercept k is an index for the annual level of muskox abundance (henceforth termed “abundance”), and the first-order autoregressive coefficient β₁ represents the degree of density dependence, with undercompensation for ß₁ < 1, perfect compensation for ß₁ = 1, and overcompensation for ß₁ > 1 (Clutton-Brock, 1997).

FIGURE 1.

Weekly counts of the number of muskoxen in the Zackenberg Valley across 17 years. Each thin line represents a separate annual value, while the black line represents the weekly mean across years.

The three coefficients in Equation 1 yielded a new data set of annual average coefficients, which, in turn, were analyzed across years by general linear regression to estimate and contrast across year variations to changes in climate, length of growth season, and sex-ratio settings. Specifically, summer temperature, spring snow cover, and the length of growth season and operational sex-ratio were selected as predictor variables with the three variables ß ¹ , a, and k as response variables, respectively (Table 1). Snow is a major controller of biotic processes in the Arctic (Meltofte et al., 2008). Temperature governs the rate of plant growth through energy input to the system (Sturm et al., 2001; Natali et al., 2012; Buchwal et al., 2013). The growth season, which is defined as the number of days with <10-cm-thick snow cover, was included as it defines the length of period with easy access to high quality forage (Meltofte et al., 2008). Abiotic data were sampled from two meteorological masts placed in the center of the valley (Schmidt et al., 2012). Finally, the adult sex ratio was integrated as an additional predictor because males compete for females during rut and might influence the selected response variables ß ₁ , a, and k (Forchhammer et al., 2005).

The effects of the selected climate and growth season variables with sex ratio on the response variables β ₁ , a, and k were analyzed using generalized linear models (Zar, 1999). Full models were reduced using stepwise backward removal of least significant variables until the most parsimonious models were achieved. Parsimony of models was estimated using the Akaike Information Criterion (AIC). All analyses were conducted in R (version 2.15.1).

TABLE 1

Explanatory variables for muskox migration, along with a description of each measurement.

FIGURE 2.

(A) Average annual observations of muskox in the Zackenberg Valley, divided into mature males and females, unspecified mature individuals, immature individuals, and calves.(B) The percentage composition between sex and calves.

Results

The data included within-year estimates of selected response variables (ß ₁ , a, and k) for each of 17 consecutive years (1996–2012). These were based on between 7 and 17 (mean ± s.d. = 11.6 ± 3.2) consecutive weekly censuses made generally between 1 July and 30 August each year (Schmidt et al., 2012). Outside this period, one census was made on 16 March 2008 and one was made on 24 October 2011. The number of muskox observed per census ranged from 1 (22 July 1996) to 372 (4 October 2011).The series mean (s.d.) was 72 (62) muskoxen (Fig. 1).

FIGURE 3.

Interannual variation in (A) intercepts estimated from annual muskox censuses, and (B) the coefficients of the 1th order autoregression, together with (C) the immigration rate. Straight line indicates trend line. The autoregressive coefficients (ß₁) and net migration rate (a) showed a temporal trend (p < 0.05), while the intercept coefficients (k) are not significant (p > 0.05).

The average annual observations of mature males, females, immature, and calves, along with the population structure, can be seen in Fig. 2. Annual temporal linear regression of each group showed that the average annual observations of adult males and females had significantly increased over the study period (p < 0.05, df = 16, α_males = 2.3, α_female = 2.6), while there is an unaltered average annual observation of immatures and calves (p > 0.05, df = 16) (Fig. 2).

The annual estimate of average level of abundance (k) increased across years although the trend was not significant (α = 0.003, p = 0.10) (Fig. 3, part A). The autoregression coeffient ß ₁ increased significantly over years, indicating an overall decrease in the strength of density dependence from 1996 to 2012 (GLM: α = 0.03, S.E. = 0.014, p < 0.05, R² = 0.25) (Fig. 3, part B). The annual net rate of immigration (a) decreased from 1996 to 2012 (GLM: α = 0.012, S.E. = 0.004, p < 0.05, R² = 0.29) (Fig. 3, part C).

The accumulated positive degrees in summer increased by 10.6 °C per year from 1996 to 2012 (GLM: α = 10.6, S.E. = 3.1, p < 0.01, R² = 0.43) (Fig. 4, part A). Spring snow cover decreased by 3% annually (GLM: α = -2.9, S.E. = 0.8, p < 0.01, R² = 0.50) (Fig. 4, part D). The sex ratio decreased from 2.14 in 1996 to 0.8 in 2010 although the trend was not significant (GLM: α = -0.04, S.E. = 0.02, p = 0.06, R² = 0.22) (Fig. 4, part C) There was no significant trend in the duration of the plant growth season (GLM: α = 1.6, S.E. = 2.3, p = 0.50, R² = 0.04) (Fig. 4, part B).

FIGURE 4.

Annual variation in the selected explanatory variables: (A) Temperature (cumulative positive degree days in June, July, and August), (B) snow-free period (time between snow cover less than 10 cm and higher than 10 cm), (C) sex ratio (female to male ratio), and (D) June snow cover (snow cover in the valley on 12 June). Straight lines are regression lines. Temperature and June snow cover all display significant changes (p < 0.05), while the snow-free period displays no significant trend (p = 0.15). Sex ratio displays an almost significant trend (p = 0.057).

We found significant parsimonious models for ß ₁ and k but not for a (Table 2). Summer temperature had a significant positive effect on ß ₁ (GLM: p <0.05, AIC = 12.0, R² = 0.27) (α = 2.2 × 10^-3, p = 0.03). Summer temperature (α = 2.4 × 10^-3, p = 0.037) and the length of growth season (α = 5.5 × 10^-3, p = 0.047) had positive effects on k (AIC = -55.6, R² =0.58). There was a strong negative relationship between the net migration rate (a) and the density dependence (ß ₁ ) (Pearson's r, ρ = -0.63).

Discussion

The results from the current study suggest that summer temperatures and length of growth season are primary drivers of intra-annual muskox population dynamics in the Zackenberg Valley The applied models demonstrated that summer temperatures both influenced the abundance of animals in the valley directly (k; Table 2), along with the strength of regulation of numbers in terms of its influence on density dependence (ß ₁ ;Table 2).

The increase in ß ₁ from 1996 to 2012 (Fig. 3, part A) indicates a decrease in the strength of density dependence (Clutton-Brock, 1997). Two conclusions follow. First, the level of intraspecific competition within the population was progressively released across the study period 1996–2012. Second, the competition release was most likely a consequence of improved foraging conditions for the muskox. Significant warming across the years (Fig. 4, part A) is likely to have resulted in a progressive increase in the annual rate of primary production (Sturm et al., 2001; Forbes et al., 2010), and consistent with this, we observed an independent positive relationship between the ß ₁ and summer temperature (Table 2). Similarly, Forchhammer et al. (2005) demonstrated that the strength of density dependence acting on female muskox was associated with variation in the annual normalized deviation vegetation index (NDVI), a proxy of biomass, which was, in turn, associated with variation in snow cover.

TABLE 2

Results from the reduction of the generalized linear models, with extracted coefficients as dependent variables, and climate and sex ratio as explanatory variables. Numbers are regression coefficients, and bold numbers indicate significant (p < 0.05) and italic almost significant regressions (p < 0.10). Bold Akaike Information Criterion (AIC) numbers give the most parsimonious model.

The annual net rate of immigration decreased (a) across the time series along with the strength of density dependence, which resulted in a negative correlation between the two variables. We interpret this in terms of a cascading effect where immigrants draw other muskox into the area, resulting in increased competition and an increase in the strength of density dependence. This result suggests that competition and immigration rate are interconnected. However, the net rate of immigration rate was not influenced by any weather variables, suggesting no coupling between annual net migration and changes in forage conditions across years.

The annual level number of muskox present in the valley (k) did not significantly change over the study period. However, our analysis suggests that the number of muskox was associated with the summer temperature as well as the length of the growing season (Table 2). As both summer temperature and growing season have previously been associated with the timing of access to forage along with forage abundance (Chapin, 1983; Billings, 1987), it implies that the annual level of abundance of muskoxen in the valley is primarily driven by current year's resource base. Specifically, an increase in summer temperature and a lengthening of the growing season would result in increasing annual number of muskoxen in the Zackenberg Valley, likely due to increased food abundance and a prolonged period of access to forage.

In summary, our study suggests that the local dynamic of muskoxen in the Zackenberg Valley is interrelated with summer temperatures and length of growth season. As the climate is predicted to change in the High Arctic (Weller, 2005; IPCC, 2014), it is likely that the increased rates of primary production and earlier access to forage will increase the abundance of muskox in the valley until checked by density dependence. The predicted climate changes could thus be favorable to the muskox population in the Zackenberg Valley, unlike the predicted mixed effects on reindeer in the High Arctic (Tyler et al., 2008), associated with spring icing of the snow-pack. However, predicted effects such as droughts and disappearance of cooling snow patches have yet to be investigated, and possible nonlinear effects of temperatures could exist as with the High Arctic reindeer (Tyler et al., 2008).