The study was conducted in the north-eastern Baltic Sea ( Figure 1 ). The sampling locations extended throughout the Estonian coastal sea with stations in the Gulf of Finland (eastern, central, and western), Gulf of Riga, Väinameri and Baltic Proper. The area is characterized by almost non-existent tides, high nutrient fluxes from large river systems and seawater freshening (Segerstråle, 1957; HELCOM, 2013). The coastline is highly heterogenic, with rocky shores at the tips of the peninsulas, mixed shores with cobbles, gravel-pebble, and sand particles on the flanks of the peninsulas and finer grained sandy and silty shores in the bays (Martin et al., 2013). Benthic communities in the area are affected mostly by depth, salinity, wave exposure, and seabed substrates. Species richness is low and spatially variable (Ojaveer et al., 2010).

The yearly catch of round goby was used as a proxy of the areal biomass of the invasive fish in each station. The catch statistics were obtained from the Veterinary and Food Board based on obligatory catch data submitted by coastal fishermen. Each round goby catch in a unique small statistical quadrat is reported as a monthly biomass (kg). The main grid of catch quadrates consists of approximately 10 x 18 km grid cells and some neighbouring grid cells near the coastline are merged. Small statistical quadrats map can be downloaded from the Estonian Government website (https://www.riigiteataja.ee/aktilisa/1311/0201/5002/VV_29102015_110m_lisa2.pdf). The annual round goby biomass was calculated as the total annual biomass (kg) divided by the catch quadrate surface area (km²) obtained with ArcMap 10.8.1. For each benthic sampling station, each sampling year-round goby biomass (kg km^-2) was assigned based on its location within a specific small statistical quadrat. For stations located near the border of two or more small catch quadrats, the biomass was calculated as the mean of the small catch quadrats. For stations located offshore, the nearest adjacent shoreline small catch quadrat was used.

Whenever possible, annual fishermen catches were compared with the census of national coastal fish monitoring (2012−2020). Data showed strong correlations in all stations between catches reported by fishermen and coastal monitoring (Spearman rank correlation: Station 125 rho=0.87, p<0.05; station Kõiguste rho=0.86, p<0.05; station Eru rho=0.83, p<0.05).

The biomass of prey species was obtained from 22 benthic sampling stations. The samples were collected between 1995 and 2020. Data from replicate samples (n=3) of each sampling station were averaged to reduce the variability due to habitat heterogeneity. Ekman and Van Veen type bottom grab samplers were used to estimate the biomass of Macoma balthica in soft sediments and scuba divers collected samples with a 0.04m² metal frame for the biomass estimations of Mytilus trossulus on hard substrates. Sampling and analysis followed the guidelines developed for the HELCOM COMBINE programme (HELCOM, 2015). All biomasses were calculated per m². Physical variables, such as depth, and the areal coverage of hard substrate and perennial species were recorded during sample collection. This habitat information was later used as covariables in prey models (see the section data analyses) as benthic habitat characteristics are known to modulate the behaviour and predation rate of round goby (Liversage et al., 2017).

Information on the physical and biogeochemical conditions (temperature, salinity, oxygen, chlorophyll a, nitrogen, and phosphorus) in the Baltic Sea were downloaded from the EU Copernicus Marine Service data portal (https://marine.copernicus.eu/). Data from the products Baltic Sea Physics Reanalysis (https://resources.marine.copernicus.eu/product-detail/BALTICSEA_REANALYSIS_PHY_003_011) and Baltic Sea Biogeochemistry Reanalysis (https://resources.marine.copernicus.eu/product-detail/BALTICSEA_REANALYSIS_BIO_003_012) were obtained to cover the time frame of 1995-2019. Unfortunately, no data from 2020 were available in these products (Copernicus Marine Service data portal visited on December 10, 2021). To match the time frame of biological data, the 2020 data were obtained from the products Baltic Sea Physics Analysis and Forecast (https://resources.marine.copernicus.eu/product-detail/BALTICSEA_ANALYSISFORECAST_PHY_003_006) and Baltic Sea Biogeochemistry Analysis and Forecast (https://resources.marine.copernicus.eu/product-detail/BALTICSEA_ANALYSISFORECAST_BGC_003_007).

The Baltic Sea Physics Reanalysis and Baltic Sea Biogeochemistry Reanalysis are produced by using the NEMO-Nordic ice-ocean model (based on NEMO-3.6, Nucleus for European Modelling of the Ocean) together with a LSEIK data assimilation scheme. The biogeochemistry product is additionally coupled with the biogeochemical model SCOBI. The Baltic Sea Physics Analysis and Forecast is produced by a Baltic Sea set up of the NEMOv4.0 ocean model. The ocean model is offline coupled with Stokes drift data from the Baltic Wave forecast product. The Baltic Sea Biogeochemistry Analysis and Forecast product is produced by the biogeochemical model ERGOM one way coupled to a Baltic Sea set up of the NEMO ocean model. The products’ documentation provides further details.

The reanalysis and forecast products were available in 4 × 4 km and 2 × 2 km resolutions, respectively. All products contained 56 vertical water column layers, but only surface and bottom layers were used in this study. Monthly means were downloaded for all products. All data were interpolated to a 1 × 1 km grid using inverse distance weighted interpolation with default settings in R package gstat (Gräler et al., 2016).

Only data from recent two years were available in the forecast products. The forecast products were used only because the reanalysis products’ timeframe did not cover the year 2020. There were some systematic differences between the values of variables of reanalysis and forecast products. To address these differences, data from 2019 was used, which was the common year in both products. Using the 2019 data, reanalysis product’s values of each variable in each grid cell were divided with those of forecast product resulting in coefficients in each grid cell. Mean values of these coefficients for each variable were then used to recalculate the values of forecast products to match the levels in the reanalysis products.

Change in the overall prey biomass was tested by calculating average annual biomass before and after round goby invasion and by Kruskal-Wallis test to test the differences in the before and after invasion biomasses in each station (Daniel, 1990).

A selection of environmental variables into final models included only these variables with documented ecological relevance to a prey species (Mytilus trossulus, Macoma balthica, Table 1 ). The Mytilus trossulus model included yearly sea bottom temperature (annual average), salinity (annual maximum), oxygen concentration (minimum of summer months, i.e. June to August average) and sea surface chlorophyll a (maximum of summer months average). Perennial plant biomass, depth and areal coverage of hard substrate were used as a proxy for habitat suitability. The study area temperature range is suitable for M. trossulus, and its variability possibly regulates biomass growth, thus annual averages were used (Marshall et al., 2021). Maximum annual salinity was used owing to M. trossulus inhabiting the edge of its distribution range and the decreasing eastward salinity gradient in the study area indicates that M. trossulus populations likely benefit from annual maximum salinity levels (Westerbom et al., 2002). Summer anoxic conditions are known to affect coastal community functioning and minimal oxygen concentrations can serve as a proxy to investigate its impact on the M. trossulus biomass (Jakubowska and Normant, 2015). Chlorophyll a served as a proxy for M. trossulus food supply, which is quite uniformly distributed in the entire study area; however, extremely high values of chlorophyll a may occur in the inner embayments, which inhibits feeding. Therefore, maximum values were used to characterize the feeding conditions of mussels (Orav-Kotta et al., 2010). The Macoma balthica model included annual sea bottom temperature, salinity, current speed (all annual averages), surface nitrate concentration (maximum of winter months i.e. December to February average), and phosphate concentration (maximum of winter months average). Depth, areal coverage of hard substrate and current speed are used as a proxy for suitable M. balthica habitats. Because extremes in temperature, salinity and current speed have not been shown to affect M. balthica biomass, annual averages were used in the models (Drent, 2002; Jansen et al., 2009). Dissolved nutrients in winter are indicative of food supply availability in summer (Piiparinen et al., 2010) and maximum events are more prone to cause extreme events which indirectly effect M. balthica biomass. The surface level chlorophyll a concentration (summer months averaged between 2000 and 2019) was used as a proxy of productivity level in each station over the entire study period.

Shape Constrained Additive Models (SCAMs) (Pya and Wood, 2015) were used to investigate the relationship between prey species (Mytilus trossulus or Macoma balthica) and key environmental variables before the round goby invasion for years prior to round goby invasion. These models were then run on the complete data to predict the biomass of the prey species. The residuals of these predictive models served as the deviation of the realized biomass and the biomass predicted after removal of environmental effects. Negative residuals indicate greater prey biomass than observed in field sampling. A series of negative residuals associated with an elevated density of round goby suggests increased predation of M. trossulus or M. balthica by the non-indigenous round goby. Prey species biomass values were square root transformed prior to modelling to reduce data skewedness and the effect of extreme values in the models. To increase normality distribution in M. trossulus dataset, depth and perennial plant biomass was log-transformed.

SCAMs are a variant of General Additive Models that allow for the general shape of the relationship between explanatory and predictive variables to be pre-selected. We chose the relationships based on visual examination of the data and theoretical expectation. M. trossulus biomass is anticipated to increase with temperature, salinity, minimum oxygen concentrations, depth, and areal coverage of hard substrate (Westerbom et al., 2002; Jakubowska and Normant, 2015). Biomass is likely to increase with the biomass of perennial plants up to a certain threshold (Kotta et al., 2015). The dataset shows that some areas in the study area exhibit overly high chlorophyll a values, so a declining relationship is expected (Orav-Kotta et al., 2010). M. balthica biomass is anticipated to be greater at moderate depths, as feeding rates are reduced in both shallow areas (too strong wave activity) and deep areas (anoxic conditions). Although a general biomass decline is expected with increased coverage of hard substrate (soft-bottom species), M. balthica is expected to accommodate both detritus feeding and filter feeding in a mixed sea bottom conditions (Olafsson, 1989), thereby suggesting a concave relationship with areal coverage of hard substrate. Low bottom current velocity is expected to improve bivalve feeding rates and thus supports biomass growth (Sobral and Widdows, 2000). M. balthica growth rate is expected to increase with salinity and temperature (Kube et al., 1996; Drent, 2002; Genelt-Yanovskiy et al., 2018). Nutrients are generally thought to enhance M. balthica growth (Seitz, 2011), with increased dissolved phosphorus concentration increasing M. balthica food availability (Beukema et al., 1977). However dissolved nitrogen is expected to promote M. balthica growth only up to a threshold, with excessive amounts increasing potential filamentous algae blooms causing increased mortality, (Tallqvist, 2001), thus a concave relationship is expected.

Spearman rank correlation test with bootstrapping (n=1000) was done between the residuals and the round goby biomass for each station. To investigate potential lags in prey biomass response to the invasion of round goby, up to 4 years of yearly lags were investigated in each station. Median and Median absolute deviation (MAD) was calculated based on the bootstrapped Spearman rho values. For each station a lag or no lag with the highest median Spearman correlation rho was chosen.

Modelling was done in the statistical software R 3.0.3 (R Core Team, 2015): package scam (Pya and Wood, 2015) and mgcv (Wood, 2011) for SCAM, and the packages boot (Davison and Hinkley, 1997; Canty and Ripley, 2021) and rcompanion (Mangiafico, 2016) for correlation and bootstrapping.

The addition of a novel demersal predator provided us with a unique opportunity to study their role in the Baltic Sea ecosystem that previously lacked such predators. Top-down trophic cascades in various ecosystems have been illustrated mainly in literature through the removal of top predators from the ecosystems (Terborgh et al., 2001; Worm and Myers, 2003; Elmhagen and Rushton, 2007; Baum and Worm, 2009). However, these experiments were performed at spatial scales that are severely restricted compared to the spatial scales in which the predator-prey interactions operate. Our study demonstrated a substantial large-scale and long-term change in the food web, with previously primarily bottom-up regulated communities currently facing strong top-down control effects. These effects have the potential to spread to upper (e.g. Liversage et al., 2017) and lower trophic levels (e.g. Kotta et al., 2018) through trophic cascades.

A true assessment of the effects of long-lasting invasion requires long-term research (Bonanno and Orlando-Bonaca, 2019; Catford et al., 2019; Strayer et al., 2019; Crystal-Ornelas and Lockwood, 2020). Previously invasive round goby effects have only been shown in experimental studies (Pagnucco et al., 2016; Henseler et al., 2021) or short-term field observations (Kornis et al., 2013; Skabeikis et al., 2019; Armoškaitė et al., 2021; van Deurs et al., 2021). Predicting long-term impact of invasion from field data is challenging as the environment can change during the invasion, which in turn can affect the results (Denley et al., 2019; Dickey et al., 2021). This study demonstrates long-term trends in benthic communities and by excluding possible effects due to the changing environment, provides more evidence that invasive species alter the functioning of coastal ecosystems over time. Invasion effects are known to be stronger when invaders occupy functions previously not fulfilled by native species (Case et al., 2016b). This makes especially fragile ecosystems that are characterized by a few species in important niches, allowing an invasive predator to have long lasting effects on benthic communities and thus on the whole ecosystem functioning.

Invasion ecology predicts that some communities are more resistant to invasion effects (González-Moreno et al., 2013; Beaury et al., 2020). Here a variety of habitats were investigated and based on the nature of the less impacted areas, we speculate areas of higher productivity are more resilient to round goby invasion. Both prey species, M. balthica and M. trossulus, rely directly or indirectly on phytoplankton primary production with greater microalgal biomass resulting in significantly greater biomass yield of the species (Beukema et al., 1977; Kotta et al., 2015). Higher growth rates enable bivalve species to avoid predation by round goby and/or recover faster. Moreover, high productivity areas are characterized by high biomasses of alternative food sources (e.g. amphipods, isopods, snails, other bivalves) further reducing the predation of M. balthica and M. trossulus by round goby. Similarly stronger effects from top-down pressure in low-productivity systems have been described among herbivores (Burkepile and Hay, 2006). Although less significant effects were observed in high productivity areas, some stations given enough time already show clear and significant impacts of the invasion. Such patterns could globally cause underestimating invasion events effects, if low productivity communities are expected to have similar reactions to invasive species and failing to look at larger patterns in community functioning.

Newly-invaded habitat often experience time lags before the effects of invasion become apparent (Thorlacius and Brodin, 2018). As such, round goby invasion in several areas has likely been too short to exhibit any effect, similar to other invaded areas (Kornis et al., 2013). In the five years (since 2016) round goby has been invading the Baltic Proper, large M. trossulus reefs have not yet been affected by round goby invasion, but lower maximum M. trossulus biomass values have already been recorded. Thus, given enough time round goby will likely significantly affect all communities it invades. Recent studies have shown that the round goby is a very mobile consumer, able to move several kilometres daily with even longer seasonal migrations (Christoffersen et al., 2019; Behrens et al., 2021). If an area contains many suitable habitats (such as in stations Aegna, Muuga, 12.5, PW), these movement patterns could prolong time lags of the effect of round goby on individual communities.

Information regarding native community recovery after invasion is limited to a few examples (Case et al., 2016a). Even after total removal of an invader, recovery to previous biomass levels has been unsuccessful (Maclean et al., 2018). With the decline of round goby abundance after invasion expansion phase, prey population recovery was observed in several stations. In areas with sufficient population of several prey species, the availability of an alternative food resource may provide a buffer allowing prey populations to recover. The observed recovery further suggests that round goby effects are not linear and given enough time prey-predator population dynamics could be observed. Apparent recovery of prey population to pre-invasion levels was observed in only a few stations (Aegna and Muuga). If large scale population biomass recovers to a substantially lower population size than before, the ability of their ecosystems to act as a natural filter to mitigate eutrophication effects are suspect. Such areas might benefit from additional eutrophication mitigation action.

Like many other ecosystems around the globe, the Baltic Sea is heavily impacted, and its communities experience multiple pressures simultaneously. Moreover, the Baltic Sea is very dynamic and benthic species inhabiting its coastal ecosystems are adapted to cope with harsh environments and are tolerant not only to variable salinity conditions, but also to other environmental variables, such as hypoxia and exposure. In such systems, physically driven fluxes often override the effects of biological interactions (e.g. Herkül et al., 2006). Thus, accurate assessment of biological interacts, such as species invasion, demands proper evaluation of confounding environmental factors.

This study provides valuable methods on how to model long term general-purpose data, to filter out environmental effects and subsequently detect possible effects of invasive species, as described in Kotta et al., 2018, where predatory crab-induced regime shift was apparent only when observed in changing environmental conditions. Similarly, our results showed that the round goby is likely the cause of the decline in benthic bivalve biomass, but this decline can be properly quantified only if other forcing factors are considered. Our “filter” model did not consider other benthic predators and therefore some of the observed decline can be attributed to interactions not covered in the current study. However, this is not likely as predation rates of the native bivalve populations by native predators are very low and all earlier evidence has suggested that the study area is primarily a bottom-up regulated system with minor small-scale top-down effects (Kotta et al., 2018). In addition, the round goby is the only fish to exploit both shallow coastal habitats with M. trossulus and deeper habitats with M. balthica (Behrens et al., 2021); thus, a simultaneous shift in the bivalve biomass in these very different habitats again point to the round goby as a causal factor. Moreover, all documented significant shifts in the biomass of benthic community have been related to shifts in nutrient loading and hence changes in the productivity of ecosystem (Bonsdorff et al., 1997; Kotta et al., 2007). Nevertheless, this aspect was considered in the modelling frame of the current study.

When dealing with marine ecosystems there are several unknowns related to sampling intensity both in space and time. Improved monitoring intensity and quality will allow researchers to apply the expected predation pressures of the round goby with a broad range of habitat characteristics beyond a single habitat and region. Moreover, the method can be used for other invasive species and regions to test hypotheses on the determinants of success of introduced species, the magnitude and dimensions of their impact, and the mechanisms sustaining successful invasions (Chabrerie et al., 2019).

To conclude, the addition of an invasive benthic predator is causing significant changes in the ecosystem with the potential for the changes to proliferate to other trophic levels by trophic cascades. We found that the round goby causes significant declines in benthic bivalve biomass in the north-eastern Baltic Sea, which is not associated with overall changes to the physical and biochemical environment. The severity of the impact is likely related to overall productivity in the area; recovery of the benthic community to pre-invasion biomass levels is infrequent. However, the ability to assess impacts of nonindigenous species is contingent on suitable monitoring of the area. The Baltic Sea provides an excellent case study to show a way forward to analyse correlative linkage of invasion and species patterns.