Macrophyte communities as bioindicator of stormwater pollution in rivers: a quantitative analysis

Study site: Psel and Bystritsa Rivers

Samples were collected on the Psel River, a left tributary of the Dnieper River (Black Sea basin) and the Bystritsa River, which flows from the left bank into the Wieprz River, a right tributary of the Vistula River-Baltic Sea basin (Fig. 1).

Both rivers are of the plain type, localized on the same geographical parallel, which suggests similar hydrological and climatic conditions. The rivers are more than 1,000 km apart and belong to different basins, which makes them independent research objects. The rivers were surveyed in sections flowing through the cities of Sumy (Ukraine) and Lublin (Poland) with approximately the same population. The studied sections on both rivers are straightened, deepened, the banks are partially turned into dams. In the studied areas, runoff from urban stormwater systems flows into the rivers. Both rivers are under the influence of channel reservoirs—a typical factor of disturbance of the natural hydrological regime of watercourses. The section on Psel river (Sumy) is located above the dam, and the section of the Bystritsa (Lublin) directly after the reservoir. This determined a significant difference in the flow velocity on the studied river sections. In the Psel River, the average flow velocity did not exceed 0.2 m/s, and in the Bystritsa River the flow velocity was in the range of 0.6–1.2 m/s. The low flow velocity in the Psel was determined by the dam located below, and in the Bystritsa high, not typical for plain rivers flow velocity was determined by water discharge from the dam and straightened canalized channel.

Sampling method

Macrophytes are a convenient object of study, but their quantitative assessment creates difficulties and does not have standard generally accepted methods and approaches. An important factor in choosing the survey length is the total length of the river or body of water that will give a representative sample. The duration of the survey is directly related to the spatial frequency of the surveys (number of surveys per water body on a given date), since the length of the survey will generally decrease as the frequency of surveys increases (Directive 2000/60/EC, 2000). The most commonly used reach length is 100 m, described by Holmes (1991). Often 100-m is simply deemed a practical size of reach to survey. For example, national environmental agencies use a length of 500 m as recommended in the Guide to Common River Monitoring Standards (JNCC, 2005). A simple takeaway from the range of approaches available is that one size cannot fit all purposes. (Pentecost, Wilby & Pitt, 2009). In our studies, we proceeded from the fact that stormwater wobbling is most pronounced 10–20 m below the runoff. For this reason, we studied sections of 10 m. Below we are described in detail the methodology which we used.

Data on the projective cover of 20 aquatic plant species were used for the analysis. In total, data from 24 stations were analyzed, six of which were located at the discharge sites of stormwater systems. The material was collected both in places where stormwater systems flowed in (before and after the flow), and in areas where there was no influence of stormwater systems.

To determine the projective cover of macrophyte species at each station, the following methodology was used. The pattern of macrophyte patches was established by measuring their boundaries along transects spaced one meter apart using a tape measure. The results of the measurements were plotted on prepared templates (Fig. 2). Connecting the points on the template gives a fairly accurate macrophyte overgrowth pattern, allowing the proportion of overgrowth by species to be assessed. An area of 100 m² was measured at each station (if the river width allowed).

The layout of macrophytes was plotted on the template. In the case of macrophytes with floating leaves, such as the yellow water-lily or waterlilies, the spot occupied by leaves was marked and the number of leaves in the spot was counted. We measured the area of 11 leaves to obtain the average area per leaf of a given species (of plant). To determine the projective area of the species, the average area was multiplied by the number of leaves within the study area. Separately, the coverage area of underwater macrophytes was estimated by plotting the configuration of their spots on a template. Single plants and small patches of vegetation were evaluated using a special floating frame divided by a wire mesh into squares of 1 dm².

Coverage may have exceeded 100%, given the fact that plants with floating leaves constituted the second (upper) stratum and those submerged underneath constituted the first (lower) stratum.

The scientific plant names have been used according to Tutin et al. (1968–1980) and Dobrochaeva, Kotov & Prokudin Yu (1987).

Statistical analyses

Ecological studies often use binary data (presence/absence) of species. Results based on abundance and presence/absence data often differ. Although we had abundance data we decided to consider the results of applying Jaccard distances to the converted to binary data. The data were processed using R version 4.0.2 (R Core Team, 2020). Hierarchical clustering was performed with hclust (function ward.D2) from core R package stats. The key to interpreting a results of hierarchical cluster analysis is to look at the point at which any given pair of cards “join together” in the tree diagram. Cards that join together sooner are more similar to each other than those that join together later (Albert & Tullis, 2008). Two types of p-values (AU and BP) for the nodes of hierarchical clustering dendrograms were calculated with R package pvclust. According to the package documentation, “AU p-value, which is computed by multiscale bootstrap resampling, is a better approximation to unbiased p-value than BP value computed by means of normal bootstrap resampling. Clusters with AU larger than 95% are highlighted by rectangles, which are strongly supported by data” (Suzuki, Terada & Shimodaira, 2019). Non-metric multidimensional scaling (NMDS) was performed with the metaMDS function from vegan (Oksanen et al., 2019). The NMDS method (Kruskal, 1964; Borg & Groenen, 2005) made it possible to obtain sets of points distributed in a 2-dimensional space in such a way that in the analyzed cases, the similar points were closest to each other, while different ones were spread far apart (Borg & Lingoes, 1987). Plots were produced using the ggplot2 (Wickham, 2016) and ggrepel packages (Slowikowski, 2020).

The Shannon–Wiener and Simpson’s indices were used to compare stations. Ordination in reduced space is often used in ecological studies so the authors considered results of nonmetric multidimensional scaling (NMDS) with different distance measures: Jaccard, Kulchinsky and Manhattan. In addition, the data were either preliminarily Hellinger transformed or not.

Species composition

The species composition of macrophytes and average values of projective cover at stations both experiencing and not experiencing the influence of the stormwater system in the Bystritsa and Psel rivers are shown in Table 1.

Projective cover

The projective cover of dominant species in the points affected by runoff on the Psel River reaches from 50% to 105%. At the same time, the share of individual dominant species varies from 15% to 35%. The total projective cover at these stations ranged from 85 to 120%.

Establishing the projective coverage of macrophytes is described in the ‘Materials and Methods’ section.

Projective cover of dominant species on the Psel River in areas not affected by runoff from stormwater systems ranged from 15% to 70%, with total cover remaining in the range of 39% to 72%. Participation of individual dominant species in these sites ranged from 10% to 30%.

The total projective cover in the Bystritsa River, which in the studied area is much narrower than the Psel River, ranged from 100% to 160% at stations affected by stormwater, while the projective cover of dominant species at these sites ranged from 66% to 117%.

In sections of the Bystritsa River not directly influenced by runoff, the total projective cover reached a maximum of 125% and a minimum of 95%. At the same time, on these sites dominant species could occupy from 31% to 73%. In areas influenced by runoff, dominant species occupied 66% to 117%, with total cover ranging from 107% to 160%.

Similarity

The similarity of species composition of macrophytes at stations influenced by stormwater and stations not influenced by them in the Psel River was 80% according to the Sørensen index, and in the Bystritsa River—85%. At the stations in the Psel River influenced by runoff, 11 macrophyte species were found, while 14 species were found at stations not influenced by runoff; at the Bystritsa River these numbers were 13 and 18 macrophyte species, respectively.

Table 2 presents species richness and biodiversity indices for non-stop discharge station and storm discharge stations.

Indicative groups

The proximity of all studied stations to each other was assessed using hierarchical clustering according to the Ward’s method. The data were previously Hellinger-transformed; for clustering, the Kulczynsky distance matrix was used (Fig. 3).

The hierarchical clustering showed that all stations with high reliability—as evidenced by the p-values based on multiscale bootstrap resampling (see Fig. 3)—split into three clusters. The stations of both rivers affected by stormwater runoff form a separate cluster. At the same time, stations not affected by runoff are distributed in the remaining two clusters according to the rivers. The fact that all stations of both rivers affected by stormwater systems were combined into one cluster suggests that the impact of runoff is the decisive factor and, according to the cluster analysis, the impact of the stormwater system was more significant than river affiliation. Indeed, if the influence of stormwater is not decisive, one should expect the stations to be divided into two groups, according to their localisation on different rivers. However, stormwater runoff was found to even out the differences between rivers, resulting in similar conditions, as reflected in the cluster analysis.

Figure 4 shows the NMDS results with Jaccard distances. No matter how many groups are selected—two, three or four—95% confidence ellipses overlap more or less (Fig. 4). When distinguishing two groups (Fig. 4A): stations of both rivers without the influence of the storm system and stations of both rivers in the area of runoff from the storm system. All stations of the Bystritsa river fall within the confidence ellipse for a group of stations in the area of runoffs. When divided into three groups: the stations unaffected by stormwater runoff on the river Psel, stations on the river Bystritsa and the stations of both rivers in the runoffs regions show that the NMDS of the matrix with Jaccard distances does not provide an opportunity to separate these groups.

Even more clearly, the impossibility of distinguishing the groups of stations based on Jaccard distances is seen in Fig. 2C when setting four groups: two for the stations of each river under the influence of storm system and without. Without exception, all the studied stations: on the river Psel, and river Bystritsa fall into a single 95% confidence ellipse for the Psel River group (Fig. 4C). Thus, the NMDS analysis based on qualitative data demonstrates only the commonality of the studied objects and does not allow speaking about reliable differences between them.

All further results were obtained using the quantitative data on the participation of individual species in the projective cover of the water body at the study stations. The division of stations into groups was the same as described above, so further on the authors will simply consider the two, three or four groups without detailing.

First, it was checked how the untransformed data work using the Kulczynski distance. The NMDS results are shown in Fig. 5. Regardless of the preassigned number of groups, the results are quite similar. When two or three groups are assigned, only the stations of both rivers under the influence of stormwater and stations on both rivers without the influence of stormwater are significantly different. An attempt to distinguish four groups (Fig. 5C) is also unsuccessful, since the points of stations of both rivers, which are under the influence of runoffs, lie mixed up on the NMDS plain and do not form two groups. Confidence ellipses for the stations without stormwater effect overlap. Thus, the NMDS, on the basis of untransformed data with the Kulczynski distances, prompts the conclusion that the objects of two types are being dealt with: rivers without the influence of stormwater systems and rivers in the area where stormwater runoff flows. This gives grounds to assert that the influence of runoffs is significant, leading to the destruction of the individuality of the river.

As it can be seen, using the Kulchinsky distances on the raw data deepened the understanding, compared to using the Jaccard distances.

It is known that Hellinger transformation of a data (Legendre & Gallagher, 2001) often contributes to better ordination results. Therefore, it was used before calculating the Kulczynski distances (Fig. 6). Here, the reliability of different groups is the same as without transformation, but the relative position of stations on the diagram is different. The position of the stations under the influence of stormwater system is especially different (Fig. 5). Without transformation, these stations of both rivers are mixed up and positioned much closer to the stations of the Psel River than the Bystritsa River. After transformation according to Helenger, the unaffected stations of both rivers are more close with their own river’s stations affected by stormwater runoffs. This allows refining the interpretation of the results. It can be stated that runoff locally transforms the composition and share of different aquatic vegetation species in the projective cover of each river, but some characteristic features of each hydro-ecosystem remain. It should be noted that when using the Jaccard distance, the stations affected by runoff did not stand out as a separate group, but in all cases had a distinct tendency to be near the stations of the corresponding river (Fig. 4).

The effect of Helenger transformation on NMDS results can also be illustrated by using Manhattan distance. Thus, without the Helenger transformation, the predefined groups on the plot overlap and do not provide a basis for a confident interpretation (Figs. 7A, 7B). After transformation according to Helenger, in both cases, a group is distinguished that unites stations affected by the runoff, and the unaffected stations from different rivers are also quite confidently separated. This coincides with the previously obtained subdivision of stations by the clustering method.

Analysis of the positions of stations on the graphs and associated plant species on the basis of hirearchical clustering (Fig. 8, right cluster) and NMDS results makes it possible to identify a group of species associated with the influence of runoff from stormwater systems. In the area of wastewater discharge from the stormwater system on the Psel River the dominated species were: Nuphar lutea, Ceratophyllum demersum, Myriophyllum spicatum, Typha latifolia, and on the Bystritsa River were: Glyceria maxima, Typha latifolia, Stuckenia pectinata and Potamogeton crispus.