Fish size spectrum as a complementary biomonitoring approach of freshwater ecosystems

Freshwater bioindicators have been developed to assess ecosystem health and responses to human-induced stressors. To date, most bioindicators primarily rely on the species identity of plant and animal communities and do not account for interactions among organisms and fluxes of energy between trophic levels. Body size is one of the most important ecological traits in aquatic ecosystems because it governs interactions among organisms and is affected by environmental conditions but, surprisingly, many of assessment approaches do not use individual body size. The community size spectrum is defined as the linear relationship between the abundance and the body size of organisms and reflects several important ecological features including ecosystem carrying capacity, predation-prey interactions, and trophic energy fluxes. In this study, we explored the potential of using the size spectrum parameters (slope, elevation and linearity) of fish community as a complementary bioindicator in 51 natural lakes and 102 reservoirs distributed across France. We determined how the fish size spectrum and other common bioindicators based on fish, macrophyte and phytoplankton communities responded to the water quality degradation, littoral habitat alterations and fish invasions. Results demonstrated that: ( i ) the size spectrum was driven by water quality degradation both in lakes and reservoirs, while size spectrum was affected by habitat loss in natural lakes and by fish invasion in reservoirs and ( ii ) the size spectrum was more sensitive to habitat loss than common bioindicators in natural lakes. This study highlights that the use of fish community size spectrum could provide additional insights into our understanding of the responses of freshwater ecosystems to global changes and could serve to improve the efficiency of management programs. This can be done at very limited additional cost because fish body size is commonly measured in biomonitoring protocols.


Introduction
Freshwater ecosystems are particularly vulnerable to multiple anthropogenic activities and are among the most threatened ecosystems on Earth (Dudgeon, 2019;Tickner et al., 2020;Vörösmarty et al., 2010).Finding quantitative, integrative, and robust metrics that summarize complex information about ecosystem health (also called "status" or "quality") is crucial to maintain their ecological integrity (O'Brien et al., 2016;Wang et al., 2021).Ecological indicators can reflect the relative long-term exposure of organisms to multiple environmental stressors (Birk et al., 2012;Dolédec and Statzner, 2010;Holt and Miller, 2011).Although bioindicators have been widely used to provide useful information on environmental pollution, limitations have been identified.For instance, most of freshwater bioindicators are based on the taxonomic composition of communities and therefore provide only limited information about biological interactions that shape community structure and ecosystem functioning (Mouchet et al., 2010;Santos et al., 2021).Therefore, bioindicators need to be reinforced by including higher levels of biological organization such as processes to provide a more comprehensive and generalizable assessment of ecosystem health (Bonilla-Valencia et al., 2022;Mouillot et al., 2013).
Body size is one of the most important ecological traits in aquatic ecosystems (Hildrew et al., 2007).Body size affects many physiological rates (e.g., respiration and growth) and determines how organisms interact with their biotic environment (Woodward et al., 2005).The empirical linear relationship between individual body size and abundance is traditionally quantified using the community size spectrum, which represents the pyramid of life from primary producers to top predators (Trebilco et al., 2013).The size spectrum integrates a myriad of bottom-up and top-down processes, including ecosystem carrying capacity, predation-prey interactions, and recruitment dynamics with three common parameters (Heneghan et al., 2019).First, the size spectrum slope indicates the efficiency of biomass transfer through body size distributions (Blanchard et al., 2009).Second, the size spectrum elevation (i.e., formally stated as the midpoint height) represents the carrying capacity of ecosystems (Murry and Farrell, 2014).Third, the size spectrum linearity describes the structure and dynamics of trophic relationships such as changes on the feeding preferences (e.g.predator-prey mass ratio, Kerr and Dickie, 2001).Taken together, these three size spectrum parameters may complement common bioindicators by including the spatial variation of the body size distributions and providing useful insights for the assessment of the ecological status of freshwater ecosystems.
Empirical studies have demonstrated that the size spectrum parameters respond to changes in environmental conditions and human impacts as increased temperature or physical habitat alteration can decrease the abundance of large-bodied individuals (Robson et al., 2005;Yvon-Durocher et al., 2011).Size spectrum elevation generally increases with increasing basal resource availability (e.g., nutrient inputs, Murry and Farrell, 2014), and deviations from linearity may arise when predator consume prey much smaller than themselves (Arranz et al., 2019).Although size spectrum studies in freshwater ecosystems have confirmed its relevance for detecting the effects of environmental stressors (Emmrich et al., 2014), the effects of multiple stressors on size spectrum variation remains undocumented (Chu et al., 2016).Multiple stressor assessment represent a general gap of knowledge, as environmental stressors may not only have a direct effect alone, but also an indirect effect through another stressors (Riseng et al., 2011;Santibáñez-Andrade et al., 2015).
The present study aims to investigate the response of lentic ecosystem health to multiple anthropogenic stressors using the size spectrum and common bioindicators.We used six water parameters quantifying agricultural or domestic effluents as a surrogate of water quality degradation, physical littoral habitat modification and biological invasion as anthropogenic stressors representing the main current stressors threatening freshwater ecosystems (Dudgeon, 2019).We first quantified the spatial variability of fish size spectrum in French natural lakes and reservoirs and identified the relationships between its natural and anthropogenic drivers.We hypothesized steeper size spectrum slopes in warm climates and artificialized littoral habitats, higher elevations with degraded water quality and a decrease of linearity with fish invasions.We then compared the fish size spectrum responses to the responses of three common bioindicators of lentic ecosystems based on three different taxa (i.e., fish, phytoplankton, and macrophyte).We then hypothesized that the fish size spectrum and the three common bioindicators would differ in their responses to anthropogenic stressors as size-based approaches can provide additional insights about community integrity than taxonomic approaches (Coccia et al., 2022).

Study area and fish data
The dataset used here was part of an extensive European Water Framework Directive survey of 153 lentic ecosystems conducted across France between 2005 and 2018 (Argillier et al., 2015(Argillier et al., , 2013)).We separated our dataset in two types of ecosystems: natural lakes (n = 51) and reservoirs (n = 102).Natural lakes were more circular (lower shoreline index) and lower connectivity than reservoirs (Appendix: Fig. A1).Fish were sampled using pelagic and benthic multi-mesh gillnets following the CEN standardized procedure (Argillier et al., 2013;CEN, 2005).Multi-mesh gillnets were set overnight, and the number of multi-mesh gillnets changed (from 3 to 81 nets) according to the lake area and depth (CEN, 2005).A total of 335,542 fish were sampled, identified to species, measured and weighed (wet mass) individually or grouped in batches (in total, 95,790) of approximately 25 individuals of the same species with similar body size.For individuals grouped in batches, we randomly assigned body size length between the minimum and maximum body sizes of the batch.When fish mass was not measured, we estimated fish body mass using a species length-weight relationship calculated from the observed data at national level.

Species assemblage variables
We calculated two integrative variables reflecting the structure of the fish species assemblage across natural lakes and reservoirs.First, community composition was defined using fish species assemblage (percentage (%) of the most common fish families: salmonids and cyprinids) and species richness as the total number of species in each lake.The three variables were reduced to a Principal Component Analysis (PCA) in the ade4 package in R (Dray and Dufour, 2007).We only used one synthetic variable (1st PCA axis, 57.4 %) showing clear differences in community composition: positive values indicate fish communities characterized by salmonid species (e.g., brown trout Salmo trutta, lake char Salvelinus umbla or rainbow trout Oncorhynchus mykiss), and low species diversity whereas negative values were characterized by high percentage of cyprinid species (e.g., roach Rutilus rutilus, bream Abramis brama or common carp Cyprinus carpio) and high species diversity (Appendix: Fig. A2a).Second, the fish invasion variable was quantified using a synthetic variable (1st PCA axis, 76.8 %) based on three different measurements of fish invasions (percentage of non-native species, percentage of non-native individuals and percentage of non-native biomass; Catford et al., 2012; Appendix: Fig. A2b): positive values indicate high levels of invasion whereas negative values indicate low levels of invasion.Species status (native and non-native) were defined following Keith et al. (2011) at national scale.

Abiotic variables
Natural lakes and reservoirs were based on their climatic condition, water quality degradation, and littoral habitat modification.First, climatic conditions were defined by the mean of 5 bioclimatic variables related to air temperature recorded by Météo France during the 5 years before fish sampling: (i) mean annual air temperature, (ii) air temperature seasonality, (iii) maximum air temperature of the warmest month, and (iv) minimum air temperature of the coldest month and (v) temperature annual range.We transformed air temperature to water temperature following Punzet et al. (2012) for temperate freshwater ecosystems.Then, we used a PCA and a synthetic variable (1st PCA axis, 58.9 %) summarizing ecosystems with "warm" climates for positive values and ecosystems with "cold" climates for negative values (Appendix: Fig. A2c).Second, water quality degradation was defined from the chemical data measured during the fish sampling season (available at https://www.naiades.eaufrance.fr).We compiled six parameters (Total phosphorus, total nitrogen, dissolved carbon, turbidity, Secchi depth and suspensive organic matter content), and performed a PCA to obtain a synthetic variable (1st PCA axis, 63.9 %) representing gradients of water quality: positive values indicate highly degraded water quality while negative values indicate minor water quality degradation (Appendix: Fig. A2d).Third, we measured an index of littoral habitat alteration by calculating the percentage of non-natural area in a 25meter buffer around the littoral corresponding to the maximum resolution of the 2012 Corine Land Cover (Büttner, 2014).The index of littoral habitat alteration ranged from 0 to 1 and a value of 1 would indicate to a total of littoral habitat artificialization (hereafter "Habitat loss").

Bioindicators
Lentic bioindicators were obtained from the Naïades database (available at https://www.naiades.eaufrance.fr) and using only values collected during the same period of the fish survey for each lentic ecosystem (maximum 5 years).Two different bioindicators were used for fish community according to ecosystem type: Indice Ichtyofaunique Lacustre (IIL) for natural lakes (Argillier et al., 2013) and Indice Ichtyofaunique pour les Retenues (IIR) for reservoirs (Argillier et al., 2015).The ILL index is based on total fish abundances and the biomass of omnivorous species with the aim to detect impacts from land cover changes and eutrophication (Argillier et al., 2013).The ILL was calibrated at the European scale using 455 natural lakes to allow comparisons between different Water Framework Directive (WFD) member states for the management of eutrophication (Argillier et al., 2013).The IIR has been developed following the same purposes (response to land cover and eutrophication level) but was calibrated on 94 French and 20 Czech reservoirs (Argillier et al., 2015).Similarly than IIL, the IIR is a multimeric index based on total fish abundances (total biomass) and feeding preference of fish (number of planktivorous and of piscivorous individuals) stated at the species level (Argillier et al., 2015).Then, we used the Indice Phytoplanctonique LACustre (IPLAC) as indicator of phytoplankton communities.IPLAC was based on species composition and total biomass metrics to reflect the trophic status of lentic ecosystems (Laplace-Treyture and Feret, 2016).Finally, we used the Indice Biologique Macrophytique Lacustre (IBML), an indicator based on macrophyte communities.The IBML is a single-metric indicator based on the sensivity/tolerance of macrophyte taxa (and their relative abundance) to high trophic level (Boutry et al., 2015).The IPLAC and IBML indicators were developed and calibrated in French lakes and reservoirs.

Fish size spectrum parameters
The fish size spectra were calculated using the normalized abundance size spectrum from individual fish body mass.We chose a binned method following other size spectrum studies in lentic ecosystems (Emmrich et al., 2014;Arranz et al., 2019).Moreover, using a linear binned approach allowed to obtain the full descriptors of size spectrum parameters (slope, elevation and linearity) using a single statistical approach.Although binning methods may perform less than maximum likelihood estimation (MLE) in determining exact slope value, normalized size spectrum is among the methods that most covary with the MLE method (Edwards et al., 2017).For size spectrum calculations, we initially classified the body mass distribution of individual fish weight into fourteen size classes starting at 1 g and following a geometric series of two.Because of the underestimation of individuals < 8 cm (equivalent to a mass of 5 g) in multimesh gillnets (CEN, 2005), we remove all individuals < 4 g (the two first classes) from size spectrum calculations to avoid any biases due to sampling inefficiency (Mehner et al., 2016;Sutton and Jones, 2020).The maximum number of size classes was defined to maintain a homogenous representativeness of the size classes across all locations (Martínez et al., 2016, Appendix: Table A1).At the end, eight size classes were used ranging from 4-8 g to 512-1024 g.Then, we calculated one size spectrum per location as the ordinary-leastsquare linear regression between the log 2 of each size class midpoint and the log 2 of the total number of individuals per size class, normalized by bin width.Size spectrum slope was extracted from the linear regression equation, size spectrum elevation was stated as the midpoint height on the y-axis after scaling as a centered intercept to avoid high correlation with slope (Sprules and Barth, 2016), and size spectrum linearity was defined as the determination coefficient (R 2 ) of the linear regression (Sprules and Barth, 2016).Details about methodology for calculating the linear size spectrum (R script) are given in the Appendix.Thus, we calculated the MLE in each lake fish community using the R package sizeSpectra (Edwards et al., 2017) and the obtained results through the MLE slope were consistent to the binning size spectrum slopes (Appendix: Table A2).

Environmental drivers of the size spectrum and bioindicators
We used Structural Equation Models (SEM) to create the causal model between each predictor of the size spectrum and main drivers by using lavaan R package (Rosseel, 2012).SEMs allowed us to test whether the potential effects of all the stressors on the size spectrum in a single statistical analysis.SEMs are composed of a set of linear regressions among inter-correlated variables of several hierarchical levels and specify an estimate of direct, indirect, and total effects (direct + indirect) and are particularly relevant for addressing causality relationships in ecological structure (Pugesek et al., 2003).SEMs were initiated by constructing an a priori causal model among all potential predictors, which we potentially adjust to maximize the fit indicator values (Capmourteres and Anand, 2016;Santibáñez-Andrade et al., 2015).We set the climate as an exogenous variable that could influence the water quality degradation or fish invasions.We based our model adjustment on 4 different SEM fit indicators: GFI, CFI, χ 2 and RMSEA (Fan et al., 2016).Habitat loss and fish invasion variables were log-transformed to fit graphically a normal distribution.According to our second hypothesis, we also tested the response of bioindicators in the final model and then extracted the normalized estimates of the total effect (direct + indirect) of water quality degradation, habitat loss and fish invasions.At the end, the structure of the final model was appropriate for our dataset with modeling indices indicating good model adjustments, even for size spectrum parameters and bioindicators except for IBML model in reservoirs where only three of the four indicators were in accordance to the threshold of good adjustment (Appendix: Table A3; Table A4).Finally, we compared their effect size values to evaluate significative differences between bioindicators and size spectrum.
Climate had a stronger effects in natural lakes than in reservoirs, both directly on the size spectrum parameters and indirectly through other stressors (Fig. 2).For example, warm climates were associated to higher water quality degradation and higher fish invasion levels (Fig. 2).By contrast, warm climates had a negative effect on community composition: more cyprinid species were caught in warmer locations.Size spectrum slopes in natural lakes were directly and negatively affected by climate and habitat loss: fish communities in warmer climatic conditions and more artificialized banks displayed steeper size spectrum slopes (Fig. 2a).On the other hand, poor water quality steeped size spectrum slope in reservoirs with a higher proportion of smallbodied individuals (Fig. 2b).Size spectrum elevation and linearity were directly driven by climate and water quality in natural lakes (Fig. 2c and e): warm climate favor higher linearity and higher carrying capacity whereas water quality degradation had a negative effect on these two size spectrum parameters.None of the size spectrum parameters were associated to community composition in natural lakes (Fig. 2a, c and e) but community composition was the unique driver of both size spectrum elevation and linearity in reservoirs: assemblages dominated by salmonids displayed lower elevation and size spectrum linearity than reservoirs containing assemblages mainly composed of cyprinids (Fig. 2d and f).

Total effects of stressors on size spectrum parameters and bioindicators
We observed a strong correlation between bioindicators in natural lakes and reservoirs (Appendix: Table A5).Similarly, the size spectrum parameters were highly correlated with each other.The size spectrum elevation was negatively correlated with the fish indicators in lakes and reservoirs (Appendix: Table A5) but in both ecosystem, elevation and fish indicators did not response to the same drivers (Fig. 3), suggesting that the covariation occurring between them was not related to the stressors tested.By contrast, the size spectrum slope was not correlated with the fish indicator or phytoplankton index in the reservoirs but they responded consistently with the negative effect of water quality degradation (Appendix: Table A5, Fig. 3).We also observed a common negative response to quality degradation in the size spectrum elevation and IPLAC in natural lakes (Fig. 3).The effect of water quality was direct on elevation while IPLAC responded indirectly through changes in fish species assemblages (Fig. 3; Appendix: Table A6).Although we did not observe a direct effect of habitat loss on the linearity of the size spectrum in natural lakes (Fig. 2e), we observed a total effect of habitat loss on linearity when indirect effects were summed (Fig. 3).Contrary to size spectrum slope and linearity in natural lakes, there was not significant association between habitat loss and bioindicators, both in natural lakes and reservoirs (Fig. 3).Finally, IIR and size spectrum linearity responded to fish invasion level in reservoirs.Direct effect was significant for IIR but the combination of indirect effect (via community composition) was required to detect an effect of fish invasion on linearity (Fig. 3; Appendix: Table A6).

Discussion
The present study is the first directly comparing community fish size spectrum to common bioindicators used in monitoring programs and demonstrates the added-value of fish size spectrum as a complementary indicator to assess the consequences of anthropogenic stressors in lentic ecosystems.As expected with the first hypothesis, we found a high spatial variability in size spectrum parameters that was driven by changes in environmental conditions in natural lakes and reservoirs.The response of the size spectrum to the three perturbations tested differed between ecosystems: the size spectrum parameters were able to detect the loss of littoral habitat and water quality degradation in natural lakes, while water quality and fish invasion were the most important drivers for assessing changes in size spectrum parameters in reservoirs.Regarding the second objective, we found that fish size spectra were sensitive to the changes in environmental conditions and these changes differed from those captured by common bioindicators: fish size spectrum responded to littoral habitat loss in natural lakes while common bioindicators were unmodified.Our results suggest that the inclusions of body size distribution can provide complementary information to existing biomonitoring of lentic ecosystems for assessing hydromorphological alterations.
Littoral habitat is a key structuring factor of aquatic communities as it provides refuge for juveniles, additional spawning habitats and resource stability for predators (MacRae and Jackson, 2001;Tokeshi and Arakaki, 2012).In our study, we observed a steeper size spectrum slope in artificialized littoral habitat in natural lakes, implying a reduction of large-bodied individuals proportion when decreasing habitat complexity.Moreover, we observed an increased linearity when increasing habitat loss across lentic ecosystems.In coral reefs, Rogers et al. (2014) found that fish communities had higher linearities in habitat with less structural complexity, partly due to the decrease in the number of small-bodied fish due to increased predation pressure.This response to habitat loss was not observed in the reservoirs and could be explained by higher water level fluctuations (e.g.irrigation, hydropower) that may naturally drive fish behavior to reduce their use of littoral habitat (Hirsch et al., 2017;Poikane et al., 2020).Thus, studying the body size distribution attests to the importance of littoral habitat for the structure of aquatic communities and would be a relevant tool to quantify directly the effectiveness of restoration measures on fish communities in natural lakes (Ritterbusch et al., 2022).
There is an urgent need to rapidly improve the methodological approaches of ecological indicators to cover different facets of community integrity in freshwater ecosystems (Wang et al., 2021).An ecological indicator has to (i) condense information about communities for rapid ecological assessment (ii) be low cost effective (iii) be easy to understand and handle for managers and (iv) be responsive independently of study scale for intercalibrate procedures (Pinna et al., 2013;Poikane et al., 2015;Pont et al., 2007).Although French bioindicators were designed to detect eutrophication level in lentic ecosystems, the fish-and macrophyte-indicators did not respond to water quality in natural lakes.This result can be explained by two main reasons: (i) the variables used to reflect trophic state in the development of bioindicators were different (generally total phosphorus concentration alone) and (ii) calibrations were performed on a wider spatial scale than the present study.For example, the IBML responded to trophic state on lakes and reservoirs combined (Boutry et al., 2015) and the IIL was calibrated on 455 natural lakes across Europe (including 40 French lakes) covering a wider eutrophic gradient than our study locations (Argillier et al., 2013).Size diversity vary less than species composition at large geographic scale but may vary more at the local scale, making it a relevant candidate metric to assess anthropogenic stressors (Basset et al., 2004;Solimini et al., 2009;Xu et al., 2016).Therefore, body size approach appeared as a pertinent tool for comparison across sites, especially to compare fish communities with heterogenous species richness and composition (Logez and Pont, 2011).
Although length-and age-based metrics have been described as relevant tools for developing ecological indicators in the Water Framework Directive (CEC, 2000), individual body size is only used for calculating total fish biomass in French bioindicators (i.e.BPUE; Argillier et al., 2013), underestimating all the ecological insights that could be obtain from this integrative trait.Other European countries already use individual body size for assessment purposes, usually through measurement of deviations from the reference sizes of the most abundant species (Ritterbusch et al., 2022).In Estonia, for example, size diversity is used directly as a metric to assess urbanization pressure (Ritterbusch et al., 2017).Individual body size is commonly measured in fish community inventory and therefore, size spectrum approach can be used, at no additional costs, in most ongoing biomonitoring program with an important added-values.Size spectrum could also be integrated for other taxa with size not systematically measured, such as planktonic and benthic invertebrate communities.Indeed, current development of high-throughput imaging techniques for measuring organisms size ( Álvarez et al., 2011;Kitahashi et al., 2018;Vandromme et al., 2012) will allow to obtain size spectrum across plant and animals communities (Petchey and Belgrano, 2010).
In conclusion, our study represents a quantitative and integrative approach demonstrating the added value of size spectrum to assess anthropogenic stressors in lentic ecosystems.Overall, our work highlights that size spectrum in biomonitoring is reliable, time-and costeffective and complementary of bioindicators.Integrating size-based approach in biomonitoring should be systematic to improve our understanding of ecosystem response to multiple stressors.

Fig. 1 .
Fig. 1.Spatial distribution of size spectrum slope (a), size spectrum elevation (b) and size spectrum linearity (c) in natural lakes (circles) and reservoirs (squares) and density plots for size spectrum slope (d), size spectrum elevation (e) and size spectrum linearity (f) in natural lakes (solid line) and reservoirs (dashed line).

Fig. 2 .
Fig. 2. Structural equation models showing the direct and indirect effects of the potential drivers on size spectrum slope (a and b), size spectrum elevation (c and d) and size spectrum linearity (e and f) in natural lakes (left) and reservoirs (right).The thickness of the arrows indicates the magnitude of the significative path coefficients (effect size) and the color indicates the direction of the response (positive in blue or negative in red).Non-significant paths (p-value ≥ 0.05) are represented by a dashed line.(For interpretation of the references to color in this figure legend, the reader is referred to the web version of this article.)