Snow mediates climatic impacts on Arctic herbivore populations

Arctic ecosystems are particularly vulnerable to impacts of climate change; however, the complex relationships between climate and ecosystems make incorporating effects of climate change into population management difficult. This study used structural equation modelling (SEM) and a 24-year multifaceted monitoring data series collected at Zackenberg, North-East Greenland, to untangle the network of climatic and local abiotic and biotic drivers, determining their direct and indirect effects on two herbivores: musk ox (Ovibos moschatus) and collared lemming (Dicrostonyx groenlandicus). Snow conditions were determined to be the central driver within the system, mediating the effects of climate on herbivore abundance. Under current climate change projections, snow is expected to decrease in the region. Snow had an indirect negative effect on musk ox, as decreased snow depth led to an earlier start to the Arctic willow growing season, shown to increase fecundity and decrease mortality. Musk ox are therefore expected to be more successful under future conditions, within a certain threshold. Snow had both positive and negative effects on lemming, with lemming expected to ultimately be less successful under climate change, as reduction in snow increases their vulnerability to predation. Through their capacity to determine effects of climatic and local drivers within a hierarchy, and the relative strength and direction of these effects, SEMs were demonstrated to have the potential to be valuable in guiding population management.


Introduction
Global climate change is accelerated in the Arctic, where surface air temperatures have risen on average 1 °C per decade, twice the global rate, since the 1980s (Post et al. 2009;Screen and Simmonds 2010;Meredith et al. 2019). Warming by up to 10 °C is projected for the year 2100 (Overland et al. 2017). Although Arctic ecosystems are recognised as extremely vulnerable to climate change, they are amongst the least studied and monitored globally due to their remote and harsh nature (Meredith et al. 2019). The vulnerability of these ecosystems is exacerbated by the relatively simple trophic structure of Arctic terrestrial systems, in which each species plays a central role (Wirta et al. 2015). Herbivores, in particular, have important influences to Arctic ecosystem structure and functioning (Legagneux et al. 2014) and may be useful indicators of climate change impacts to these systems more broadly. However, improved understanding of the relationships between climatic drivers, local weather, and trophic interactions is necessary for effective management and protection of Arctic ecosystems (Wirta et al. 2015). This study applies a systems-level perspective to trace the influences of regional climate drivers-broad-scale atmospheric circulation patterns and regional sea ice extent-on herbivore populations via intermediary effects of local snow depth, air temperature, extreme temperature events, and vegetation quality and availability (Fig. 1). The relative importance of predation pressure is also considered. Each of these drivers is first briefly reviewed.

Climate: Arctic oscillation
Climate drives many ecological patterns and processes in the Arctic (Forchhammer et al. 2008a;Post et al. 2009). Arctic climate variability is captured by teleconnection indices such as the Arctic Oscillation (AO), North Atlantic Oscillation (NAO) and Atlantic Multidecadal Oscillation (AMO). The AO winter index has been shown to be significantly correlated with seasonal temperature variability and snow cover in the winter and spring throughout the Arctic (Thompson and Wallace 1998;Rigor et al. 2002;Bamzai 2003). Furthermore, seasonal indices of the AO and NAO may be better predictors of ecological responses than local weather variables (Hallett et al. 2004;Forchhammer et al. 2008b), as they capture the complexity of local weather, including cumulative forcings, which can have the greatest effect on species (Picton 1984;Hallett et al. 2004).
The AO characterises the air pressure exchange between the Icelandic low system and weaker pressure centres in the North Pacific (Thompson and Wallace 1998). When the AO is in a positive phase, there is an intense Icelandic low creating strong westerly winds, confining cold air to the Arctic. Conversely, the Arctic is warmer when the AO is negative: a weaker pressure gradient between the Arctic and the North Atlantic and North Pacific and weaker westerly winds in the Arctic allow cold air to leave the Arctic moving southward (Thompson and Wallace 1998;NOAA 2019a). There is considerable interannual variation in AO and the processes linking teleconnection indices, such as AO, and sea ice extent are not well understood (Overland and Wang 2010;Jaiser et al. 2012). Positive AO phases have previously been associated with reduced sea ice extent due to increases in wind-driven transport (Rigor et al. 2002); however, in the past few years this relationship has not held as reduced summer sea ice extent has followed extreme negative AO winter phases (Liu et al. 2004;Stroeve et al. 2011;NSIDC 2019a). It is possible that a continuing decline in Arctic sea ice extent will result in more frequent negative phases, as the reduced sea ice concentration in the summer results in stronger heat release to the atmosphere, which influences large scale atmospheric circulation favouring more negative AO phases, whilst also feeding back to further reduce sea ice concentration (Overland and Wang 2010;Jaiser et al. 2012;Nakamura et al. 2015). Through their interconnectedness, annual indices of the AO and Arctic sea ice extent encapsulate an interactive system of global climate drivers affecting Arctic ecology.

Climate: Arctic sea ice extent
Arctic sea ice extent is an indicator of climatic conditions (Maslanik et al. 1996) with considerable biological and environmental implications . Sea ice extent has persistently decreased by ~ 45,000 km 2 per year since 1981 (Maslanik et al. 2007;NSIDC 2019b). This decline is likely to accelerate over the coming decades, due to rising temperatures and the sea ice-albedo feedback mechanism, in which an increase in open water area leads to greater heat absorption, further accelerating sea ice decline (Kaplan and New 2006;Serreze and Francis 2006). Similarly, sea ice decline itself amplifies Arctic warming as ice-free waters result in large heat fluxes to the atmosphere from the warming ocean (Kaplan and New 2006;Screen and Simmonds 2010;Serreze and Barry 2011), with effects on local weather conditions including temperature variability and precipitation (Liu et al. 2012;Meredith et al. 2019). Reduced sea ice is likely to Fig. 1 Conceptual model identifying the key classes of drivers within the system and the expected relationships between them increase precipitation due to increased available atmospheric water content. However, whether this precipitation falls as snow (Callaghan et al. 2011;Liu et al. 2012;Bintanja and Selten 2014) or rain (Post et al. 2009;Overland et al. 2017) depends on air temperature and geographical location.
Through these influences on local weather, sea ice decline influences the seasonality and productivity of associated ecosystems. For example, localised warming of terrestrial environments adjacent to areas of sea ice loss has been linked to earlier soil thaw and spring vegetation green-up, which increases the potential for trophic mismatch, in which the beginning of the vegetation growing season no longer aligns with the beginning of the herbivore breeding season (Kerby and Post 2013;Post et al. 2013).

Local weather: snow depth
Snow is one of the most important local weather drivers of Arctic ecosystems. Snow affects the distribution, biodiversity and productivity of vegetation, which influences higher trophic levels (Callaghan et al. 2011). Snow insulates the soil surface during winter (Hinkler et al. 2008;Bokhorst et al. 2009), with ground temperatures under snow shown to be maintained at 0 °C whilst air temperatures reach − 20 °C (e.g. Pomeroy and Brun 2001). Snow depth, therefore, affects the overwinter survival of vegetation and productivity in the following spring (Walsh et al. 2003) and the timing of the spring snowmelt determines the beginning of the vegetation growing season (Ellebjerg et al. 2008;Callaghan et al. 2011). Snow also provides winter habitat and nesting sites for smaller mammals and invertebrates (Stenseth and Ims 1993;Bokhorst et al. 2012). Conversely, it obstructs larger herbivores in accessing forage ).

Local weather: air temperature and extreme temperature events
Air temperature has been shown to strongly influence primary productivity and the structure and intensity of trophic interactions in the Arctic (see below sections 'Vegetation: Productivity and Phenology' and 'Arctic Herbivores'; Walker et al. 2006;Elmendorf et al. 2012;Legagneux et al. 2014). Summer temperature variability has also been shown to influence ecosystem dynamics. Summer temperature extremes have been revealed as a key determinant of the duration and timing of the vegetation growing season, and to influence plant productivity (Ellebjerg et al. 2008;Callaghan et al. 2011;Piao et al. 2015). They have also been shown to impact the reproductive success of ungulates including reindeer and caribou (Post et al. 2009). In extreme midwinter warming events, positive air temperatures thaw the insulating snow layer, exposing vegetation upon the return to harsh winter temperatures (Callaghan et al. 2011). These events are therefore linked to increased plant mortality and reductions in productivity and fitness (Bokhorst et al. 2009). Rainfall during a warming event, known as rain on snow (ROS) events, creates ice layers within the snowpack, which can render the vegetation inaccessible to herbivores. ROS events have been implicated in population crashes of musk ox (Ovibos moschatus) (Forchhammer and Boertman 1993;Rennert et al. 2009), reindeer (Rangifer tarandus platyrhynchus) (Kohler and Aanes 2004) and caribou (Rangifer tarandus) throughout the Arctic, and have been shown to influence the amplitude of collared lemming (Dicrostonyx groenlandicus) population fluctuations (Domine et al. 2018). These mid-winter warming events and their impacts are becoming increasingly common under climate change (Callaghan et al. 2011;Bokhorst et al. 2012;Hoegh-Guldberg et al. 2018).

Vegetation: productivity and phenology
The quality and availability of vegetation, which determine the resources available to herbivores, are heavily influenced by broad-scale climate and local weather conditions (Forchhammer et al. 2002;Kerby and Post 2013;Schmidt et al. 2019). Increasing productivity and Arctic greening has been associated with increasing summer air temperatures (Blok et al. 2011), and earlier spring snowmelt and extended growing season length (Euskirchen et al. 2006). The productivity of vegetation in the current year also influences productivity of the following year, as pre-formed buds and carbohydrate stores enable rapid bud-burst after the next spring snowmelt (Ellebjerg et al. 2008).
Vegetation productivity has been linked to herbivore reproductive success, abundance and migration patterns (Forchhammer et al. 2005;Berg et al. 2008;Wipf and Rixen 2010). Growing season timing is also important due to tight coupling of the beginning of herbivores' reproductive season (Forchhammer et al. 2005;Berg et al. 2008). Because the Arctic growing season generally lasts 60-100 days, the potential for phenological mismatch is high (Ellebjerg et al. 2008;Post et al. 2009) and consequent reductions in offspring survival have been reported (Post and Forchhammer 2008;Plard et al. 2014).

Arctic herbivores
Musk ox (O. moschatus) and collared lemming (D. groenlandicus) are endemic herbivores that are dominant in many Arctic systems (Gunn and Forchhammer 2008;Cassola 2016). Although neither is considered threatened (IUCN Red List), they are both non-migratory and therefore may be at risk from the rapid climate change in the Arctic.

Musk ox
Musk ox are large ungulates with a long, thick woollen coat and horns. They typically live in mixed-sex herds of 20-30 individuals (Gunn and Forchhammer 2008). Musk ox are generalist herbivores, relying primarily on Arctic willow (Salix arctica) and graminoids (Thing et al. 1987;Berg et al. 2008) and requiring access to forage year-round (Forchhammer et al. 2002). Musk ox breed in August, with single calves born following 250 days gestation in May-June. The availability of good forage during these months is essential to breeding success (Forchhammer et al. 2002;Post and Forchhammer 2008). Musk ox hold cultural and economic importance for indigenous communities who hunt them for their meat, horns and fur (Nuttall et al. 2004).

Collared lemming
The collared lemming is a microtine rodent inhabiting Arctic tundra of North America and Greenland (Stenseth and Ims 1993). They generally spend most of their time belowground; when aboveground, they divide their time evenly between foraging and being vigilant Schmidt et al. 2008). During winter, lemming build nests within the snowpack and after an approximately 20 day gestation period they can produce a litter of 1-11 young (Stenseth and Ims 1993;Gilg et al. 2009). Lemming have a number of predators including the long-tailed skua (Stercorarius longicaudus), the Arctic fox (Alopex lagopus), the snowy owl (Nyctea scandiaca) and the stoat (Mustela erminea) Gilg et al. 2009;Barraquand et al. 2014).

Objectives and hypotheses
This study took advantage of a long-term ecological monitoring programme in Zackenberg, East Greenland (GEM 2020) to investigate the network of climatic, local weather and trophic drivers ( Fig. 1) and their influence on herbivore populations from a systems perspective. The hierarchy of interacting drivers highlights that ecological systems are too complex to be fully understood with pairwise analysis; therefore, this project used structural equation modelling (SEM; Shipley 2000; Pugesek et al. 2004) to analyse the time series. SEMs have been used to evaluate hypothesised networks of causal relationships in ecological systems (e.g. Byrnes et al. 2011;Duffy et al. 2015;Jing et al. 2015), but have rarely been applied to Arctic ecology (but see Juhasz et al. 2020). Specifically, the objective of this research was to determine the direct and indirect effects of climatic and local biotic and abiotic drivers on musk ox and lemming abundance in Zackenberg. The results are interpreted in the context of climate change projections and inferred impacts to the focal herbivore species, and the management implications of the causal mechanisms supported by the SEMs are discussed.
Climate change impacts to herbivore populations are uncertain and may be idiosyncratic. We hypothesised musk ox would be more abundant in years and climate conditions with (i) less local spring snow depth, increasing access to forage, and (ii) an earlier spring snowmelt, making vegetation more readily available and productive during their calving season (Forchhammer et al. 2002;Berg et al. 2008;Post and Forchhammer 2008). As musk ox do not have any major predators, they were expected to be 'bottom-up' controlled, with vegetation availability and weather conditions expected to have significant influences on their population dynamics (Hunter and Price 1992). In contrast, lemming were expected to be most abundant in years and climate conditions with high snow depth, as they rely on deep snow to create winter nests and for safety from predation Gilg et al. 2009). Lemming have been shown to exhibit strong cyclical population variation throughout time, which is potentially due to 'top-down' control by their predators (Hunter and Price 1992;Stenseth 1999;Schmidt et al. 2008). Both herbivores were expected to be affected by density-dependent effects, in which the population's abundance is regulated by the previous year's population size (Stenseth 1999;Forchhammer et al. 2002).

Study site: Zackenberg
The Zackenberg Research Station is located in the Zackenbergdalen valley on the coast of Wollaston Foreland peninsula in northeast Greenland (74° 30′ N, 20° 30′ W) (Fig. 2). It has a high Arctic climate with average temperatures of − 19.5 °C and 5 °C in the winter and summer respectively (measured over 1996-2019; GEM 2020). The Zackenbergdalen valley is ~ 2.5 km wide and includes the drainage basin for the river Zackenbergelven, which has a catchment area of 514km 2 that is ~ 20% glacier-covered (Meltofte and Rasch 2008). It is a mountainous area with diverse topography and peaks reaching 1400 m. Vegetation is tundra, with species such as Arctic willow, white Arctic bell heather (Cassiope tetragona), Arctic poppy (Papaver radicatum), mountain avens (Dryas octopetala) and graminoids such as dark cottongrass (Eriophorum triste) . Wildlife species include collared lemming, musk ox, Arctic fox, Arctic hare (Lepus arcticus), polar bear (Ursus maritimus) and diverse arthropods and migratory birds . Comprehensive environmental monitoring has been conducted annually at this site since 1995 (GEM 2020). The majority of the monitoring and research takes place within 20 km 2 surrounding the climate station (Fig. 2).

Data
This study used publicly available monitoring data from the National Oceanic and Atmospheric Administration (NOAA; climate variables) and the Zackenberg Basic monitoring programme (GEM 2020; local weather and biotic variables) for each year between 1996 and 2019, totalling 24 years of observations (Table 1). All variables were represented with a single value per year; Table 1 describes how measurements made at greater frequency were aggregated to annual or seasonal estimates.
Broad-scale climate was represented by winter AO and Summer (June-August) sea ice extent (Table 1). AO was chosen because, unlike the correlated NAO, its primary centre of action is over the Arctic. AO is best expressed during winter when the largest variability is captured (Thompson and Wallace 1998;Forchhammer et al. 2008b). Summer sea ice extent measurements were taken to correspond to the growing season (Kerby and Post 2013;Post et al. 2013).
In situ measurements of local weather and biotic variables were collated from the Greenland Ecosystem Monitoring (GEM) database subgroups GeoBasis, ClimateBasis and BioBasis (GEM 2020). Most in situ measurements were collected during the June-August field season; the climate station undertakes automated year-round monitoring. The variables were measured across slightly different spatial and temporal footprints (Fig. 2, and see details in Table 1).
Measures of local weather (Table 1) included average hourly snow depth over the summer and the winter seasons, average seasonal air temperature to capture typical  The variance of air temperature throughout the season was determined as a proxy for the frequency of extreme temperature events (Mearns et al. 1984). Air temperature was recorded hourly at 2 m at the Zackenberg Climate Station (Fig. 2) and then the variance across June, July and August was calculated. Temperature Variability Winter The process was as above; however, the variance was calculated across December (year t − 1), January and February (year t). Vegetation Productivity and Phenology Source: GEM BioBasis NDVI (Normalised Difference Vegetation Index) NDVI measures the photosynthetic capacity of vegetation; the normalised difference of reflectance of light in the near-infrared (NIR) to light in the RED wavelengths (NDVI = (NIR-RED)/(NIR + RED)). Values of NDVI range from − 1 to + 1, with dense, green vegetation typically between 0.5 and 0.9. The NDVI was calculated from a Landsat satellite image (30 m pixel resolution) acquired as close to July 31 as possible, considering imaging complications. All images were corrected for atmospheric and terrain effects. For musk ox models, the NDVI was averaged across the 47 km 2 musk ox census area, and for lemming models, the NDVI was measured within the 1.06 km 2 lemming census area (Fig. 2).
conditions, and air temperature variability in winter and summer as proxies of the occurrence of extreme temperature events (Mearns et al. 1984). Both measures of air temperature were included as increases in extreme temperature events under climate change (Hoegh-Guldberg et al. 2018) may occur due to either a shift (i.e. an increase in the average) or a broadening (i.e. an increase in the variability) of the whole temperature distribution (Meehl et al. 2000). Vegetation was represented by growing season productivity and flowering phenology of Arctic willow, the dominant shrub species (Table 1). Monitoring recorded population counts of musk ox, including demographic variables of age, sex, and records of carcasses, from which the demographic rates of fecundity and mortality were derived (Table 1). Lemming abundance was proxied by counts of winter nests (Gilg et al. Each local variable was measured annually at Zackenberg, Greenland between 1996 and 2019. Variables were either measured automatically throughout the year, or manually during the annual field season from June-August  (Thing et al. 1987;Stenseth and Ims 1993). At weekly intervals throughout the field season, all 7 phenology plots ( Fig. 2) were visited and at least 100 buds were counted in each plot. Once the ratio of open flowers to buds exceeded 50/100, the day of year was recorded. Musk Ox Abundance and Vital Rates (in year t and t − 1) Source: GEM BioBasis

Population
The average number of musk ox counted during a weekly field census throughout July-August, within the 47 km 2 musk ox census area (Fig. 2). The annual raw total of musk ox counted was divided by the total number of censuses performed during the season, to account for the variability in the total number of censuses performed each year and to account for any potential double counting. Mortality The number of adult (≥ 2 years) fresh musk ox carcasses counted during weekly field censuses, divided by the number of adult live musk ox counted, corrected for the total number of censuses performed during the season. Fecundity The number of calves observed, divided by average number of females of child-bearing age (3-4 + years) observed during weekly field censuses, corrected for the total number of censuses performed during the season. Lemming Abundance (in year t and t + 1) Source: GEM BioBasis

Lemming Nests
The total number of new lemming winter nests recorded over 60 transects. Undertaken during the field season within the 1.06 km 2 lemming census area (Fig. 2), which is the main area used by the lemming population. Upon examination of the nests, they were destroyed to avoid double counts. This is not harmful as lemming build new nests annually and the surveys were conducted after the nests had been discarded for the season (Stenseth and Ims 1993). Lemming Predator Abundance (in year t) Source: GEM BioBasis

Fox and Stoat Predation Rate
The number of lemming winter nests showing evidence of predation from stoat or fox, expressed as the percentage of the number of new winter nests. Nest predation is recognised by large amounts of lemming fur inside the nests used for insulation when nest is taken over as a stoat residence, and/ or the presence of fox and stoat faeces nearby.

Skua Breeding Pairs
The number of skua breeding pairs counted within the bird census area (Fig. 2) during the annual bird census undertaken in mid to late June.
1 3 , 2009Krebs et al. 2012). As the nest counts reflect winter abundance and are expected to represent the integrated effect of ecological processes in the previous year, analyses paired nest count measurements made in year t + 1 with measurements of all driving variables for year t. Winter snow depth for year t + 1 was also included in the lemming models as a driver to lemming abundance, allowing investigation of the direct effect of the snow depth during the period when lemming nests are constructed (as both variables reflect conditions from the same period (December of year t-February of year t + 1). Lemming predator abundance was represented by the number of breeding pairs of long-tailed skua present in the summer and by the predation rate on lemming winter nests by stoat and Arctic fox ( Table 1). The variables included in these analyses are a subset of variables available within the GEM dataset and were selected following expectations based on the current literature. Some alternative indicator variables, including lagged effects of climate and the local environment in the previous year, were considered but were not supported following preliminary analyses. To reduce the number of potential variables and avoid over-specifying the models, lagged effects of climate and environmental variables were assumed to be captured through density dependence, with the inclusion of the previous year herbivore abundance as a driver to the current year abundance. See Table 1, Skov et al. (2015) and Schmidt et al. (2017) for further details regarding the sampling procedures and estimation of the included variables.

Statistical analyses
Ecological systems are highly complex, comprising of a network of interdependent relationships across scales. Consequently, pairwise analyses can be problematic and inefficient. A useful multivariate alternative is structural equation models (SEMs; Shipley 2000), which translate hypothesised direct and indirect relationships between variables and across trophic levels into a composite hypothesis based on statistical dependencies (Shipley 2000;Grace et al. 2010). SEMs are a confirmatory approach; they begin with a hypothesised model based on existing knowledge which is refined iteratively to produce a final model that fits the observed data. The capacity to determine the indirect and direct effects of independent variables on the dependent variables makes SEMs more powerful than generalised linear models (GLMs) and other pairwise analyses when studying ecological systems (Shipley 2000;Pugesek et al. 2004). They can also include multiple dependent variables and allow for correlations between independent variables (Pugesek et al. 2004).

Model development
First, a priori models were developed, specifying the network of relationships based on existing knowledge from the literature. These are graphical models, known as specification diagrams, illustrating expected relationships between conceptual drivers. In this study, four different models are presented; a musk ox abundance (MOA) model (Fig. 3a), a musk ox vital rates (MOVR) model (Fig. 3c), a lemming abundance model excluding predation as a driver (LAEP) (Fig. 3b) and a lemming abundance model including predation effects (LAIP) (Fig. 3d). All models assessed all climate, local weather and vegetation drivers, and included the previous year's abundance of the focal herbivore to account for density-dependent effects (Fig. 3).
Musk ox abundance was the key response variable in both MOA and MOVR, but in MOVR the effects of environmental variables on musk ox abundance were mediated via effects on the vital rates, fecundity and mortality (Fig. 3c). This therefore investigates whether mortality or reproduction is more important for driving abundance and how the environmental conditions differentially influence these demographic processes. Both lemming models used lemming winter nest count as the key response variable (Fig. 3). However, LAEP only included vegetation variables as biotic factors, whereas LAIP also contained skua summer abundance and fox and stoat winter predation rate. Comparison between these models can suggest the relative importance of 'bottom-up' control from vegetation or 'top-down' control by predators.
Specification diagrams were translated into path diagrams, illustrating the hypothesised relationships between measured variables (Online Resource 1-4). Each relationship in the path diagram was evaluated with preliminary pairwise GLMs to test the assumptions of the maximum likelihood estimation method used in the SEMs (Schumaker and Lomax 2011). The GLM residuals were tested for normality, using the Shapiro-Wilks test (Shapiro and Wilk 1965), and homogeneity of variances, from visual inspection of the plot of residuals vs. fitted values. The musk ox total population, lemming nest count and skua breeding pair numbers required log transformation to ensure normal distribution (Shipley 2000;Ives 2015). All similar variables, such as measures of air temperature average and variability, were tested for correlation (all r found to be < 0.4) since collinearity in independent variables can increase the Type II error rate in SEMs (Grewal et al. 2004). The units of analysis for all SEMs were years. Any yearly measurements for a variable that were missing were left blank. The SEM deals with missing data through listwise deletion, in which the entire record containing a missing value is excluded from the analysis. Despite the exclusions, the musk ox models had a final n = 22 years and the lemming models had a final Fig. 3 The top panel shows the consistent basis of each specification diagram for each structural equation model. These describe the network of hypothesised relationships between conceptual drivers in the system (circles) and identifying the indicator variables (rectangles) used for each. The bottom four inserts outline the differences between each of the four models. The base a Musk Ox Abundance (MOA) and b Lemming Abundance Excluding Predation (LAEP) Models add only measurements of the focal species abundances to the general model structure. Expanded models incorporate, c fecundity and mortality processes, in the Musk Ox Vital Rates (MOVR) Model, and d the effects of predation, in the Lemming Abundance Including Predation (LAIP) Model. All variables were measured in the focal year t, unless otherwise indicated n = 21 years. Specifically, we used the piecewise structural equation modelling implementation (Shipley 2009;Lefcheck 2016). All SEMs were performed in 'R' (version 3.4.1) (R Core Team, 2018) with the package 'piecewiseSEM' using maximum likelihood parameter estimation (Lefcheck et al. 2020, package version 2.1.1).

Model refinement
Models initially included all possible hypothesised relationships between variables (Online Resource 1-4) and were refined by removing the relationship returning the most nonsignificant p value. This process was continued iteratively (following Grace et al. 2010) until only marginal (p < 0.1) or significant (p < 0.05) paths remained in the model.
At each step of model refinement, model fit was evaluated using the Fisher's C statistic, and the hypothesised model was determined to fit the observed dataset when the Fisher's C p value > 0.05 (Lefcheck 2016). Akaike's Information Criterion (AIC) was used to select between competing final models (Schumaker and Lomax 2011), providing their original specification was the same. A smaller AIC signifies a better balance of model fit and complexity, with a reduction of two AIC points or more considered to be a better model (Akaike 1974).
Results present final path diagrams revealing the set of relationships supported by the data following model refinement. Standardised path coefficients, ranging between − 1 and 1, were used to compare the relative importance of the direct and indirect effects remaining in the models (Shipley 2000). All indirect coefficients were estimated as the product of all coefficients along a path, summing alternate paths connecting two variables. Model performance was also assessed with the R 2 values of all endogenous variables, including internal response variables and the abundance of the focal herbivore, allowing comparison of the alternative formulations MOA and MOVR, and LAEP and LEIP.

Results
Initial models were large, containing many non-significant relationships. Following removal of non-significant relationships, the four presented models (Figs. 4 and 5) successfully fit the data (Fisher's C, p > 0.05; Table 2). Variables retained in all models were AO, sea ice extent, snow depth summer , temperature average summer and Arctic willow phenology. Snow depth winter was retained in the lemming models only, which kept the snow depth winter measurements from both years, t and t + 1. Both fecundity and mortality were retained in the MOVR model. The only measure of predation retained in the LAIP model was the fox and stoat predation rate. The effect of density dependence was retained in both musk ox models, whilst lemming abundance in the previous year was removed from both lemming models. NDVI, temperature variability summer , temperature variability winter and temperature average winter were not retained in any final models.
In all models, sea ice extent explained a moderate amount of the variation in average summer temperatures (R 2 for the endogenous variable was 19% for both musk ox models, and 14% for both lemming models; Table 3) and had an indirect positive effect on summer snow depth (cumulative effect coefficient for sea ice extent on summer snow depth was 0.23 and 0.12 for musk ox and lemming models, respectively; Figs. 4 and 5). Likewise, AO had consistent positive effects on snow depth, with a direct positive effect on summer snow depth found for musk ox models (Fig. 4) and a direct positive effect on winter The significance of each relationship is described by the weight of the arrow (see legend) and the standardised path coefficients are presented next to each arrow. The coefficients are based on a correlation matrix thus they are standardised and can be compared. Unless specified, all variables are in the year t All models provided a good fit to the data (i.e. p value > 0.05). C value represents the Fisher's C statistic, Df degrees of freedom. Model names are MOA musk ox abundance, MOVR musk ox vital rates, LAEP lemming abundance excluding predation, LAIP lemming abundance including predation C value refers to Fisher's C statistic, which is a χ 2 distribution with 2k degrees of freedom (Shipley 2013), in which k degrees of freedom is the number of independence claims derived from the basis set (i.e. the number of direct and 2nd order direct links which are not specified in the original model) ( snow depth, which itself had a direct positive effect on summer snow depth, in lemming models (Fig. 5). These snow depth variables explained over 70% of the variation in Arctic willow phenology in all models (Table 3), which was a function of summer snow depth alone in musk ox models (Fig. 4), and both summer and winter snow depth in lemming models (Fig. 5).

Musk ox
Overall, the timing of the Arctic willow growing season had the strongest total effect on musk ox abundance (total effects: − 0.51, − 0.28; Table 4). Within both the MOA and MOVR models, summer snow depth had a direct, positive effect on Arctic willow phenology (direct effects: 0.85; Fig. 4). Therefore, summer snow depth had indirect negative effects on musk ox abundance (Table 3), as decreased snow depth led to an earlier growing season for Arctic willow, increasing musk ox abundance (Fig. 4). All of these paths were highly significant. As summer snow depth itself was positively affected by AO (direct effects: 0.34; Fig. 4) and indirectly positively affected by sea ice extent (Fig. 4), these broad-scale climate drivers had indirect negative effects on musk ox abundance (total effects: − 0.15, − 0.08 and − 0.10 and − 0.05 respectively, Table 4). Summer temperature average had an indirect positive effect on musk ox abundance (total effects: 0.23, 0.12; Table 4) via its direct negative effect on summer snow depth and indirect negative effect on Arctic willow phenology (Fig. 4). A strong, positive effect of density dependence was observed in both musk ox models (direct effects: 0.63, 0.57; Table 3).
The inclusion of vital rates increased the explanatory power of musk ox models (MOA R 2 = 0.58; MOVR R 2 = 0.65; Table 2). The MOVR model showed that mortality had a slightly stronger effect on populations than fecundity (direct effects: − 0.37, 0.31; Fig. 4), but that fecundity was under stronger environmental control (R 2 : 0.25, 0.12; Table 2). Both fecundity and mortality were affected by the growing season timing of Arctic willow, which had a positive effect on mortality (direct effect: 0.34, Fig. 4) and stronger negative effect on fecundity (direct effect: − 0.50, Fig. 4).

Lemming
Summer snow depth had a total negative effect on lemming abundance (total effect: − 0.54, − 0.32; Table 4) which, in the LAEP model, was the sum of a direct negative effect (direct effect: − 0.79, Fig. 5) and an indirect positive effect mediated by Arctic willow phenology (total effect: 0.25, Fig. 5). The direct negative effect of summer snow depth was weaker and the indirect link through Arctic willow phenology was lost when predation was included in the model (snow depth summer direct effect: − 0.79, − 0.32, Fig. 5). Conversely, winter snow depth (year t + 1) had a positive direct effect on lemming abundance (direct effect: 0.34, 0.37; Fig. 5), which was slightly stronger when predation was included in the model. In the LAIP model, the positive effect of winter snow depth (year t + 1) was slightly greater than the negative effect of summer snow depth (direct effects: 0.37 and − 0.32 respectively, Fig. 5). Winter snow depth (year t) was found to have a direct positive effect on the following summer's snow depth (direct effect: 0.62, 0.61; Fig. 5). Summer and winter snow depth (year t) both had a positive effect on Arctic willow phenology (direct effects: 0.40, 0.53 respectively, Fig. 5), yet Arctic willow phenology only had a direct effect on lemming abundance in the LAEP model (direct effect: 0.62, Fig. 5). The inclusion of predator terms removed the direct effect of Arctic willow phenology on lemming abundance and increased the explanatory power of the model (LAEP R 2 = 0.3; LAIP R 2 = 0.57; Table 3). In the LAIP model, the fox and stoat predation rate had a direct negative effect on lemming abundance (direct effect: 0.44; Fig. 5). Summer temperature average had an indirect positive effect on lemming abundance (total effects: 0.17, 0.10; Table 3) through its direct negative effect on summer snow depth (Fig. 5) and sea ice extent had an weak negative total effect (total effects: − 0.06, − 0.04; Table 4) on lemming abundance, through its negative effect on summer temperature average (direct effect: − 0.44; Fig. 5). AO had a positive direct effect on winter snow depth (direct effect: 0.34; Fig. 5) but only a weak negative total effect on lemming abundance in the LAIP model (Table 4).

Central role of snow
SEMs of herbivore abundance found that snow depth plays the central role within the hierarchy of drivers in this Arctic system, by mediating the effects of broad-scale climatic conditions, and with both direct and indirect influences of its own on herbivore populations. Snow depth primarily influenced musk ox negatively and indirectly, through effects on vegetation phenology, whilst snow depth influenced lemming directly, through its role in reducing predation pressure, and possibly indirectly, through its influence on vegetation phenology. Despite the initial inclusion of a range of other local weather variables in the analyses, snow depth was the only variable to play a consistently significant, mediating role between the both broad-scale climate drivers and the ecological responses, although the effects of sea ice extent were related to processes involving average summer temperature, which was linked to the ecological responses by snow depth.

Broad-scale climatic drivers of snow depth
Snow conditions in the winter and subsequent summer months have been shown to be highly correlated with AO throughout the Arctic (Bamzai 2003;Brown et al. 2010) and at Zackenberg (Forchhammer et al. 2008b). The results of these analyses agree with those findings, with positive links between AO and summer snow depth in the musk ox models and with winter snow depth in the lemming models. The effect of sea ice extent on snow has been variably reported throughout the Arctic (Park et al. 2013), with some areas reporting snowier winters in years with low sea ice extent (Deser et al. 2010;Liu et al. 2012), whilst the majority of the Arctic is reporting reduced snow, an increase in rainfall and an earlier start to the growing season under these broad-scale forcings (Høye et al. 2007;Brown et al. 2010;Liston and Hiemstra 2011;Bintjana and Selten 2014). The positive relationship observed within this study between sea ice extent and snow depth, mediated by the negative effect of sea ice extent on summer temperature average (Screen and Simmonds 2010), therefore supports the latter, that a continued reduction in sea ice extent may cause a decline in snow depth at Zackenberg. AO and sea ice extent are closely linked (Rigor et al. 2002), with negative AO phases being more commonly linked to periods of low sea ice extent in recent years (Liu et al. 2004;Stroeve et al. 2011), potentially creating a feedback with one another. Thus it is unsurprising that they were both found to exert a positive effect on snow depth, meaning that both a reduction in sea ice extent and the increasing frequency of negative AO phases would be likely to affect a reduction in snow at Zackenberg. It is worth noting that both AO and sea ice extent have previously been correlated with seasonal temperatures (Thompson and Wallace 1998;Rigor et al. 2002;Screen and Simmonds 2010); however, in this study, there was no correlation with average winter temperatures or temperature variability. Sea ice extent was only correlated with average summer temperature, which then directly affected summer snow depth. The two climatic drivers were therefore shown to hold an important and connected role in this hierarchy of drivers through their effects on snow.

Snow as a driver of vegetation conditions
Snow emerged as the central driver within this system as, not only was snow depth correlated with the climatic drivers, snow also had direct effects on vegetation phenology, and both direct and indirect effects on herbivore populations. A strong effect in all models was the positive relationship observed between summer snow depth and the phenology of Arctic willow, which is to be expected as snow depth determines the timing of the snowmelt, and the timing of the growing season is dependent on the snowmelt (Ellebjerg et al. 2008;Callaghan et al. 2011). Snow effects on musk ox were mediated by these effects on Arctic willow phenology, with reduced summer snow depth associated with an earlier Arctic willow growing season and increased musk ox abundance. Within the lemming models, winter snow depth was also found to have a positive direct effect on Arctic willow phenology. Winter snow also had an indirect effect on Arctic willow phenology, by its positive effect in determining the following summer's snow depth. Reduced summer and winter snow depth were both found to be associated with an earlier growing season in the lemming models. However, unlike for the musk ox, there was only equivocal support for an indirect effect of snow on the lemming population, mediated by Arctic willow phenology. This indirect effect was positive in the lemming models that excluded predation as a driver-i.e. earlier growing seasons were associated with reduced lemming abundance, but was not supported in the lemming models that did incorporate predation terms.

Snow as a habitat
Deep snow has previously been related to lemming success, as this species relies on snow for nesting and to escape predation, especially during the winter months (Gilg et al. 2009;Duchesne et al. 2011;Reid et al. 2012;Bilodeau et al. 2013a;Fauteux et al. 2015) when they are preyed upon by stoat and Arctic fox (Gilg et al. 2006;Schmidt et al. 2008Schmidt et al. , 2012. This is supported by the present study, as a positive direct effect of winter snow depth (corresponding to the nesting period for which the nest counts were performed; year t + 1) on lemming abundance was observed. Moreover, a negative relationship was observed between fox and stoat predation rate and lemming abundance. Thus, reduced winter snow depth was associated with a decreased lemming abundance, likely due to diminished protection from winter predation by the stoat and Arctic fox (although it is important to note that a direct effect of snow depth on Arctic fox and stoat predation rate, which would further strengthen this interpretation, was not found).
The preceding summer's snow depth, however, was found to have a direct negative effect on lemming abundance. This effect cannot be explained by forage availability, as NDVI was not retained in the model and Arctic willow phenology had a separate path. It also cannot be explained by a link between snow and summer temperature, with high summer snow depth indicating harsh summer conditions, as no effects of temperature on lemming were detected. This direct effect of summer snow depth may therefore reflect that minimal or absent summer snow could lessen summer predation, which is dominated by birds; long-tailed skua and snowy owl (Schmidt et al. 2012), to which the contrasting background of snow would be a hunting advantage as the collared lemming summer morph exhibit a brown-grey pelage (Nagy et al. 1995). Whilst long-tailed skua abundance has previously been found to be closely coupled with lemming abundance (Schmidt et al. 2012;Barraquand et al. 2014), skua abundance was not retained in the LAIP model. This may reflect that their diet is not specialised to lemming, as they are also known to prey heavily on wader chicks at Zackenberg ). The snowy owl was not included in the analysis due to a lack of data. As summer predation alone has previously been shown to drive significant collared lemming population decline (Gilg 2002), further examination into the role of summer predation on lemming abundance at Zackenberg, including the influence of snow on predation rates, is necessary.
There is ongoing debate whether lemming are under top-down (Stenseth and Ims 1993;Gilg et al. 2006;Krebs 2011;Fauteux et al. 2016) or bottom-up (Turchin et al. 2000;Oksanen et al. 2008) control, with greater population regulation by predators or forage, respectively. These analyses suggest that top-down control is the dominant mechanism of control operating within this lemming population, as there was a direct effect of predation and the effect of vegetation on lemming abundance was removed when predation was added to the model. In addition, the positive effect of winter snow depth, which is known to be relied upon to escape predation, increased when predation was added to the model, and the relationship between summer snow depth on lemming may also be related to protective camouflaging from predators.

Expectations of herbivore success under a changing climate
Continued, pronounced climate change is projected for the Arctic (Meredith et al. 2019). The results of these analyses can be used to suggest potential climate change impacts on musk ox and lemming, which may inform management strategies for climate change adaptation of these species. Current projected trends towards reduced sea ice extent and more frequently negative AO phases are associated with an increase in average summer temperature and decrease in snow at Zackenberg (although AO and sea ice extent do not account for all of the variability of average temperature and snow depth). This is consistent with other areas of the Arctic which are experiencing increasing summer temperatures, reduced snowfall, earlier spring snowmelt, and an earlier start to the vegetation growing season (Høye et al. 2007;Liston and Hiemstra 2011).
As snow depth was found to have a strong positive effect on the timing of the Arctic willow growing season in all models, climate change may continue to bring the growing season earlier in Zackenberg. This supports recent findings that Zackenberg has generally been experiencing an advancement of spring since monitoring began (Høye et al. 2007;Westergaard-Nielsen et al. 2017), despite some exceptions following winters with extreme snowfall (Schmidt et al. 2019). The timing of the growing season is particularly important for musk ox, as good access to forage in May/June at the beginning of the calving season is integral for reproductive success (Forchhammer et al. 2002;Post and Forchhammer 2008). The consequences of a mismatch between calving and forage availability can be dire Post and Forchhammer 2008;Schmidt et al. 2019). Our models suggest that projections for earlier springs are expected to be beneficial for musk ox abundance. However, this is likely to only be true so long as spring advancement does not overshoot the calving season entirely. Over the time series, the median day of year when 50% of Arctic willow flowers were open was 23rd June ± 7.6 days. There is therefore a threshold of several weeks where the spring could advance before a mismatch with the May/June calving season is anticipated.
An earlier growing season may also lead to a mismatch with the lemming population, as collared lemming are known to have winter and summer morphs in which they use photoperiod as a cue for an annual increase in body mass (Nagy et al. 1993(Nagy et al. , 1995. Lemming increase their body mass in August, when day length is long but decreasing rapidly, and maintain this larger morph until mid-winter (Nagy et al. 1995). Their maximum growth rate has therefore developed to coincide with peak vegetation productivity and biomass (Nagy et al. 1995;Ellebjerg et al. 2008). As such, with an earlier growing season, the timing of the photoperiod-driven change in body size may mismatch with the timing of the availability of the most productive vegetation, meaning a decline in summer snow depth may be anticipated to have an adverse effect on lemming abundance via its effect on Arctic willow phenology. However, the direct connection between Arctic willow phenology and lemming abundance was only included in the LAEP model, thus our results do not provide strong support for climate change impacts to lemming due to phenological disruption. As well, in addition to effects on vegetation phenology, reduced summer snow depth was also shown to relate to increased lemming abundance via a direct pathway. Our results suggest that the overall impact of a reduction in summer snow depth across these two paths is positive, but the contrasting effects of summer snow on lemming warrant further investigation.
The effect of a reduction in winter snow depth on lemming abundance is more straightforward, as it is expected to have an overall reducing effect on the lemming population due to their heavy reliance on snow for protection from predation (Gilg et al. 2009;Duchesne et al. 2011;Bilodeau et al. 2013a). It can be expected that continually declining snow conditions would have increasingly negative impacts on lemming's ability to escape winter predation and build nests to reproduce (although it is important to note that how the predators will independently be affected by the changing climate beyond the effect of winter snow depth, has not been considered in the models). Although our study found lemming to have opposite responses to summer and winter snow depth, on balance the standardised path coefficient was slightly stronger for winter than for summer snow depth in the better-performing model of lemming abundance (LAIP), suggesting that the overall effect of climate change on lemming may be negative.
It is clear that care should be exercised when inferring climate change impacts from these results. The studied relationships between climate, local weather and ecosystem responses are almost certainly non-linear over longer gradients and may exhibit tipping points (May 1977;Groffman et al. 2006). For example, it is expected that the relationship between musk ox and Arctic willow phenology is not infinitely linear, as if the snowmelt is consistently shifted earlier, the most nutritious new growth period will eventually mismatch with the musk ox calving season (Post et al. 2009), as ungulate calving seasons generally display limited interannual variability and would not be expected to keep up with the advancing spring (Kerby and Post 2013;Paoli et al. 2018). Therefore, this mismatch could be expected at Zackenberg, similar to the mismatch observed between caribou and vegetation productivity in response to rising summer temperatures in West Greenland resulting in increasing mortality (Post and Forchhammer 2008;Kerby and Post 2013). It is therefore important to avoid extrapolating interpretations from statistical analyses, such as these, beyond the observed range of values in the time series. This is unlikely to be prohibitive, however; due to the high variability in the Arctic, the range of values in long-term monitoring records is likely to be broad, as was the case with the Zackenberg time series.

Management implications
Musk ox and collared lemming both play important roles in Arctic ecosystems and at Zackenberg (Legagneux et al. 2014). Lemming are an important prey item of predator communities (Schmidt et al. 2012). Musk ox promote plant species diversity through effects of grazing (Post and Pedersen 2008) and hold substantial cultural and economic value for indigenous populations of Greenland (Nuttall et al. 2004). Therefore, where climate change is anticipated to have negative impacts, as these analyses suggest may be the case via some mechanisms for lemming (and may occur beyond certain tipping points of spring advancement for musk ox), it may be desirable to proactively develop management strategies to mitigate the risks.
Climate impacts are challenging to manage. Most management strategies are effective against direct human impacts (Heller and Zavaleta 2009). In contrast, climate changes cannot be directly mitigated, so management should ensure that both the species and the systems that sustain them have adaptive capacity (Gallopín 2006), through active management of cumulative impacts or the facilitation of climate change responses such as establishing connected networks of protected areas to accommodate distribution shifts (Heller and Zavaleta 2009). However, this is difficult in high latitudes, where direct impacts are low, climate change poses the greatest threat (Post et al. 2009), and systems are already at the ends of climate gradients. Protection of heterogeneous landscapes may provide climate change refugia for vulnerable species (Morelli et al. 2016) and accommodate divergent ecological requirements and species-specific responses. In particular, this study suggests that sites with heterogeneous snow cover may support populations of both musk ox and lemming. This particularly applies to the lemming population which had contrasting responses to snow, but relies on snow for protection in the winter. Under reduced snow conditions, developing sites where winter snow cover can most effectively accumulate, such as with snow fencing, may assist in supporting the lemming population Bilodeau et al. 2013b).
Compensatory conservation activities, in which adaptive management of alternative processes is performed to balance climate change impacts (Bauduin et al. 2018), may also be viable and understanding from our SEMs may be used to suggest appropriate approaches. As these SEMs initially included both summer and winter conditions, yet winter conditions were largely removed due to non-significance within the musk ox models, these results suggest that the musk ox population is largely driven by the summer conditions. Whilst the current climate change projections indicate potentially favourable conditions for musk ox (providing the changes exist within a certain threshold), it is valuable to be aware that the potential impacts to musk ox are likely to arise during summer and be mediated by forage availability. For example, years with high snow cover and late spring green-up would likely result in population declines due to phenological asynchrony. This population impact may be reduced or avoided with supplemental feeding (Boutin 1990;Kowalczyk et al. 2011). These analyses indicate that these management strategies are best aimed at reducing mortality which had a stronger effect on musk ox abundance.

Monitoring implications
Long-term monitoring is essential for adaptive, evidencebased management and for obtaining adequate knowledge of complex climate change impacts (Meredith et al. 2019). The twenty-four year time series used here is a long and valuable dataset from an ecosystem management perspective, especially within the Arctic. However, it is a relatively small sample from a statistical perspective (Grace 2006) and the addition of just a few years of data may change the results and interpretations. Continuation of this monitoring is therefore crucial; the longer a time series becomes, the more effectively it will capture environmental gradients and support more complicated models, with more relationships, including interactions, and greater power to detect an effect (Shipley 2000;Pugesek et al. 2004). Within these analyses, the number of identified pathways was likely limited by the relatively small sample size and the relatively modest levels of significance observed are expected to improve with larger sample sizes.
However, efficient monitoring should focus on variables that are the most effective indicators, closely related to the driving ecological processes. The identified role of snow depth highlights the importance of ongoing, accurate snow measurements at Zackenberg. Whilst our analyses used snow depth measured at a single point, a more spatially distributed measure of snow would be preferred. However, though snow cover measurements were made across the study area (but restricted to the summer field season), due to equipment malfunctions resulting in several years of missing data, these data could not be included in the models. This variable may have been particularly useful and, as such, we also recommend extending snow cover measurement to include springtime (March-May) to align with the herbivore reproductive season, to the extent possible (acknowledging that Zackenberg is not inhabited outside of the field season), and supplementing these measurements with snow modelling studies. This would facilitate evaluating if snow cover is a more effective indicator variable linking weather variation to herbivore populations. In addition, the density of the snowpack has previously been determined to have an effect on the vulnerability of lemming populations to predation (Bilodeau et al. 2013c) and whilst snow density is a measured parameter at Zackenberg within the GEM database, there are several missing years of data and the data is taken only during summer. If possible, determining snow pack density during winter, when the snow pack is most crucial to protection against predation, would provide a more complete picture of the environment and allow its impacts to be incorporated into models of lemming abundance. Finally, we were surprised that our proxy of extreme weather events (temperature variability) was not retained in any of the models. This may be due to the relatively indirect indicator variable that was used. More useful indictors may be direct measurements of ROS events at Zackenberg, if these could be detected from automated data series, or the development of alterative indicators of extreme events from the temperature time series.
Population measures related to ecological processes, such as fecundity and mortality, are particularly valuable in understanding population dynamics and developing management strategies. Within this study, the inclusion of fecundity and mortality substantially improved the musk ox abundance model (Table 2). Thus, access to a population count and measures of vital rates for the lemming population could be a very useful addition to the Zackenberg monitoring programme.

Conclusions
Climate impacts should be incorporated into population monitoring and management, especially in the Arctic where climate change is accelerated. This study demonstrated how SEMs can be a useful tool for analysing time series data and studying the direct and indirect effects of climate in a network of abiotic and biotic ecosystem drivers. These analyses suggest that snow conditions play the central role within the system, mediating climate impacts on musk ox and lemming populations. With projected climate trends, an earlier growing season and reduced snow depth could be expected at Zackenberg, which our results suggest will have positive effects on musk ox abundance (where the changes occur within a threshold of conditions) and both negative and positive effects on lemming abundance. These opposing effects highlight that there is no 'one size fits all' approach to incorporate climate effects into population management. Although climate impacts cannot be directly mitigated, SEMs can evaluate the support for hypothesised networks of relationships between variables based on data, providing a broader, system-level understanding which can be used to suggest appropriate management actions to compensate for these climate change impacts.
Acknowledgements Data were provided from the Greenland Ecosystem Monitoring (GEM) Programme and the National Oceanic and Atmospheric Administration (NOAA). The authors would like to thank GEM for their dedicated efforts collecting data at Zackenberg for 25 years, and for the open access to the BioBasis, ClimateBasis and GeoBasis databases. The authors would especially like to thank Dr. Kristy Langley (ClimateBasis) and Dr. Niels Schmidt (BioBasis) for their insightful feedback and revisions on the manuscript prior to submission. The authors are grateful to Dr. Tuomas Kankaanpää and Prof. John-William Shipley for their comprehensive and valuable reviews which greatly improved the manuscript.
Author contributions RJD and MCF developed the project, RJD and MEA designed the analyses, RJD performed the data analyses and interpreted the results-to which MCF and MEA contributed valuable discussions, RJD wrote the manuscript, and MEA provided critical feedback on the drafts.
Funding Open access funding provided by The University Centre in Svalbard. The authors received no funding for this study.
Data availability All GEM databases are available free of charge from http:// data.g-e-m. dk and all NOAA databases are available free of charge from https:// www. ncdc. noaa. gov/.

Conflict of interest
The authors declare that they have no conflict of interest.

Ethical approval
No permits or ethics approvals were required for this study.
Open Access This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if changes were made. The images or other third party material in this article are included in the article's Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article's Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this licence, visit http:// creat iveco mmons. org/ licen ses/ by/4. 0/.