Spatial heterogeneity in climate change effects decouples the long-term dynamics of wild reindeer populations in the high Arctic.

The 'Moran effect' predicts that dynamics of populations of a species are synchronized over similar distances as their environmental drivers. Strong population synchrony reduces species viability, but spatial heterogeneity in density dependence, the environment, or its ecological responses may decouple dynamics in space, preventing extinctions. How such heterogeneity buffers impacts of global change on large-scale population dynamics is not well studied. Here, we show that spatially autocorrelated fluctuations in annual winter weather synchronize wild reindeer dynamics across high-Arctic Svalbard, while, paradoxically, spatial variation in winter climate trends contribute to diverging local population trajectories. Warmer summers have improved the carrying capacity and apparently led to increased total reindeer abundance. However, fluctuations in population size seem mainly driven by negative effects of stochastic winter rain-on-snow (ROS) events causing icing, with strongest effects at high densities. Count data for 10 reindeer populations 8-324 km apart suggested that density-dependent ROS effects contributed to synchrony in population dynamics, mainly through spatially autocorrelated mortality. By comparing one coastal and one 'continental' reindeer population over four decades, we show that locally contrasting abundance trends can arise from spatial differences in climate change and responses to weather. The coastal population experienced a larger increase in ROS, and a stronger density-dependent ROS effect on population growth rates, than the continental population. In contrast, the latter experienced stronger summer warming and showed the strongest positive response to summer temperatures. Accordingly, contrasting net effects of a recent climate regime shift-with increased ROS and harsher winters, yet higher summer temperatures and improved carrying capacity-led to negative and positive abundance trends in the coastal and continental population respectively. Thus, synchronized population fluctuations by climatic drivers can be buffered by spatial heterogeneity in the same drivers, as well as in the ecological responses, averaging out climate change effects at larger spatial scales.

INTRODUCTION

Fluctuations in the abundance of wildlife populations are often strongly influenced by stochastic weather effects and, in the longer run, climate (Sæther, 1997). Because climate and weather patterns are typically autocorrelated in space (Koenig, 1999), Moran (1953) suggested that fluctuations of spatially distinct populations should be synchronized similarly if their dynamics are driven by these weather patterns, and when assuming equal (log‐linear) density dependence (i.e. the ‘Moran effect’, Royama, 1992). Spatial population synchrony has indeed been observed across large areas and in a wide range of taxa (e.g. Bjørnstad, Ims, & Lambin, 1999; Koenig, 1999; Liebhold, Koenig, & Bjørnstad, 2004; Sæther et al., 2007; Walter et al., 2017), yet it has often proved difficult to disentangle the synchronizing role of weather from those of dispersal and trophic interaction effects (Engen, Lande, & Sæther, 2002; Engen, Lande, Sæther, & Bregnballe, 2005; Ims & Andreassen, 2000; Lande, Engen, & Sæther, 1999; Liebhold et al., 2004; Moran, 1953; Royama, 1992; Ydenberg, 1987). Nevertheless, several studies have now demonstrated that the Moran effect is often the main cause of spatial co‐fluctuations in population abundances (e.g. in feral sheep Ovies aries, Grenfell et al., 1998; roe deer Capreolus capreolus, Grøtan et al., 2005; caribou and reindeer Rangifer tarandus, Post & Forchhammer, 2004; passerine birds, Sæther et al., 2007; plants, i.e. an analogue to the Moran effect, Defriez & Reuman, 2017; Koenig & Knops, 1998, 2013; marine plankton, Defriez, Sheppard, Reid, & Reuman, 2016).

While fluctuations in asynchronous local populations will be cancelled out in large‐scale population sizes, a high level of population synchrony may lead to mass extinctions, for instance under extreme events (Heino, Kaitala, Ranta, & Lindstrom, 1997). Thus, strong, weather‐driven population synchrony is expected to increase extinction risk at the ‘meta‐population’ or species level (Engen et al., 2002; Heino et al., 1997; Ranta, Kaitala, Lindstrom, & Linden, 1995) under widespread environmental change, such as global warming (Post & Forchhammer, 2002). However, the observed strength and spatial scale of synchrony in population dynamics is usually much lower than that of climate itself (Sæther et al., 2007, 2003; Stenseth et al., 1999; Walter et al., 2017). A multitude of other drivers, which vary spatially, may also act on population growth rates or the ecological carrying capacity, desynchronizing population dynamics (Engen & Sæther, 2005). Even the extent to which different populations respond to a specific weather driver often varies locally due, for instance, to habitat heterogeneity (Anders & Post, 2006; Engen & Sæther, 2005; Post, Brodie, et al., 2009; Sæther, 1997). This may lead to uncorrelated dynamics in space that could moderate the effects of common environmental drivers, such as weather and climate, at the level of meta‐population or species (Hilborn, Quinn, Schindler, & Rogers, 2003), thereby enhancing species persistence (Anders & Post, 2006). Accordingly, to understand climate change effects at large spatial scales, we should assess how spatial heterogeneity in the environment or in the ecological responses averages out the effects of shared climate forcing acting across populations (Post, Forchhammer, et al., 2009; Sheppard, Bell, Harrington, & Reuman, 2016).

The relative contribution of weather fluctuations to animal population dynamics seems to be particularly large in marginal habitats, such as at high latitudes (Ims & Ehrich, 2013), so Arctic tundra food webs tend to display relatively strong population synchrony within and across species (Hansen et al., 2013; Ims & Fuglei, 2005; Post & Forchhammer, 2002, 2004). The climate is changing more rapidly in the Arctic than in any other biome on earth, especially during winter (AMAP, 2017; Bintanja & Andry, 2017). As a consequence, the frequency of extreme warm spells and associated heavy rain‐on‐snow (ROS) events in winter is increasing, causing dramatic changes to the snow‐pack characteristics (Hansen et al., 2014; Kohler & Aanes, 2004; Peeters et al., 2019; Putkonen & Roe, 2003) and, hence, the feeding conditions for overwintering herbivores. Especially in permafrost conditions, heavy ROS percolating through the snow may lead to the formation of basal ice (‘ground‐ice’), encapsulating the vegetation. Such ice‐locked tundra has potentially ecosystem‐wide consequences (Hansen et al., 2013). For instance, icing events can cause overwinter body mass loss, reduced skeleton growth, and reduced survival and fecundity in muskox (Ovibos moschatus) and reindeer (Albon et al., 2017; Berger, Hartway, Gruzdev, & Johnson, 2018; Stien et al., 2012). These effects of fluctuations in ROS, which can interact with density‐dependent mechanisms (Hansen et al., 2019), may result in large annual fluctuations in animal densities (Forbes et al., 2016; Forchhammer & Boertmann, 1993; Hansen, Aanes, Herfindal, Kohler, & Sæther, 2011; Kohler & Aanes, 2004; Miller & Gunn, 2003). However, such negative influence of warmer and wetter winters may—at least locally—be counteracted by the positive impacts of a shorter winter and a longer and warmer plant growth season (Albon et al., 2017; Tews, Ferguson, & Fahrig, 2007), causing an increase in ecological carrying capacity due to increased primary production (Van der Wal & Stien, 2014) and ‘Arctic greening’ (Jia, Epstein, & Walker, 2003; Vickers et al., 2016).

Disentangling such contrasting climate change effects on population dynamics is difficult, regardless of spatial scale (for a simulation study of caribou, see Tews et al., 2007). Arctic ungulate populations are typically influenced by a variety of other factors, such as hunting, insect harassment, infrastructure development, and predators, and these impacts often vary spatially (Gunn, Russell, White, & Kofinas, 2009). Poor data quality and discontinuous time series also characterize many monitored populations. Few long‐term, continuous time series have therefore been available with sufficient spatial replication to investigate spatial population synchrony (but see Post & Forchhammer, 2004, based on harvest records), and the buffering role of spatial heterogeneity, either in the environment or ecological responses. Here, we take advantage of such a rare, spatially replicated dataset, to analyse spatiotemporal patterns of population dynamics in wild Svalbard reindeer (Rangifer tarandus platyrhynchus) in the high Arctic, a hot‐spot of climate change. We ask: what is the role of the Moran effect in the large‐scale net outcome of recent climate change for this northernmost Rangifer subspecies? We assess this by examining (a) local and total population size trends for 10 reindeer populations; (b) how annual fluctuations in weather (i.e. ROS), and their pattern of spatial autocorrelation, synchronize reindeer dynamics in space (Moran, 1953); and (c) how the effects of this synchronization are moderated by spatial heterogeneity in climate, climate trends (i.e. long‐term changes in weather), and density‐dependent climate responses. In particular, we account for local variation in strength of density dependence and how it may interact with ROS effects on fecundity, mortality and population growth rates (Hansen et al., 2019).

MATERIALS AND METHODS

Study system

The climate in Svalbard is mild for the latitude (74–81°N, 10–35°E), with mean annual, summer (July to August) and winter (November to April) temperatures at Svalbard Airport of −6.7, 5.3 and −13.9°C respectively. Conditions are generally oceanic, yet dry (mean annual precipitation = 190 mm), and become increasingly dry with distance from the coastline (Van Pelt et al., 2016). Furthermore, weather patterns differ considerably between western parts, which are strongly influenced by the warm Gulf Stream, and the east, which is more influenced by the Arctic Ocean and sea‐ice cover (Johansen, Karlsen, & Tømmervik, 2012; Sakshaug, Johnsen, & Kovacs, 2009).

The Svalbard tundra is characterized by a simple food web, with few species. The wild Svalbard reindeer is the only mammalian herbivore, besides a small local population of sibling voles Microtus levis (Ims & Ehrich, 2013). The reindeer are largely solitary and virtually free from predation, although rare attacks by polar bears (Ursus maritimus) on sick or weak reindeer have been observed (Derocher, Wiig, & Bangjord, 2000). Thus, direct density dependence and large annual variations in weather conditions—notably the amount of ROS, but also the length and warmth of the summer—shape the reindeer's body condition and vital rates (Albon et al., 2017; Solberg et al., 2001; Stien et al., 2012). This, in turn, causes large fluctuations in population abundances from year to year (Albon et al., 2017; Hansen et al., 2011, 2013; Kohler & Aanes, 2004; Solberg et al., 2001). There is some evidence of spatial population synchrony (Aanes et al., 2003), most likely explained by spatially autocorrelated fluctuations in weather rather than dispersal, which is strongly restricted by natural barriers such as open sea, glaciers and steep mountains. Genetic analyses also indicate significant population differentiation across very small distances (Côté et al., 2002), largely confirming the observation that Svalbard reindeer are resident within annual home ranges of only a few km2 (Tyler & Øritsland, 1989). Overall, individuals do not move between neighbouring populations, except under special circumstances, such as following overgrazing or extreme ROS events (Hansen, Aanes, & Sæther, 2010; Loe et al., 2016). Combined with the strong climatic signals observed in vital rates (Albon et al., 2017), the characteristics of this study system provide a rare opportunity to investigate the impacts of climate‐induced environmental change on the population dynamics and trends of a high‐Arctic ungulate, both locally and at a meta‐population scale.

Data

We obtained daily time series of mean temperature and total precipitation data from five manned weather stations 14–410 km apart (Figure 1). Data from Barentsburg were downloaded from Tutiempo Network, S.L. (https://en.tutiempo.net/climate/ws-201070.html), while data from Hornsund were provided by the Institute of Geophysics of the Polish Academy of Sciences. For the three remaining weather stations (i.e. Svalbard Airport, Ny‐Ålesund and Hopen), we downloaded data from the Norwegian Meteorological Institute (http://eklima.met.no).

Figure 1

Map of Svalbard (inset) and the study area. Crosses show the locations of weather stations, coloured polygons delineate the monitored reindeer populations (see also Figure 2) and black dots show the position of the Ny‐Ålesund (Brøggerhalvøya) and Svalbard Airport (Adventdalen) weather stations

Figure 2

Time‐series data and trend estimates. Annual fluctuations and temporal trends in rain‐on‐snow (ROS) and reindeer population parameters in Svalbard are shown for the period used for population synchrony analyses (1997–2015). (a) Annual ROS (mm) recorded at five weather stations. (b–e) Annual reindeer fecundity (proxy), mortality (proxy), population sizes and growth rates based on summer ground and helicopter counts. (f, g) Linear trend estimates of ROS and reindeer population sizes over the study period. Whiskers indicate 95% confidence interval. In (a–e), solid black lines show ‘across Svalbard‐scale’ estimates from linear mixed regression models

Reindeer population monitoring data were obtained using multiple census methods across 10 reindeer populations. Total counts of live reindeer, calves and carcasses (easily visible as patches of white fur in the terrain) were carried out annually by helicopter in early August 1997–2015 across six populations (Figure 1). Five of these six populations are hunted in the fall with a fairly constant annual offtake corresponding to <5% of the local population size. In addition, total population counts were performed annually on foot in late June/early July 1979–2015 in Adventdalen (Hansen et al., 2013; Tyler, Forchhammer, & Øritsland, 2008) and included the number of calves and carcasses. Finally, following the reintroduction of 12 reindeer to Ny‐Ålesund, Brøggerhalvøya, in 1978 and the subsequent population irruption and spatial expansion southwards (Aanes, Sæther, & Øritsland, 2000), three semi‐isolated populations in this northwestern area have been monitored by annual total ground counts (including number of calves). Counts were conducted by snowmobile in early spring (i.e. late March to April 1978–2015 on Brøggerhalvøya) and on foot in summer (i.e. July; Brøggerhalvøya 1999–2015; Sarsøyra 2000–2015, except 2011; Kaffiøyra 2002–2015, except 2011–2012) (Hansen et al., 2011; Le Moullec et al., 2017). The two seasonal measures of population size on Brøggerhalvøya were highly correlated (r = .92, Figure S1). For Brøggerhalvøya, we used (a) the summer count time series in the analysis of Svalbard‐scale population dynamics and patterns of population synchrony (see below), to maximize methodological consistency; and (b) the spring count time series for the comparison of population dynamics with Adventdalen, to maximize time‐series length.

Svalbard‐scale reindeer parameters

Based on all 10 reindeer populations, we estimated overall annual abundances and vital rates, and their temporal trends. For consistency, and because most population time series were established in 1997, we restricted the study period to 1997–2015. Note that, because of missing years in some populations, simple averaging or adding of population sizes was not possible. For each population, we obtained annual proxies of ‘fecundity’ by dividing the number of calves by the total number of reindeer, and ‘mortality’ by dividing the number of carcasses by the total number of reindeer the previous year (see Peeters et al., 2017). Note that carcass counts are likely to give underestimates of mortality rates, although presumably consistent between years. Annual estimates of the ‘across Svalbard‐scale’ fecundity and mortality (i.e. across all populations) were then obtained through linear mixed effects regression models of the local fecundity or mortality estimates (on a logit scale) as a function of year (as a fixed factor) and population (random intercept effect). These regressions were weighted by the population‐specific mean population size (over each time series) to allow larger populations to contribute more to the across Svalbard proxies than smaller ones. Estimates of across Svalbard‐scale proxies of relative abundance (‘population size’ N Svalbard) and, thereby, population growth rates (R Svalbard = ln(N Svalbard,/N Svalbard,t)), were also obtained using a similar linear mixed effects regression model of local population sizes N.

Estimating the contribution of Svalbard‐scale climate to spatial synchrony

We analysed spatial synchrony in fecundity, mortality and population growth rates R between the 10 local populations by implementing a nonparametric covariance function, which uses smoothing splines to estimate spatial covariance as a function of distance (Bjørnstad & Falck, 2001). This method estimates spatial synchrony based on pairwise correlations and distances among populations. We set the maximum degrees of freedom for the smoothing spline to 4 (2 for mortality due to low sample size), and confidence intervals around the nonparametric curve were based on bootstrapping with 1,000 iterations. We used ≥5 years of data overlap as a criterion to calculate pairwise correlations, and the available population time series resulted in 45 (8–324 km), 15 (8–64 km) and 45 (8–324 km) pairwise correlations for reindeer fecundity, mortality (data lacking for Brøggerhalvøya, Sarsøyra, Kaffiøyra and Svenskøya; Table S2), and population growth rates respectively.

Rain‐on‐snow and population density are clearly the main drivers of annual fluctuations in vital rates and abundance (Albon et al., 2017; Hansen et al., 2011, 2013, 2019; Kohler & Aanes, 2004; Solberg et al., 2001; Stien et al., 2012), and ROS and density effects may also interact (Hansen et al., 2019; see also Coulson et al., 2001). More specifically, the ROS effect diminishes at very low population densities because available ice‐free patches per capita will anyways be sufficient, while at high densities, ‘all’ reindeer will suffer from a multiplicative effect of ROS on resource competition (Hansen et al., 2019). We therefore tested the extent to which density‐dependent ROS effects contributed to the synchrony in fecundity, mortality and population growth rates by also estimating the correlations between populations after accounting for this climate–density interaction effect, i.e. the correlations in model residuals. ROS was calculated as the annual amount of rain in winter (November − April +1), i.e. ln(mm precipitation + 1, to avoid log of zero) falling at temperatures ≥1°C. For each local population, we fitted generalized linear regression models of fecundity and mortality (quasi‐binomial family, logit link), and linear regression models of population growth rates, as functions of ROS in interaction with detrended ln population size. As nearby weather station data were not available for most populations (Figure 1), we used a common ‘across Svalbard‐scale’ ROS estimate obtained from a linear mixed effects regression model of ROS at the five weather stations as a function of year (as fixed factor), and with weather station as random intercept effect. The residuals from the population‐specific models of fecundity, mortality and population growth rates were then analysed for pairwise correlations and spatial synchrony across the range of distances. The contribution of density‐dependent ROS effects to the observed synchrony in fecundity, mortality and population growth was estimated based on the difference in nonparametric bootstrapped replicates of average synchrony (n = 1,000) before and after accounting for the density‐dependent ROS effect.

Although any (temporal variation in) spatially correlated monitoring efforts could ‘artificially’ increase the measured synchrony, the different monitoring methods between populations reduce this risk. For the counts performed on the ground, the accuracy seems relatively high (i.e. less biased and more precise; Le Moullec et al., 2017; Tyler et al., 2008). For instance, the correlation between early spring and summer population size on Brøggerhalvøya, counted on snow mobiles and on foot, respectively, was very high (r = .92). Likewise, the correlation in population size between our Adventdalen count data and independent counts performed by Tyler et al. (2008) in the same valley in 2001–2007 was also very high (r = .97), indicating high precision of the counts.

In‐depth analysis of climate effects on the dynamics of two focal populations

To obtain a more detailed understanding of local population dynamics and trends, and the causes and consequences of spatial heterogeneity, we analysed the population dynamics of the two longer reindeer abundance time series from Adventdalen (1979–2015) and Brøggerhalvøya (1978–2015). These two populations are each near a manned weather station (Figure 1), enabling the inclusion of local weather conditions in the analyses. We first looked for linear trends and sudden break points in the mean (i.e. ‘regime shifts’) weather conditions, by implementing the binseg function (with segment length restricted to minimum 5 years) in the changepoint package (Killick, Haynes, & Eckley, 2016) for R (R Core Team, 2018).

For the ‘continental’ population in Adventdalen, central Spitsbergen (Figure 1), we fitted a linear regression model of population growth rates against the following local‐scale predictors: density N, ROS, the interaction N × ROS, annual amount of snowfall (in mm, calculated as total precipitation (mm) falling at temperatures <1°C; Solberg et al., 2001), mean summer temperature (July to August; Hansen et al., 2013), and winter length (see also Albon et al., 2017). Winter length was estimated as the number of days between the time when the 10‐day running mean of daily temperature in the fall fell below freezing (0°C) and stayed below for a minimum of 10 consecutive days, and the time in spring +1 when the 10‐day running mean of daily temperature was above freezing and stayed above for a minimum of 10 consecutive days (Le Moullec, Buchwal, Wal, Sandal, & Hansen, 2019).

For comparison of climate effects between the two populations, a similar model was then applied to the coastal Brøggerhalvøya population, adding a two‐level factor ‘period’ (irruptive phase and post‐irruptive phase, see Results) in the N × ROS interaction (i.e. a three‐way interaction, period × N × ROS) to account for a possible change in density dependence in this irruptive population from the crash year in 1994 and onwards (Kohler & Aanes, 2004). All analyses were performed in R (R Core Team, 2018), implementing the packages lme4 for mixed‐effects models (Bates, Maechler, Bolker, & Walker, 2015) and ncf for the nonparametric covariance function (Bjørnstad, 2016).

RESULTS

Spatiotemporal patterns in weather and reindeer dynamics across Svalbard

During 1997–2015, data from weather stations across Svalbard showed an overall positive trend in both ROS (Figure 2a,f) and summer temperature (Figure S2). This indicates worsening winter feeding conditions due to more frequent icing (shown through in situ ice measurements, Figure S3a), yet, at the same time, improved food abundance and carrying capacity due to higher plant productivity in the warmer summers (inferred by the Normalized Difference Vegetation Index, Figure S3b). Although reindeer population size trajectories (Figure 2d) showed locally contrasting linear trends, varying from negative to positive (Figure 2g), there was evidence of a statistically significant increase over time of the ‘across Svalbard‐scale' reindeer abundance index (Figure 2g), weighted for local mean population sizes. Fecundity, mortality and population growth rates were all characterized by large annual fluctuations (Figure 2b,c,e), resulting in saw‐shaped patterns indicating direct density dependence. For instance, population growth rates R were particularly small in 2005, 2007 and 2009—i.e. years with fairly high levels of ROS at most weather stations, as well as high estimated mortality and low fecundity. Conversely, population growth rates were relatively large in 2006, 2008 and 2010. Thus, at the Svalbard scale, correlations with annual ROS were negative for fecundity (Pearson's r = −.63, p < .01), positive for mortality (r = .56, p < .05) and negative for population growth rates (r = −.44, p = .07) (Figure 3). Linear regression models accounting for potential density dependence in this ROS effect confirmed the overall presence of direct density dependence but also indicated inconsistent patterns of climate–density interactions across populations (Table S1).

Figure 3

Rain‐on‐snow (ROS) effects on Svalbard reindeer. Annual (a) fecundity (proxy), (b) mortality (proxy) and (c) population growth rates of the ‘across Svalbard‐scale’ reindeer population during 1997–2015 as a function of rain‐on‐snow. Note that the y‐axis in (b) is on natural log‐scale

Rain‐on‐snow showed a high mean pairwise correlation between weather stations (ρ ROS = 0.62 [95% confidence interval = 0.48, 0.76]), and the correlation decreased with distance (Figure S4). Similarly, pairwise correlations between the reindeer populations in annual fecundity (mean pairwise correlation ρ Fecundity = 0.52 [0.39, 0.67]), mortality (ρ Mortality = 0.83 [0.77, 0.89]; note the shorter extent of distances) and population growth rates (ρ R = 0.23 [0.02, 0.43]) were generally positive and decreased with increasing distance between populations (Figure 4a–c; Table S2). When accounting for only an additive effect of ROS in population‐specific linear models, the average pairwise population correlations (in residuals) remained similar for mortality (ρ resid = 0.83 [0.76, 0.91]) and population growth rates (ρ resid = 0.21 [0.01, 0.39]), but were lower for fecundity (ρ resid = 0.44 [0.28, 0.62]). However, when accounting for a density‐dependent ROS effect (N × ROS, Table S1), correlations between populations were lower for population growth rates (ρ resid = 0.15 [−0.18, 0.48]) and, in particular, mortality (ρ resid = 0.67 [0.49, 0.88]), while less so for fecundity (ρ resid = 0.47 [0.26, 0.64]) (Figure 4d–f). Thus, based on whether the lower 5% percentile (i.e. an a priori expected change, ‘one‐tailed’) of the change in pairwise correlation obtained from bootstrapping was above or below zero (Figure 4g–i), this change in correlation due to the N × ROS effect was far from ‘significant’ for fecundity (mean difference = 0.06 [−0.12, 0.26]), significant for mortality (0.44 [0.17, 0.79]) and marginally non‐significant for population growth rates (0.12 [−0.02, 0.25]). Note that synchrony models with other weather variables (analyses not presented) or with only density N included as covariate caused relatively little or no reduction in average regional synchrony (fecundity: ρ resid = 0.52 [0.31, 0.71]; mortality: ρ resid = 0.73 [0.60, 0.89]; population growth rates: ρ resid = 0.23 [0.01, 0.47]). Thus, we conclude that density‐dependent ROS effects were the main detected contributor to the observed synchrony in population dynamics, yet marginally non‐significant, and mainly operating through synchronization of mortality.

Figure 4

Spatial synchrony in reindeer population parameters across Svalbard, and the impact of rain‐on‐snow (ROS). (a–c) Pairwise correlations between populations in (a) fecundity (proxy), (b) mortality (proxy) and (c) population growth rates during 1997–2015, as function of distance (see also Table S2). Note the different axes scales for the different population parameters. The thick solid line represents the smoothed spline function, while the average regional synchrony (i.e. mean correlation) and bootstrapped 90% confidence intervals are shown with dashed and dotted blue lines respectively. (d–f) Pairwise correlations between populations of residuals (when accounting for density‐dependent ROS effects; Table S1) of (d) fecundity, (e) mortality and (f) population growth rates, as functions of distance. (g–i) The contributions of the density‐dependent ROS effect to synchrony in (g) fecundity, (h) mortality and (i) population growth rates, shown as histograms of the difference in nonparametric bootstrap replicates of average synchrony (n = 1,000) before and after accounting for the density‐dependent ROS effect. Thick (50%) and thin (5% and 95%) blue dashed lines illustrate the percentiles

Causes and consequences of spatial heterogeneity

Can the same climate or weather parameters cause synchronous annual population dynamics while simultaneously causing diverging long‐term trends? Population growth rates were indeed positively correlated (r = .46, Table S2) between our two focal reindeer populations (i.e. Adventdalen [‘continental’] and Brøggerhalvøya [coastal], 126 km apart, Figure 1), from which we had both local weather data and long‐term, ground‐based count data. However, some small yet important differences in climate trends between these two populations have become evident since the 1990s (Figure 5a,b). Firstly, the linear regressions indicated a significant increase in the annual amount of ROS (on natural log‐scale) during 1978–2015 at the weather station on Brøggerhalvøya (i.e. Ny‐Ålesund; β = .055, t = 2.44, p < .05), but not at the station in Adventdalen (i.e. Svalbard Airport; β = .030, t = 1.63, p = .11). A regime shift with overall higher annual ROS since around year 2000 was detected in Ny‐Ålesund (Figure 5a), but not at Svalbard Airport (see Peeters et al., 2019 for a longer time span). Secondly, while linear regression models indicated a summer temperature increase over time for both Ny‐Ålesund (β = .028, t = 3.09, p < .01) and Svalbard Airport (β = .053, t = 4.85, p < .001), the slope was twice as steep at the latter station. Thus, a sudden regime shift with elevated summer temperature levels from around 1998 was detected at both weather stations, but in this case the increase was greater at Svalbard Airport (ca. 1.3°C) than in Ny‐Ålesund (ca. 0.9°C, Figure 5b). We found no evidence for other differences between Adventdalen and Brøggerhalvøya in temporal trends of potential drivers of reindeer population growth (Figure S2), i.e. winter length (equally strong linear decrease in both areas) and total snowfall (no strong trends).

Figure 5

Time‐series data and population trends in two focal reindeer populations. Annual fluctuations in (a) rain‐on‐snow (ROS), (b) summer temperatures and (c) standardized reindeer population sizes, and the estimates of (d) population size trends (with varying trend starting year; whiskers showing 95% CI) over the period 1978–2015 in Brøggerhalvøya (blue; Ny‐Ålesund weather station) and Adventdalen (red; Svalbard Airport weather station). Horizontal dashed lines in (a) and (b) denote detected ‘regime shifts’, i.e. change points in mean. Population size trends in (c) are shown as solid lines, for the period after the irruptive phase in Brøggerhalvøya

The differences in the relative magnitude of temporal change in ROS and summer temperature between populations contributed to the long‐term reindeer population trajectories in these two focal populations. Overall, the effects of local weather variables on population growth rates were broadly similar (i.e. similar signs, but with different strengths) in Adventdalen and on Brøggerhalvøya (Figure 6). Thus, linear regression models (for a cross‐validation procedure, see Figure S5), which were run separately for each population due to the period/irruption factor on Brøggerhalvøya, suggested an overall negative ROS effect on population growth rates. The ROS effect was stronger at high reindeer densities, and, importantly, this interaction effect between ROS and density was stronger on Brøggerhalvøya than in Adventdalen (Figure 6). In contrast, the positive effect of summer temperature was stronger in Adventdalen than on Brøggerhalvøya. Snowfall amount tended to have an overall negative effect, while there was no evidence of an influence of winter length in either population.

Figure 6

Climate–density effects on two focal reindeer populations. Standardized parameter estimates are shown for effects of weather and population size (N) on reindeer population growth rates on Brøggerhalvøya (blue; estimates from the post‐irruptive phase are shown) and in Adventdalen (red). Whiskers show 95% confidence intervals

As expected from the slightly different effects of local weather on population growth and the local differences in climate trends, we found contrasting abundance trends between the two populations. More specifically, the population size trend was strongly positive in Adventdalen (irrespective of trend starting year; Figure 5d) where the local weather station data indicated no significant increase in ROS (Figure 5a) but a strong regime shift in summer temperatures from the late 1990s (Figure 5b). In contrast, the trend slopes for the post‐irruptive Brøggerhalvøya population—which experienced weaker summer warming but stronger ROS increase—tended to be negative when the trend starting year was before ~2001 (Figure 5d). Since around that time, the population size first declined for then to stabilize around relatively low densities (Figure 5c). Accordingly, there was a consistent difference in estimated trends of abundance between the two focal populations (i.e. generally no overlap between confidence intervals and the mean estimate for the other population), irrespective of choice of trend starting year (Figure 5e).

DISCUSSION

This study has demonstrated that, although weather fluctuations tend to synchronize short‐term (i.e. annual) dynamics, the longer term dynamics of wild reindeer populations can be decoupled by spatially heterogeneous climate change effects. Annual fluctuations in the amount of ROS during winter contributed to synchronized mortality rates and, thereby, population dynamics across large distances in Svalbard (Figure 4), while paradoxically, the very same environmental driver contributed to opposite trends in local reindeer abundances (Figures 2g and 5d). We document how such diverging population trends can occur due to local variation in both the relative strength of winter versus summer climate change trends (Figure 5) and the density‐dependent responses to warming (Figure 6). However, the net outcome of improved carrying capacity versus worsening winter‐feeding conditions appears to be a growing total population of Svalbard reindeer (Figure 2d,g), in contrast to the recent patterns of decline observed in some other Rangifer populations (Mallory & Boyce, 2017; Vors & Boyce, 2009).

Post, Forchhammer, et al. (2009) pointed out that spatial environmental heterogeneity and scale‐dependent ecological responses are key to understanding and predicting net implications of the various climate change effects reported in Arctic tundra study sites. In caribou and reindeer, some observed population declines have indeed been linked to global environmental change, but the proposed drivers seem to differ locally and include habitat fragmentation (Festa‐Bianchet, Ray, Boutin, Côté, & Gunn, 2011), Arctic ‘shrubification’ (Fauchald, Park, Tømmervik, Myneni, & Hausner, 2017), increased insect harassment (Mallory & Boyce, 2017) and worsening winter‐feeding conditions (Forbes et al., 2016; Hansen et al., 2011). Also, the drivers of change may vary in time. For instance, the early population size trajectories on Brøggerhalvøya since the reindeer reintroduction in 1978 followed irruption‐like dynamics (Aanes et al., 2000), with habitat degradation in the early 1990s (Hansen, Henriksen, Aanes, & Sæther, 2007) that may have contributed to the 1994 population crash and, thereby, our trend estimates. The more recent negative trends, which appear to be a consequence of negative net effects of climate change, were in sharp contrast to the Adventdalen population increase during the same period. Accordingly, because of the spatial and temporal diversity of factors influencing Rangifer population growth rates and ecological carrying capacity, Gunn et al. (2009) warned against generalizing observed local abundance trends and impacts of climate change. The concern that there is a ‘global’ decline in Rangifer abundance due to increased anthropogenic impact (Vors & Boyce, 2009) also lacks quantitative support. Most Eurasian reindeer populations show no sign of decline, but rather the opposite (Uboni et al., 2016), which also appears to be the case for Svalbard reindeer. Importantly, however, the substantial local variation in abundance trends calls for great caution when extrapolating results from single‐population studies of long‐term climate change effects, even on such small spatial scales (in a Rangifer context).

The extinction risk of a meta‐population or (sub)species depends in part on the degree of population synchrony, and thereby the strength of climate as a synchronizing agent (Engen et al., 2002; Heino et al., 1997; Ranta et al., 1995). The estimated regional level of synchrony of Svalbard reindeer dynamics was high for single vital rates, but this was less apparent in the population growth rates, which are more strongly influenced by sampling error. However, as expected based on the observed pattern of autocorrelation in ROS (Figure S4), synchrony declined with distance for all these population parameters. Uboni et al. (2016) found no such pattern of spatial scaling in their pairwise population correlations in population growth rates across large distances (up to ~7,000 km) in Eurasian reindeer, but for our range of distances (up to ~324 km) their average regional synchrony (ca. 0.24) was very similar to ours (ρ R = 0.23). Their lack of evidence for climate effects on population dynamics and synchrony could, at least in part, be a result of not accounting for local variation in density dependence and climate–density interactions, the importance of which has been clearly demonstrated here. Thus, even in ‘simple’ trophic systems such as Svalbard, where the influence of human presence is low and where few external factors other than weather fluctuations are expected to influence annual reindeer population fluctuations (Ims & Ehrich, 2013; Tyler et al., 2008), it has proven crucial to account for density dependence in studies of population dynamics (Aanes et al., 2000; Albon et al., 2017; Hansen et al., 2011; Tyler et al., 2008). In this study, direct density dependence was generally present across populations and population parameters, but there was little evidence of a consistent pattern of ROS–density interaction effects (as found in Hansen et al., 2019) across the 10 populations. The contribution of this interaction effect to population synchrony was therefore not huge; its synchronizing effect was significant for mortality, yet negligible for fecundity (which seemed better explained by an additive ROS effect), and, thereby, only close to significant for population growth rates. This lack of significance, and the local variation across Svalbard in estimated effects of ROS, and its interaction with density, could in part also be due to the slightly ‘noisy’ and short datasets used for this analysis. Furthermore, the Svalbard‐scale ROS variable probably reflected local conditions better for some populations (close to one or more weather stations) than others, e.g. the remote Svenskøya population in eastern Svalbard (Figure 1). Nevertheless, the in‐depth analysis of the two high‐quality datasets confirmed that important local differences in density regulation indeed occur. Accordingly, when density dependence and its potential interaction with climate effects differ locally, the assumption of the Moran effect that there is a similar (log‐linear) density dependence among populations does not hold. Contrasting density dependence desynchronizes population dynamics (Engen et al., 2005) and should clearly be accounted for in the analysis of climate effects on population synchrony.

One major reason for the overall poor understanding of density‐dependent and ‐independent drivers of Arctic ungulate population dynamics, and how they vary in space, is that time series of population monitoring data are notoriously discontinuous (Mallory, Campbell, & Boyce, 2018; Uboni et al., 2016; Vors & Boyce, 2009). This complicates any attempt to quantify how climate change shapes the dynamics and trends at the level of population/herd, meta‐population, and subspecies, as well as patterns of spatial population synchrony. The accuracy of ‘total population count’ data is also rarely assessed, but aerial counts are often assumed to be uncertain and underestimates of the actual population sizes (Poole, Cuyler, & Nymand, 2013). We have no reason to believe that the helicopter counts used here differ in that sense; however, ground counts appear highly accurate (Le Moullec et al., 2017). Uncertain aerial counts contribute to noise in the dataset, reducing the measured level of population synchrony and the ability to detect effects of climate on population growth rates and their spatial synchrony. However, because of the solitary, non‐migratory behaviour of Svalbard reindeer, and the rather small, well‐defined study areas in open tundra landscapes, we believe that these helicopter counts generally captured major changes occurring from year to year, as well as the long‐term trends. This is especially the case for the vital rates, which were highly correlated across relatively nearby populations, as also indicated in other subspecies of Rangifer (Hegel, Verbyla, Huettmann, & Barboza, 2012).

Climate projections suggest that the rapid regime shifts observed in winter (Peeters et al., 2019) and summer climate in Svalbard represent a bellwether for future changes in other parts of the Arctic (IPCC, 2013). Based on our results (see also Albon et al., 2017), it seems likely that continued warming leads to a changing relative importance of seasons for reindeer population dynamics, locally and at the meta‐population level. ROS events and the phenomenon of ice‐locked pastures are anticipated to increase further in frequency, magnitude and spatial extent, and not only in Svalbard (Bintanja & Andry, 2017; Rennert, Roe, Putkonen, & Bitz, 2009). Because the local impact of such events may depend on the population state, often with a multiplicative effect at high density due to the strong resource competition, a change in event frequency will influence population stability and local extinction probability in Rangifer (Hansen et al., 2019). Given the local variation in density–climate interactions found here, altered patterns of large‐scale population dynamics and synchrony can therefore also be expected.

Other snow pack changes than basal ice formation, notably snow depth (Tyler, 2010), may modify the effects of ROS events on winter‐feeding conditions (Peeters et al., 2019), and snowfall amount per se was also found to have a negative effect on both the Adventdalen and Brøggerhalvøya populations. Furthermore, a progressively longer snow‐free season (Albon et al., 2017) and improved carrying capacity due to summer warming and ‘greening’ (which is likely highly autocorrelated in space; Milner, Stien, & Van der Wal, 2018) can have positive effects on population growth, both in the short and long run, as observed in Adventdalen (but see Fauchald et al., 2017 for the low Arctic). As demonstrated in this study, the net effect of changes in environmental drivers linked with summer or winter warming could be expected to vary spatially due to local variation in climate trends and density‐dependent responses, decoupling Rangifer population dynamics in space and, potentially, reducing extinction risk (e.g. Heino et al., 1997).

The results from this study have ecological implications far beyond the levels of population, meta‐population and subspecies. Like many other deer species (Côté, Rooney, Tremblay, Dussault, & Waller, 2004), Svalbard reindeer play a key role in the ecosystem through top‐down and bottom‐up effects of trophic interactions. Impacts of grazing, trampling and high nutrient turnover rates resulting from increased reindeer densities may cause vegetation state transitions from lichen‐ to moss‐ and graminoid‐dominated vegetation (Hansen et al., 2007; Van der Wal, 2006). Thus, similar or contrasting long‐term trends in abundance between reindeer populations, distributed across large areas, may influence the spatial autocorrelation in vegetation structure, with implications for other herbivores and ecosystem components (Ims & Ehrich, 2013). Importantly, reindeer carcasses represent a major food source for the Svalbard population of Arctic fox (Vulpes lagopus), and stochastic ROS events therefore cause lagged synchrony in the two species’ population growth rates (Hansen et al., 2013). The Arctic fox does not kill reindeer but is the most important predator of the other overwintering herbivore, the Svalbard rock ptarmigan (Lagopus muta hyperborea), as well as migratory birds. This key ecosystem role of Rangifer through trophic interaction effects is not unique to Svalbard (Ims & Ehrich, 2013; Legagneux et al., 2014). Thus, how future climate change affects the level and spatial scaling of abundance synchrony in caribou and reindeer and, in turn, their food plants and consumers, are likely to impact spatiotemporal ecological dynamics occurring at the entire tundra community level.

Click here for additional data file.