Biological invasions and eutrophication reshape the spatial patterns of stream fish size spectra in France

The global patterns of body size distributions are affected by global environmental changes (GEC), but our knowledge of the interaction effects between GEC and natural drivers is still limited. In this study, we aimed to test the effects of these interactions on fish community size spectra, that is, the variation in a community property across the body size range of individuals in the community.


| INTRODUC TI ON
Global environmental changes (GEC) are reshaping the spatial patterns of community structure (Piggott et al., 2015), and understanding these responses requires approaches that emerge from the individual-to community-level processes. Body size provides reliable information on the biological rates of individuals (e.g. metabolism, reproduction, growth), and allometric body-size scaling relationships can help to better represent the structure and dynamics of communities (Woodward & Warren, 2007). The community size spectrum, that is, the variation in a community property across the body size range of individuals in the community (Rice & Gislason, 1996) that can take different forms like the biomass spectrum (Boudreau & Dickie, 1992), is one of the universal body-size patterns because it reflects community structure and stability (Trebilco et al., 2013), predator-prey dynamics (Hatton et al., 2015), and ecosystem functions (Bartrons et al., 2020;Mehner et al., 2018). The community size spectrum can therefore represent a powerful indicator of community responses to GEC (Petchey & Belgrano, 2010).
Despite the importance of the size spectrum as a key framework to simplify the community structure (Andersen, 2019;Kerr & Dickie, 2001), surprisingly few attempts have empirically determined if GEC interacted with natural factors to shape the community size spectrum. The shape of the community size spectrum is consistent under steady-state conditions due to the energetic constraints of food assimilation and acquisition from smaller to larger organisms (Andersen, 2019; Kerr & Dickie, 2001). However, substantial variations of the community size spectrum have been reported empirically among contrasted environmental contexts (Benejam et al., 2018;Pomeranz et al., 2022). For instance, warmer conditions can steepen the community size spectrum slope (i.e. an increase of the relative proportion of small-bodied organisms; Daufresne et al., 2009;Pomeranz et al., 2022) because increased temperature generally benefits smaller organisms through increased metabolic rates (Atkinson, 1994;Bergmann, 1847;Daufresne et al., 2009). In stream fish communities, Benejam et al. (2018) found consistent altitudinal changes, with a flatter size spectrum towards the upstream direction corresponding to systematic changes in life-history traits and species richness (Santoul et al., 2005;Vannote et al., 1980). Global environmental changes is a complex and multifaceted phenomenon with a potential interplay between multiple stressors likely to generate unpredictable effects (i.e. ecological surprises; Jackson et al., 2016) by amplifying or alleviating the independent effects of natural factors (Piggott et al., 2015). Freshwater ecosystems face GEC that affects individuals' behaviour and metabolism and that cascade into changes in community size spectra (Woodward et al., 2010). The community size spectrum is therefore central in understanding and predicting ecological disturbances in response to GEC, such as climate warming, biological invasions, and eutrophication. Climate warming steepens the size spectrum under experimental conditions (Dossena et al., 2012;Yvon-Durocher et al., 2011) although this response was not universal in natural ecosystems (O'Gorman et al., 2017). Steeper size spectra were responsive to an increase in total phosphorus (TP) in the water for fish (Mor et al., 2021) and planktonic (Atkinson et al., 2021) communities because high levels of eutrophication can reduce the number of pollution-sensitive predators and resulted in energy limitation to the higher trophic levels (Mor et al., 2021). Biological invasions can blur community body size patterns since the body size of non-native species is often larger than native species (Blanchet et al., 2010). Biological invasions could therefore disrupt ecological communities towards a topheavier pyramid structure, equivalent to a flatter size spectrum (Arranz et al., 2021;Kopf et al., 2019). Taken together, one may expect that the natural patterns of community size spectra can be modulated by GEC factors in ecological communities facing multiple anthropogenic pressures.
In this study, we used individual body sizes of stream fish across a large environmental gradient in France to (i) quantify the spatial variability of size spectra in stream fish communities and (ii) test the interacting effects of natural factors and GEC in driving fish community size spectra. We focused our investigation on climate warming, biological invasions and eutrophication because of their critical impacts on the body size distributions of freshwater fish (Brucet et al., 2013;Côte et al., 2019;Leprieur et al., 2009). We hypothesise that GEC drivers modulate the effects of natural factors on spatial patterns of size spectra in stream fish communities through synergistic or antagonistic interaction responses (Jackson et al., 2016).

| Data collection
We compiled a database of fish species abundance and individual body size in stream locations from collections made by the Office Français de la Biodiversité (OFB; Irz et al., 2022). A standardized methodology was employed in all locations using an electrofishing sampling protocol (CEN, 2002;Poulet et al., 2011). Sampled fish were individually identified to species, counted, measured to the nearest millimetre (fork length), and weighted before released. When individual fish mass was not weighted during sampling, it was estimated using specific-species length-mass relationships. From this database, we implemented a data filtering process to select stream locations that can accurately reflect the spatial distribution of body size structure in stream fish communities. We focused on stream locations sampled once between 2015 and 2018 to limit temporal variability. We selected stream locations sampled during low-flow periods (June-October). To avoid potential biases in size spectrum calculations, we also selected stream locations where the total number of fish sampled was >100 individuals (mean = 621.8 ± 653.07 SD). This database contained 1095 stream locations and 680,832 individuals belonging to 61 fish species.

| Natural factors
We considered two hydromorphological variables to represent the upstream-downstream gradient in the stream network: surface area of the drainage basin upstream of the sampling site (km 2 ) and the distance from the headwater source to the sampling point along the stream network (km) at the spatial resolution of 0.5 km (Table S1). The hydromorphological data were obtained from BD CarTHAgE database (Table S1). We used a principal component analysis (PCA) to summarize these two variables. The first axis of the PCA, accounting for 97.5% of the total variability, was used as a synthetic variable describing the upstream-downstream gradient in the stream network (Buisson et al., 2010). The negative values of the first PCA axis corresponded to the most downstream locations and positive values to the most upstream locations. We also estimated climatic conditions (°C) as the average temperature between 1960 and 1970 for the summer (June-August) period (Table S1) to obtain baseline conditions. We used daily air temperatures provided by Météo France and extracted from the highresolution (8 km by 8 km grid) SAFRAN atmospheric analysis over France (Le Moigne, 2002). To account for the sigmoidal-shaped relationship between air and water temperature (Mohseni & Stefan, 1999), we transformed air temperature to water temperature following Punzet et al. (2012). Finally, species richness was determined as the number of species sampled in each stream fish community (Irz et al., 2022; Table S1).

| Global environmental changes
For each stream location, we collected publicly available data of TP concentration (Table S1) and assessed the eutrophication level by calculating the mean annual concentration of TP (mg L −1 ) measured during the same year as the fish sampling (Table S1). Climate warming was computed by the ordinary-least-square (OLS) linear regression trend starting from the average water temperature between 1960 and 1970 until the average water temperature in the year of the fish sampling for the summer (June-August) period ( Table S1).
The linear regression estimate was a surrogate of climate warming (°C year −1 ), with positive coefficients indicating an increase in water temperature over time while negative coefficients indicated a decrease. Finally, biological invasions were quantified using the percentage of individuals belonging to non-native species in each stream location (Table S1).

| Community size spectrum
We calculated the size spectrum by the rate of decline in abundance across body size classes (Kerr & Dickie, 2001) for each location. First, we classified the body size distribution of individual fish mass into 10 size classes following a geometric series of two where size intervals were narrow for small body sizes but became progressively wider with increasing body size (Table S2; Sprules & Barth, 2016).
We used fish body mass rather than body length to calculate size spectrum because it reflects the amount of energy within organisms and follows early bioenergetic size spectrum models (Platt & Denman, 1977;Silvert & Platt, 1978). This, however, would not modify the results since the size spectrum slope remains unchanged depending on the body size type, such as mass, length or volume (Guiet et al., 2016;Sprules, 2022). Because electrofishing can lead to biases by underrepresenting the smallest individuals (<4 g; Hense et al., 2010), we grouped the smallest fish into the first size class (1st size class, <4 g, Table S2). However, accumulating the smallest fish into the first size class did not lead to biased results because further calculations of the size spectrum slopes showed similar results when the smallest fish (<0.5 g, representing 6.12% of the total fish data set) were removed ( Figure S1). Second, we normalized the fish abundance by dividing it by the width of each size class (Table S2) because the width of the size classes progressively increases with body size, having a potential influence on size spectrum calculations (Sprules & Barth, 2016). The normalization avoids these effects and allows comparisons between studies (Sprules & Barth, 2016). The size spectrum slope was computed in each fish community from an OLS regression model between the normalized fish abundance (y-axis) and the midpoint of the size intervals (x-axis) at a log-log scale. Flatter slopes (closer to 0) suggest a relatively large abundance of large-bodied individuals. In contrast, steeper slopes (more negative values) indicate a higher relative proportion of small-bodied individuals (Kerr & Dickie, 2001;Sprules & Barth, 2016). Coefficients of determination (adjusted R 2 from OLS regressions) for the spatial variability of size spectrum slopes were uniformly high (mean = 0.89 ± 0.10 SD), validating the quality of the OLS regression model.
In addition, we also calculated size spectrum slopes using alternative methods, including binning with different bin widths and maximum likelihood methods using modifications of the 'sizeSpectra' package (Edwards et al., 2017; Figure S1). The size spectrum slopes differed among methods, but they were all strongly correlated (Pearson's rank correlations ≥ .80; Figure S1). We selected the OLS regression approach because Xiao et al. (2011) showed that OLS regressions on binned data work best under the multiplicative error structure that almost certainly holds for our data (86% of the stream locations followed a multiplicative error structure).

| Statistical analyses
Multicollinearity among the explanatory variables was checked using variance inflation factors (VIF; Dormann et al., 2013), and VIF values among predictors were low (all VIF < 3). The mean annual concentration of TP and percentage of individuals belonging to nonnative species were log-transformed to minimize skewness and kurtosis. All predictors were standardized and converted into a z-score by subtracting the sample mean from each variable and dividing by its standard deviation. The size spectrum slopes calculated using fish body mass were used as response variables. We performed a model-averaging approach to cope with model uncertainty and analysed all possible candidate models using the 'glmulti' package (Calcagno & Mazancourt, 2010). To prevent overfitting and facilitate computation tasks, we discarded models with more than five explanatory variables and interactions among natural factors (in total, 27,896 candidate models). We ranked the models based on corrected Akaike's information criterion (AIC c ; Akaike, 1974). The set of models with a difference in AIC c (ΔAIC c ) < 2 from the best model were considered to have equivalently strong empirical support and equal plausibility (Burnham & Anderson, 2002). The AIC c weight of the best model represents the probability that a particular model is selected as the best fitting model if the data are collected again under identical circumstances (Lukacs et al., 2010;Whittingham et al., 2005). We then calculated model-averaged coefficients across the selected models only when the predictor appeared once or more. We did not control for spatial autocorrelation in the modelling approach because there was no significant spatial autocorrelation in size spectrum slopes among sites (Mantel correlogram analysis; global Pearson's r = .009; p-value = .270). All analyses were conducted with the software R 3.5.2 (R Core Team, 2021).
Finally, results obtained with the other size spectrum methods were qualitatively similar (Table S3), highlighting the robustness of our findings.

| RE SULTS
The community size spectra were highly varying among sites and showed no clear geographical pattern (Figure 1a,b). Model-averaging approach selected three main models that best explained the spatial variability of size spectrum slopes (Table 1; Figure S2). Natural factors were involved in four significant interactions with GEC (Table 1).
Specifically, we observed significant and positive interaction effects (antagonistic effects) between the percentage of non-native individuals and water temperature, and between eutrophication and water temperature (Table 1). Fish communities with higher levels of biological invasions experienced an increase in size spectrum slopes (i.e. flatter slopes) with increasing water temperature, as opposed to fish communities with lower levels of biological invasions (Table 1; Figure 2a). This result indicated that, in warmer conditions, a high proportion of large-bodied non-native individuals disrupted the negative effect of water temperature on fish size spectra (Table 1; Figure 2a). In addition, stream locations with higher levels of eutrophication also experienced flatter size spectrum slopes with increasing water temperature, as opposed to fish communities with lower levels of eutrophication (Table 1; Figure 2b).
The GEC factors also interacted with the upstream-downstream gradient in modifying fish size spectra (Table 1), with a significant and negative interaction effect (synergistic response) between the upstream-downstream gradient and eutrophication (Table 1; Figure 2c). The size spectrum slopes were generally steeper (i.e. large proportions of small-bodied fish) in upstream than downstream locations, but these slopes were even steeper in situations with higher levels of eutrophication (Table 1; Figure 2c). Finally, contrasting effects of species richness on size spectrum slopes were observed depending on the levels of eutrophication, with flatter slopes observed for species-rich communities, but only in stream locations with higher nutrient concentrations (Table 1; Figure 2d).

| DISCUSS ION
By testing for the interactions between natural and GEC drivers in shaping size spectra of stream fish communities, this study empirically F I G U R E 1 (a) Spatial distribution of the studied stream sites (n = 1095) and (b) density plot of community size spectrum slopes across France.
demonstrates that the effects of GEC are highly context-dependent and likely to modify natural community patterns. Although GEC factors did not interact with each other, we detected interaction effects between environmental factors (climate and upstreamdownstream gradient) and several GEC factors likely to reverse size spectrum responses along natural environmental gradients (Jackson et al., 2016). Specifically, our results confirm those from previous experimental studies demonstrating that the negative effects of temperature are difficult to detect in communities impacted by multiple stressors (Bouraï et al., 2020;Morris et al., 2022). Moreover, our results suggest that antagonistic interactions are more pervasive in shaping body size patterns than synergetic or additive interactions, highlighting the need to unravel the mechanistic basis of these interactive effects (Brown et al., 2013;Côté et al., 2016). Overall, these findings provide a compelling argument confirming previous studies emphasizing the need to include interaction effects in global change research (Jackson et al., 2016;Piggott et al., 2015). TA B L E 1 Averaged estimates for the fixed and interaction effects from the best three generalized linear models (GLM) explaining the variability in community size spectrum slopes. Note: Significant results (p-value < .05) are provided in the last column.

F I G U R E 2 (a-d) Significant
interactions between global environmental changes and natural drivers in predicting the size spectrum slopes. Colour lines and points represent different levels of each global environmental change from the modelaveraging approach.
Our findings give empirical evidence about how the effects of temperature on community size spectra are not systematically observed because the non-native species reshape community body size structure. This effect of biological invasions is consistent with recent size spectrum studies in freshwater fish (Arranz et al., 2021;Kopf et al., 2019) and has been explained by the high presence of large-bodied individuals in non-native fish species (Kopf et al., 2019).
Biological invasions can modify food web structure by lengthening food chains or changing trophic interactions (Britton et al., 2010;Kopf et al., 2019). Non-native species can modify global body size patterns (Blanchet et al., 2010) by alleviating temperature-size rules (Bergmann, 1847), implying those interacting mechanisms that involve temperature and non-native species may be important determinants of community size structure. In our study, the mean body mass of non-native species was significantly larger than native species ( Figure S3). Moreover, non-native species tend to have higher physiological optima than native species, thus favouring their establishment in warmer locations (i.e. thermal tolerance hypothesis ;Cucherousset et al., 2007;Kelley, 2014). Apart from the competitive and reproduction advantages of non-native species over native species, this study demonstrated for the first time a reconfiguration of the size spectrum patterns in highly invaded communities in warmer locations.
Although we did not find significant interaction effects among GEC factors in our study, they may occur in the future, as climate warming may make freshwater ecosystems more vulnerable to eutrophication by increasing nutrient inputs to ecosystems (Jeppesen et al., 2010).
Fish size spectrum slopes were generally steeper upstream than downstream, but this effect was stronger in locations with more nutrients. Thus, eutrophication mediated the effects of natural factors on the spatial variability of the community size spectra, confirming that eutrophication can be a good predictor of shifts in stream food webs (Mor et al., 2021;Roman et al., 2019). Eutrophication can affect the biological performance of fish through high stress levels (Atkinson et al., 2021;Roman et al., 2019), especially in large-bodied individuals whose relative abundances are likely to be reduced (Mor et al., 2021). In stream locations with low levels of eutrophication, our results were opposite to those of Benejam et al. (2018) (Santoul et al., 2005).
Our findings highlight the need for further research to examine the downstream-upstream gradient in stream fish size spectrum.
In conclusion, our study integrates multiple GEC to test how they modulate community size spectrum responses to natural factors.
This sheds new light on the growing evidence from other studies that GEC can exacerbate or mitigate community-level responses to perturbations (Bouraï et al., 2020;Morris et al., 2022). Although some meta-analyses (Crain et al., 2008;Wahl et al., 2011), and experimentations  have been conducted, there is still a paucity of studies aiming to quantify and capture the interaction effects between natural and GEC factors at a large spatial scale. Focusing on single GEC such as climate warming has produced highly variable results (Morris et al., 2022), and this variability may be caused by interaction effects among environmental drivers (Tylianakis et al., 2008). Research aiming at identifying the main interaction effects will be particularly important now since freshwater ecosystems are commonly threatened by multiple GEC factors occurring simultaneously (Côté et al., 2016;Jackson et al., 2016;Piggott et al., 2015). Our results outline that fish community size spectra are likely sensitive to freshwater GEC factors resulting from a re-organization of the natural patterns of community size structure. This may result in changes in the energy and nutrient fluxes at the ecosystem scales because food web structure is particularly sensitive to body-size distributions (Woodward et al., 2005). As future GEC factors will likely alter body size composition in natural communities, we emphasize that integrating spatial contingency in size spectrum research is crucial to lessen the uncertainties in global change predictions and better guide future conservation and management strategies.

ACK N OWLED G M ENTS
This work was funded by Office Français de la Biodiversité (OFB) as part of the SPECTRA project. We are indebted to the OFB for providing fish data.

CO N FLI C T O F I NTE R E S T
All authors have no conflicts of interest to disclose.

DATA AVA I L A B I L I T Y S TAT E M E N T
Data for review is available in the public repository Zenodo (https:// zenodo.org/recor d/75533 94#.Y8o-AOyZPJw).

PEER R E V I E W
The peer review history for this article is available at https://publo ns.com/publo n/10.1111/ddi.13681.