Spatiotemporal variability in Swedish lake ecosystems

Studying ecosystem dynamics is critical to monitoring and managing linked systems of humans and nature. Due to the growth of tools and techniques for collecting data, information on the condition of these systems is more widely available. While there are a variety of approaches for mining and assessing data, there is a need for methods to detect latent characteristics in ecosystems linked to temporal and spatial patterns of change. Resilience-based approaches have been effective at not only identifying environmental change but also providing warning in advance of critical transitions in social-ecological systems (SES). In this study, we examine the usefulness of one such method, Fisher Information (FI) for spatiotemporal analysis. FI is used to assess patterns in data and has been established as an effective tool for capturing complex system dynamics to include regimes and regime shifts. We employed FI to assess the biophysical condition of eighty-five Swedish lakes from 1996–2018. Results showed that FI captured spatiotemporal changes in the Swedish lakes and identified distinct spatial patterns above and below the Limes Norrlandicus, a hard ecotone boundary which separates northern and southern ecoregions in Sweden. Further, it revealed that spatial variance changed approaching this boundary. Our results demonstrate the utility of this resilience-based approach for spatiotemporal and spatial regimes analyses linked to monitoring and managing critical watersheds and waterbodies impacted by accelerating environmental change.


Introduction
Shifting ecosystem dynamics have broad implications for coupled human and natural systems. Because these systems are complex and integrated, no action exists in isolation and any disturbance may start a chain of events propagating across spatial and temporal scales [1]. Accordingly, the dynamic reality of these linked systems and potential impacts (e.g., increasing a1111111111 a1111111111 a1111111111 a1111111111 a1111111111 pollution levels, loss of ecosystem services, declining biodiversity), highlight the need to better understand patterns of change in social-ecological systems (SES). In the past, data availability was a significant limiting factor in these efforts, however, advancements in data collection tools and techniques have enhanced access to large-scale datasets (e.g., environmental, public health, economic). In addition, geospatial data and analytics have grown in their usefulness over the years, and data from remote sensing, drones, artificial intelligence, etc., have continued to be integrated into a variety of research areas. While access to such rich geotagged data streams provides significant opportunities for discovery and development, these data also present a challenge for researchers to capture useful information and emergent characteristics in the midst of a data deluge [2,3]. Further, the implications of known and emerging global issues (e.g., climate change, novel viruses) necessitates the need to continue developing methods that help monitor and manage SESs in the face of accelerating environmental change.
Resilience provides a means for studying environmental change and relates to the capacity of a system to withstand change and remain within a regime with the same processes and structures [4]. While the theory has a well-documented history of development and application, the term "resilience" has many different meanings, even in ecology. Allen et al. [5] identified three main definitions of resilience for SESs to include resilience as a: rate, process or an emergent property. Resilience as a rate, also known as recovery, engineering resilience or resiliency [6], focuses on a system's rate of return to equilibrium. Resilience as a process emphasizes building, maintaining, or enhancing system resilience of the current regime, e.g., community resilience, disaster resilience [5]. Resilience as an emergent property (e.g., socialecological resilience) explicitly accounts for the possibility that SESs can exist in different configurations or regimes [4,7]. Social-ecological resilience provides the foundation for a set of theories, tools and approaches for studying environmental change that examine the complex nature of ecosystems, including nonlinear change in system dynamics, tipping points and the possibility of alternate regimes in SESs [4,8]. Here, we employ social-ecological resilience as the overarching approach to explain the dynamics of SESs, as it encompasses recovery, stability, variability, adaptation, transformation and alternate regimes, and accounts for the interconnectedness of SES on multiple spatial and temporal scales, as well as cross-scale interactions [7].
SESs have been assessed using a variety of quantitative approaches (e.g., discontinuity analyses, early warning indicators, information theory, multivariate time series modeling) which have been effective in assessing a variety of model and real, human and natural systems [3,[9][10][11][12][13][14][15][16][17]. These approaches examine change in system dynamics which is a critical aspect of socialecological resilience [4,7]. Prior research has shown that these resilience-based approaches not only detect ecosystem change but may provide warning in advance of transitions (i.e., regime shifts), possibly allowing time for management interventions [3,11,12,14,17]. While these studies primarily focused on temporal dynamics, it has been recognized that quantitative resilience studies need to be expanded to assess change in spatial dynamics at the landscape scale [18]. Research on spatial boundaries has typically resulted in the development of ecoregion maps based on ecological potential, but this method does not consider the impact of external stressors (e.g., land use, climate) on underlying system conditions [17]. These stressors may have broad implications for human and natural systems through temporal and spatial scales. Hence, further research is need on the application of techniques that evaluate ecosystem dynamics in a spatially explicit way [3,17,19].
Researchers have employed spatial correlates of classic early warning indicators (e.g., variance, autocorrelation, skewness) with limited success [20][21][22][23][24][25][26] due to a number of concerns over inconsistent results and detection accuracy [27][28][29], although artificial intelligence approaches may help overcome some of these problems [30]. In addition, many of these methods require temporal data to document and verify changes in results from spatial approaches [31][32][33]. Newer methods like recovery length and network-based indicators lack utility for assessing different system types and heterogeneous patterns, and require large, highresolution datasets, highlighting the need for multiple methods to detect change in system dynamics [3,[34][35][36]. Further, many of these studies tended to focus on examining single or a few selected system variables. While this may be appropriate for models or systems that are well understood, the dynamics of real ecosystems are complex and linked, and often drivers (e.g., fast and slow variables such as flash floods or climate change) are unknown or unmeasured. Accordingly, there is a great need for methods that can track trends in multiple variables, thereby reflecting the inherent complexity in ecosystems [11,37,38].
A variety of multivariate methods have been used for spatial analysis including spatial variance, sequential t-test of regime shifts (STARS) and constrained hierarchical clustering to monitor spatial regimes along a rangeland transect [26,39]. For example, Angeler et al. [40] employed multivariate spatial modelling (canonical ordination) to examine spatial pattens and functional redundancies in invertebrate community structure collected across a Swedish boreal lake landscape. Each of these approaches provides useful insights on spatial patterns, regimes and resilience that may facilitate monitoring and management. However, these methods typically require the use of multiple techniques to analyze the data (e.g., principal components analysis, redundancy analysis, variation partitioning) [26,40].
Fisher information (FI) is a statistical measure of dynamic order in systems [41] and has been used to detect patterns in complex system behavior by examining trends in variables that characterize the state (condition) of the system [42]. The method involves simultaneously assessing the dynamics of multiple, underlying system variables by collapsing them into an index that can be used to track changes in system condition [3, 9-12, 14, 17, 43, 44]. FI has been established as a useful measure of resilience as it assesses change in complex system behavior to include the identification of ecological stability and variability, regimes and regime shifts [3,11,12,14,17,43]. The method does not need additional analytical techniques, has no strict data requirements, does not require prior knowledge of system dynamics (3) and can handle multiple variables and a variety of data types [3,[10][11][12]17]. This information theory approach has been used to examine patterns and trends in SESs [9-12, 14, 17, 44-48]; however, its application to spatial data has been limited [3,17]. FI has previously been used to successfully detect spatial regimes in bird and zooplankton communities along transects at sub-continental and oceanic scales [17]. To transform the approach into a general method for spatial analysis, Eason et al. [3] enhanced the method to detect distinctive spatial patterns over a landscape, however the application was largely limited to simulated data at discrete time periods.
In this study, our aim is to further test the utility of FI by using it to assess spatiotemporal patterns in the biophysical condition of multiple Swedish lakes that are embedded in and form a network in a terrestrial matrix (i.e., "lake landscape") [49]. Using readily available data supplied by the Swedish University of Agricultural Sciences, we used FI to explore spatial and temporal patterns in the physicochemical environment and phytoplankton community structure data gathered in eighty-five Swedish Lakes from 1996-2018. This study assesses whether FI can detect changes in the biophysical environment of the Swedish lakes through time and explores the implications for social-ecological resilience. Further, we examine whether there is a relationship between spatial regimes and increased variance as found in other ecosystems that are spatially explicit [17,26] and non-explicit (summarized in Angeler et al. [50]). Specifically, we consider and analyze a northern and a southern ecoregion as two independent spatial regimes, which are separated by a "hard boundary", the Limes Norrlandicus (information below). A recent study, based on canonical ordination, has shown that despite environmental factors steering ecological community dynamics distinctly in these two spatial regimes, the community dynamics per se were stable in a time span of approximately two decades [51]. Our study aims at complementing these previous results not only by assessing ecological patterns of system variability within both regimes, but also more specifically whether ecological dynamics change with approximation to regime boundaries, a characteristic inherent in spatial resilience and regime theory [17,18].

Study area and data
In the late 1980s, Sweden initiated a long-term monitoring program of multiple habitats and trophic levels of lakes. Eighty-five lakes (Fig 1) with the longest continuous time series were selected from the monitoring data for this study. The lakes were situated in two ecoregions in Sweden, separated by the Limes Norrlandicus (LN), a strong biogeographical and climatic divide between northern and southern Sweden. These ecoregions differ in terms of vegetation (e.g., boreal/alpine in the north vs. hemiboreal in the south), precipitation (duration of snow cover), soil type and air temperature. Coinciding with different biotic structures, the pressures affecting the ecoregions above and below the LN are different due to higher human population density, greater agricultural land use, higher storm intensity, and historically more acidification in the southern ecoregion [51].
Despite different biophysical and anthropogenic factors characterizing the northern and southern regime at the ecoregion level (i.e., macroscale), all Swedish lakes were situated in least disturbed catchments (mesoscale). The lakes span gradients in abiotic characteristics (Secchi depth 1.1 to 10.7 m; water color <0.01 to 7.98 ±mg Pt L-1; pH 4.6 to 7.3, alkalinity <0.01 to 0.30 meq L-1; and lake size 0.11 to 5.43 km 2 ). For this study, physicochemical and phytoplankton community assemblage data were compiled from data collected by the Swedish University of Agricultural Sciences between 1996-2018 (S1 and S2 Data) and used to characterize the condition of the lakes. Physicochemical variables include watercolor (absorption), water temperature (˚C), pH, alkalinity (mg/L), SO 4 2 -concentration, electrical conductivity (mSiemens/cm) and nutrients (NH 4 -N concentration (μg/L), PO 4 -P concentration (μg/L), total organic carbon (mg/L) and total phosphorus (μg/L). Data from surface water samples were selected from the periods when phytoplankton biovolume was sampled, usually in August, at a mid-lake station in each lake. Water was collected with a plexiglass sampler and kept cool during transport to the laboratory. Phytoplankton was sampled by taking a water sample from the epilimnion using a 2-m-long plexiglass tube sampler and preserving it with Lugol's iodine-solution. Phytoplankton counts were made using an inverted light microscope and the modified Utermöhl technique commonly used in Nordic countries [52]. Specifically, taxa were identified to the finest taxonomic unit possible (usually species), and species-specific biovolume measures were calculated based on three-dimensional geometric shapes following methods employed by Blomqvist et al. [53]. All physical and chemical analyses were conducted at the Department of Aquatic Sciences and Assessment, Swedish University of Agricultural Sciences following international (ISO) or European (EN) standards, when available [54]. Additional information is available at the website of the Department of Aquatic Sciences and Assessment, Swedish University of Agricultural Sciences (http://www.slu.se/vatten-miljo) and Fölster et al. [54]. The data for this study was managed, analyzed and visualized using Microsoft Excel 365, Matlab (R2021) and JMP 14.

Fisher information (FI)
The form of FI used for this study (Eq 1) is based on the probability of observing a particular state (condition) of the system (p(s)) [55,56]. Given that FI is proportional to the change in probability of a particular state (dp(s)) over the change in state (ds), the higher the probability of a particular state, the more orderly (stable) the system (i.e., consistent patterns) and the higher the FI.

PLOS ONE
Conversely, when a system is more disordered (highly variable) and exhibits no distinctive patterns or bias toward a particular state, FI is low and tends toward zero [57]. System states are based on the condition of the system as characterized by a set of measurable variables (s: x 1, x 2. . .. x n ) that may be assessed over time and/or location. Because of internal dynamics or external stressors, system condition may vary such that it experiences small perturbations or could undergo substantial change. To distinguish minor fluctuations from significant variation, a size of states (sost) is defined for each variable such that if the measured value varies within a finite range (e.g., |x 1 (t 1 )-x 1 (t 2 )| � sost x1 ), the two points are indistinguishable and counted as observations of the same state. The p(s) is calculated by counting the number of points that fit within a particular state [42]. Using a series of derivation steps, Eq 1 was adapted to handle empirical data in line with this grouping (binning) approach, such that FI is calculated as shown in Eq 2: where q(s) � p p(s) [42]. In accordance with our discussion on resilience and alternate regimes, the basis of interpreting FI is understanding that different system regimes are characterized by distinct processes and structures [4,58]. Accordingly, tracking trends in FI affords the ability to detect change in system condition [3]. While rising FI signifies increasing dynamic order and greater stability, declining FI indicates loss of order and increasing variance which is indicative of degraded resilience in a particular regime [11,12,42,44]. To compare different systems, time periods or geographical regions, basic statistical measures like the mean (μ), standard deviation (σ) and coefficient of variation (cv = s m ) of FI are used to help classify variation in the condition to include identifying regimes and regime shifts [3,10,17,59,60]. System regimes are periods or areas in which relatively stable patterns are observed. Stable regimes exhibit more consistent patterns (greater dynamic order) and display less variability; hence, they are defined by relatively high mean FI ("μFI), low standard deviation in FI (#σFI) and low coefficient of variation of FI (#cvFI). Transitional dynamics are periods where the patterns and conditions in a particular regime begin to degrade and may warn of an impending regime shift. With instability serving as the hallmark, these regime shift precursors are identified as periods with "σFI, "cvFI [3,17]. More details on the derivation, calculation and interpretation of FI may be found in [3,17,42,56,61,62].

FI for spatial analysis
Using FI to assess spatial dynamics begins with compiling measurable data from each survey station. The data should include station identifiers (unique ID), location (latitude and longitude) and measured values for variables that characterize system condition (e.g., physicochemical and plankton data) over the study period: While time is a natural ordering principle for temporal studies, spatial analysis requires using a space for time substitution [17]; accordingly, Eason et al. [3] tested multiple approaches and selected Haversine distance for ordering the stations over the landscape as it most efficiently accounts for the curvature of the Earth. To order the stations, the Haversine distance is computed as [3,63]: ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffiffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi where, λ = longitude, φ = latitude, r = radius of the earth (r~6.371 km) [63].
The distance is computed from a reference location (minimum latitude and longitude of the dataset). After which, the stations are sorted by distance and the data collected at the sampling stations are used to compute FI in a moving window. Each window encompasses small geographical areas over the lake landscape, such that the first window will contain the data from stations that are closest to the reference site and the second window will advance forward to the next closest station. Then, in line with the original (temporal) FI algorithm [3,12,17]: • In each window, group (bin) points into states of the system using the level of uncertainty (sost) around the variables. If the measurement uncertainty for each variable is unknown, then uncertainty data may be gathered from a similar system or by using a relatively "stable" period within the dataset as a proxy. This stable period may be found by separating the data into moving windows, computing a standard deviation for each variable and using the mean standard deviation of the window with the least variation to set the sost [61].
• Count the number of points grouped into each state and divide this value by the total number of points in the window to produce p(s).
• Compute q(s) = p p(s) and calculate FI using Eq 2.
The steps above are repeated for each moving window producing a FI value that is plotted at the corresponding sampling stations (latitude and longitude) over the geospatial area. FI code is readily accessible at https://github.com/csunlab/fisher-information [2].

Comparative analysis
To explore possible linkages between phytoplankton community structure and changes in the physicochemical environment, we computed FI for both the phytoplankton (PHYTO FI) and physicochemical data (PC FI) over the lake landscape and compared the patterns and trends in the FI results, spatially and temporally. GIS plots and Kendall's Tau-b correlations provided a means of visualizing the results and assessing whether the phytoplankton community dynamics are correlated with changes in the physicochemical environment. That is, we examined correlations between PHYTO FI and PC FI over the lake landscape to assess the congruence in spatiotemporal patterns between abiotic and biotic factors and explored whether distinct patterns tended to emerge above and below the LN. In addition, we assessed whether the lake ecosystem dynamics changed monotonically over time by correlating FI with the factor "year" (from 1996 to 2018) to infer whether the spatial variability of the lake landscape changed for the period of the study. Further, to help streamline the assessment, summary statistics (mean FI: μFI, standard deviation in FI: σFI and coefficient of variation in FI: cvFI) were used to examine changing condition over the landscape and through time. When identifying "stable" regimes during spatiotemporal analysis, the metrics must meet all three criteria = > "μFI: μFI year > μFI 1996-2018 , #σFI: σFI year < mean σFI 1996-2018 and #cvFI: cvFI year > mean cvFI 1996-

Results
The Swedish University of Agricultural Sciences dataset consists of physicochemical (PC) and phytoplankton (PHYTO) community data on 623 phytoplankton species collected at the 85 lakes between 1996-2018. These data characterize the biophysical environment of the lake landscape and were used to compute the physicochemical FI (PC FI) and phytoplankton FI (PHYTO FI). The PC and PHYTO datasets were sorted by Haversine distance in each year and stored in separate workbooks (S1 and S2 Data).
Using data from the first year of the study (1996) and a window size of 10 (i.e., 10 lakes), we calculated the PC FI and PHYTO FI which were plotted at the lake locations (latitude and longitude) corresponding to each window. Fig 2 provides GIS plots of the results using color and marker size gradients to differentiate low (blue) and high (red) FI values and help visualize the changing biophysical condition over the lake landscape. While FI values tended be higher above and decline approaching the LN for both datasets, the LN served as a much clearer demarcation for lower PC FI values. In addition, we found that the PC FI and PHYTO FI values were positively correlated (τ: 0.32, p-value <0.05) with the PC FI values which displayed more variability (higher σFI (PC: 1.01 > PHYTO: 0.86) and higher cvFI (PC: 0.24 > PHYTO: 0.17)). GIS plots of the PC FI and PHYTO FI from 1996-2018 are provided in the (S1-S6 Figs). To examine how these spatial patterns emerged and changed over time, we used Kendall's Tau correlations and summary statistics (μFI, σFI and cvFI) to streamline the comparisons and facilitate the assessment.

Assessing spatiotemporal change
Scatterplots in conjunction with Kendall's Tau correlations showed that the PC FI and PHYTO FI values at least moderately corresponded over time (Fig 3). The strength of the association ranged from weak (τ 2018 = 0.15, p-value<0.05) to strong (τ 2010 = 0.48, p-value<0.05) and declined during the study period (Fig 4). Figs 5 and 6 provide the summary statistics which help capture changes in spatial variation for PC and PHTYO dynamics over time. The most stable years are found in the lower right-hand side of the μFI vs σFI plot ("μFI, #σFI; Figs 5A and 6A) and in Figs 5B and 6B, the highlighted values are those meeting the "stable"    (Figs 5B and 6B). Fig 7. provides a plot of the cvFI for PC FI and PHYTO FI over time. Since cvFI is a measure of dispersion around the mean, the cvFI alone is a useful aggregate for capturing the spatial variation over the lake landscape, such that higher cvFI values indicate greater variation. Results show statistically significant declines in cvFI for both PC (τ PC : -0.39, p-value < 0.05) and PHYTO FI (τ PHYTO = -0.36, p-value < 0.05) and lower cvFI values for PHYTO over the entire period indicating greater stability.
To examine the difference in spatial variability between Sweden's northern and southern ecoregions, we partitioned the data and computed the cvFI on the lakes residing above (N = 31) and below (N = 54) the LN. Results showed that the cvFI values were higher below the LN indicating greater variability (Fig 8). Further, we found evidence of declining trends in cvFI for PC and PHYTO above the LN and PHYTO below the LN (t � -0.30, p-value <0.05).

Discussion
Lake ecosystems vary temporally, spatially and across scales. Although they are typically able to maintain processes, structures and functions, accelerating environmental change can erode lake resilience, reducing their capacity to withstand disturbances, and resulting in undesirable consequences (e.g., regime shifts to eutrophic lakes). Social-ecological resilience presents a means of understanding SES dynamics and accounts for the fact that natural systems rarely change monotonically, instead they are often characterized by ecological thresholds, nonlinear change and alternate regimes [4]. Resilience-based approaches have proven to be powerful tools for examining coupled human and natural systems. One such method is FI, however, there has been limited application of this approach to assessing spatial dynamics [3, 9-12, 14, 17, 44, 47]. This work serves as the first time that FI was used to study spatiotemporal dynamics at the landscape scale. Using the approach, we were able to assess variation in the biophysical conditions of the Swedish lakes landscape, which was useful for gaining better insight into the dynamics of these lakes from the perspective of spatial resilience and spatial regimes.
To evaluate patterns of stability and variation within the Swedish lakes landscape from the viewpoint of spatial regimes [18], this study involved compiling data on the physicochemical environment and phytoplankton community structure of 85 Swedish Lakes from 1996-2018

PLOS ONE
and determining whether FI was able to detect spatiotemporal change in biophysical conditions across the Swedish lake landscape. Results showed that FI values tended to be higher with less variation at sampling stations in the northern spatial regime relative to the southern spatial regime which suggests that the condition of lakes in the northern spatial regime is more stable both spatially and temporally (Figs 2, 8, S1-S6 Figs). This finding is in agreement with a previous study based on benthic macroinvertebrates, which showed that differences in stability were influenced by different environmental variables [51]. The PC and PHYTO data displayed distinct levels of variation with more spatial stability found in phytoplankton community structure, suggesting greater spatial resilience. However, there was more congruence and less variability in biophysical condition near the end of the study period (Figs 5-7). This observed congruence between results of phytoplankton and physicochemical conditions may be useful for management. There is a rich body of ecological literature which explored the congruence of patterns between taxonomic groups as indicators of environmental change and ecological status (e.g. [64]), which suggested that weak congruence between results from different ecological communities requires multiple lines of evidence to help strengthen the inference [65][66][67]. Given that monitoring lakes to the extent done in the Swedish lakes monitoring program is extremely costly, our results suggest that measuring either phytoplankton or the physicochemical environment may lead to similar results about spatial and temporal variability in the lake landscape; thereby, substantially reducing monitoring costs. However, we emphasize that such substitutions may not be generalizable and need to be considered with caution and more research.
Our findings showed that while the biophysical conditions of the lakes generally increased in stability over time. Physicochemical conditions and phytoplankton community structure varied from year to year and experienced a lengthy period of instability before both appeared to find more consistent patterns in 2016 (Figs 5-7). Such a monotonic trend fits the pattern observed in previous research which found that spatial phytoplankton community turnover (i.e., the spatial differentiation of community structure) increased while nestedness (a surrogate for species extinctions) decreased linearly over time [68]. Community stability despite higher community turnover has also been observed in other studies (e.g., [69]). Other research focusing on the abiotic environment has shown linear changes in decreasing total phosphorus concentrations across lakes due to climate change and recovery from acidification [70]. Whether the increasing phytoplankton spatiotemporal resilience observed in this study is linked to changes in regional nutrient patterns or other factors remains to be explored.
Focusing in on the spatial regimes [17] aspect of our study, our analysis of the spatial dynamics above and below the Limes Norrlandicus (LN) identified patterns of decreased spatial variation (greater stability) in biophysical conditions (both PC and PHYTO) in the vicinity above the LN and below the LN in PHYTO dynamics (Fig 8). The analysis revealed a spatial transition zone ("σFI and "cvFI) which coincides with the ecotone boundary that separates the southern and northern biogeographical areas of Sweden (Fig 2; S1-S6 Figs). The results of our analysis are in agreement with changes in community structure previously observed in lake invertebrate communities [71], and delineating the northern and southern biogeographical areas as spatial regimes [17], based on their biophysical and climatic status (see Methods). These results are also supported by analyses finding differences in the correlations between ecological stability patterns and environmental variables [51]. Lower variation in dynamic order as captured by cvFI seems to indicate spatial resilience as a function of transitions between spatial regimes in the Swedish lakes landscape. These changes differ between ecoregions that are separated by a clearly delineated stable boundary in the Swedish landscape [51], the Limes Norrlandicus [72]. Therefore, from a spatial regimes perspective, our results provide empirical support of non-random patterns of variability at the edge of regimes which has been discussed in ecological resilience [73]. Researchers have found that variables like animal fitness, nomadism and species invasions are variable at ecological transition zones [50]. Our study is, to the best of our knowledge, among the first assessing such non-randomness in an explicit spatial context by providing evidence for the assumption that higher variability manifests at spatial transition zones, and that spatial resilience might decrease towards spatial regime boundaries [18]. While our study assessed spatial variability in relation to a stable biogeographical frontier, a recent study has shown how spatial regime boundaries can be assessed when they change dynamically as a function of land-scape level changes in vegetation cover and animal communities [39].
Results also revealed a decline in the correlation between PC and PHTYO dynamics indicating an increased decoupling over time. Phytoplankton biodiversity and community structure is largely shaped by water chemistry (particularly nutrients) but can also be influenced by biotic interactions and spatial effects at the inter-basin scale and the resulting spatially structured abiotic lake environment (e.g. [67]). Therefore, there is high complexity mediating community dynamics that is difficult to scrutinize; accordingly, at this stage, we can only speculate about likely effects of biotic interactions influencing phytoplankton in a way to potentially explain the increasing decoupling over time. It is also speculative to relate an increasing occurrence of blooms and the spatial expansion of the raphidophycean flagellate, Gonyostomum semen, to heightened variability in community dynamics [74] affecting the observed correlation patterns. Our resilience-based study builds upon the evidence provided by more traditional analyses used in ecology. Specifically, the successful application of FI is supported by a high congruence between the identified FI patterns and the spatial structure revealed in traditional [71] and resilience-based analyses of littoral macroinvertebrate assemblages [40].
In line with our results, researchers have found evidence that despite the long-term changes in the abiotic and biotic environment observed, Swedish lakes have been quite stable thus far [75], meaning that there is no evidence that the entire lake landscape is at eminent risk of a broad-scale, regional regime shift. Investigations about the contribution of individual lakes to the stability of the lake landscape have started recently and are thus preliminary [51]. Additional studies and other methods can complement our work by helping to determine the specific conditions that changed in the SES, and underlying causes. As more spatial data accumulate, future work could aim to better link resilience-based approaches to real world applications of spatial resilience and spatial regimes detection by further examining index results to uncover trends in underlying ecological variables and identifying key drivers useful for monitoring and management. Such efforts are warranted given that accelerating environmental change substantially alters the configuration of and shifts between spatial regimes [76,77].
We conclude by highlighting that the monitoring data only provided us with the opportunity to analyze data from a single month (usually August). While providing important insight into patterns of spatial variability, we acknowledge that phytoplankton community structure can undergo substantial seasonal change. Our results are therefore only relevant for the summer and should not be generalized for other seasons. Notwithstanding, our analyses support the value and power of a resilience-based approach (FI) for assessing spatiotemporal change in SES. Further development and application of such approaches would provide more focused and relevant information to support systems-based management and policy decisions and facilitate efforts to manage social-ecological resilience in the face of rapidly accelerating environmental change.