We used species inventories from the three largest peri-Alpine French lakes: Lake Annecy (45°51′ 41.489″ N, 6°10′ 2.364″ E), Lake Bourget (45°43′ 46.842″ N, 5°52′ 10.484″ E) and Lake Geneva (46°26′ 27.213″ N, 6°30′ 38.177″ E) originating from annual monitoring surveys (i.e., recording of environmental parameters and biodiversity samplings) and scientific reports (©SOERE OLA-IS, INRAE Thonon-les-Bains, SILA, CISALB, CIPEL; Rimet et al., 2020). As these lakes are highly similar in terms of biodiversity (Jacquet et al., 2014), species co-occurring in the three lakes (n = 118; from primary producers to large vertebrates) were retained to represent typical species of these ecosystems and were characterized by their taxonomy (i.e., subphylum, group, family, genus, species), their average body size (μm) commonly used in community ecology studies, a habitat trait (i.e., littoral, pelagic/littoral or pelagic) and a feeding trait (i.e., carnivorous, omnivorous, herbivorous or primary producer) referenced in species inventories. Animal species were then clustered to the family level (S = 48) and vegetal species were clustered to the class level (S = 10). This taxonomic aggregation led to 58 species clusters (hereafter called SC; Supplementary Table 1) gathered according to similar ecological functions/requirements and could hence be considered as functional nodes in the reconstructed food webs (Allesina and Pascual, 2009). This clustering procedure allowed to avoid unnecessary food web complexity by taking into account the main SC that may be directly or indirectly impacted by catfish in line with recommendations emerging from recent ecosystem modeling studies (Geary et al., 2020), to promote the computational efficiency and to facilitate ecological interpretations of the processes involved in abundance changes of SC due to S. glanis.

We used the aNM (Vagnon et al., 2021) to infer trophic interactions between SC included in the typical peri-Alpine lake food web and to reconstruct the “non-invaded” (i.e., without S. glanis) and “invaded” food webs (i.e., with S. glanis). This model relies on the niche model principles (Williams and Martinez, 2000) stating that the niche position of consumer j is given by its average body size bs_j and that its resources fall within a body size range bs_r_j centered on bs_c_j. The range bounds were estimated using quantile regressions (i.e., bs_r_jmin = QR at 5% and bs_r_jmax = QR at 95%) as suggested by Gravel et al. (2013) and are specifically fitted whether consumers are vertebrate or invertebrate (Vagnon et al., 2021).

The aNM allowed obtaining a binary squared matrix (M_b) of trophic link occurrences of the whole food webs. These links were then weighted for each consumer j considering a Gaussian probability density function, similarly to Williams et al. (2010), with μ_j = bs_c_j and σ_j = standard deviation of 100 points evenly spaced over bs_r_j (i.e., scaling of the normal distribution to bs_r_j). Weighted links were then normalized by the maximum value of the normal distribution to obtain a maximum weighting (0.95) at bs_c_j and a minimal weighting (0.1) corresponding to prey SC with body sizes at bs_r_j bounds. The resulting weighted links were finally converted as the proportion of resources i in the diet of consumer j (ω_ji, eq. 1), so that the total proportion of species in the diet of consumer j sums to 1.

In a first step, trophic interactions between SC were first inferred without S. glanis to reconstruct the food web before invasion. Trophic positions (TP) of resident species were estimated according to the method of Levine (1980) with TP_{Primary producers} = 1 and TP_Consumers = 1 + mean (TP_Resources).

In a second step, we independently introduced 40 nodes representing S. glanis with increasing body sizes (i.e., 5–200 cm by 5 cm) in the SC inventory to study changes in the food web topology as S. glanis could cause different impacts due to its size-dependent diet (i.e., ontogenetic diet shift from invertebrates to fish; Carol et al., 2009; Copp et al., 2009; Alp, 2017). The directed connectance (i.e., number of actual links over the number of possible links) of each invaded food web was calculated to provide an estimate of interaction variations among resident species and S. glanis of different body sizes (Supplementary Figure 1; Bersier et al., 2002). Three body sizes associated with the highest variations of the directed connectance were selected to simulate three invasion scenarios: small body size 40 cm (S40), medium body size 85 cm (S85) and large body size 150 cm (S150). The trophic interactions inferred for the corresponding three invaded food webs were included in the following steps of the analysis.

The simulations of population dynamics were based on an updated version of the allometric population dynamic model proposed by Yodzis and Innes (1992) and extended by Brose et al. (2006). The dynamics of each primary producer i are given by:

The parameters are the intrinsic growth rate r_i and biotic capacity K_i. Both are allometrically scaled (Supplementary Table 2) differently from the ones proposed by Brose et al. The functional response (i.e., the per capita consumption rate of consumer j on resource i), F_ji(N→) follows a Holling Type II and is given by:

The parameters are x_j, the metabolic rate body mass dependent, and y_j the maximum consumption rate relative to the metabolic rate of consumer j. The former is newly allometrically scaled, while the latter depends on characteristic of the species (Supplementary Table 2). ω_ji is the proportion of resource i in the diet of consumer j (eq. 1) obtained from the weighting procedure previously described. The handling time h_jk of consumer j (i.e., time for handling and consuming the resource k) is also allometrically scaled (Supplementary Table 2), differently from the parameterization of Brose et al. for providing realistic estimations of this parameter. Finally, the dynamics of consumer j is given by:

The parameters are the morality rate m_j (allometrically scaled) and the consumption efficiency ε_j. (depends on species characteristics; Supplementary Table 2). Note that we choose N_i and N_j as the SC abundances (number of individuals in the SC) and not as the biomass for more stable calculations. Therefore, all allometrically scaled parameters are per individual (see Supplementary Table 2).

We used 200 simulations of 20,000 time steps for all scenarios (no. S. glanis, S40, S85, and S150) and final results were obtained by averaging simulation outputs at each time step for each SC. Variability in metabolic rates was included in calculations by generating white noise following 𝒩 (μ = 0, σ² = 0.1% of metabolic rate). For each simulation, random initial abundances (N₀) were assumed to be uniformly distributed (Brose et al., 2006) on the interval [0.15–1] and were ranked according to average body sizes of SC (i.e., SC with the largest body size had the lower N0; Peters and Wassenberg, 1983). N₀ of S. glanis was fixed at 0.1 in all simulations to be lower than N₀ of the resident apex predator SC (Esocidae), to avoid bias resulting from the initial abundance variations of S. glanis.

The impacts of S. glanis on resident species were investigated through two main ecological processes that are SC extinctions and changes in non-extinct SC abundance. Extinctions were assumed to be effective when the abundance of a SC fell below 1.10^–6 (Ovaskainen and Hanski, 2003) and extinct SC were not permitted to reintegrate the system. Extinctions were explored both qualitatively (i.e., SC taxonomic category) and quantitatively (i.e., number of extinctions and time lags between time at extinction with and without S. glanis).

The impacts on the abundances of non-extinct SC were investigated by calculating the relative abundance difference (RAD; Eq. 5) at each time step to describe increase/decrease in abundances with S. glanis compared to abundances without S. glanis similarly to Zhang et al. (2019).

The median RAD (MRAD) for the 200 simulations was used to quantify the extent (i.e., amplitude of changes) as well as the type (positive, MRAD > 0, or negative, MRAD < 0) of S. glanis impact on SC abundances (i.e., increase or decrease in abundances).

The model robustness was assessed by including an increasing variability in the metabolic rate using white noise (i.e., 𝒩 (μ = 0, σ² = 0.1, 0.5, 1, 5, 10, and 20%; Supplementary Figures 3–9) using MRAD as the response variable.

Impacts of S. glanis on SC (extinct or with modified abundances) were classified according to three main ecological interactions based on the food web structure inferred with the aNM:

(i) Predation when a SC was a prey for S. glanis;

(ii) Competition when a SC shared common resources but did not directly interact with S. glanis. In this case, Schoener’s overlap index was complementarily calculated to quantify diet similarities between both competitors (Vera-Duarte and Landaeta, 2017) as follows:

where P_xi and P_yi are the proportions of resource i in the inferred diet of competitor x and y;

(iii) Trophic cascades were considered when a SC was not a prey of S. glanis and corresponded to repercussions across multiple trophic levels (Carpenter et al., 2010).

Interactions not corresponding to these three ecological interactions were classified as “Others.”

All statistical and graphical displays were performed using R.3.5.1 (R Core Team, 2018) with the packages ade4 (Dray and Dufour, 2007), deSolve (Soetaert et al., 2010), cheddar (Hudson et al., 2013), igraph (Csárdi, 2019), NetIndices (Kones et al., 2009), foodweb (Perdomo, 2015), ggplot2 (Wickham, 2016). R Markdown files with data and R codes summarizing the conducted analyses are available at https://github.com/chloevagnon/aNM-and-population-dynamics to provide reproducible examples.

Extinctions did not relate to S. glanis body sizes as the same extinctions occurred for the three different invasion scenarios (Figure 1). Extinct species clusters (SC) corresponded to different trophic levels (from primary producers to secondary consumers) and were mainly characterized by a high pelagic reliance (Supplementary Figure 11). In contrast, S. glanis body size was determinant in the amplitude of abundance changes of the non-extinct SC (i.e., MRAD from –0.07 to +0.15). The strongest negative impacts were found for SC with TP (Trophic Position) close to the invader’s TP in the food webs while positive impacts mainly concerned SC with low TP (Figure 1). The lowest and the highest MRAD were observed for S85 and suggested positive impacts on primary littoral and pelagic consumers, and negative ones on SC with TP > 2.5 (Figure 1). Littoral secondary consumers were mainly impacted for S40 and the lowest MRAD were found for S150.

Among the 58 SC, 14 were extinct at the end of the simulations in each scenario (Figure 2A). Three major extinction phases were noticeable and concerned in the first place zooplankton (up to 1,043.3 ± 0.5 time steps), followed by invertebrates and phytoplankton (up to 10,475.8 ± 4.3 time steps), and finally one invertebrate (at 17,915.5 ± 35.2 time steps). The same extinctions patterns occurred for all scenarios while time lags between extinctions differed, particularly for the latest extinctions (Figure 2B). Extinctions tented to occur earlier for S85 and S150 while those could appear later for S40 compared to extinctions without S. glanis.

RAD (Relative Abundance Differences) varied among SC taxonomic categories and for the different S. glanis body sizes. Fish SC were impacted for all scenarios while presenting variable changes as shown by the boxplots representing RAD at each time step (Figure 3A). The three smallest fish SC and the largest one were especially negatively impacted by S85 (lowest RAD and highest variability), followed by S150 especially for the largest fish SC. In contrast, negative effects of S40 were more pronounced for medium-size fish SC and even appeared positive for two fish SC.

Half of invertebrate SC had higher abundances in presence of all S. glanis and presented a high variability along time steps, particularly with S85 (i.e., highest RAD; Figure 3B). The other half of invertebrate SC was negatively impacted by S40 for all time steps (i.e., all RAD values < 0) and the three largest invertebrates were negatively affected by S85 and S150.

RAD for phytoplankton, phytobenthos and zooplankton SC were similar among scenarios and more variable in time for phytobenthos SC (Figure 3C). Abundances of zooplankton and phytobenthos SC slightly increased for S40 while they decreased in the scenarios with the two larger S. glanis (S85 and S150). Abundances of half of phytoplankton SC increased with S. glanis body size while the reverse was found for the other half of phytoplankton SC.

Extinctions may have been caused by trophic cascades regarding the low trophic position of extinct SC compared to that of the three S. glanis, even for one invertebrate SC predated by S40 (Figure 4A and Supplementary Table 3), but did not seem directly triggered by S. glanis as they also occurred without the predator in the food web.

Ecological interactions for non-extinct SC were much more diverse among scenarios and corresponded to both negative and positive impacts mediated by direct and indirect interactions, even for three interactions with fish SC for S40 and S85 (Figure 4B and Supplementary Table 3).

Direct negative impacts corresponded to predation with 13 SC for S40 (91% of invertebrate SC), seven SC for S85 (large invertebrates and medium-size fish SC) and 10 SC for S150 (crayfish and all fish SC).

Indirect negative impacts were emphasized through competition identified for different fish SC when considering the invasion of S40 (e.g., Lotidae, Cyprinidae, Ictaluridae; Schoener’s index of 0.78, 0.54 and 0.53, respectively), S85 (e.g., Esocidae; Schoener’s index = 0.91) and S150 (Esocidae; Schoener’s index = 0.71). Negative impacts were also due to trophic cascades, concerning three SC for S40 and increasing with S. glanis body size (i.e., 10 SC for S85 and 14 SC for S150).

Among the different ecological interactions, Indirect positive interactions were preponderant with 21 trophic cascades identified for S40 (composed at 67% of invertebrate SC, 5% of zooplankton and 28% of primary producer SC), 22 trophic cascades for S85 (composed of 82% of invertebrate SC, 9% of zooplankton and 9% of primary producer SC) and 19 trophic cascades for S150 (100 % of invertebrate SC).

The sensitivity analysis revealed a gradual increase in MRAD variability in response to the increasing metabolic noise. However, MRAD distribution remained consistent among scenarios and SC categories for metabolic noise up to 5% and then started to flatten, traducing a decrease in the consistency of the results with higher metabolic noises (Supplementary Figures 3–9). Especially fish SC seemed more sensitive than other SC to metabolic noise. Therefore, our simulations could be sensitive to high metabolic rate variability while patterns in responses remained robust according to taxonomic categories and among the invasion scenarios.

The food web topology was only slightly modified by the invasion of S. glanis, which exerted moderate impacts on connectance and abundance modifications (i.e., low connectance changes and restrained amplitude of abundance modifications). This result may be explained by the rather high resolution of the food web, as associated to omnivory that both promote weak interaction strengths (Supplementary Table 4) and stability in the persistence of communities (Emmerson and Yearsley, 2004; Dunne et al., 2005; Landi et al., 2018; Kawatsu et al., 2021). Ecological interactions and abundance dynamics were nevertheless influenced by S. glanis and were modulated according to its body size.

The highest negative impacts of S. glanis were related to direct interactions for resources acquisition. Predation was indeed the main ecological interaction associated with the decrease in SC abundances, in lines with numerous impact studies on invasive fish (Gozlan et al., 2010; Van der Veer and Nentwig, 2015; David et al., 2017). Specifically, the decrease in abundance of invertebrate SC (e.g., Ephemeridae) and small fish SC (e.g., Blenniidae) would be caused by S. glanis of 40 cm while the decrease of larger fish SC (e.g., Percidae) and large invertebrate SC (e.g., crayfish) would mostly result from the invasion of larger S. glanis (85 cm and 150 cm). Percidae and crayfish represented consistent fraction of S. glanis diet in empirical studies (Carol et al., 2009; Ferreira et al., 2019) and the presence of S. glanis was suspected to alter their abundance in few ecosystems (Copp et al., 2009; Guillerault et al., 2015; Vagnon et al., 2022), suggesting plausible predictions from our approach. Competition between S. glanis and large fish SC was also suggested to negatively alter fish abundances and particularly dampened the abundance of the resident apex predator (Esocidae; high diet overlaps for S85 and S150) while fish SC with weak diet overlaps (Schoener index < 0.1) were poorly or not impacted. These results supported both plausible inferences of trophic interactions for predators and species abundance alterations due to both predation and competition. Indeed, invasive predators frequently cause decreases in abundances of their prey (Mills et al., 2004; David et al., 2017), they are expected to have a strong competitive effect in aquatic ecosystems and they are known to induce changes in species assemblages (Gozlan et al., 2010; Allesina and Tang, 2012; David et al., 2017).

Interestingly our results highlighted that considering two trophic levels (i.e., consumer/resource) or two competitors does not allow to identify major impacts of invasive species and their repercussions at the whole food web scale as we underlined main abundance changes for SC through trophic cascades, and particularly we found higher positive impacts than negative ones. A succession of negative and positive impacts was indeed noticeable along the food web and emphasized typical patterns of top-down cascades found in aquatic ecosystems (Carpenter et al., 2010; Heath et al., 2014; Su et al., 2021). In these processes, the negative impacts on abundances of SC such as fish result in a relaxed predation on lower trophic levels (e.g., grazers) and thus in a reduction of primary producers (Carpenter et al., 2008; Heath et al., 2014; Koning and McIntyre, 2021). Consequently, positive and negative impacts were higher for S85 than for S40 and S150, as fish and large invertebrate SC were particularly impacted in this invasion scenario (i.e., competition with fish SC and predation on both fish and large invertebrate SC). Complementarily, trophic cascades were suggested to modulate time at extinctions (earlier for S85 and S150 and later for S40), yet the presence of S. glanis did not qualitatively influence extinction processes (i.e., same extinctions with/without S. glanis). These results underline the ability of S. glanis to cause top-down trophic cascades by regulating mesopredator abundances while suggesting that it may not be a major source of species extinctions, similarly to empirical studies conducted in reservoirs, lakes and rivers (Copp et al., 2009; Vejřík et al., 2017). In fact, the opportunistic feeding behavior of S. glanis (Copp et al., 2009; Cucherousset et al., 2018; Vagnon et al., 2022) could foster its reliance on a diversified prey set limiting strong interaction strengths usually known to induce stronger impacts on resident species populations than weak interactions (Terraube et al., 2011; Wootton and Stouffer, 2016).

Our study based on the combination of an allometric niche model and an alternative version of the population dynamic model of Brose et al. (2006) succeeded at supporting a comprehensive assessment of resident species abundance modifications and the involved trophic interactions considering a multi-trophic system, often challenging to evaluate only based on experiments or on traditional empirical methods (Crystal-Ornelas and Lockwood, 2020), notwithstanding process simplifications inherent to the elaborated method and to our study objectives.

Firstly, we considered an average body size to represent species nodes in food webs and we integrated independently a unique S. glanis body size node in simulation scenarios. While we recognize that diets of both resident and S. glanis populations can be ontogenetic-dependent, we conserved this approach commonly approved in community ecology to preserve the real significance of each node in the food web (i.e., one S. glanis node = one invertebrate node = one phytoplankton node). Secondly, we assumed environmental and anthropic drivers constant to capture specifically the invasive species impacts in the recipient ecosystem. Although these factors are known to be dynamic in real systems (Brose and Hillebrand, 2016), our approach appeared relevant considering the different time scales involved in species abundances modifications following exposure to various external pressures. Indeed, the catfish impacts could be observed more rapidly than external factors such as climate change (i.e., decades vs. several decades). Our dynamic model could obviously be completed in future studies to account for possible other drivers of species trophic interactions and their dynamics (e.g., metabolism modifications following temperature increase due to climate change) but remained out of the scope of this study.

However, the combination of allometric models would be appropriate for an application to a broad aquatic ecosystem array where body size governs trophic interactions, mainly thanks to the minimal required data inputs often available from monitoring surveys or literature (i.e., species inventory and the average species body sizes). Here we used species clusters for convenience and for limiting computation concerns such as singularities but different taxonomic resolutions can be considered depending on the initial study scope. Moreover, allometric parameterizations in the dynamic model could also be replaced by empirical data (e.g., body mass, metabolic rate) and/or can be completed by other calculation methods, for instance, to infer interaction strengths between consumers and resources (Calizza et al., 2021). When these data are not available, our initial simulation parameterization appears valuable regarding the convergence between our results and literature on predator impacts in food webs and freshwater ecosystems (Cucherousset and Olden, 2011; Cucherousset et al., 2012; Jackson et al., 2017), although empirical validations of our simulations still remain out of reach. Indeed, management plans mostly focus on a part of species compared to the totality of resident species to survey, mainly due to the restricted resources available for management actions (Vander Zanden and Olden, 2008), thus limiting the consideration of all taxa responses at long term that could be underlined in our study.

Overall, our study framework addresses the impact of invasive species and may be relevant regarding the increasing rate of species introductions representing a major threat to ecosystems (Lockwood et al., 2013). We underlined the importance of considering interspecific interactions at the whole food web scale for the assessment of invasive species impacts. The introduction of invaders indeed frequently involves a wide diversity of new ecological interactions, resulting in modifications of species abundances through both direct and indirect interactions and thus causing non-negligible impacts on resident communities. The balance between negative and positive aspects of invasions is also a significant factor to consider as positive effects of invasions can appear non-negligible, and may attenuate a priori expectations (Gozlan, 2008; Tablado et al., 2010; Schlaepfer et al., 2011). Our study can thus participate to the growing corpus of methodologies trying to reach comprehensive assessments and predictions of invader impacts, considering their direct/indirect and positive/negative effects in freshwater ecosystems.