Linear and nonlinear effects of nutrient enrichments on the diversity of macrobenthos in lowland watercourses

The study concerns the relationships between taxonomic, functional and phylogenetic diversity of benthic invertebrates inhabiting watercourses and abiotic parameters associated with excessive nutrients load (concentration of Kjejdahl nitrogen, nitrates, phosphorus, organic carbon and dissolved oxygen, values of BOD5 and electrolytic conductivity). The research used data on the species composition of leeches, molluscs and larval forms of odonates and chironomid dipterans. Their description using mathematical functions allowed to determine whether the diversity reaches maximal values at extreme or moderate values of nutrients enrichment. In most cases, statistically significant relationships were unimodal—the highest diversity was observed at intermediate values of nutrients content and associated parameters, however the different patterns of relationships, monotonic and inverse quadratic, were also observed. Indirect impacts of nutrients enrichment on diversity were found as the most significant relationships. Significant responses of functional diversity were clearer and stronger than responses of taxonomic and phylogenetic diversity. The identification of fauna to the species level allowed for obtaining precise results that could enable selection of appropriate parameters for effective assessment of environmental degradation.


Introduction
Water pollution, including excessive nutrients enrichment is one of the most important global threats to freshwater biodiversity (Dudgeon et al. 2006;Quadri and Bhat 2020). Its direct and indirect effects are one of the most important factors determining ecological quality of watercourses thus their impact on the diversity of bottom invertebrates was the subject of many studies (Johnson and Hering 2009;Johnson and Angeler 2014;Camargo 2019). Since benthic invertebrates are a group with high indicative value, their usefulness in biological monitoring of organically polluted watercourses has been tested intensively for years (Wright et al. 1995;Fore et al. 1996). The extensive literature on biological monitoring of watercourses, using macro-invertebrates is therefore, for the most part, a study of the relationships between the composition and structure of faunal assemblages and the nutrients enrichment (but effects of other types of alteration have been also presented-e.g. Friberg et al. 2009;Dunbar et al. 2010). This consideration concerns also the diversity of invertebrate assemblages, but among them there is relatively little research based on the comparisons between different components of biodiversity: taxonomic, phylogenetic and functional.
The inflow of nitrogen compounds is one of the most nature-degrading anthropogenic factors (Good and Beatty 2011) and the inflow of orthophosphates is the main reason for eutrophication processes in both standing and flowing waters (Häkanson 1999). High values of nutrients load may be the result of their excessive enrichment by human activity, but both terms should be clearly separated because a significant part of the nutrients content in the watercourse is not anthropogenic (as the results of runoff from the catchment area and the natural processes of decomposition of organic matter). An increase in nutrients content above a certain level usually causes an increase in biomass and a decrease in the beta diversity of various groups of stream organisms, including macrobenthos (Bini et al. 2014;Cook et al. 2018). Such changes have been documented on numerous occasions at sites subjected to strong pressure from nutrient-rich pollutants (Wright et al. 1995;Mykrä et al. 2011;Peralta et al. 2020). The published results of previous studies on their effects on the diversity of benthic assemblages are however far from unambiguous conclusions. It has been shown both a lack of response of macrobenthos diversity to contrasting nutrients concentration (Lewis and McCutchan 2010), a decrease in macrobenthos diversity with an increase in nutrients concentration (Hering et al. 2004), an increase in diversity along the decrease in nutrients load (Camargo 2019). A decrease in the diversity of certain benthos groups, taxonomic or functional, with the simultaneous increase of others was also often stated (Cross et at2006;Fleituch 2013). The results of studies on the effects of differences in the N:P ratio on the diversity of benthic organisms are also ambiguous and indicate the dominant role of local differences in the processes of degradation (e.g. Evans-White 2009).
The causes and mechanisms of the usually observed unimodal relationship between watercourse productivity due to nutrients enrichments and species richness have been the subject of many theoretical considerations and the contribution to formulate concepts such as IPH (Intermediate Productivity Hypothesis), MPD (Multivariate Productivity-Diversity Hypothesis) and DEM (Dynamic Equilibrium Model) (Townsend et al. 1997;Cardinale et al. 2009;Huston 2014) and more. Maximum richness of benthos in watercourses observed at intermediate values of the parameters describing the inflow of nutrients was presented by, among others Adamek and Jurajda (2001), Rozenzweig (1995) and Tonkin et al. (2013) and such effects were also noted in the laboratory (Kassen et al. 2000). Contrary to these observations, the assumption of a monotonic decrease in diversity measured at the level of larger taxonomic units, usually families, along with increasing gradients of nutrients load and its indirect effects underlies i.a. procedures for the development and implementation of multimetric indices (Zhang et al. 2010) and has some empirical support (e.g. Czerniawska-Kusza 2015).
The supply of excessive loads of nutrients causes eutrophication of the watercourse, but the mechanisms of changes in the ecological niches of invertebrates in watercourses are very complex and depend on many factors. Such effects can be easily divided into clearly negative (toxic) and those that stimulate the production and diversity of certain elements of the biocenosis, but the boundaries between these two types of effects are difficult to establish. Excessive nutrients load typically results in an increase in the productivity of microorganisms and the most frequent symptoms are local changes in the values of commonly measured chemical and physical parameters: an increase in algal chlorophyll content, electrolytic conductivity, biological and chemical oxygen demand and organic carbon content, and decrease in dissolved oxygen content in water (Allan and Castillo 2007). These parameters considered indirect effects of nutrients enrichment are linked to excessive productivity of microorganisms (Baldy et al. 2007;Dunck et al. 2015). The limiting effects of abiotic parameters clearly indicated an excessive increase in the productivity of microorganisms on the biological diversity of various groups of freshwater organisms, including benthic invertebrates in watercourses have been reported (e.g . Fleituch 2013;Curtean-Bȃnȃduc 2016;Ogbuagu et al. 2019), on the other hand, some studies show no such effect (Stojanovic et al. 2017). However, excessive nutrients content in the water and parameters that can be classified as its indirect effects should be considered separately, despite the fact that they were caused by the same root causes. Attempts to find specific reasons for presented regularities are further complicated by: • Physiological and ecological specificity of analyzed groups of benthos-functional and taxonomic groups react differently to the same changes. • Analysed range of the nutrients gradient-noticing changes in diversity only in a small part of the occurring gradient can lead to erroneous generalizations (Koperski 2009). • The level of taxonomic accuracy used for the identification (Wright et al. 1995). Most often, biodiversity of invertebrates has been unfortunately determined at the family level and its estimation is often transferred to the whole of macrobenthos diversity, using questionable assumptions of the approach called taxonomic surrogacy, as it was critically discussed by, e.g. Jones (2008), Koperski (2011), and Š iling and Urbanič (2016). • Distinct differences in reactions of stream biocenoses to nutrients load between the upper and lower section (Fleituch et al. 2002) and between different abiotic types of watercourses (Koperski and Meronka 2017) • The type of diversity component taken under analysis (Koperski and Meronka 2017) Because different responses of taxonomic groups selected for the study along the gradients of environmental factors related to nutrients enrichment were expected, it was considered that, unlike most studies of this type, the most accurate possible level of identification of all individuals -species, should be used.
The main aims of the study were: • To determine at which level of nutrients enrichment and related abiotic parameters (e.g. electrolytic conductivity) biodiversity of macrobenthos (taxonomic, functional and phylogenetic) achieve the highest values. The analyzes were used to test alternative hypotheses: H1. Maximum biodiversity occurs with minimal/maximal values of the gradients of nutrients enrichment and related abiotic parameters H2. Maximum biodiversity occurs with intermediate values of these gradients • To find, to test and to describe as mathematical functions relationships between diversity of different groups of benthos, concentrations of nutrients and values of related abiotic parameters.
It has been assumed that patterns of different measures of diversity along the gradients will differ between taxonomic groups of benthos and between smaller and larger watercourses (streams and rivers). Obtaining answers to these questions using data based on invertebrates identified to the species level seems to be important for the practice of biological monitoring in watercourses, allowing the selection of appropriate measures for effective assessment of the degree of water degradation.

Samples
The results are based on selected archival and previously unpublished data obtained from the analysis of benthic samples, collected in the years 2003-2015. Before 2012, samples were collected and processed according to protocols BMWPpl (described in Fleituch et al. 2002, Galas et al. 2014, later, according to RIVECO protocols. Such archival databases appear to be a very valuable source of information on the environmental effects of water pollution. Inferring from them does not imply killing further aquatic animals for analysis. This argument is being raised more and more frequently in the contemporary literature on animal research (Koperski in press). Both sampling methods used are fully quantitative, the same samplers (hand-net, Surber sampler and Guenther sampler) are recommended to collect fauna from different types of bottom substrate (multihabitat sampling) and from the similar surface area (not-fixed in BMWPpl, typically 1-2 m 2 and 1.25 m 2 in RIVECO). Both methods produce similar results in terms of taxonomic composition and taxonomic richness. The differences in the mean abundance of samples and number of taxa found between the analyzed years (Kruskall-Wallis test, p [ 0.05) and between the sampling methods used (Mann-Whitney test, p [ 0.05) are not statistically significant. The presented study is partly meta-analysis because some of the results are newly analyzed raw data already presented in published articles by Koperski (2005Koperski ( , 2006Koperski ( , 2009Koperski ( , 2010Koperski ( , 2017Koperski ( , 2019; Karasek and Koperski (2015), Koperski and Meronka (2017). Some results based on the previous analysis are included in above articles. The analysis included data from samples collected between mid-May and mid-July between 2003 and 2015, for which it was possible to obtain accompanying data on physical and chemical parameters, carried out at an interval of not more than one month from the date of benthos sampling.
Samples collected in lowland watercourses belonging to different biotic types, classified as streams (catchment area \ 100 km 2 ) and small and mediumsized rivers (catchment area ranged between 100 and 10,000 km 2 ) were selected for analysis. It means that data from large lowland rivers (catchment area [ 10,000 km 2 ) were excluded as heterogeneous and hindering the subsequent interpretation of the results. Preliminary analyzes have shown that the division of habitats into categories below and above 100 km 2 better reflects their differences in terms of the taxonomic composition of selected groups taken into analysis than the division according to the biotic type or the type of bottom substrate. The analysis covered 127 sampling sites (63 on streams, 64 on rivers) on 91 watercourses (57 streams and 34 rivers). At some sites, the samples were taken in two or even three consecutive years. Sampling sites covered a significant part of the northern, lowland area of Poland, but were most densely distributed in the northeastern area of the country (Fig. 1).
The analysis included 241 samples of four common and ecologically important groups of benthic invertebrate: larvae of the Chironomidae (Diptera), molluscs (Mollusca), leeches (Euhirudinea) and odonate larvae (Odonata). These groups are considered less sensitive to direct organic pollutants than EPT (larval forms of Ephemeroptera, Plecoptera and Trichoptera). However, the analyzed database lacked sufficiently extensive data on EPT identified to the species level. Samples were collected by various teams during the implementation of various research projects: The remaining 14 samples were collected and prepared by the Author between 2005 and 2015 and the results of their analysis have not been published before. The list of sampling sites with geographical coordinates and the number of samples assigned to them is presented in Additional Data A. All animals were determined under an optical microscope and stereo-microscope with the maximum accuracy obtainable by methods based on morphological characters. Leeches were determined to the level of species on the basis of the key by Nesemann and Neubert (1999), molluscs were determined to the level  (2016), odonate larvae were determined to the level of species on the basis of the key by Heidemann and Seidenbusch (2002). Chironomidae larvae were identified to the highest level of taxonomic resolution, following taxonomic keys: Cranston (1982), Wiederholm (1983), Klink and Moller-Pillot (2003), and Vallenduuk and Morozova 2005.

Diversity of macrobenthos
Four, widely used indices describing biological diversity at different levels were calculated for each sample: • Rarefied species richness, typically used for nonbiased comparing taxonomic diversity in samples of different sizes was estimated with EstimateS 9 software (Colwell 2013). • Shannon diversity index.
• Taxonomic distinctness, proposed and recommended by Warwick and Clarke (1995) as more sensitive index of community perturbation than species diversity (i.e. species richness or Shannon index). Taxonomic distinctness was expressed by the Warwick and Clarke's diversity index calculated with PAST 3 software (Hammer et al. 2001). Taxonomic distance between individuals in the sample was determined using: A five-point scale for Odonata: suborder, group of families, family, genus, species (Heidemann and Seidenbusch 2002;Saux et al. 2003).
Six-point scale in the case of Chironomidae: subfamily, tribe, genus, subgenus, group of species, species (Cranston 1982).
• Functional diversity based on RAO function (quadratic entropy) being a measure of diversity of ecological communities based on the relative abundance of species and dissimilarity among them. The dissimilarity is based on a set of specified functional traits. It was calculated using software package FunctDiv (Lepš et al. 2006).
For this purpose selected biological traits were determined for each species. 15 traits were determined for Mollusca, 8 for Hirudinea, 6 for larvae of Odonata and Chironomidae (see Additional Data B for details). These groups of animals differ greatly in biological characteristics-modes of feeding, respiration and reproduction, mobility, behavior and diet composition-hence different sets of biological traits were used for each group. First of all, qualitative differences were used to determine their categories, and in the case of quantitative differences, the basis for categorization were discontinuities in the distribution of the variables.

Abiotic parameters
The results regarding chemical and physical parameters of water, listed below, were taken from the website of Head Inspectorates of Environmental Protection, Republic of Poland (www.gios.gov.pl) and the description of analytical methods and procedures used to obtain them is given there. The exceptions are the results of water analyzes, which accompanied the collection of 59 samples from 2002 to 2003, described in subsection Samples-they were obtained with Merck Spectroquant Nova 60 Photometer and Corning Checkmate equipment (Gołub and Koperski 2006).
To estimate the level of nutrients load and the intensity of related abiotic parameters, being its indirect effects at each of the sites in the time close to the moment of sampling the invertebrates, the values of seven basic parameters, measured during routine monitoring, were used. It must be emphasized that only these parameters related to nutrients enrichment, available at this website have complete data sets. To estimate the nutrients load, three parameters were used: concentration of total Kjejdahl nitrogen (TKN, mg/l), concentration of nitrates (NNitr, mg/l) and concentration of total phosphorus (TP, mg/l). To estimate the intensity of indirect effects of nutrients enrichment, the following parameters were used: DisOx (concentration of dissolved oxygen, mg/l), BOD5 (biological oxygen demand in 5 days, mg/l), TOC (concentration of total organic carbon, mg/l) and electrolytic conductivity (Conductivity, lS/cm).
Linear relationships between abiotic parameters and above four indices describing biodiversity at taxonomic, phylogenetic and functional level were analyzed using the MANOVA procedure. This procedure independently tested the statistical significance of dependent variables, independent variables and the significance of 256 relationships between their values-8 abiotic parameters as independent variables, 4 measures of diversity as dependent variables, in 4 taxonomic groups of fauna and in 2 types of environments (streams and rivers).
It was considered that the most effective approach would be looking for significant relationships between the components of biodiversity of benthos groups and the values of complex indices describing the level of nutrients load and its indirect effects in the same environment, based on the values of measured abiotic parameters. For this purpose, the values of all parameters were standardized and normalized within the type of a watercourse (separately for streams and rivers). On their basis, three complex indices were developed. They describe, using single score, estimated nutrients load (NutrientsSt), estimated indirect effects of nutrients load linked to excessive productivity (ProductivitySt) and relative content of the most important nutrients (N:P) at each sampling site: NutrientsSt is the average of the standardized values of TKNSt, NNitrSt and TPSt. The values of this index ranged between -0.644 and 4.901 in streams and between -0.969 and 2.319 in rivers. ProductivitySt is the average of the standardized values of DisOxSt, BOD5St, TOCSt, ConductSt. The values of this index ranged between -1.131 and 2.61 in streams and between -0.816 and 1.357 in rivers. N:P is the ratio of total dissolved nitrogen (TN) content to total dissolved phosphorus (TP) contentthe first value was obtained by adding Kjejdahl nitrogen concentration (TKN) and nitrate concentration (NNitr). The values of this index ranged between 1.057 and 69.00 in streams and between 1.053 and 54.146 in rivers.

Validity of patterns
Linear and nonlinear relationships between the variables assessed were considered and compared. The models of linear or nonlinear regression were calculated with PAST 3 software (Hammer et al. 2001) and developed to describe relationships between indices explaining nutrients load or those explaining its indirect effects and different components of the diversity of benthic invertebrates. For each of the 24 relationships modelled (3 indices, 4 groups of benthos, 2 types of environments), a comparison was made between the best fit linear model and the best fit nonlinear ones: logarithmic, exponential, quadratic, power or logistic. For each of them p values (probability that an index and a measure of diversity are not correlated-calculated with T test) were calculated as well as three parameters used to determine the fit accuracy of the model: extent the diversity is determined by the values of the complex index. • Residual coefficient of variation RCV-high values of RCV indicate a loose and low value indicates a tight fit of the model to the data • Akaike criterion AIC-it estimates the relative quality of statistical models and deals with the trade-off between the goodness of fit and the simplicity of the model The best fit model was considered to be one that met at least two of the three criteria: it had a higher R 2 value, had a lower AIC value, and had a lower RCV value. The relationships for which p was less than 0.05 were considered significant.

Benthic animals
A total of 22,240 individuals were identified (5340 Chironomidae, 10,502 Mollusca, 1954 odonates, 4444 leeches). Because no significant correlations were unexpectedly found during the initial analyzes between the abiotic parameters studied and the diversity of Mollusca, it was decided to remove all molluscs breathing atmospheric air (order Pulmonata) from the analysis as not informative. Thus, further analysis included only 7169 individuals of Mollusca that respirate with oxygen dissolved in water using gills-presented below as RDO Mollusca (Mollusca Respirating with Dissolved Oxygen) found in 48 samples. The total number of benthic animals included in the analysis was then 18,907, found in 241 samples (see Additional Data C).

Abiotic parameters
In streams (Table 1), the values of nitrate nitrogen concentration (ranged between 0.02 and 20.54) and N:P ratio (0.16-69) differed the most among sites, while in rivers concentrations of nitrate nitrogen (0.06-13.37) and Kjejdahl nitrogen (0.53-12.77) differed the most between sampling sites. Significant linear relationships were found between diversity components of selected groups of invertebrates and the studied abiotic parameters (Additional Data D.). Using the MANOVA procedure, it was found that the strength of such relationships in only 20 out of 256 modelled relationships was statistically significant at the level of p = 0.05 but to recognize their actual significance in accordance with requirements of the procedure, the contribution of dependent and independent variables in explaining the variability according to MANOVA must also be considered a-priori significant. Taking this condition into account, only 3 out of 256 relationships were considered statistically significant: • Functional diversity of Chironomidae were negatively related to BOD5 in rivers (R 2 = 0.394; y = -0.0264x ? 0.5206; MANOVA p = 0.0134) • Taxonomic distinctness of Chironomidae were also negatively related to BOD5 in rivers (R 2 = 0.261, y = -0.109x ? 2.831; MANOVA p = 0.0096) • Functional diversity of RDO Mollusca was negatively related to concentration of total nitrogen in streams (R 2 = 0.376; y = -0.031x ? 0.439; MANOVA p = 0.00226).

Diversity patterns
Considerably more statistically significant correlations between biodiversity and three complex indices (NutrientsSt, ProductivitySt, N:P ratio), calculated on the basis of abiotic parameters were found. Apart from 50 cases of the lack of any significant relationship (p [ 0.05, R 2 very low, RCV very high) in remaining 46 cases five types of diversity pattern along a gradients of NutrientsSt, ProductivitySt and N:P ratio were observed: • 46 significant models of regression were found among 96 analyzed relationships between Productiv-itySt, NutrientsSt, N:P and different components of diversity of four taxonomic groups of invertebrates. Most often it was an unimodal quadratic relation (29 cases), in 10 cases a monotonically decreasing linear function and in 7 cases an inverse quadratic function. The patterns were distinctly different: between streams and rivers; between taxonomic groups; between reactions to different complex indices and between different components of diversity. Relationships between complex indices and the diversity of each of the taxonomic groups in as many as 19 cases were significant for Chironomidae, for RDO Mollusca in 10 cases, Odonata in 11 cases and in Hirudinea in 6 cases.
The values of all diversity measures of all taxonomic groups in streams monotonically decreased according to a linear function with the increase in the value of NutrientsSt. In the case of Chironomidae these were statistically significant relationships (Fig. 2). In the rivers, on the contrary, the relationship patterns between different measures of diversity and NutrientsSt values were very different in distinct taxonomic groups (Fig. 3). Statistically significant dependency models were observed in Chironomidae (unimodal, quadratic function for Rarefied Richness and Shannon index, UQD pattern in the case of taxonomic distinctness and linear inverse relationship in the case of functional diversity), RDO Mollusca (UQ relation with Rarefied richness and taxonomic distinctness, MLD with functional diversity) and Odonata (MLD pattern of relationship with functional diversity).
The values of all diversity measures of Chironomidae and Odonata in streams along the gradient of ProductivitySt are consistent with the unimodal quadratic function and this relationship is highly statistically significant. However measures of the diversity of RDO Mollusca and Hirudinea decrease monotonically in accordance with linear function, and this relationship is only significant in the case of rarefied richness (RDO Mollusca) and functional diversity (RDO Mollusca and Hirudinea) (Fig. 4). In rivers, the values of all measures of diversity of all taxonomic groups depended on ProductivitySt according to the unimodal quadratic function, and in most cases it was a statistically significant relationship (Fig. 5).
Patterns of the diversity of four taxonomic groups along the N:P ratio were completely different in streams (inverse quadratic relationship: RQ or RQI, Fig. 6) and in rivers (it was always an unimodal quadratic relation: UQ or UQD, Fig. 7). In streams, dependency models built for Odonata (all measures of diversity), RDO Mollusca (Shannon index and functional diversity) and Hirudinea (Shannon index)  (Fig. 4) and the lowest recorded NutrientsSt values (Fig. 4). In rivers, maximum diversity of all analyzed benthic groups was observed at intermediate values of  (Fig. 5). In rivers, the patterns of differences in measures of diversity in the gradient of nutrient load were the most diverse (Fig. 3) and their maxima were observed at the lowest (e.g. Shannon diversity and taxonomic distinctness of Hirudinea and Odonata) or intermediate (e.g. rarefied species richness and distinctness of Chironomidae and RDO Mollusca) values of NutrientsSt. Based on the values of complex indices (NutrientsSt, Productivi-tySt, N:P ratio) corresponding to the maxima of diversity (Fig. 8), it is possible to estimate the range of basic abiotic parameters (Additional Data E). For example, values of ProductivitySt ranged between 0 and 0.22, at which the maxima of diversity of Chironomidae and Odonata were found in streams correspond, taking into account 95% confidence intervals to values 6.6-6.8 mg DOX, 4-4.5 mg BOD5, 13-14.5 mg TOC and 510-560 lS cm -1 Conductivity. In turn, values of NutrientsSt equal to -0.97, in which the maximum of functional diversity of Mollusca RDO is observed in rivers, correspond to the values of 5-15 mg TKN, 6.2-8.5 mg NNitr and 0.34-0.42 mg TP (Additional Data E).

Discussion
The obtained results are complex and their interpretation is not simple, but several generalizations can be made on their basis:  . Many more significant models among the relationships studied were found along the Productiv-itySt gradient than along the NutrientsSt gradient. It seems that differences in macrobenthos diversity, especially Odonata and Hirudinea, are more strongly determined by differences in indirect effects of nutrients enrichment than by differences in nutrients concentrations alone.
5. The influence of very high values of Productivi-tySt on the decrease in Shannon index, distinctness and functional diversity of all groups in rivers is particularly strong. 6. The increase in Odonata diversity is particularly pronounced at very high N:P ratio in streams.
The enrichment of the watercourse with nutrients causes a cascade of changes in the circulation and transfer of various forms of organic matter through the individual elements of the trophic network (Allan and Castillo 2007; Mor et al. 2019). Many models have   (Biggs 2000;Hoyle et al. 2014) and the taxonomic structure of periphyton on the diversity of stream invertebrates (Tonkin et al. 2014) has been shown. Thus, algivores and detritivores such as Mollusca and most species of Chironomidae seem to be a particularly good model for studying the effects of nutrient enrichment and its effects in the form of excessive productivity on the taxonomic structure and diversity of lotic fauna. Near the top of the trophic network in the watercourse are odonate larvae, predatory Chironomidae and leeches, which on the one hand regulate the composition and diversity of algae and detritivores, on the other hand their biomass and diversity is limited by the prey populations (Allan and Castillo 2007).
Many studies have shown that macrobenthos diversity is strongly dependent on habitat heterogeneity-a complex of ecological factors that is very  Black diamonds-Chironomidae, white circles-RDO Mollusca, grey triangles-Hirudinea, grey circles-Odonata difficult to reliably assess (e.g. Death and Winterbourne 1995;Brown 2007). This parameter usually includes information on the structure and diversity of the bottom substrate, aquatic and riparian vegetation as well as the diversity of hydrodynamic parameters. The set of environments used for the analysis certainly present a wide range of habitat variability, which was not included in the analysis as a factor because the databases used do not provide information enabling reliable assessment of this parameter. The results presented above clearly show the differences between small and larger watercourses in the patterns of diversity along the gradients of measured parameters. In watercourses with a catchment area less than 100 km 2 , a monotonic decrease with an increase in nutrient content and an increase in the symptoms of excessive productivity is more common. In watercourses with a larger catchment area, unimodal distribution of benthos diversity was more often observed, similarly to the results presented by Göthe et al. (2015).
Studies analyzing the relationships between nutrients enrichment or its indirect effects and the diversity of invertebrates in flowing waters are typically based on data in which animals are identified with low accuracy, e.g. at a family or generic level, using the concept of ''taxonomic surrogacy'' (e.g. Heino and Soininen 2007;Törnblom et al. 2011). Based on numerous analyzes (Bertrand et al. 2006;Koperski 2011;Bevilacqua et al. 2012;Heino 2014) it can be clearly stated that taxonomic surrogacy, according to which diversity at the level of larger taxonomic units brings similar and sufficient information as the diversity obtained the basis of species identification is based on incorrect assumptions and leads to erroneous conclusions. The assessment of the taxonomic diversity obtained in accordance with the taxonomic surrogacy procedure is the more unreliable the larger and more varied the number of species is included in the individual larger taxa being analyzed. Therefore, the greater explanatory value should be assigned to models regarding the relationship between nutrient enrichment and its effects on benthic diversity in flowing waters, based on data in which all invertebrates are determined with maximum accuracy (e.g. Johnson and Angeler 2014;Koperski 2017Koperski , 2019 than those based on, e.g. on family richness or Shannon index values calculated at family level (Koperski and Meronka 2017;Bis and Makulec 2014).
Gill breathing Mollusca are presented in the literature as typical positive bioindicators but some common species also occur in highly eutrophicated environments (Jakubik and Lewandowski 2014). In the paper by Jakubik et al. (2014) it was found that the species richness of molluscs in a lowland river significantly correlates negatively with the concentration of nitrogen compounds and the content of organic matter but not with the concentration of phosphorus. It seems that the very large bioindicative potential is hidden in the species-rich family Sphaeriidae, very rarely identified in publications on a biological assessment with maximum accuracy. In this work, a significant part of the species richness of RDO Mollusca in the studied environments accounts Sphaeriidae (on average 49% in streams and 47% in rivers). Therefore, the contribution of species from this family to significant differences in Mollusca diversity along the gradients of the studied parameters is important. It seems that the abandonment of detailed identification and thus not taking into account the different preferences of species and their response to gradients of abiotic parameters in almost all biotic indices significantly reduces the bioindicative value of Mollusca in biological monitoring. In results presented above, the diversity of Mollusca showed very different, statistically significant patterns along the studied gradients: both a monotonic decrease in the gradient of changes in nutrients enrichment in streams and effects of productivity in rivers, unimodal distribution in the case of a productivity gradient in rivers, and also an increase along the gradient of N:P ratio in streams.
Important role of Chironomidae in biological monitoring results from their great species richness, outstanding ecological diversity and diverse lifehistories, but their specificity results from the difficulties in identifying species. Environmental quality, including organic pollution and other ecological gradients linked with nutrients load seems to be an important factor influencing the diversity of lotic chironomid assemblages (Orendt 2018;Głowacki et al. 2011) but some authors (e.g., Rabeni and Wang 2001) suggested removing them from the protocols to make bioassessment methods more efficient, because of difficulties in taxonomic identification. Serra et al. (2017) present results clearly showing that identification of Chironomidae to the level of family or subfamily masks species sensitivities. The reduction of the taxonomical precision from species of Chironomidae to genus level lowers the statistical significance of assessed environmental preferences (Orendt 2018). Changes in the taxonomic and functional structure of the lotic assemblages of larval Chironomidae along the increase in organic pollution or nutrients concentration have been found many times (e.g. Rae 1989;Shieh et al. 2003;Koperski 2009), but there are only few convincing reports of the accompanying decrease in their diversity (Moore and Palmer 2005). A clear reduction in the diversity of this group along the increase in nutrient content, described by Koperski (2009) in small lowland watercourses is probably associated with a small range of nutrient levels in the analyzed environments-including a wider range of nutrient diversity in a larger database (Koperski 2019) indicates a typical unimodal distribution, as demonstrated by Lenat (1983). In the presented research, among all studied taxonomic groups, the diversity of Chironomidae most often showed statistically significant correlations with the studied abiotic factors. These were estimated most accurately by unimodal quadratic function except for the monotonically decreasing relationship with nutrients load in streams.
Hirudinea are presented as organisms preferring organically polluted environments (e.g. Kazancı et al. 2015) and in most biomonitoring protocols they are treated as negative bioindicators. However, the response of their diversity to pollution and the increase in nutrients load is not easy to interprete, because particular leech species distinctly differ in environmental preferences (Koperski 2005(Koperski , 2006(Koperski , 2017Kubova et al. 2013;Cichocka et al. 2015). Basic abiotic parameters seem to be clearly correlated with the taxonomic structure of lotic leech assemblages (Kubová and Schenkova 2014;Jabłońska-Barna 2017). Hirudinea being specialized (Glossiphonidae) or generalistic predators (Erpobdellidae) may show a functional response to changes in the abundance of their prey-including Chironomidae and Mollusca. Koperski (2010) presented that the diversity of Hirudinea in small lowland watercourses shows a significant positive relationship with the level of organic pollution. In presented results, among all the taxonomic groups studied, the Hirudinea diversity least often showed statistically significant relationships with the abiotic factors studied. These were mainly unimodal relationships (with ProductivitySt and N:P ratio in rivers) and only in the case of functional diversity in streams its maximum was observed at the lowest values of ProductivitySt.
A decrease in diversity of Odonata along an increase in nutrients in the watercourses has been demonstrated by Catling (2005). In results presented by Koperski (2010), Odonata diversity in lowland watercourses was much more dependent on hydromorphological parameters than on organic pollution. Honkanen et al. (2011) showed a significant impact of nutrients content on the diversity of Odonata, which, however, was clearly weaker than the presence and diversity of macrophytes and pH. In presented results, Odonata diversity relatively rarely showed statistically significant relationships with the abiotic factors studied. They were mainly unimodal, quadratic relationships with the values of ProductivitySt and inverse quadratic relationships in streams with the gradient of N: P ratio. It is worth adding that in rivers only the functional diversity reached its maximum at the lowest NutrientsSt values.
In the cases of Mollusca and Chironomidae, a monotonic decrease in diversity (MLD and UQD) was observed along with the NutrientsSt and Productivi-tySt gradients more often than in the cases of Hirudinea and Odonata-it leads to the conclusion that the first two groups actually have a higher bioindicative value in biological monitoring. Unimodal (hump-like) type of response makes it difficult to distinguish the values of diversity indices at high values of ecological quality from those at low values (De Pauw and Vanhooren 1983;Barbour et al. 1999) and consequently it makes this group not suitable for the purposes of biological assessment (Bini et al. 2014). Of the four biodiversity measures used, the most statistically significant models of regression were found for functional diversity expressed by RAO index (15) while for the remaining measures the number of significant models ranged between 10 and 11. Functional diversity seems to be an effective and promising tool for biological monitoring (Hulot et al. 2000). Its weakness seems to be arbitrariness-the values of functional diversity indices calculated for different groups of fauna are strongly dependent on the applied biological traits and scores assigned to their levels. Higher effectiveness of prediction of ecological parameters based on various benthos groups, expressed by functional or phylogenetic diversity than by traditional taxonomic indices have already been demonstrated earlier (Crozier et al. 2005;Orendt et al. 2012;Koperski and Meronka 2017).
The presented results cover only a small part of about 1700 sites designated on Polish watercourses and being a subject to routine monitoring, based on the composition of macrobenthos. However, they seem to constitute a random, representative sub-sample, and the obtained regularities can be considered typical for these types of environments. Only some of the abiotic factors, considered relevant as ecological determinants of the composition and diversity of lotic macrobenthos were used in the analysis. Routine monitoring of watercourses in Poland also covers data on many other important parameters related to nutrients enrichment: concentration of chorophyll a, concentration of suspended particles, concentration of orthophosphates or chemical oxygen demand, but the databases used in years covered by the analysis are incomplete or not available. Unfortunately, this made it impossible to use any advanced, modern and very effective method of selecting the most important abiotic parameters. Such can be considered e.g. selection of predictors from an abundant set of abiotic variables, e.g. by the BST method (Boosted Regression Trees) to create an empirical model (e.g. Species Distribution Models, SDMs) describing the relationship between the fauna of watercourses and the environment (Irving et al. 2020).
Not too much emphasis should be placed on specific values of abiotic parameters, corresponding to the maximal values of diversity of different benthic groups. They are only approximate values, they result from local, instantaneous values of abiotic parameters and they would certainly differ if regression models were based on data from other watercourses. However, it should be noted that the differences in the values of the abiotic parameters extrapolated for the diversity maxima between taxonomic groups are really significant. For example, the values of the taxonomic distinctness of RDO Mollusca and Chironomidae in rivers correspond to the values of ProductivitySt, approx. 0.2-0.4 and -0.4 to -0.2, respectively, which in turn corresponds to the extrapolated BOD5 values 1-1.5 mg O 2 /dm 3 for Chironomidae and 4-4.5 for Mollusca RDO. Such differences also suggest that when assessing the ecological status of a watercourse using, as in many cases, the diversity of benthos, we can classify completely different sites to the highest category than on the basis of physical and chemical analysis of water.
The obtained results clearly show that it makes no sense to analyze the response of the diversity of all benthic invertebrates as a one group to any environmental variables-individual macrobenthos groups, such as those presented in this paper show different or even completely different responses to the analyzed parameters. This is generally consistent with the conclusions of Johnson and Hering (2009), who observed that ''response trajectories differ between taxonomic groups and stressor, and even with stream type''. However, most of the monitoring protocols used to assess ecological status of freshwater environments so far used in practice are based on the taxonomic composition and taxonomic diversity, usually at a higher than the species level. The obtained results indicate that the departure from this practice towards a more intensive use of functional and phylogenetic diversity of selected groups of benthos, identified to the species level, will improve the assessment of the effects of excessive nutrients inflow in lowland watercourses.