Spatiotemporal variation in climatic conditions across ecosystems

: Environmental variation in time and space affects biological processes such as extinction risk and speed of adaptation to environmental change. The spatial structure of environmental variation may vary among ecosystems, for instance due to differences in the flow of nu -trients, genes and individuals. However, inferences about ecosystem spatial scale should also include spatial autocorrelation in environmental stochasticity, such as fluctuations in weather or climate. We used spatially structured time series (19−36 yr) on temperature from 4 different ecosystems (terrestrial, limnic, coastal sea and open ocean) to assess the spatiotemporal patterns of environmental variation over large geographical scales (up to 1900 km) during summer and winter. The distance of positive spatial autocorrelation in mean temperature was greatest for the terrestrial system (range: 592−622 km), and shorter for the open ocean (range: 472−414 km), coastal sea (range: 155−814 km) and the limnic systems (range: 51−324 km), suggesting a stronger spatial structure in environmental variation in the terrestrial system. The terrestrial system had high spatial synchrony in temperature (mean correlation: winter = 0.82, summer = 0.66) with a great spatial scaling (>650 km). Consequently, populations of terrestrial species experience similar environmental fluctuations even at distances up to 1000 km, compared to species in the aquatic systems (<500 km). There were clear seasonal differences in environmental synchrony in the terrestrial and limnic systems, but less so in the other systems. Our results suggest that biological processes affected by environmental stochasticity occur at the largest spatial scale in terrestrial systems, but their magnitude depends on whether the process is affected by winter or summer conditions. on bootstrapping of spline models all of between and the suggesting strong evidence for seasonal differences the


INTRODUCTION
Environmental variation in time and space is the main driver affecting processes of ecological and evolutionary change, including population dynamics, ecosystem change, evolution and speciation. Ac -correspond closely to those of the environment. Accordingly, a thorough understanding of such spatiotemporal environmental dynamics is pivotal in order to predict spatial patterns of ecological and evolutionary processes and dynamics. This is particularly true in the face of the increasing human impact on ecosystems where the need for management and conservation actions at appropriate spatial scales is challenged by decreasing access to natural habitats (Moore et al. 2010, Knapp et al. 2017, Brennan et al. 2019. Ecosystems vary in the temporal and spatial scale of environmental and biological processes (Cole et al. 1991, Powell & Steele 1995. The recognition that open aquatic systems have a higher potential for dispersal and migration, and flow of nutrients or other components, than most terrestrial systems (Myers et al. 1997, Carr et al. 2003, has led to the perception that 'blue' marine systems have a greater temporal and spatial scaling of biological patterns and processes than 'green' terrestrial systems (Mayr 1954, Steele 1991, Carr et al. 2003, Vasseur & Yodzis 2004. These intrinsic physical factors are unquestionably important for processes affecting both population dynamics and gene flow. However, both theory (Frank 2005) and a large body of empirical studies (e.g. Steen et al. 1996, Myers et al. 1997, Stenseth et al. 1999, Jones et al. 2003, Tack et al. 2015 suggest that the spatial scale of environmental variation, such as weather or climate, is at least as important as flow of genes or individuals for stochastic ecological and evolutionary processes. It is therefore also necessary to consider the spatial structure of environmental variation when assessing the spatial scale of biological processes. For instance, regional extinction rates are related to how much dynamics of local populations are synchronised (Heino et al. 1997, Liebhold et al. 2004). If environmental stochasticity affecting population growth is synchronised over great spatial scales, the fluctuations in population abundances will be synchronised accordingly (Moran 1953), which may have cascading ecosystem effects (Elton 1924). Still, a thorough evaluation of spatial structure of such environmental stochasticity across systems is rarely done in comparative studies (Hansen et al. 2020).
Sensitivity to environmental variation varies among species (Bjørkvoll et al. 2012, and accordingly, the synchronising effect of the environment on population dynamics can differ among species (Marquez et al. 2019). Often this variation can be related to the fast−slow continuum of life histories (Oli 2004). Because the dynamics of species with fast and slow life histories are affected mainly by repro-duction and adult mortality, respectively (Oli 2004), the spatial structure of the environmental variation affecting population dynamics may differ. This occurs because reproduction and adult mortality often occur in different parts of the year, which in temperate and boreal regions is typically summer and winter, respectively. Accordingly, any seasonal difference in the spatial scaling of environmental conditions can lead to species-specific patterns in their spatial structure of ecological and evolutionary processes. Such differences may also occur if the same life history event (e.g. breeding) occurs at different times, even if the different species inhabit the same area, which is the case for many marine fish species (Olsen et al. 2010). Consequently, the spatial structure of environmental variation must be related to the critical period of a species (Saether et al. 1996), which further emphasises the need for detailed em pirical assessment of system-specific environmental scaling (Mokany et al. 2010).
Here we describe the spatial properties (patterns of spatial autocorrelation and scale of synchrony; Walter et al. 2017) of environmental variation in 4 systems: terrestrial, limnic, coastal sea and open ocean. These systems differ considerably in biological and physical characteristics (Steele 1985, Cole et al. 1991, Steele & Henderson 1994, Carr et al. 2003; however, they experience much of the same challenges with respect to biodiversity threats due to human activity (Carr et al. 2003, Brondizio et al. 2019. We used spatially structured long-term data on temperature to compare the spatial structure of environmental variation between these systems. Temperature is one of the most important environmental variables affecting individual growth, reproduction and survival in species across all taxa (Clarke 2017). Accordingly, all species have lower and upper thermal bounds (the thermal niche), and conditions outside these thermally viable envelopes will reduce fitness and population growth . Knowing the patterns of spatiotemporal variation in temperature across ecosystems may thus be a first step to assess general patterns of the spatial scale of biological processes.

Data
Daily mean temperature from weather stations in Norway and Sweden were downloaded from the European Climate Assessment database (Klein Tank et al. 2002). We excluded stations on remote islands. The daily temperatures were aggregated to mean temperature during winter (February−March) and summer (July−August) per station and year. The total number of stations with at least 10 yr of data were 156 for both summer and winter seasons (Fig. 1).
Data on water temperature from Norwegian lakes were obtained from the Norwegian Water Resources and Energy Directorate (NVE) and from Swedish lakes from the Swedish University of Agricultural Sciences database for Swedish lakes and watercourses (http://miljodata.slu.se/mvm/, in Swedish). These data were collected at different depths which were classified into 1 m (surface to 1 m), 2 m (1−2 m), 5 m (2−5 m) and 15 m (5−15 m) depth bins. The data were often sampled at short intervals, mainly during the months February, March and July−August. For each station, we calculated the seasonal mean value per year (same definition of seasons as above). The total number of stations with sufficient length of the time series (at least 10 yr) for summer and winter were 134 and 80, respectively, with fewer locations for the 15 m depth class.
For the coastal areas, sea temperature was ob tained from 9 stations along the Norwegian coast from the Institute of Marine Research in Norway (Aure & Østensen 1993, Saetre et al. 2003. Daily sea temperature was available from the following depths: 0.5, 50, 100 and 300 m. Data were aggregated to seasonal (see above) averages per year per depth. All stations had data from both summer and winter at all depth classes, except the 300 m depth (8 stations).
Since the 1970s, the Institute of Marine Research in Norway has conducted different annual surveys in the Barents sea where temperatures have been measured systematically at different depths using CTD casts (see e.g. Jakobsen et al. 1997, Stiansen & Filin 2007. The samples are approximately uniformly distributed in space, although the positions of these stations have not been constant between years. Sea temperature was categorised in depth groups as follows: surface (0−5 m); 50 (30−75 m); 100 (75−150 m); 200 (150−250); 300 (250−350 m). Because the locations were not exactly the same from year to year, we aggregated locations based on a stratified grid of equal area hexagons (size = 8100 km 2 , see Fig. 1 and Marquez et al. 2019). Temperature was then aggregated by season (see above) and depth class within each hexagon. The number of hexagons with temperature data (minimum 10 yr) during summer and winter was 72 and 60, respectively, but somewhat lower for the deepest depth class.
Hereafter, we refer to the spatial location of the temperature data as 'location' in all 4 systems. We restricted data to the years 1980−2015 for all systems in order to get comparable study periods. Moreover, we used only time series with a minimum of 10 yr, and when calculating correlation in time series be tween pairs of locations we excluded pairs with less than 10 yr of overlapping data. This was done to ob tain robust

Estimating spatial autocorrelation
The spatial autocorrelation of seasonal temperatures for a system was described by correlograms based on the spatial pattern of Moran's I based on mean seasonal temperatures over the study period for each location. I was then described as a function of distance between pairs of locations, modelled by a spline function with maximum number of knots = 4 (Bjørnstad & Falck 2001). In such correlograms, 0 represents the overall similarity in temperature across all locations within a system. Distances where I > 0 mean that the locations at these distances are more similar than what would be expected by chance, whereas values < 0 mean that locations are more dissimilar than what would be expected by a random distribution of the locations in space. The typical spatial pattern is that close lo cations are more similar than by chance (i.e. Moran's I > 0). The curve that describes the relationship between distance and Moran's I will cross 0 at a distance D I,0 , representing the distance for which there is no spatial autocorrelation in the temperature.

Estimating spatial synchrony in temperature
The spatial synchrony in temperature was assessed by calculating the pairwise correlation in the standardised first-order differential of the time series from a system and season (Hansen et al. 2020). For each pair of locations, we get a correlation, ρ, and a distance between the locations. Following the ap proach of Bjørnstad & Falck (2001), we used a spline function with maximum number of knots = 4 to describe how the correlation in the dynamics of temperature was related to distance, expecting a de crease in the correlation with increasing distance between locations (Bjørnstad et al. 1999). We calculated the mean correlation among all pairs of stations, ρ -, and based on the spline model, the distance at which the estimated curve crossed ρ -, λ ρ -.
For both the spatial autocorrelation and the synchrony, we ran a bootstrapping procedure (n = 1000) by drawing pairs of locations at random and refitting the spline (Bjørnstad & Falck 2001). We then used the 95% credible interval (CI) of differences in the bootstrap of 2 systems or seasons to assess the strength of evidence for system-specific or season-specific differences in spatial properties of environmental variation.

RESULTS
Estimates of ρ -(and hence λ ρ -), and to some extent also D I,0 , will be affected by the distribution of distances between locations within a given system. For instance, limnic (winter) and open ocean mainly have locations closer than 800 km, i.e. half the range as for terrestrial systems. Accordingly, in limnic and open ocean, ρwill be based on only close locations compared to terrestrial, and thus be biased towards higher ρand shorter λ ρ -. This can be accounted for by only including pairs of locations that were < 650 km from each other (see Fig. 3, dotted lines). Below we report the estimates for the complete range for a system; however, the 95% CIs to assess evidence for differences between systems or seasons are based on pairs of locations closer than 650 km (e.g. Table 1 and see Table S1 in the Supplement at www.int-res.com/ articles/suppl/c086p009_supp.pdf).

Spatial autocorrelation
There was a strong spatial autocorrelation in terrestrial temperature (Fig. 2), with closer locations being more similar than the average similarity across all stations, and locations far from each other being more dissimilar than expected by chance. This pattern was quite similar for summer and winter (summer D I,0 = 592 km, winter D I,0 = 622 km, 95% CI overlapped 0, Table 1). The limnic and coastal sea systems had a much weaker spatial structured autocorrelation in temperature compared to the terrestrial system ( Fig. 2; Table S1). The coastal sea had a stronger spatial structure during summer, with D I,0 > 750 km at several depths; however, the uncertainty in these estimates is high due to the rather few spatial locations ( Figs. 1 & 2, Table 1). The open ocean system had a spatial autocorrelation structure similar to the terrestrial, with high positive autocorrelation at short distances and negative autocorrelation at large distances. However, D I,0 was shorter than for the terrestrial system both in summer and winter and for all depths ( Fig. 2; Table S1).

Synchrony in the dynamics of temperature
In winter, the spatial synchrony in temperature was considerably higher for the terrestrial system (ρ -= 0.82) than for the 3 aquatic systems (ρ -< 0.64, Fig. 3, all 95% CI > 0; Table S1). The limnic systems had the lowest synchrony of all during winter, with ρ -< 0.17 at all depths. The coastal sea had winter ρfrom 0.32 to 0.64, and the open ocean was quite similar (winter ρfrom 0.33 to 0.46). The terrestrial system had considerably lower ρduring summer (ρ -= 0.66) compared to winter (Table 1), but it was still higher than in all the other systems ( Fig. 3; Table S1). Of the 4 systems, only the limnic system had consistently higher ρduring summer than in winter, particular at the shallowest depths (ρ -= 0.39, Fig. 3, Table 1). Coastal sea had somewhat lower ρduring summer (0.16−0.38), as did the open ocean (0.28−0.51, Fig. 3, Table 1).

Correlation between spatial autocorrelation and spatial synchrony
The correlation between the scaling of the spatial autocorrelation and the pattern of synchrony was analysed based on pairs with < 650 km distance to avoid bias due to difference in the extent of the study system. For the terrestrial system during winter, the estimated D I,0 was not estimable because of positive autocorrelation up to 650 km. Accordingly, we set this to 650 km for the terrestrial system during winter. The correlation between D I,0 and λ ρ -(both variables ln-transformed) was significantly positive when analysed across all seasons and systems (r 2 = 0.615, t = 3.82, df = 24, p = 0.001, Fig. 4). Omitting the coastal system, which had few points and great uncertainties in the measures of spatial structure (e.g. Fig. 2), showed an even stronger positive correlation between D I,0 and λ ρ -(r 2 = 0.870, t = 7.07, df = 16, p < 0.001). The correlation between D I,0 and ρwas significantly positive both when including coastal sea (r 2 = 0.549, t = 3.22, df = 24, p = 0.004), and excluding it (r 2 = 0.803, t = 5.39, df = 16, p < 0.001, Fig. 4). ρand λ ρ -were positively correlated when excluding coastal sea (r 2 = 0.733, t = 4.31, df = 16, p = 0.001), but not if coastal sea was included (r 2 = 0.329, t = 1.70, df = 24, p = 0.101, Fig. 4). These positive correlations suggest that the spatial autocorrelation in temperature to a great extent also reflects the spatial scale of synchrony in temperature.

DISCUSSION
Ecosystems worldwide and in all biomes face challenges related to human activity (IPBES 2019). This calls for conservation actions in order to maintain  Table 1. Seasonal differences in measures of the scale of spatial autocorrelation (D I,0 ), the mean synchrony ρand the spatial scaling of the synchrony λ ρ -in temperature in the 4 systems. Differences are presented as summer−winter, i.e. positive values mean higher spatial autocorrelation and synchrony in summer than in winter, and are based on bootstrapping (n = 1000) of spline models from pairwise distances (see Bjørnstad & Falck 2001) for locations <650 km for all systems, in order to allow for comparison when the range of distances differs between systems and seasons. Parentheses give the 95% credible intervals. Bold font highlights intervals that do not contain 0, suggesting strong evidence for seasonal differences in the spatial property of the environmental variation in system biodiversity and ecosystem functioning. The physical differences between many systems has raised the question whether it is possible to generalise conservation practices and strategies across systems, for instance between terrestrial and marine systems, which vary greatly with respect to processes such as dispersal distance and movement patterns (Carr et al. 2003, Mokany et al. 2010). Here, we show that the spatial structure of environmental variation also differs among systems. However, our results suggest that the scale of the spatial structure and synchrony is not necessarily greatest in marine systems (e.g. Carr et al. 2003). Temperature in the terrestrial systems had a high autocorrelation over large distances, which to a lesser extent was the case in the ocean systems, but not in the limnic system (Fig. 3). The generality of our finding relies to some extent on comparable data between systems. For instance, ρwas calculated as the mean correlation in the dynamics of temperature of all pairs of locations. For the ter-14 Fig. 2. Spatial autocorrelation in temperature during winter (dashed lines) and summer (solid lines) in 4 systems and at different depths in the aquatic systems. The vertical lines indicate the distance at which the spatial autocorrelation does not deviate from the overall spatial autocorrelation of temperature in the specific system and season, D I,0 . The shaded areas give the 95% confidence interval based on non-parametric bootstrapping (N = 1000, Bjørnstad & Falck 2001) restrial system, the distribution of distance between pairs of locations reached >1900 km, with a high number of pairs across the range of distances. For the open ocean and limnic system during winter, most locations were < 700 km from each other. Because environmental correlation is typically higher at short distances, the estimate of ρfor systems where we had data over larger areas, such as the terrestrial and limnic systems during summer, will be underesti-mated compared to systems with shorter distances between pairs (Hansen et al. 2020). However, accounting for such a bias by calculating ρand λ ρ -by using a similar maximum distance (e.g. 650 km) between pairs for all systems still showed the largest spatial scale and highest synchrony in the variation in temperature in the terrestrial system. The variation in temperature differed greatly between systems. In some terrestrial locations, the 15 Fig. 3. Spatial synchrony, ρ, in temperature in terrestrial, limnic, coastal sea and marine systems (Fig. 1) during winter (February−March) and summer (July−August). In the 3 aquatic systems, temperature was measured at 4 depths. The horizontal dashed lines show the overall mean correlation, ρ -, in the temporal dynamics of temperature across all locations in a system. The vertical dashed lines show λ ρ -, the distance at which the predicted relationship between ρ and distance is below ρ -. The light gray dotted lines show ρand λ ρ -based on pairs of locations that are maximum 650 km from each other Fig. 3 continued on next page range in winter temperature was >14°C during the study period, which will never occur in the aquatic system during any season, at least not at greater depths. However, organisms are adapted to temperature ranges occurring in their environments (Clarke 2017, Gvoždík 2018, which is more narrow in marine and limnic compared to most terrestrial systems. In the context of temporal variation in temperatures and how this is synchronised among locations, we believe that the standardised variation around the mean is most appropriate for understanding the spatiotemporal properties of environmental variation across sys-tems, because these values to a large extent determine the thermally viable envelopes for species . Spatial heterogeneity in environmental fluctuations is central for local and regional extinction probability both from classic meta-population theory (Hanski & Gilpin 1996) as well as more recent models of spatial ecology (Heino et al. 1997, Engen & Saether 2005. The basic concept is that if environmental conditions synchronise populations over large distances, the likelihood that all populations experience periods of critical low population sizes at the same time increases (Heino et al. 1997). In periods of small population sizes, the negative influence of demographic stochasticity on population growth increases, and so do the loss of genetic variation and dispersal rates. Such negative feedback may drive the species into a negative vortex ending with local and regional extinction (Amarasekare 1998, Tack et al. 2015. High environmental synchrony may cause such processes to occur over large areas simultaneously, which reduces the likelihood for a local rescue effect by immigration from neighbouring populations. Such processes can be scaled up to communities (Hansen et al. 2013, Koenig & Liebhold 2016, where spatiotemporal environmental variation affects the spatial patterns of fluctuations in species composition or richness (Mutshinda et al. 2009, Bellier et al. 2014. Species differ in their sensitivity to environmental conditions, and this has been used to explain variation in population synchrony among species with contrasting life histories (Marquez et al. 2019). Another striking consequence of life history variation, for instance along the fast−slow continuum (Oli 2004, Bjørkvoll et al. 2012, is that the sensitivity of population dynamics to different life history stages varies among species, and these life history stages often occur during different periods of the year. Our results suggest that seasonal variation in environmental synchrony can cause patterns of spatial scaling properties of species dynamics because the Moran effect is stronger when the environment is more synchronised (Moran 1953, Royama 1992, Hansen et al. 2020). In the terrestrial system, this occurs during winter, whereas, for instance, the limnic system has a higher synchrony in summer (Fig. 3, Table 1). However, few studies have assessed to what extent seasonal variation in environmental synchrony can explain speciesspecific patterns of population synchrony (but see Herfindal et al. 2020 for an example on synchrony in life history traits). Here, comparative studies between systems that differ in the seasonal patterns of synchrony (e.g. marine vs. terrestrial systems, Fig. 3) may shed light onto the relative importance of species demography and environmental variation in time and space for the spatial scaling of population dynamics.
We believe that our results give some indications of the spatial scale at which the conservation of populations and ecosystems may occur given the environmental condition in time and space affecting stochastic processes in populations and communities. This does not mean that other characteristics, such as dispersal capacity or nutrient flow, are not important (Mokany et al. 2010). However, stochastic processes must also be considered in practical conservation planning, even if modelling such processes in time and space on species or communities may be challenging. Assessing the patterns of the major driver of such processes, i.e. the environmental variation in time and space, may be one step to understand the relative importance of these processes across ecosystems. This is particularly true in the face of the ongoing global changes in climate, habitats or other important drivers. Although there are quite clear predictions about the temporal patterns of climate change, less is known about how climate change affects the spatial structure of the environment (but 17 Fig. 4. Correlation plots between system-specific measures of spatial patterns of autocorrelation. D I,0 is the distance at which the spatial autocorrelation in temperature becomes similar to the overall mean similarity based on Moran's I. ρis the overall synchrony in the temperature dynamics in the system, and λ ρ -is the distance at which the predicted relationship between distance and synchrony, ρ, decreases below ρ -. Colours represent systems: orange = terrestrial, purple = limnic, green = coastal sea, blue = open ocean. Symbols with black outline are winter, symbols with no outline are summer. Note that D I,0 and λ ρ -are ln-transformed. To allow comparisons across systems with different spatial extent, only pairs of locations closer than 650 km were used in the estimation of the 3 parameters see Koenig & Liebhold 2016). If the spatial structure of the environment changes, so must our conservation and management of nature. It is thus pivotal to acknowledge the spatial properties of environmental variation.