Spatial regime shifts over time and space
To assess how the dominance of large predatory fish (Eurasian perch Perca fluviatilis, northern pike Esox lucius) and three-spined stickleback (Gasterosteus aculeatus) changed over time and space, we collated juvenile fish abundance data from 13073 samplings conducted during 39 years (1979–2017) in 486 bays along a 1200 km stretch of the Swedish Baltic Sea coast. We used juvenile fish surveys because they (i) include stickleback, (ii) were conducted along the entire coast, including the outer archipelago, and (iii) capture patterns of recruitment failure; a proposed driver of local perch and pike stock decline34,35,37,45. The samplings were conducted by various monitoring programs and research projects to quantify fish recruitment. Much of the data was extracted from the Swedish national database for coastal fish (http://www.slu.se/kul). Other fish species also occurred in the data, but we—like others25,34,45—focused on the most common and strongly interacting species.
The timing and placement of most of the fish surveys was not chosen with this study in mind. We therefore included data from as many surveys as possible, ensuring that the dataset covered (i) gradients in distance to open sea and wave exposure, (ii) the entire Swedish east coast, and (iii) the longest time period possible (Supplementary Fig. 7). Nearly all sampling (94%) was conducted during July-September. To achieve the best possible spatial coverage we also included some bays only sampled in October and (for one bay) June. Initial exploration of data from bays sampled monthly from June to October suggested there were no large differences in the October and June data. Nine of the 486 bays occurred much further into the archipelago (49–67 km) than the rest (<41 km), resulting in very poor spatial coverage of the innermost half of the archipelago gradient. Moreover, these nine bays were all sampled after 2011, resulting in that space and time were confounded. Consequently, we excluded the nine bays (which were all predator-dominated) from further analyses.
The surveys were conducted in shallow coastal bays at 1.7 ± 0.6 depth (mean ± SD, range: 0.35–4.5 m, 99.9% ≤3.5 m, N = 606). Juvenile fish were sampled using low-impact pressure waves; a standard method in the area59. In short, the ignition of a small, underwater explosive charge generates a pressure wave that stuns or kills all small fish with a swim bladder within the blast radius. Using the current Swedish standard (10 g Pentex explosive ignited by a 1 g non-electric charge), 2–20 cm fish are sampled within approximately a 5 m radius (ca. 80 m2). The fish are then collected using swing nets and snorkeling, identified and counted. All sampling was conducted by certified personnel and with required ethical permits.
The type and amount of explosives varied between surveys (1–25 g), and some surveys noted only floating (not sunken) fish. To account for these differences, we recalculated the total number of individual fish per species to the current detonation standard using experimentally derived conversion factors59. For surveys only reporting floating fish, we also calculated the expected total number of fish using the sunken:floating fish ratio from other surveys (N = 47–152 per species).
Since the method samples highly mobile organisms within a small surface area and very short time span, it is prone to high variability and false absences. We therefore averaged the species abundance of perch, pike and stickleback per bay and year (based on 5–177 samplings [detonations] per bay-year: mean ± SD = 15.6 ± 12.2). Most (75%) of the 477 included bays were sampled once, but 25% during multiple (2–34) years, generating 833 unique bay-year combinations. Since fish assemblages within a single bay can shift from perch and pike to stickleback dominance (for example, see Fig. 3), the 833 bay-years were treated as individual replicates. For each bay-year we then calculated the relative predator dominance (the summed abundance of perch and pike divided by the summed abundance of perch, pike and stickleback, ranging from 1 to 0 where 1 = complete perch and pike dominance, and 0 = complete stickleback dominance). We primarily used relative abundance because it captures the community state well60,61,62 and reduces the “noise” in absolute abundance caused by the notoriously high year-to-year variability in fish recruitment63. Perch constituted 96% of the pooled abundance of perch and pike, but we included pike as well because they were occasionally more common than perch.
We also included data on four covariates known to influence juvenile fish abundance36, when they were sampled: water depth (nearest 0.1 m, N = 606 bay-years), water surface temperature (nearest 0.1 °C, N = 598), salinity (nearest 0.1 psu, N = 300) and water visibility (Secchi depth to the nearest 0.1 m, estimated from turbidity, N = 325). Finally, we for each bay used GIS to calculate the distance to open sea (m shortest water distance from the baseline) and the relative wave exposure (m2 s−1), based on a surface wave model (SWM) utilizing fetch and wind data averaged across 16 compass directions. The model mimics diffraction using empirically derived algorithms (for details, see ref. 64).
Temporal regime shift at Forsmark
To test for a more classic temporal regime shift at the local (non-spatial) scale, we used the longest time-series from a single area; a 34-year data set from an annual monitoring program (1981–2017, no sampling in 1982, 1989, and 2003) conducted outside Forsmark (60.4°N, 18.2°E). The data was collected in an area used as a reference for estimations of effects of the release of heated cooling water from the Forsmark nuclear power plant. Only perch and stickleback (no pike) occurred.
Importance of predator–prey reversal for fish recruitment
To assess the relative importance of predation and predator–prey reversal for juvenile fish recruitment, we combined an ecosystem field survey with causal path analysis; a powerful approach to tease apart the role of direct and indirect interactions in ecosystems36,49. In 2014, we surveyed 32 shallow bays (<3.5 m depth) situated >10 km apart (or separated by natural fish movement barriers like deep water or land) along a 360 km stretch of the central Swedish Baltic Sea coast (Supplementary Fig. 5). Together, the bays formed gradients in distance to open sea and wave exposure, but also in perch, pike and stickleback abundance. This enabled us to statistically separate the influence of biotic interactions (e.g. predation) from abiotic conditions (e.g. wave exposure). In spring (May), after adult fish had migrated into the bays to spawn, we sampled the abundance of adult fish using 3–5 standard 30 m Nordic survey gill nets set overnight36. In 6–8 stations per bay we also quantified two environmental covariates known to influence fish: total % bottom cover of habitat-forming aquatic vegetation (visually estimated by snorkeling, separating rooted from non-rooted macrophytes; for details, see36) and density of zooplankton (juvenile fish food37, estimated as the total plankton density per liter water). Zooplankton were sampled using a 25 cm diameter Epstein net (80um mesh), slowly pulled three times vertically from 0.7 m above the seabed to the water surface. Zooplankton were fixated in 5% formalin and then counted in the lab. In late summer (August) we at 6–8 stations per bay estimated fish recruitment (density of young-of-the-year fish, using low-impact pressure waves; see above) and % cover of habitat-forming benthic vegetation. The fish sampling was evaluated and approved by the ethical board on animal experiments of the County court of Uppsala, Sweden, permit C 139/13.
Statistics and reproducibility
Spatial regime shifts over time and space: we hypothesized that a shift from perch and pike to stickleback dominance started in wave-exposed bays in the outer archipelago, but then propagated towards the inner archipelago over time (Fig. 1a)—i.e. a statistical effect of space (distance to open sea), time (year) and, possibly, their interaction. Regime shifts have often been identified as breakpoints in time-series using non-parametric change-point detection methods. However, spatial regime shifts in heterogeneous systems should in theory be gradual (the mean of many small-scale shifts) at the whole-system scale12,13. Moreover, while multiple drivers could in theory influence whether bay assemblages are perch and pike- or stickleback-dominated, change-point detection methods cannot handle multiple predictors and interactions between them. Finally, our response variable (relative predator dominance) is bounded between 0 and 1 and was bimodal, based on Hartigan’s dip test of unimodality65 with 10,000 permutations (Fig. 1b). Consequently, we used a binomial generalized linear model with a logit link function66 to explore the effects of time (year), distance to open sea and their interaction on relatively predator dominance, while controlling for the influence of wave exposure (log-transformed) and latitude. After assessing assumptions of homoscedasticity by plotting deviance residuals vs. observed values for each predictor (Supplementary Fig. 8), and multicollinearity using the variance inflation factor (VIF; most <2, all <5), we identified the most parsimonious model by (i) comparing candidate models using Akaike’s Information Criterion (AIC) and (ii) stepwise removal of non-significant terms (at α = 0.05). All statistical analyses were conducted using R v. 3.6.067.
Since 75% of the bays were sampled only once we could not include “bay” as a random factor. To assess the robustness of our results given the repeated sampling in 25% of the bays, we (i) generated a smaller dataset including one randomly chosen observation (year) per bay (N = 477), (ii) refitted the best binomial model using this new dataset, (iii) extracted the estimates, standard deviations and standard errors for each parameter, (iv) repeated the whole procedure 500 times, and iv) summarized the average results (Supplementary Table 4). The additive effects of time, distance to open sea, wave exposure and latitude remained, meaning that the conclusion of a spatial regime shift were robust. The interaction between time and distance was only significant in 29% of the runs; a discrepancy most likely explained by the lower power of a test based on a smaller sample size (N = 477 vs. 833 in the full dataset), as well as the influence of randomly including more perch- and pike- or stickleback-dominated years from the 25% (131) bays sampled >1 year.
To test whether the maximum spatial extent of stickleback dominance along the archipelago gradient increased over time, we first selected all stickleback-dominated bay-years (relative predator dominance ≤0.1, i.e. ≥90% stickleback). For each sampling year we then extracted the maximum distance from open sea, separating wave-sheltered vs. -exposed bays (using a cutoff of log10(m2 s−1) >4). Since no bays were sampled >20 km from the open sea prior to 1995, we only used data from 1995–2017. Finally, we used a general linear multiple regression model to test how sampling year, wave exposure (two levels) and their interaction explained the maximum distance.
The bimodality in relative predator dominance seen over time and space (Figs. 1–3) could in theory be caused by variability in local abiotic conditions, and not by predator–prey interactions. We therefore tested whether any of four abiotic conditions estimated locally—water surface temperature, salinity, turbidity, water depth—could explain the variability in deviance residuals from the binomial glm (Supplementary Fig. 8), using a regular linear model (after assessing model assumptions and ruling out multicollinearity; see above). Salinity and turbidity (but not depth and temperature) had statistically significant but weak influences (R2 = 0.11), and the bimodality clearly remained (Supplementary Fig. 9).
Finally, we explored how time, distance to open sea, wave exposure, latitude and their two- and three-way interactions influenced log-transformed abundances of (i) predators (perch and pike pooled) and (ii) stickleback, using general linear models. We identified the most parsimonious models as outlined above.
Temporal regime shift at Forsmark: using the 34-year Forsmark time series, we first tested for temporal breakpoints in logit-transformed relative predator dominance data using change-point detection (strucchange) for linear models47. This method estimates the optimal number and (if identified) position of breakpoints using the Bayesian Information Criterion (BIC). Second, we explored what temporal change(s) in perch and stickleback abundances that preceded the shift. Because of highly non-linear patterns we modeled the temporal changes using generalized additive models (GAM) as implemented in the mgvc package68.
Importance of predator–prey reversal for fish recruitment: to assess the relative importance of predator–prey reversal for perch, pike and stickleback recruitment, we used statistical model selection based on path analyses; a form of structural equation modeling that can be used to tease apart direct vs. indirect (mediated) relationships between multiple (>2) variables, and thereby assess the relative importance of direct vs. indirect relationships in systems48,69. Initial data exploration using multiple regression showed that one bay was a clear outlier due to 0 juvenile perch and pike, generating (i) too high leverage (influence on statistical relationships), (ii) heteroscedasticity and (iii) non-normally distributed errors. Since we suspected that juveniles had already migrated out of this bay, the bay was excluded (resulting in N = 31). Removing this statistical outlier resulted in that the model fulfilled test assumptions and the overall fit more than doubled (adjusted R2 increased from 0.17 to 0.37).
Based on ecological knowledge of the study system, we expressed 14 multivariate hypotheses of the direct and indirect drivers of perch and stickleback recruitment, as graphical network models of interacting paths36,49. Due to the relatively low sample size we restricted the number of paths to 7. The two simplest models assumed that perch+pike and stickleback juvenile abundance in summer (i.e. recruitment) was influenced by adult abundance in spring (stock recruitment) and cumulative cover of rooted vegetation, while adult abundance was explained by spring cumulative cover of all vegetation species36,46, and distance to open sea or wave exposure46. The more complex models included combinations of known predator–prey interactions: perch and pike controlling adult stickleback in spring through predation36, stickleback feeding on juvenile perch and pike45, and stickleback competing with juvenile perch for zooplankton prey37. We then analyzed each model using piecewise path analysis as implemented in the piecewiseSEM package48. First, we tested the goodness of fit of each model to the data using Shipley’s test of directional separation (D-sep)70. If missing paths were identified, they were included in a new model. For the models that fitted the data (p > 0.05) we used the Akaike’s Information Criterion corrected for small samples (AICc) (calculated using Shipley’s general approach to calculate AIC for path analysis71) to compare relative model fit. A summary of all candidate models and details of the best-fitting model are presented in Supplementary Tables 1 and 2, respectively. The strength of paths in the best-fitting model are presented using standardized path coefficients, which (based on the best-fitting model, see Fig. 3) mean that 1 SD increase in pooled adult perch and pike abundance reduces adult stickleback abundance by 0.53 SD. We also calculated the amount of variation (R2) in adult stickleback abundance, pooled juvenile perch and pike, and juvenile stickleback abundance, that was explained by the paths. Finally, we tested whether the relative predator dominance of adult and juvenile fish in this smaller dataset (N = 31) was also uni- or bimodal, using Hartigan’s dip test (for details, see above).
Reporting summary
Further information on research design is available in the Nature Research Reporting Summary linked to this article.
Source: Ecology - nature.com
