Location of the study
The research was conducted in the summer of 2023 in the waters of the Dnipro-Orilskyi Nature Reserve (Dnipro region, Ukraine) (Fig. S1). The Dnipro-Orilskyi Nature Reserve was established in 1990. It covers an area of 3,766 ha. The reserve is located within the Dnipro Reservoir, which is part of the Dnipro River cascade of reservoirs. As a result of the construction of dams on the Dnipro, the average flow velocity during spring floods at the mouth of the river into the Black Sea has decreased by 1.9 times. In the reservoirs, the flow rate decreased from 0.6–0.8 to 0.02–0.3 m/s (Timchenko 2010). The Dnipro reservoir was created after the construction of the Dnipro hydroelectric dam in 1932. After the construction of the dam, the water level of the Dnipro was raised by 1.5–2 m, which corresponds to an average level of 49.7 m above sea level. During the Second World War, in 1941, the dam was destroyed and the water level returned to its previous level. After the dam was rebuilt in 1950, the water level was restored and after the commissioning of the second unit of the Dnipro HPP in the 1960s and the construction of the Dniprodzerzhynska (Kamyanska) HPP, the water level was raised to 51.4 m above sea level. Thus, after the construction of the cascade of Dnipro reservoirs, the total rise of the Dnipro river level in the Dniprovsko-Orilskyi Reserve was 3.0–3.5 m above the natural level, which resulted in the flooding of a part of the floodplain, changes in the configuration of the banks and the area of water barriers. A sharp increase in the use of fertilisers in agriculture led to an almost tenfold increase in nutrient inputs to the Dnipro from the 1960s to the 1990s (Borysova et al. 2005; Strokal and Kroeze 2013). The water flow through the Kamianske hydroelectric station accounts for 98% of the inflow to the Dnipro reservoir, which averages 46.6 km3 per year. Precipitation on the water surface is 0.1–0.3 km3 per year. Evaporation is 0.2–0.3 km3 per year. Industrial cities consume about 3 km3 per year, of which about 1.8–2.5 km3 per year is returned to the river as domestic and industrial discharges (Denisova et al. 1989). The water level in the reservoir depends on the dynamics of water releases from the Kamyansk Dam, as the dynamics of the Zaporizhzhya Dam follows the rhythm of the Kamyansk Dam due to the small working volume of the Dnipro Reservoir. In spring and summer, the average water level in the reservoir is 51.4 ± 0.0048 m and varies from 51.1 to 51.9 m. In spring and summer, the water level in the reservoir tends to decrease. The rate of drawdown during this period is 0.7–2.7 mm per day. The average daily water level fluctuations range from 0 to 0.60 of the water level range during the season. In addition to the daily water level fluctuations, the rhythm of the reservoir water level includes oscillatory processes with a frequency of ≈ 20–25 and 40–80 days. (Zhukov et al., 2022).
Morphology of water bodies
The water bodies of the Reserve are divided into the following water body systems: water bodies of the Dnipro riverbed, water bodies of the Oril Canal (artificial formation), water bodies of the Mykolaivka ladge system, Obukhiv floodplaine system and the Taromske ladge system (Fig. S1). The database of registered water bodies of the Reserve consists of 302 objects. The total area of the floodplain and the reservoirs close to it is 24.2 km2, and the area of the reservoirs is 6.03 km2. The current area of the Reserve is considered with the projected areas that could be added to the Reserve due to their high conservation value. The water bodies were surveyed, during which the depths of the water bodies were measured using a Humminbird HELIX 9 CHIRP MEGA SI + GPS G4N echo sounder. The area of the studied water bodies was relatively evenly covered with measurement points. The coordinates of the measurement points were recorded using a GPS device. The configuration of the water bodies was recreated using detailed satellite images from the Bing Maps service (https://www.bing.com/maps), and was also refined during field surveys. The Dnipro riverbed depth map was based on Navionics SonarChart™ data (www.navionics.com). The data from the Advanced Land Observation Satellite – ALOS (http://www.eorc.jaxa.jp/ALOS/en/index.htm) were used to create the digital elevation model. The spatial resolution for this model is 12.5 meters. The DEM was resampled to a resolution of 1 m using the kriging procedure (Susetyo 2016). The digital elevation model was combined with a digital model of the reservoir bottom to create an integrated digital elevation model of the surface and underwater terrain. The morphometric characteristics of the water bodies were calculated using the lakemorpho library (Hollister and Stachelek 2017). The following morphometric characteristics of water bodies were calculated: shoreline length, surface area, shoreline development, maximum depth, mean depth, volume, maximum length, maximum and mean width, major and minor axis, axis ratio, fetch. For further analysis, the morphometric indicators were subjected to principal component analysis, which resulted in the extraction of three principal components with eigenvalues greater than unity (Zhukov et al., 2024). Principal component 1 was interpreted as an integral measure of the size of the reservoirs. Principal Component 2 was interpreted as an integral indicator of the elongation aspect of the shape of the reservoirs. Principal component 3 was interpreted as an integral measure of shoreline development. The scores of these principal components were used as predictors in the following statistical analyses.
Measurement of water properties
The dissolved oxygen content of water, temperature, electrical conductivity and transparency are important indicators of the quality of the water environment (Csábrági et al. 2019). The measurements of electrical conductivity, temperature, water clarity and dissolved oxygen content were carried out at 291 monitoring points on 5, 6 and 7 August 2023. This period corresponds to the seasonal peak in temperature when hypoxic crisis phenomena are most likely to occur in water bodies (Fedonenko et al. 2022). The measurements were made between 11 and 15 am, when the daily temperature variation was the smallest, and the measured water temperatures can be considered as close as possible to the maximum daily temperatures. The electrical conductivity of the water was measured using a Portable Water Conductivity MeterHI 993310. The dissolved oxygen content was measured using a Hanna HI 98193 Oximeter. Water transparency was determined using a Secchi disc (Secchi 1866). The depth of the Secchi disc is largely determined by optically active particulate matter, namely total suspended solids (TSS) or phytoplankton, which is usually measured using chlorophyll-a (Chl-a) (Feng et al. 2019).
Constrained community ordination and community inertia partitioning
Canonical correspondence analysis is a community ordination procedure suitable for measuring the dependence of numerous selected taxa on environmental factors. We performed a canonical correspondence analysis with water temperature as a constrain predictor. Calculations were performed using the vegan package (Oksanen et al., 2018). The optimum and tolerance of species along the axis were estimated using the scores and tolerance function. Inertia partitioning analysis was used to evaluate the independent and combined contribution of environmental variables to explain the variability in temperature estimates (Peres-Neto et al. 2006).
Plant records
Submerged, floating and emergent vegetation of all water bodies was surveyed in August 2023. The field survey procedure included the use of the tape transect method (Kolada et al. 2014). The method consists in observing aquatic vegetation along 20 m wide strips laid perpendicular to the shoreline and covering the entire vegetation zone, from the upper eulitorium to the outer limit of macrophyte growth. Observations were made in the wade and by boat using a rake and a bathyscope. At each site, the maximum depth of occupation and average plant cover were determined, and all submerged, floating and surface plants were identified. The hydrobotanical relevés were organised into a community matrix that includes 79 taxa observed at 291 monitoring stations. The projected plant cover data were subjected to an arcsine transformation:
p tr = asin((p/100%)0.5),
where ptr is the transformed data, p is the projected plant coverage (%). The community matrix was subjected to the Hellinger transformation (Legendre and Gallagher 2001).
Plant taxonomy was based on Euro + Med Plantbase (http://ww2.bgbm.org/EuroPlusMed).
Spatial distribution of plant species within the European continent
The response of aquatic macrophyte species in habitat and geographical terms was compared, and the response of species within the geographical range within Europe was determined. Information on the points of occurrence of plant species was extracted from the GBIF database (http://www.gbif.org/) using the rgbif library (Chamberlain and Boettiger 2017). Prior to developing the model, we applied spatial filtering to the presence point data in R to avoid overfitting and improve predictive performance. In particular, we applied a raster thinning technique designed to reduce duplication of presence information in individual grid cells (Bai et al. 2022). Pseudo-absence points were randomly generated across the European continent. The response of species to the temperature gradient of the warmest month of the year (bioclimatic variable 5) was investigated. Data on the bioclimatic variable were obtained using the getData function of the raster library (Hijmans 2023).
Assessment of species optimum and tolerance using the weighted average method
The species optimum is calculated using the weighted averaging method. To estimate the optimal value of a species along an environmental gradient, the weighted averaging method can be used when dealing with symmetrical bell-shaped response curves (ter Braak and Looman, 1986). The weighted averaging method is effective when samples cover the whole range of a species and the response of the species to the environmental factor can be described by a bell-shaped curve (ter Braak and Gremmen, 1987). The projected coverage of plant species can be used as a weight in the calculation of the average value of the environmental factor (Šmilauer and Lepš 2014):
$$WA=\frac{\sum _{i=1}^{n}{Env}_{i}\times {Pr}_{i}}{\sum _{i=1}^{n}{Pr}_{i}}$$
,
where Envi is the value of the environmental variable in the i-th sample, and Pri is the projected coverage of the species in the i-th sample.
The species tolerance, represented by the width of the bell-shaped curve, can be calculated as the square root of the weighted average of the squared differences between the species optimum and the actual value in the sample. The value is analogous to the standard deviation (Šmilauer and Lepš 2014):
$$SD= \sqrt{\frac{\sum _{i=1}^{n}{({Env}_{i}-WA)}^{2}\times {Pr}_{i}}{\sum _{i=1}^{n}{Pr}_{i}}}$$
.
If the entire range of the species is covered and the response of the species is symmetrical and bell-shaped, then the weighted average is correct. However, if only part of the range is covered, the estimate is biased towards the tail, which is not truncated, relative to the real value. The number of species with truncated ranges will increase if the covered part of the gradient is short, and as a result, the estimates will be biased. Therefore, it is important to cover the entire ecological gradient to obtain optimal estimates for all species (Šmilauer and Lepš 2014).
Estimation of species optimum and tolerance by regression method
Regressions can be used to model the symmetrical bell-shaped response of a species to an environmental gradient. Species-environment relationships can be modelled using generalised linear mixed models. This approach is particularly useful for unimodal symmetric species responses (Jamil and ter Braak 2013). The model can be fitted in R using the following code (Oksanen 2004):
mod <- glm(y ~ x + I(x^2), family = poisson)
b <- coef(mod)
u <- -b[2]/2/b[3]
t <- sqrt(-1/2/b[3])
where y is the species response, x is the ecological gradient, u is the optimal value of the species, and t is the species tolerance.
Huisman, Olff and Fresco hierarchical models (HOF)
The use of symmetric Gaussian response functions in gradient analysis is not a universal approach due to the systematic deviation of real data from symmetric response (Austin, 1976; Austin, 1999, 2013). Huisman, Olff and Fresco (Huisman et al. 1993) proposed hierarchical models (HOF) that include both symmetric and skewed responses. In addition to the five HOF models, two bimodal response shapes were included to account for species that are limited to gradient extremes due to competition (Jansen 2013; Jansen and Oksanen 2013; Michaelis and Diekmann 2017). The Huisman-Olff-Fresco models expanded by Jansen-Oksanen (HOFJO) are ranked according to the increasing complexity of biological information they contain (Supplement 1). The stability of the model form was evaluated using the Index of Qualitative Variation (IQV). This index takes a value of zero if all repeated runs result in the same model shape, and a value of one if all model types are chosen equally often (Mueller and Schuessler 1962). The calculation of the index is as follows:
$$IQV=\frac{1-{\sum }_{i=n}^{}{p}_{i}^{2}}{\frac{1}{n}\times (n-1)},$$
where n is the number of model types, p is the proportion for each models (Michaelis and Diekmann 2017).
The ecological niche parameters were calculated based on response models, including the species optimum, central borders, and outer borders. The species optimum represents the highest probability of species occurrence along the environmental gradient. The outer and central borders represent the rate of decrease in the response in both directions from the optimum independent of each other. These borders are determined by the distance from the optimum required for the response curve to drop a certain amount. The central borders are calculated as the gradient values at which the response reaches exp(–1/2) = 0.61 of the maximum. The outer borders of the species niche are calculated as the gradient values, where the response reaches exp(–2) = 0.14 of the top (Heegaard 2002).
Temperature ecological indicator values
There are two types of ecological scales: point scales (Diekmann 2003; Landolt 2010; Chytry et al. 2018) and range scales (Ramenskiy et al. 1956; Didukh 2011). For the first type of habitat assessment, the point of location of the species is used directly on the factor scale, which represents the ecological optimum of the species. In the second case, special formulas are used (Buzuk 2017). The Ellenberg-like scale was used for the temperature evaluation (Dengler et al. 2023). The weighted average of the indicator values of the species by their projected coverage was calculated for each community. The ideal indicator method was used to estimate the temperature according to Didukh's (Didukh 2011) range scales (Buzuk, 2017; Zhukov et al., 2018).
Phytoindication scores using range scales are converted into point scores in the process of calculations using the weighted average method based on the assumption of a symmetrical response of the species to the impact of environmental factors (Didukh 2012). For example, the data presented on the JUICE version 7.1 website on Didukh's phytoindicative scales (https://www.sci.muni.cz/botany/juice) have already been converted to a point format. In this case, all the advantages of range scales disappear, and the disadvantages of point scales are added. The main disadvantage is an increase in the shift of scores to the centre of the distribution when approaching marginal conditions. Averaging on point scales will almost never give the maximum or minimum value of a factor, and range scales will never give such estimates, even if the environmental factor is really close to its maximum or minimum value. The Buzuk ideal indicator method solves this problem (Buzuk 2017). By analogy with statistics, both approaches can be called "parametric" as they are based on the assumption of a Gaussian (normal) distribution of species in the space of factors. The ideal indicator method is "non-parametric". This approach is indifferent to the shape of the species response curve, but the properties of the "ideal indicator species" are estimated based on the regression relationship between cardinal points and tolerance. In addition to insensitivity to the properties of the species response curve (although this is only a formal drawback, like the drawbacks of non-parametric statistics compared to parametric ones), the disadvantages of the approach include the fact that it does not take into account quantitative measures of species presence in the community. The ideal indicator method can be applied only to data on the occurrence of plant species and does not take into account information on the projected coverage of plants. We have proposed an improvement to this approach (Zhukov et al., 2018), which is described below.
The information on changes in species presence or abundance along environmental gradients is used to determine species response curves. Bioindication aims to establish the value of environmental factor gradients based on quantitative characteristics of species response curves. Asymmetric species distributions are common, and their patterns of asymmetry are natural. The distribution tails typically skew towards more favourable environmental factors (Austin et al., 1990). The asymmetric pattern of species distribution requires adjustments to the procedure for determining the species optimum as an indicator value for assessing environmental properties. To model the response function of species in a gradient of environmental conditions, the β-function was proposed (Austin, 1976):
V = k × (x – a)α ×·(b – x)β,
where V is a measure of species abundance; k is a constant; a and b define the lower and upper limits of the species in the gradient of the variable x; α and γ are parameters of the distribution shape.
The location of the optimum and the skewness cannot be estimated independently in the case of beta functions, and the correlation between the response shape and the location of the optimum is an artifact of the response function. For this reason, the beta function cannot be used to estimate the asymmetry of a species' response (Oksanen, 1997), but this property allows the beta function to be used to model different response shapes of species when the actual shape is unknown. In the range of values of a and b, the β-function can demonstrate a wide variety of possible distribution shapes: from a close approximation of the Gaussian distribution to an extremely asymmetric one. The shape parameters α and γ determine three intuitive characteristics of the response function. The position of the maximum of the function, which can be interpreted as the optimum of the species in the gradient of the environmental factor, taking into account the endpoints a and b, can be determined as follows:
$$Opt=\frac{\alpha b-\beta a}{\alpha +\beta }.$$
The optimum (maximum of the function) is in the middle of the range when α = β. Asymmetry of the response, which is determined by the ratio of the parameters α and β. When α = β, the response is symmetrical, and the asymmetry increases with the difference in the values of the parameters. The kurtosis (or flatness) of the response is described by the absolute values of the parameters α and β. The parameters a and b of the beta function, which determine the smallest and largest species limits in the environmental factor gradient, fully correspond to the indicator values of the range scales. The shape of the response function of a plant species for all the factors that are indicated is unknown, as there is no information on the shape parameters α and β. However, the estimation of indicator values of species based on range scales significantly affects the results of phytoindication and depends on the peculiarities of the distribution of species responses in the gradient of the indicated factors. The view that this distribution can be approximated by Gauss' law is a significant simplification of the real pattern, so the transition from range scales to point scales by averaging marginal values is unjustified to the extent that the actual distribution of species in the factor gradient differs from the symmetrical one.
Finding out the actual distributions of species in the gradient of factor values is of considerable interest and is the task of a separate study. Most likely, species whose optimum is close to the global optimum will be characterized by a symmetrical bell-shaped distribution. Such a distribution can be described by Gauss' law or another symmetric distribution. For example, the β-function is symmetric if the shape parameters α and β are equal. These features are characteristic of species whose range indicator values are symmetrical with respect to the average value of the global range of the environmental factor. For such species, the optimum and point estimate of the indicator value will be equal to the global average value of the factor. Species whose optimum is located at the marginal position of the factor should be characterized by a distribution whose maximum is located at the corresponding limit of the environmental factor. This category of species can include species for which one of the range indicator values corresponds to the smallest or largest global value of the environmental factor. Accordingly, for such species, the optimum and point estimate of the indicator value will be equal to the global marginal value of the factor. For the purposes of simplification, we can assume that transitional variants of species optima location from the most marginal to centrally symmetric are characterized by intermediate asymmetric unimodal distributions with distribution tails directed to more favorable values of environmental factors (Austin et al., 1990; Austin, 2013).
The β-function is a more flexible distribution than the Gaussian distribution, as it can model both symmetrical and asymmetrical distributions. However, unlike the Gaussian distribution, which is characterized by only two parameters (mean and variance), the β-function has four parameters. The range indicator values quantitatively characterize the range of a species in the gradient of an environmental factor. By incorporating information on the parameters α and γ into the known range indicator scales, we can determine the nature of the species dependence in the gradient of the corresponding factor. To simplify the analytical representation of the β-function, we can impose the dependence of the unknown shape indicators α and β on the known parameters a and b, such that α = a and β = 1 - b. This results in the beta function taking the form required to solve the problem:
V = k × (x – a)a×exp5*(b−a)) × (b – x) (1−b)×exp(5*(b−a)),
where x is normalized to the range 0–1.
When a = 1 – b (meaning that the range scales are equidistant from the range boundaries), the β-function will be symmetrical (as shown by the purple line in Fig. 1). If a = 0 or b = 1, the maximum of the function (the optimum of the species) will be at either x = 0 or x = 1, respectively, resulting in a maximally asymmetric response curve of the species (as shown by the blue line in Fig. 1). For transitional values of a and b, the beta-function will be asymmetric and will gradually shift between the two extreme distributions shown in Fig. 1 (red and green lines).
The assessment of a species' environmental conditions in a given ecosystem can be made by analysing the probable response pattern of a species in a gradient of environmental factors and the projected range of the species. The response pattern is intersected by a horizontal line representing the projected coverage, and the resulting intersection points indicate the 'local' indicator values of the scale range. When the projected coverage of a species is 100%, the range scale becomes a point scale. Indicator values, traditionally used to assess environmental factors, neglect the asymmetric distribution of species response patterns. In our model, however, we consider symmetrical responses as one of the options for species response. In the central part of the gradient, both symmetrical and asymmetrical species responses to a factor are possible. However, as the gradient approaches, the probability of an asymmetrical response increases. Our model takes this into account.
The maximum projected coverage of a species indicates that the conditions in a given ecosystem are closest to optimal for that species. In contrast, a low presence in the community suggests that ecological conditions may fall within the entire range of indicator values. An increase in the specific abundance of a species narrows the range of conditions where a species with the appropriate projective coverage can be found.
The ideal indicator method was used for the phytoindicative assessment of environmental factors. This approach is based on the assumption that an ideal indicator species, with the same minimum and maximum indicator points, can most accurately indicate environmental conditions. Although such a species does not exist, the properties of an ideal indicator can be calculated by establishing the linear regression dependence of the maximum and minimum indicator values on the tolerance, which is the difference between these values. The minimum and maximum indicator values were estimated using the 'tabular' values of the indicator scales, which were adjusted to account for the projected coverage of a specific ecosystem. The constants of the linear regression dependencies of the upper and lower tolerance limits on tolerance are equal and indicate the optimum of the ideal (hypothetical) species, in which tolerance is zero. Such an ideal species can inhabit only under conditions that correspond to a given biotope, and any deviations in the conditions of existence will lead to its extinction due to its zero tolerance.