Mayfly response to different stress types in small and mid-sized lowland rivers

Abstract Freshwater ecosystems are endangered worldwide by various human pressures, resulting in dramatic habitat and species loss. Many aquatic invertebrates respond to disturbances in their habitat, and mayflies are among the most sensitive ones. Therefore, we investigated mayfly response to anthropogenic disturbances at 46 study sites encompassing slightly to heavily modified small and mid-sized lowland streams and rivers. Mayfly nymphs were sampled between April and September 2016 using a benthos hand net. A total of 21 species was recorded, with Cloeon dipterum (Linnaeus, 1761) being the most frequently recorded one. Nevertheless, the taxa richness was rather low per site, i.e., between zero and nine. Assemblage structure had a high share of lower reaches and lentic (potamic and littoral) elements, and detritivores (gatherers/collectors and active filter feeders). This indicates that hydromorphological alterations lead to assemblage “potamisation” in small and mid-sized rivers. More mayfly species were related to higher oxygen concentration and lower water temperature, abundance of aquatic vegetation and total organic carbon. Additionally, the assemblage diversity and abundance were negatively associated with increasing intensive agriculture area at the catchment scale. This study confirms mayfly bio-indicative properties, i.e., their sensitivity to alterations of their habitat and pollution, but also provides new data related to mayfly response to the impacted environment. Those data can be used for management and protection activities of lowland rivers and their biota according to the requirements of the European Water Framework Directive.


introduction
Freshwater ecosystems represent an indispensable resource of water supplies for humans (Carpenter et al. 2011), but they also have a crucial role in biodiversity maintenance and conservation (Previšić et al. 2009;Ivković and Plant 2015). Therefore, it is essential they remain in good ecological status (Dudgeon et al. 2006;Vörösmarty et al. 2010). Nevertheless, the status of many aquatic systems is far from good worldwide (Carpenter et al. 2011). Various anthropogenic impacts represent major threats to aquatic biodiversity and make lotic habitats among the most endangered ones (Malmqvist and Rundle 2002;Hering et al. 2006;Stoddard et al. 2006). Human population growth, increased urbanisation and industrialisation have led to increased demands for land use for purposes of agriculture, forestry, irrigation activities and wetland drainage, resulting in alterations of habitat morphology, hydrological regime and causing degraded water quality, pollution and increased sediment erosion into lotic systems (Waters 1995;Dudgeon et al. 2006;Woodward et al. 2012). By altering their natural condition, such activities largely downgrade the habitat integrity, which results in reduced ecological function and biodiversity (Steffen et al. 2015), including native species loss (Carpenter et al. 2011). The habitat characteristics change dramatically: formation of macrophyte assemblages is disturbed (Jones et al. 2014;Turunen et al. 2017), habitat heterogeneity and availability for macroinvertebrates is reduced (Jones et al. 2012), while primary production ) and decomposition of organic matter (Lecerf and Richardson 2010) are highly altered.
As freshwater organisms live almost continuously in the aquatic environment, they clearly respond to all those environmental stresses (Morse et al. 2007;Vilenica et al. 2019;. The aquatic assemblages can respond to alterations of their habitats with their structure differing from a reference state, i.e., they can show characteristics of "rhithralisation" (e.g., caused by channel straightening) or "potamisation" (e.g., caused by the impounding) (Jungwirth et al. 2000;Moog and Chovanec 2000;Kokavec et al. 2018;Vilenica et al. 2016;2019), or there is a change in the trophic structure (Brasil et al. 2013). By observing the assemblages' structural alterations, we can conclude that the lotic system has been altered, which in the end indicates a certain level of ecological disturbance (Moog 2002;Vilenica et al. 2016). Mayflies are able to colonise all kinds of freshwater habitats but are found to be the most diverse in lotic ones. They are among particularly sensitive aquatic macroinvertebrates, mainly disappearing when faced even with small-scale disturbance in their habitat (Firmiano et al. 2017;Vilenica et al. 2019). Previous studies demonstrated that the majority of species can tolerate a rather narrow range of environmental factors, being highly sensitive to oxygen depletion, acidification, and various contaminants such as metals, ammonia, nitrogen, phosphorous (Moog et al. 1997;Vilenica et al. 2017Vilenica et al. , 2019. Therefore, the absence/ presence of a particular species can tell us a lot about the quality of the environment it inhabits. Ecological assessments in different regions worldwide, as well as at habitats of various ecological status are necessary for effective conservation and management of freshwater habitats and their biota (Hughes et al. 1986;Stoddard et al. 2008). There-fore, in order to obtain additional data on mayfly response to anthropogenic disturbances in their habitat, we investigated mayfly assemblages and their relationship with environmental factors at 46 slightly to heavily modified lotic habitats.

Study area
The study encompassed 46 lotic slow-flowing study sites (Tables 1, 2, Fig. 1), including heavily modified streams and rivers (by, for instance, channelling and/or modification of the water flow or riverbed, removal of the riparian vegetation and pollution). The majority of the study sites are located in the vicinity of agricultural areas or cattle farms. Sampling was conducted between April and September 2016. Within the research, it was not possible to include a reference site. True reference sites are not available due to long-lasting and strong anthropogenic influence. The relatively high ratio of urban areas and even more agricultural ones are present in their catchment. The majority of the rivers have been channelled for agricultural land use purposes, or have limited lateral movement because of dykes protecting urban areas and settlements. During RFI (River Fauna Index) and assessment system development, the best available sites were chosen. The reference RFI and metrics value was calculated by adding 20% of the metric range to the high/good boundary. Study sites are part of the national monitoring program. From 25 m (small streams) to 50 m (mid-sized rivers) long sampling area was selected to cover the greatest possible diversity of microhabitats representative of the reach.
The study area is located in the Croatian part of the Pannonian lowland ecoregion (ER11) (Illies 1978). The area is characterised with temperate humid climate with warm summer (Cfb, Köppen classification) where the average temperature of the warmest month is below 22 °C (Šegota and Filipčić 2003). The average annual air temperature is around 12 °C and average annual rainfall is between 800 and 1100 mm (Zaninović et al. 2008).

Sampling protocol
Mayfly nymphs were collected together with other macroinvertebrates (AQEM protocol-AQEM expert consortium 2002). At each site, 20 subsamples were collected proportionally according to available microhabitat presence, using a benthos hand net (25 × 25 cm; mesh size = 500 μm) and pooled into one composite sample. The substrates were mainly composed of fine sediment (sand, silt, mud), lithal (stones, gravel), and aquatic vegetation (submerged and emergent). Samples were stored in 96% alcohol and analysed in the lab.
In the laboratory, subsampling was done to reduce the effort for sorting and identification. At least 1/6 of the sample was sorted until the minimum targeted number of 700 individuals was reached. The rest of the sample was also inspected searching for table 1. List of the 46 degraded lowland streams and rivers investigated in Croatia, with environmental parameters measured at the time of macroinvertebrate sampling. Codes of the study sites are as in Fig. 1. Legend: River size -S -small rivers (catchment area less than 100 km 2 ), M -medium-sized rivers (catchment area less than 1000 km 2 ). Channel width and water depth are expressed in meters. HYMO Group in SIMPER analysis -according to RFI EQR (1 -good and high; 2 -moderate; 3 -poor and bad). Tw -water temperature (°C), Oxy -dissolved oxygen content (mg/L), Con -conductivity (μS/cm), pH -pH, dominant substrates -lithal -stones, gravel; fine sediment -silt, mud, sand; phytal -aquatic vegetation. List of the 46 degraded lowland streams and rivers investigated in Croatia, with environmental parameters presented as mean value of 12 composite samples collected over a one-year period (January-December 2016) (including standard deviation, SD). Codes of the study sites are as in Fig. 1. Legend: NH 4 + -ammonium (mgN/L), NO 3 --nitrates (mgN/L), TN -total nitrogen (mgN/L), PO 4 3− -orthophosphates (mgP/L), TOC -total organic carbon (mg/L), BOD 5 -biological oxygen demand (mgO 2 /L), COD Mn -chemical oxygen demand (mgO 2 /L). macroinvertebrates which are not part of subsample analysed. Mayflies were identified to the lowest possible taxonomical level (very juvenile and/or damaged individuals were identified only to the genus or family level) using Müller-Liebenau (1969), Malzacher (1984) and Bauernfeind and Humpesch (2001). All voucher specimens are deposited at the Department of Biology, Faculty of Science, University of Zagreb, Croatia.

Environmental factors
At each study site, the following environmental parameters were measured at the time of macroinvertebrate sampling: water temperature, dissolved oxygen concentration (using the oximeter WTW Oxi 330/SET), conductivity (with the conductivity meter WTW LF 330), pH (using the pH-meter WTW ph 330), mean channel width and maximum water depth (using a hand meter on approximately 100 meter long reach of specific site) ( Table 1). The remaining environmental parameters are presented as the mean value of 12 composite samples collected over a one-year period (January -December 2016) (  (2014) gives a score of HYMO alternation based on the response of macroinvertebrate assemblages. The scores are then normalised with regard to reference states in the form of the WFD (Water framework directive) recommended EQRs (ecological quality ratios) and range from 0 (the worst HYMO conditions) to 1 (reflecting reference states). The HYMO evaluation of rivers has been performed by European Standards EN 14614 and EN 15843. Type specific RFI was used as a relative measure of HYMO alternation because HYMO evaluations for all of the investigated rivers are not available.

Data analysis
Mayfly assemblages from sites classified as high and good by the RFI EQR (EQR > 0.6) represented Group 1, from sites classified as moderate (0.4 < EQR <0.6) represented Group 2 and from sites classified as poor and bad (EQR < 0.4) represented Group 3 in the analysis of similarity percentages (SIMPER) of the (Bray-Curtis) similarity (Clarke 1993) between mayfly assemblages. This was done in order to determine how mayfly assemblages differ among sites of different degrees in HYMO alternation in terms of species composition and abundance contribution.
The composition of mayfly assemblages in terms of the trophic structure and longitudinal zonal associations of species at each study site was analysed using the classification given by Buffagni et al. (2009;, while the methodology was described in Vilenica et al. (2018). Study sites without mayfly records, and sites with one taxon where we could not identify the specimens to the species level (i.e., sites 17 and 18) were excluded from the analysis.
In order to ordinate mayfly occurrence with respect to environmental variables, the Canonical Correspondence Analysis (CCA) was used. The analysis was performed using data for 21 taxa (rare species were downweighed) and 14 environmental variables. The Monte Carlo permutation test with 499 permutations was used to test the statistical significance of the relationship between all taxa and all variables.
Mayfly taxa abundances were correlated against agricultural land cover data, using the Spearman coefficient, in order to determine if and to what extent does this type of land cover in the catchment area influence specific taxa occurrence. Mayfly species richness, abundance and local diversity (Shannon index) were plotted against the ratio of intensive agriculture in the catchment in order to determine the "general" mayfly response in relation to increased agricultural pressures.
The Bray-Curtis similarity index, Shannon diversity index and SIMPER analyses were conducted in Primer 6 (Clarke and Gorley 2006). The CCA analysis was performed using CANOCO 5.00 (ter Braak and Šmilauer 2012). Mayfly/intensive agriculture graphs were plotted, and regression equations were calculated and tested for significance using Statistica 13.0 (TIBCO Software Inc. 2017). The species data were log-transformed prior to analyses. All figures were processed with Adobe Illustrator CS6.
The SIMPER group similarity analysis (Table 4) showed that all groups of sites were dominated by juvenile instars of Baetis sp. and had significant abundances of Cloeon dipterum present at most sites. Baetis fuscatus (Linnaeus, 1761) and Baetis buceratus Eaton, 1870 were associated with sites of both ends of the HYMO gradient (Group 1 and Group 3). Furthermore, Baetis vernus Curtis, 1834 individuals were associated with sites that had a lower degree of HYMO degradation (Group 1 and Group 2). Juvenile instars of Caenis sp. were usually associated with more degraded sites (Group 3), whereas Serratella ignita and Caenis luctuosa (Burmeister, 1839) were associated only with sites of good and high ecological status following the RFI.
Generally, a high share of lower reaches and lentic elements (potamic and littoral elements) was recorded: it was dominant (> 50 %) at 13 study sites, eight sites had an equal share of lower reaches/lentic and upper reaches elements (crenal and rhithral) (50:50 %), while16 study sites were dominated by upper reaches elements (> 50 %) (Fig. 2a). We also recorded a high share of detritivores (gatherers/collectors and active filter feeders): they were dominant at 21 study sites and equally represented as grazers/ scrapers at the rest of the sites (Fig 2b).

Mayflies and environmental variables
The results of the ordination of species and environmental data of the CCA are presented on the F1 × F2 ordination plot (Fig. 3). The eigenvalues for the first two CCA axes were 0.40 and 0.25 and explained 50.9 % of the species-environment relations. The Monte Carlo permutation test showed that the species-environment ordination was significant (first axis: F-ratio = 4.23, p = 0.002; overall: trace = 1.28, F = 1.54, p = 0.006) indicating that mayfly assemblages were significantly related to the tested set of environmental variables. Axis 1 was related to total organic carbon (R = 0.49)   and dissolved oxygen (R = -0.46), and axis 2 to aquatic vegetation (R = -0.37) and water temperature (R = -0.36), indicating that these were the most important parameters in explaining patterns of mayfly assemblages (Fig. 3). Mayfly species richness, abundance and consequently also local diversity, were found to significantly decrease with increased ratios of intensive agriculture areas in the catchment area (Fig. 4).

Discussion
Our results indicate that a relatively high number of mayfly species can be found in anthropogenically impacted freshwater habitats. Nevertheless, at a large part of the study sites (i.e., 72 %) taxa richness was low, i.e., between zero and four taxa, corroborating previous studies (Vilenica et al. 2016;2019). Mayflies inhabit both lotic and lentic habitats, although upper and middle reaches of fast-flowing streams, and ecologically intact large rivers harbour the highest mayfly diversity (Bauernfeind and Soldán 2012;Vilenica et al. 2016;2018). Therefore, such low species richness, not typical for a lotic habitat (Bauernfeind and Moog 2000;Zedková et al. 2014;Vilenica et al. 2018), could be a consequence of various disturbances present at those sites, such as channelling, eutrophication, pollution, and microhabitat homogeneity (Axelsson et al. 2011;Carvalho et al. 2013;Ligeiro et al. 2013). In many cases, we observed shoreline erosion, as the emergent vegetation along the habitat edges, together with surrounding vegetation was mowed. This could have resulted in an increased input of sediments into the habitats, which could have influenced the habitat physico-chemical characteristics and hydrological cycle, resulting in reduced water quality and habitat heterogeneity (Mendes et al. 2017 and references herein). Consequently, these habitats showed to be less favourable for a high number of mayfly species. The majority of study sites were inhabited by widespread and generalist species (Popielarz and Neal 2007;Bauernfeind and Soldán 2012), yet sites with more microhabitat heterogeneity and higher water velocity, had also several microhabitat specialists, such as Baetis lutheri and Ecdyonurus torrentis for mesolithal, and Centroptilum luteolum as specialists for macrophytes (Buffagni et al. 2009;. The Zelina stream in Božjakovina (site 22) and Toplica River upstream from Daruvar town (site 36) showed somewhat higher species richness, yet their assemblages mainly consisted of species inhabiting a wide range of habitats, such as Baetis rhodani, Centroptilum luteolum, Serratella ignita and Caenis luctuosa (Buffagni et al. 2009;Bauernfeind and Soldán 2012). The most interesting finding was a record of Oligoneuriella rhenana at Toplica River, which is considered rare in Croatia (Vilenica et al. 2015;2018). Although the species can tolerate some variations of environmental factors, its presence indicates that the ecological condition of Toplica River upstream from  Table 2. Legend: Environmental variables (red arrow symbols): Tw -water temperature (°C), Oxy -dissolved oxygen content (mg/L), Con -conductivity (μS/cm), pH -pH, NH 4 + -ammonium (mgN/L), NO 3 --nitrates (mgN/L), TN -total nitrogen (mgN/L), PO 4 3− -orthophosphates (mgP/L), TOC -total organic carbon (mg/L), BOD 5 -biological oxygen demand (mgO 2 /L), COD Mn -chemical oxygen demand (mgO 2 /L), vegetation -aquatic vegetation/phytal, fine sediment -silt, mud and sand, lithal -stones and gravel.
Daruvar town is not as poor as at the majority of other sites (Găldean, 1999;Petrovici and Tudorancea 2000). Another interesting species was the rarest in our study, a riverine Heptagenia flava, uncommon in Croatian waters (Vilenica et al. 2015). Although the species was reported to have rather high ecological plasticity, usually it does not inhabit heavily polluted rivers (Vidinova and Rusev 1997). Therefore, the species record at Česma River in Narta (site 28) could be considered as an accidental finding, as shown by Vidinova and Rusev (1997). On the other hand, two eurytopic and euryvalent species (i.e., with wide tolerance towards the environmental conditions and habitat type), Cloeon dipterum and Serratella ignita, were recorded as the most common and the most numerous, respectively (Buffagni et al. 2009;Bauernfeind and Soldán 2012;Vilenica et al. 2019). Nevertheless, while discussing the total species richness at a particular site, we need to keep in mind that standardised sampling methods generally do not include sampling of underrepresented microhabitats, which could be important for some rare species (Haase et al. 2008). Therefore, in order to obtain a more complete species list, it might be beneficial to complement standardised quantitative sampling with a qualitative one.
Stream channelling is a widely used engineering practice designed for flood control and wetland draining, which affects the majority of hydrogeomorphological characteristics and processes at the channelled habitat. Due to these changes, the biota is also severely affected (Hupp 1992), i.e., the community structure and composition are changed and poorer (Waters 1995). Our results showed that mayfly assemblages have mainly consisted of taxa of potamic (lower reaches) and lentic preferences (e.g., Baetis buceratus, Caenis horaria) or wide range (e.g., Cloeon dipterum, Centroptilum luteolum, Serratella ignita) habitat type preferences (Buffagni et al. 2009;Bauernfeind and Soldán 2012). Moreover, Baetis vernus, Caenis luctuosa and Serratella ignita, species with relatively strong rhithral affinity (Biss et al. 2002) were predominantly associated with hydromorphologically less degraded sites, while species with more prominent potamic preference, such as Baetis buceratus and Baetis fuscatus (Schöll et al. 2005) were present at sites both with low and high degree of hydromorphological degradation. Some study sites showed a higher share of rhithral elements, yet that was mainly due to the dominance of eurytopic Cloeon dipterum (Buffagni et al. 2009;. As the majority of sites are characterised by low microhabitat diversity, a high level of sedimentation and nutrients, assemblages were dominated by detritivores (Buffagni et al. 2009(Buffagni et al. , 2020. Previous researches showed that mayflies are highly dependent on specific environmental cues, and many species rapidly disappear when faced with anthropogenic disturbances in their habitat (Bauernfeind and Moog 2000;Goulart and Callisto 2005;Stepanian et al. 2020). Our results corroborate previous studies that showed negative responses of mayflies to high water temperature (e.g., Chadwick and Feminella 2001;Alhejoj et al. 2014) and low oxygen concentrations (e.g., Nebeker 1972;Lock and Goethals 2011). Sites that were characterised by high water temperatures were also often accompanied by low oxygen content and dense aquatic vegetation. High levels of nutrients in the water support such dense growth of vegetation, leading to a decrease of oxygen level (Boeykens et al. 2017). Moreover, the decay of organic matter (especially aquatic vegetation), together with bacterial growth, animal/human metabolic activity and various synthetic sources (such as pesticides, fertilisers, pharmaceuticals, detergents) lead to elevated concentration of total organic carbon (TOC) in water (e.g., Volk et al. 2002). A part of the TOC can be explained by the increased shoreline erosion due to management and clearing of vegetation in the shoreland zone, which probably also negatively affected mayflies in this study. Riparian buffers, especially undisturbed vegetated riparian zones situated adjacent to river and streams, can greatly mitigate nutrients, sediment from surface and groundwater flow through the processes of deposition, absorption and denitrification (e.g., Peterjohn and Correll 1984). Finally, the strong negative association of mayfly assemblages with intensive agriculture in the catchment area corroborates results of previous studies that showed high mayfly sensitivity to agricultural pollution (Siegloch et al. 2014;Zedková et al. 2015). Here, as especially sensitive showed Alainites muticus, Baetis lutheri and Oligoneuriella rhenana, species with low and moderate tolerance to water pollution (mainly occurring in oligosaphrobic and beta-mesosaphrobic waters) (Bauernfeind et al. 2002;Mihaljević 2011). In addition, another species was distinguished as sensitive to such kind of pollution, Baetis rhodani. Those results could come as a surprise, as this eurytopic mayfly has a wide ecological tolerance, and generally contributes as a major part of the macroinvertebrate biomass in many European streams and rivers (Elliott et al. 1988). Nevertheless, as Baetis rhodani is a species complex (Williams et al. 2006), those results should be inspected in more details, using molecular analyses. Our results confirm that water pollution is one of the largest limitation factors for the majority of mayflies (Van Dijk et al. 2013;Zedková et al. 2015).

Conclusions
This study contributes to our knowledge of mayfly relationship with environmental conditions in heavily modified and anthropogenic habitats. Various anthropogenic pressures resulted in changes in mayfly assemblage composition and structure, whereas species richness decreased. For instance, the assemblages consisted mainly of a relatively low number of widespread generalists and species characteristic for lower reaches and lentic habitats. This indicates that hydromorphological alterations could have resulted in assemblage's "potamisation". Moreover, highly polluted sites, with high temperatures and low oxygen content, were inhabited almost exclusively with the euryvalent Cloeon dipterum, or were completely unsuitable for any mayfly species, confirming the high sensitivity of mayflies to disturbances in their habitats. Our results can enable planning of management and conservation activities of lowland rivers and their biota according to the requirements of the European Water Framework Directive.