Ecological Preferences and Indication Potential of Freshwater Bryophytes–Insights from Croatian Watercourses

A comprehensive survey of Croatian watercourses covering the whole of the national territory and investigating inherent watercourse heterogeneity was conducted to explore the ecological responses of the most frequent freshwater bryophytes with respect to water chemistry variables and land use within the catchment area. Direct multivariate ordination (CCA) of vegetation data paired with 18 environmental variables revealed that freshwater bryophytes and their assemblages were segregated along the gradients of water chemistry and the proportion of natural and urban area within the catchment. Generalized additive models (GAM) were employed to explore the ecological responses of individual species. The results showed that most of the investigated species preferred natural, clean, well-oxygenated watercourses, with low nutrient and organic matter content, as well as with low electrical conductivity. Species such as Palustriella falcata, Eucladium vertcillatum, Dichodontium flavescens and Jungermannia atrovirens had narrow ecological niches and were restricted to pristine watercourses, while the most frequent and widely distributed species, such as Fontinalis antipyretica, Rhynchostegium riparioides, Cratoneuron filicinum, Fissidens crassipes, Cinclidotus fontinaloides and C. riparius, had a wide ecological tolerance. Riccia fluitans and Leptodyctium riparium had wide ecological ranges, but with optima in hypereutrophic waters with high nutrient and organic content, as well as high electrical conductivity. Furthermore, these two species were frequently associated with a high share of intensive agriculture and a low share of natural land within the catchment.


Introduction
Freshwater bryophytes are a common but unassuming and frequently overlooked component of freshwater ecosystems. They constitute an important part of macrophyte vegetation, especially in headwater streams, mountain and upland watercourses, and highly seasonal and intermittent rivers [1][2][3][4][5], and they are the most prominent part of the vegetation of waterfalls and cascades [6][7][8][9]. Bryophytes dominate the vegetation of such lotic habitats due to a wide variety of structural and physiological adaptations [6,10,11]. These make them resilient to seasonal desiccation, high water velocity and associated mechanical stress, which other macrophyte representatives cannot withstand. Freshwater bryophytes play a significant role in these harsh environments as dominant primary producers and influence the overall nutrient and trophic dynamics. Bryophyte beds make an excellent shelter and habitat for various invertebrate assemblages and provide a surface for the growth of epiphytic algae [12,13]. Freshwater bryophytes occur in middle and lower river reaches as well, but in much lower abundance, with the exception of some of the Mediterranean rivers. In middle and lower reaches, they are in general represented by a lower number of rheophyte taxa restricted to larger and thus more stable substrates, or by a somewhat more diverse set of semi-aquatic species inhabiting periodically flooded river margins [14,15].
While the presence and coverage of freshwater bryophytes are primarily determined by riverbed stability and substrate size [16][17][18], physiographic, geological, climatic and 1.
Determine how water physicochemical factors and land use influence species occurrences and segregation; 2.
Explore ecological responses of freshwater bryophyte species and augment data on the autecology and ecological preferences of freshwater bryophytes; 3.
Infer the bioindication potential of selected species from their optima and ecological tolerance.

Results
A total of 21 freshwater bryophyte species (8 rheophytes, 1 hydrophyte, 2 amphyphytes and 10 hygrophytes) were collected from at least five localities, out of 648 localities surveyed, and were thus included in further analysis. Finally, 182 localities were retained by this criterion, 124 within the Dinaric and 58 within the Pannonian Ecoregion of Croatia ( Figure 1). This encompassed diverse watercourse types, from oligotrophic small streams to eutrophic large rivers and artificial canals (Table A1, Appendix A), therefore covering wide gradients suitable for the study of the autecology of selected species. In the Pannonian Ecoregion, the highest bryophyte richness was recorded in small lowland watercourses (16 species). In the Dinaric-Continental Subecoregion, this was true for montane and mid-altitude medium and large watercourses (20 species), followed by montane and mid-altitude small watercourses (18 species). In the Dinaric-Mediterranean Subecoregion, the highest bryophyte richness was recorded in lowland and mid-altitude small watercourses (16 species) (Table A2, Appendix A). As for lowland medium and large watercourses, they harboured a similar set of a total of 11 species and similar overall bryophyte frequency both in the Pannonian Ecoregion and the Dinaric-Mediterranean Subecoregion, while these were substantially higher in the Dinaric-Continental Subecoregion. Species number, as well as their overall frequencies and abundance, in both Dinaric subecoregions exceeded those in the Pannonian Ecoregion (Table A2, Appendix A, Table S1).
Among the most frequent species in our study were Fontinalis antipyretica, Rhynchostegium riparioides, Cratoneuron filicnum and Leptodyctium riparium, with the latter being more frequent in the Pannonian Ecoregion, i.e., its lowland watercourses. Species such as Didymodon tophaceus, Eucladium verticillatum and Jungermannia atrovirens were exclusive to the Dinaric Ecoregion, while Cinclidotus aquaticus, Apopellia endiviifola, Brachythecium rivulare, Dichodontium flavescens and D. pellucidum were more frequent in the Dinaric Ecoregion and within the Pannonian Ecoregion, restricted to only lowland small watercourses.
Species medians revealed the differences in species preferences across the gradients of selected variables, with Riccia fluitans and Leptodyctium riparium being most obviously separated from the rest of the species in terms of water chemistry and land use variables. These species showed preferences for hypereutrophic water with high nutrient loads, biochemical (BOD) and chemical (COD) oxygen demand and high electrical conductivity, and were more frequently associated with a higher share of intensive agriculture and a low share of natural land within the catchment area ( Figure 2). Most other species had their

Ordination Results
Results of the CCA revealed similar patterns in the ecological preferences of the studied species, except for Riccia fluitans, which, as an outlier, was excluded from this analysis. The first axis of the CCA explained 51.06% and the second 16.01% of the variation in the relationship between vegetation data and environmental factors, while eigenvalues of the first and the second axis equalled 0.4 and 0.1, respectively. Forward selection of environmental variables revealed a set of eight most-contributing and non-redundant variables in determining the freshwater bryophyte distribution (Appendix B). Total nitrogen concentration made the largest contribution to explaining the observed variation in bryophyte composition, followed by the share of the natural area within the catchment and chemical oxygen demand. Overall analysis was statistically significant (p < 0.002), which was confirmed by the Monte Carlo test (999 permutations).
The CCA analysis revealed the main compositional gradient representing the water quality along axis 1-from sites with well-oxygenated water, situated within unchanged catchment areas under little anthropogenic influence to sites with high nutrient loads, i.e., high concentration of total nitrogen and total phosphorous, as well as high chemical oxygen  (Figures 3 and 4). Additionally, the share of the natural area within the catchment was strongly negatively correlated with agricultural land-extensive (AGE, r s = −0.72, p < 0.001) and intensive agriculture (AGI, r s = −0.81, p < 0.001). Regarding the nutrients, total phosphorus was highly positively correlated with orthophosphates (PO 4 3− , r s = 0.82, p < 0.001), total suspended solids (TSS, r s = 0.73, p < 0.001) and biochemical oxygen demand (BOD, r s = 0.73, p < 0.001), and total nitrogen with nitrates (NO 3 − , r s = 0.85, p < 0.001), while COD correlated with biochemical oxygen demand (BOD, r s = 0.76, p < 0.001). Sites with higher nutrient loads, as well as COD and higher shares of agricultural land, were mostly those on the watercourses situated in the Pannonian Ecoregion and were quite well separated from the sites situated in the Dinaric Ecoregion on the CCA ordination biplot, especially those of its Mediterranean Subecoregion ( Figure 3). The Dinaric Ecoregion included sites with intermediate values of the studied parameters, as well as sites within natural catchments, with oligotrophic, well-oxygenated water and somewhat higher pH.  Fissidens crassipes, Cinclidotus riparius and Rhynchostegium riparioides were associated with a higher share of urban area within the catchment and higher electrical conductivity.

Species Responses to Environmental Variables
The GAM response curves of species abundances for the most influential environmental variables from the CCA additionally corroborated observed patterns ( Figure 5). Response curves of GAMs fitted for Riccia fluitans and Leptodyctium riparium against total nitrogen and total phosphorus showed the preference of both species for hypereutrophic water with a high nutrient load. Regarding total nitrogen, the means of both species were quite high, 1.93 mgN/L for Riccia fluitans and 1.45 mgN/L for Leptodyctium riparium (Appendix C). Leptodyctium riparium displayed unimodal response with an optimum at 4.90 mgN/L, while Riccia fluitans had a monotonically increasing response curve and a maximum at 8.9 mgN/L. Other species had their optima at low levels of total nitrogen (<0.5 mgN/L). A steep decline in abundance with increasing concentrations of total nitrogen was characteristic of most of these species ( Figure 5). Species such as Palustriella falcata, Jungermannia atrovirens, Dichodotium flavescens, Eucladium vericillatum, Didymodon tophaceus, Chilosyphus polyanthos and Apopellia endiviifolia had quite low maxima and narrow niches, while the maxima of Cratoneuron filicinum, Brachythecium rivulare, Cinclidotus fontinaloides, C. riparius, Fissidens crassipes, Fontinalis antipyretica and Marchantia polymorpha were between 2 and 3 mgN/L, and the maximum of Rhynchostegium riparioides reached 5.18 mgN/L. These sets of species showed similar behavior with respect to nitrate and The CCA analysis demonstrated the affinity of particular freshwater bryophyte species to different water chemistry and land use ( Figure 4). Mosses such as Cinclidotus aquaticus, Eucladium verticllatum and Palustriella falcata, as well as the liverworts Apopellia endiviifolia, Chiloscyphus polyanthos and Jungermannia atrovirens, preferred clean, well-oxygenated water and were more strongly correlated with a higher share of natural area within the catchment. On the other hand, Leptodyctium riparium was of all species investigated the most prominently associated with eutrophic water, while Cratoneuron filicinum, Cinclidotus fontinaloides, C. riparius, Fontinalis antipyretica, Fissidens crassipes and Rhynchostegium riparioides displayed intermediate behavior along the first axis. Among these, Fissidens crassipes, Cinclidotus riparius and Rhynchostegium riparioides were associated with a higher share of urban area within the catchment and higher electrical conductivity.

Species Responses to Environmental Variables
The GAM response curves of species abundances for the most influential environmental variables from the CCA additionally corroborated observed patterns ( Figure 5). Response curves of GAMs fitted for Riccia fluitans and Leptodyctium riparium against total nitrogen and total phosphorus showed the preference of both species for hypereutrophic water with a high nutrient load. Regarding total nitrogen, the means of both species were quite high, 1.93 mgN/L for Riccia fluitans and 1.45 mgN/L for Leptodyctium riparium (Appendix C). Leptodyctium riparium displayed unimodal response with an optimum at 4.90 mgN/L, while Riccia fluitans had a monotonically increasing response curve and a maximum at 8.9 mgN/L. Other species had their optima at low levels of total nitrogen (<0.5 mgN/L). A steep decline in abundance with increasing concentrations of total nitrogen was characteristic of most of these species ( Figure 5). Species such as Palustriella falcata, Jungermannia atrovirens, Dichodotium flavescens, Eucladium vericillatum, Didymodon tophaceus, Chilosyphus polyanthos and Apopellia endiviifolia had quite low maxima and narrow niches, while the maxima of Cratoneuron filicinum, Brachythecium rivulare, Cinclidotus fontinaloides, C. riparius, Fissidens crassipes, Fontinalis antipyretica and Marchantia polymorpha were between 2 and 3 mgN/L, and the maximum of Rhynchostegium riparioides reached 5.18 mgN/L. These sets of species showed similar behavior with respect to nitrate and ammonium concentration, as well as other water chemistry parameters associated with water quality-orthophosphates, total phosphorus, the amount of organic compounds in water (COD, BOD) and electrical conductivity (Figures 1 and 5).
Regarding the total phosphorus concentration, most of the studied species displayed mean values below 0.04 mgP/L, while the mean values of Riccia fluitans and Leptodyctium riparium were 0.177 and 0.094 mgP/L, respectively, which corresponded to eutrophic water. The optima of both species were in hypereutrophic water; the optimum of Leptodyctium riparium amounted to 0.21 mgP/L, and of Riccia fluitans to 0.699 mgP/L ( Figure 5). All other species had their optima in oligotrophic water with respect to phosphorus, with Palustriella falcata, Eucladium verticillatum, Jungermannia atrovirens, Apopellia endiviifolia, Chiloscyphus polyanthos, Didymodon tophaceus and Dichodontium flavescens having steep monotonically decreasing response curves and low maxima, and Cinclidotus aquaticus and Brachythecum rivulare tolerating somewhat higher concentrations and disappearing at 0.1 and 0.165 mgP/L, respectively. The abundance of rheophytes such as Fontinalis antipyretica, Cinclidotus fontinaloides, C. riparius, Fissidens crassipes and Rhynchostegium riparioides, as well as of the amphyphyte Cratoneuron filicinum and the hygrophytes Chiloscyphus pallescens and Ptychostomum pseudotriquetrum, decreased with increasing concentration of total phosphorus, but these species displayed tolerance to eutrophic water. Total phosphorus concentration was highly correlated with total suspended solids, and species preferring, or displaying a high degree of tolerance to, eutrophic water were tolerant of more turbid water as well. These were Leptodyctium riparium, Fontinalis antipyretica, Rhynchostegium riparioides and Riccia fluitans, persisting at TSS values over 40 mg/L. Nevertheless, the TSS means of all species were below 15 mg/L (Appendix C), i.e., in clear water.
Species response curves for chemical oxygen demand showed that the majority of the species had low optima (<1 mgO 2 /L), while Leptodyctium riparium and Riccia fluitans peaked at high levels, 8.2 and 7.9 mgO 2 /L, respectively, and Dichodontium pellucidum at 4.9 mgO 2 /L ( Figure 5). Again, the majority of rheophytes, as well as the amphiphyte Cratoneuron filicinum, showed tolerance to high levels of COD, although their abundance decreased with increasing COD. This was also the case with the hygrophyte Chiloscyphus pallescens, which was recorded at a maximum COD level of 10.67 mgO 2 /L. GAMs were not successfully fitted for Didymodon tophaceus and Ptychostomum pseudotriquetrum, but these species tolerated very high COD levels as well, with maxima around 10.6 mgO 2 /L. Species restricted to low levels were Eucladium verticillatum, Jungermannia atrovirens, Apopellia endiviifolia, Cinclidotus aquaticus, Brachythecium rivulare and Chiloscyphus polyanthos. Similar patterns in species preferences can be observed from the descriptive statistics regarding biochemical oxygen demand as well (Appendix C). Response curves fitted against dissolved oxygen revealed that the majority of species investigated had their optima in well-oxygenated water ( Figure 5). Mean values of this parameter ranged from 10.13 to 11.58 mgO 2 /L, except for Riccia fluitans. This species peaked at 6.5 mgO 2 /L, but was present within a quite wide range, from 5.25 to 9.34 mgO 2 /L. The optimum of Leptodyctium riparium was at 4 mgO 2 /L, with a prominent decrease in abundance with increasing oxygen concentration. This coincided with the rise in the abundance of other rheophyte species. Nevertheless, Leptodyctium riparium was still present in well-oxygenated watercourses, displaying wide tolerance. The set of species restricted to high dissolved oxygen concentration was similar to that characteristic of low nutrient and COD situations. These were Brachytehcium rivulare, Palustriella falcata, Jungermannia atrovirens, Dichodontium flavescens, Chiloscyphus polyanthos, Apopellia endiviifolia and Cinclidotus aquaticus, as well as Dichodontium pellucidum and Chiloscyphus pallescens. Other species, such as Rhynchostegium riparioides, Fontinalis antipyretica, Fissidens crassipes, Cratoneuron filicinum, Cinclidotus fontinaloides, C. riparius and Ptychostomum pseudotriquetrum had wider ranges and were present in less oxygenated watercourses, although in lower abundance.
The response curve was not successfully fitted for Didymodon tophaceus, but the mean value was at 10.13 mgO 2 /L and the minimum as low as 7.51 mgO 2 /L (Appendix C).
Regarding electrical conductivity, most of the species had optima below 400 µS/cm but displayed different tolerances to higher values. Hygrophyte species had quite narrow niches, with Apopellia endiviifolia and Chiloscyphus polyanthos and Ptychostomum pseudotriquetrum being more tolerant to elevated electrical conductivity than other hygrophytes. Response curves were not fitted for Dydimodon tophaceus and Marchantia polymorpha but descriptive statistics indicated wide tolerance with respect to electrical conductivity (Appendix C). Similar behavior can be seen in the case of the majority of rheophyte species, except for Brachythecium rivulare and Palustriella falcata, which were restricted to waters with lower conductivity, while Rhynchostegium riparioides, Cratoneuron filicinum, Cinclidotus aqauticus, C. riparius and C. fontinaloides showed tolerance to high values of electrical conductivity. The optimum of Leptodyctium riparium was around 600 µS/cm, while the abundance of Fontinalis antipyretica and Riccia fluitans displayed a continuous rise until the maxima above 900 µS/cm ( Figure 5). Concerning water pH values, most of the species preferred basic water and showed an increase in abundance above pH 7.9. Leptodyctium riparium peaked at 7.87, which was followed by a relatively steep decrease in abundance, which coincided with the increase in the abundance of basophilous species ( Figure 5). Furthermore, the abundance of Riccia fluitans decreased along the pH gradient and the species disappeared when pH exceeded 7.87. Similarly, Chiloscyphus pallescens showed a sharp decrease along this gradient, being most abundant around pH 7, but still present at 8.2.
GAM response curves fitted against the share of the natural land within the catchment area revealed that the majority of the freshwater bryophytes studied had their optima at values over 80%. Exceptions were again Riccia fluitans and Leptodyctium riparium. Riccia fluitans displayed a steep monotonically decreasing curve with an optimum in the least natural catchments ( Figure 5). The optimum of Leptodyctium riparium was at 33.57% of the natural area within the catchment but it was present in more natural catchments as well, although in lower abundance. Mean values of this parameter were over 80% for the set of species that were restricted to waters with low nutrient content, low COD, BOD and electrical conductivity, and high dissolved oxygen (Eucladium verticillatum, Apopellia endiviifolia, Cinclidotus aquaticus, Jungermannia atrovirens, Palustriella falcata, Chiloscyphus polyanthos, Dichodontium flavescens and D. pellucidum), and ranged between 60 and 80% for all other species with wider niches, except for Riccia fluitans with a mean value equal to 34.91% (Appendix C). This species was highly associated with intensive agriculture, with the mean value of this parameter being 52.20% and a maximum of 97.69% ( Figure 1, Appendix C). On the other hand, species such as Eucladium verticillatum, Brachythecium rivulare, Palustriella falcata, Dichodontium flavescens, D. pellucidum, Chiloscyphus polyanthos, Apopellia endiviifolia and Jungermannia atrovirens had mean values of intensive agriculture within the catchment area lower than 5% and maxima lower than 15%. Considering the urban area within the catchment, most species preferred low levels, with rheophytes such as Rhynchostegium riparioides, Cinclidotus fontinaloides, C. aquaticus, C. riparioides, Fontinalis antipyretica and Cratoneuron filicinum peaking below 12% and then decreasing in abundance and finally disappearing at around 30% of urban area within the catchment.

Discussion
Our results are in line with previous studies, confirming that different water chemistry factors associated with water quality influence the distribution and segregation of freshwater bryophytes [37,38]. Similarly, changes in land use within the catchment area, such as an increase in the share of intensive agriculture or urban area, which are most often associated with water pollution and eutrophication, also influence freshwater bryophytes [20].
The results of our study reveal that most of the studied freshwater bryophytes prefer natural catchments with clear, well-oxygenated water that has low nutrient levels, both nitrogen and phosphorus, as well as a low amount of organic matter. Similar results were obtained from the study of the freshwater bryophytes in the Tiber River Basin (Italy), where the majority of the species showed a general preference for fast-flowing, clear, cold and oxygenated water with a low nutrient load, especially low ammonia and orthophosphates [21].
Although most species in our study had their optima in water of good quality, differences in their responses to investigated ecological gradients revealed a set of species with very narrow ecological niches concerning this group of parameters. Species such as Palustriella falcata, Eucladium verticillatum, Jungermannia atrovirens and Dichodontium flavescens were mostly characteristic of the Dinaric Ecoregion and were recorded in pristine karstic watercourses with low nutrient content and low organic matter. Similar results were reported earlier for these species. In a study dealing with Natura 2000 petrifying sources in Belgium, the occurrence of Eucladium verticillatum, a tufa-forming species, was negatively correlated with ammonium and phosphate concentration, and the species preferred open habitats with lots of light [41]. Palustriella falcata, known as calcicole, is mostly linked to base-rich, neutro-alkaline, cold and turbulent oligotrophic streams with low conductivity [3,18,39], while Jungermannia atrovirens was reported in both oligotrophic and meso-oligotrophic waters [55]. Dichodontium flavescens was in our study almost exclusive to the Dinaric-Continental Subecoregion and restricted to clean, cold and turbulent karstic watercourses with almost completely natural catchment areas. Similarly, according to Dierßen [56], Dichodontium flavescens inhabits areas under no or only weak human impact.
Additionally, species such as Cinclidotus aquaticus, Chiloscyphus polyanthos, Apopellia endiviifolia and Didymodon tophaceus had somewhat wider niches, with optima in oligotrophic waters but tolerating elevated nutrient levels to some extent. However, their presence in higher abundance may be indicative of good water quality. Cinclidotus aquaticus, was previously reported as a species of clear, cold and turbulent oligotrophic waters and highlighted as a valid indicator of water quality [21], unlike Apopellia endiviifolia which was found across a broader spectrum of water quality, from oligotrophic [38,57] to eutrophic [21,41,58], tolerating high nitrate, ammonium and orthophosphate concentration.
All the above-mentioned species are basophilous and are either exclusive to or dominantly occur in karstic rivers of the Dinaric Ecoregion. These rivers flow over carbonate bedrock which influences water pH and alkalinity, as well as species assemblages where basophilous, i.e., acid-sensitive, taxa dominate. An important influence of geology, water pH and alkalinity on aquatic bryophytes has been already demonstrated on the European level [1,2,4,5,30,39,57] and beyond [6,17,40,59], and clear segregation of aquatic bryophytes along the alkalinity and pH gradient was recently demonstrated for Croatian bryophytedominated watercourses as well [22]. The latter study identified three communities characterized by different basophilous species which were mostly associated with karstic rivers of the Dinaric Ecoregion, and two communities in small rivers situated in the Pannonian Ecoregion, which were dominated by a high share of hygrophyte taxa inhabiting periodically flooded river margins. The influence of the water pH on aquatic bryophytes was evident in our study as well, although the investigated pH gradient was quite short. The majority of the species included in our study preferred near-neutral to basic water, with the mosses Eucladium verticillatum, Didymodon tophaceues, Palustriella falcata, Brachythecium rivulare, Fissidens crassipes and Cinclidotus aquaticus, as well as the liverworts Apopellia endiviifolia and Chiloscyphus polyanthos, being most strongly associated with higher water pH and having high frequencies and abundance in the watercourses of the Dinaric Ecoregion. This ecoregion is known to harbor greater diversity regarding freshwater bryophytes [15] and their communities [22] since its fast, cold montane and semi-montane karstic rivers with larger and more stable substrates provide more suitable habitats than the lowland rivers of the Pannonian Ecoregion, which are usually slow and warmer, with dominantly sandy and gravelly substrates. This was demonstrated in our study as well, with the highest number of species recorded in montane and mid-altitude small watercourses, followed by montane and mid-altitude medium and large watercourses of the Continental Subecoregion, which are both characterized by the dominance of large and medium substrates.
Widely distributed and frequent species such as Fontinalis antipyretica, Rhynchostegium riparioides, Cratoneuron filicinum, Fissdens crassipes, Brachythecium rivulare, Cinclidotus riparius and C. fontinaloides displayed in our research broad ecological tolerance regarding the investigated environmental gradients, with their maxima of ammonium, nitrate and total nitrogen concentrations exceeding the thresholds for good water chemistry status set for all water body types included in this analysis according to the Croatian Regulation on Water Quality Standard [60]. Furthermore, the same was observed in the case of their maxima for total phosphorus, except for Brachythecium rivulare and Fontinalis antipyretica. The maxima of the latter two species did exceed the values set for good status considering total phosphorous for all water body types of the Dinaric Ecoregion, which are naturally more oligotrophic. However, for all investigated Pannonian water body types, they exceeded the values set for very good status. Regarding the organic matter content, the COD maxima of Fontinalis antipyretica, Rhynchostegium riparioides, Cratoneuron filicinum, Fissdens crassipes, Cinclidotus riparius and C. fontinaloides again exceeded thresholds set for good status for all types in the study, while that of Brachythecium rivulare exceeded thresholds set for all investigated Dinaric types, as well as thresholds for very good status regarding this parameter for investigated Panonnian types. Nevertheless, all of these species had their optima in clean, oligotrophic, oxygenated and slightly basic water of low electrical conductivity, except for Fontinalis antipyretica, which preferred high values of the latter parameter. Fontinalis antipyretica and Rhynchostegoium riparioides, which are both widespread and the most common species in Croatia [15], are often found together occupying a wide range of freshwater habitats [18,39,61,62]. Their broad trophic range, and in general wide ecological behavior, have been emphasized by several authors [3,18,24,25,38], with Fontinalis antipyretica being more tolerant of eutrophication and elevated electrical conductivity [21,38]. However, its frequency was reported to increase with decreasing concentrations of nitrates and phosphates [21], as was the case with its abundance in our study. Interestingly, Fontinalis antipyretica, as well as Fissides crassipes, displayed distinctly different response curves in two hydrographic networks in France and Belgium [37], having a maximal frequency in oligotrophic water in one and tolerating the most polluted waters in the other, while their overall optima within the study were in eutrophic waters when data from both hydrographic networks were considered simultaneously. The authors suggested that such species may include several ecotypes with different trophic requirements within different hydrographic networks, which is possibly a result of microevolution, favoured by the fact that river basins are rarely interconnected. The observed differences in autecology between populations of the same species complicate the use of certain aquatic bryophytes as bioindicators of water quality on a large scale, and thus further research on their distribution and ecological responses, as well as their taxonomy, microevolution processes and ecophysiology, is still welcomed to elucidate the influence of environmental factors on the species that have so far shown contradictory behavior.
Furthermore, studies encompassing larger geographic areas and gradients of water quality parameters could improve the knowledge of the autecology and ecological tolerance of species that have shown contradictory behavior in different studies. An example of contradictory findings being reported for the same species can be seen in the case of the common basophilous species, Cratoneuron filicinum. It was reported as a valid bioindicator of water quality with narrow ecological tolerance, preferring clear, turbulent waters, with temperature below 12 • C, low conductivity (below 300 µS/cm), and low concentration of nutrients (phosphates about 0.01 mg/L, a maximum concentration of ammonium 0.10 mg/L and nitrates 0.90 mg/L), in a study covering the Tiber River Basin [21], while tolerance to light to moderate eutrophication has been reported in several other studies as well [20,41,58]. In our study, which encompassed a larger and more diverse geographic area with many different hydrographic networks, a broader ecological behavior of this species was observed. It was more frequent in colder watercourses, but occurred in much warmer waters as well, with a maximum of 18.75 • C. Similarly, it preferred medium conductivity of about 420 µS/cm, while tolerating levels as high as 941 µS/cm. Regarding the nutrients, its optimum was in clean water of low trophy level, but the species persisted in waters with a great nutrient load as well (e.g., maximum concentration of total phosphorous 0.28 mgP/L, orthophosphates 0.09 mgP/L, total nitrogen 2.00 mgN/L and ammonium 0.86 mgN/L).
Riccia fluitans and Leptodyctium riparium showed markedly different behavior from all the other species in our study, preferring neutral, warmer, hypereutrophic water with high electrical conductivity and organic matter content, which is in line with previous findings [18,21,23,26,37,63]. In our study, although these species displayed quite wide water quality niches, their abundance was the highest in eutrophic situations, with a prominent fall in abundance with an increase in water quality. Similarly, the frequency of these species was positively correlated with the concentration of phosphates and ammonia, as well as electrical conductivity, in a study conducted by Ceschin et al. [21] in Italy. While these authors referred to both species as indicators of eutrophication, others found that Leptodyctium riparium exhibited a broad ecological range [37], from oligotrophic streams to hypertrophic rivers, and did not appear as a reliable indicator, although its frequency increased under eutrophic conditions [38]. Noteworthy is that in our study, the optima of Riccia fluitans for total nitrogen and total phosphorous (8.9 mgN/L and 0.7 mgP/L, respectively) were considerably higher than those of L. riparium (4.9 mgN/L and 0.21 mgP/L, respectively), and the species favoured the sites with the least natural catchment areas with a large proportion of intensive agriculture. Leptodyctium riparium was also associated with a low proportion of natural area within the catchment, with an optimum of 33.57%. The case was similar in highly seasonal rivers in Bulgaria, where the species was characteristic for sites located in regions with increased intensive agriculture and watercourses with reduced flow and pronounced silting, as well as elevated total nitrogen concentration [20]. Both species had a higher frequency in the Pannonian than in the Dinaric Ecoregion of Croatia, which is clearly related to the characteristics of its water bodies; these are more frequently slow and eutrophic lowland streams, rivers and canals with unstable sediment [64], which are known as less suitable habitats for bryophytes and support only modest diversity [10,15,65]. Here, the aquatic form of Riccia fluitans was recorded floating mostly in stagnant waters of hypereutrophic artificial canals, while Leptodyctium riparium was growing on rarely present large stable rocks, dead wood, periodically submerged tree bases and margins of the watercourses, having the highest frequency in lowland small watercourses, followed by lowland medium and large watercourses. These watercourses naturally have higher trophic status and the vast majority of them are additionally subjected to substantial changes in land use associated with high nutrient input, as well as hydromorphological degradation [66], which are known to reduce habitat quality for bryophytes, resulting in reduced cover, diversity and changes in community structure [14,17,20,22]. Thus, as expected, freshwater bryophytes which occur in higher frequency within this region are those which can tolerate poorer water-quality and inhabit less-natural catchment areas, such as Leptodictyum riparium, Cratoneuron filicinum, Cinclidotus riparius, Rhynchostegium riparioides, Fontinalis antipyretica and Fissidens crassipes.
As the data used in this study were gathered in the course of the national macrophyte monitoring conducted for the purpose of assessing the ecological status of waterbodies as required by the WFD, we want to emphasize the importance of its implementation, as it encouraged research into freshwater bryophytes and their ecology on a national [15,22,[67][68][69][70][71][72] and European level [1,2,4,5,14,28,29,50,51,53] by including this group as a part of macrophyte vegetation. Our research is the first into the ecology of the aquatic bryophytes in Croatia, exploring the ecological responses of the most frequent species and determining the influence of different environmental variables on their occurrence. These results make a solid base from which the bioindication potential of these particular species can be inferred, based on their optima and niche width, and a starting point crucial for the improvement and adjustment of the national methodology regarding the bryophytes as an integrative part of macrophyte vegetation. It should be noted that WFD takes a more holistic approach than traditional monitoring practices and requires the assessment of the ecological status, which must be determined type-specifically. Namely, for each type of water body recognized by the national typology, reference conditions should be identified, and degradation has to be quantified as the deviation in species composition and abundance from those that would be present at reference conditions [73]. Having this in mind, our results are a good starting point, because they add information on species currently not included and additional information on ecological responses for a few species already included in the Croatian methodology. However, particular species scores and the inclusion of new species into this methodology should be derived considering the type-specific reference conditions to meet the requirements of the WFD.

Study Area
A total of 648 sampling sites on 382 different watercourses were surveyed during vegetation seasons from 2016 to 2021. Surveys were carried out within the national surface water monitoring scheme, which is conducted to assess the ecological status of water bodies as required by the Water Framework Directive (WFD). The sampling sites were preselected to encompass the heterogeneity of different water body types recognized by the typology developed as a basis for the monitoring of surface waters [60], with all water body types represented proportionally and fulfilling the requirements of the stratified sampling. This typology recognizes two hydrological and biogeographical regions in Croatia-the Pannonian and the Dinaric Ecoregion, the latter being subdivided into Continental and Mediterranean subecoregions (Figure 1).
The Pannonian Ecoregion refers to the continental, lowland part of the country, largely converted into agricultural areas. The geological bedrock is dominantly siliceous, while the climate is temperate, without a dry season, with warm summers (Cfb), becoming hotter towards the east (Cfa) [74]. The Dinaric Ecoregion refers to the central and western part of the country, with a dominant karstic landscape, developed on limestone and dolomite bedrock. It is divided into the Continental Subecoregion, characterized by a temperate climate (Cfb), and the Mediterranean Subecoregion with mostly Mediterranean climate, with dry and hot summer months (Csa) [74]. The Pannonian watercourses and the majority of the watercourses of the Dinaric-Continental Subecoregion belong to the Black Sea Basin, while the watercourses of the Dinaric-Mediterranean Subecoregion belong to the Adriatic Sea Basin.

Vegetation Data Sampling
Macrophyte vegetation was surveyed from June to September, during the main vegetation period and the lowest water discharge levels. Following the national methodology for macrophyte sampling [60], watercourses were surveyed along 100 m-long transects from the banks and by zigzagging across the channel if the water depth was low enough. The vegetation survey included all macrophyte representatives (bryophytes, vascular plants and macroalgae), and the cover and abundance of each species was assessed using the standard Central European phytocoenological methodology, i.e., extended Braun-Blanquet scale (r = one individual, + = up to 5 individuals, 1 = up to 50 individuals, 2m = over 50 individuals but coverage < 5%, 2a = coverage 5-15%, 2b = coverage 15-25%, 3 = 25-50%; 4 = coverage 50-75%; 5 = coverage over 75%) [75][76][77], which was further transformed to the van der Maarel scale from 1 to 9 [78] (Appendix D). To investigate the ecological preferences and autecology of freshwater bryophytes, further analysis included only bryophytes with ≥5 occurrences that fall into categories of greater water affinity according to Dierßen [56]. These were collected from various substrates (e.g., rocks, boulders, pebbles, xylal) within the riverbed, as well as from the periodically flooded river margins. Voucher specimens were deposited at the Herbarium collection ZA [79]. The nomenclature follows Hodgetts et al. [80].

Environmental Data Sampling and Acquisition
All localities were also sampled for basic water physicochemical and chemical analysis once a month throughout the year. Water temperature, electrical conductivity, pH and dissolved oxygen were measured in situ with a Hach HQ40D Portable Multi Meter under standard conditions. Furthermore, water samples were collected and analyzed in an accredited laboratory (Central Water Management Laboratory, Zagreb) for total alkalinity, total suspended solids, biochemical oxygen demand and chemical oxygen demand, as well as for nitrogen and phosphorus compounds (ammonium, nitrites, nitrates, total nitrogen, orthophosphates and total phosphorus) ( Table 1). The land use in the catchment area of each sampling site was obtained from the database of Hrvatske vode-the legal entity for water management. Here, four distinct categories are recognized and calculated from the CORINE land cover dataset [81]-natural area, urban area, and intensive and extensive agricultural land (Table 1).

Data Analysis
Data analysis included rheophyte, hydrophyte, amphyphyte and hygrophyte species [56] occurring in at least five of the surveyed localities (a total of 21 species from 182 localities) (Table S1) matched with 18 environmental variables (Table S2). This data set included thirteen different types of watercourses according to the current national typology (Appendix A, Table A1) and was used to compile a frequency table from the species occurrence within each type (Appendix A, Table A2). Basic descriptive statistics (min, max, mean, SE, SD and median) of all environmental variables were calculated for the species (Appendix C) in Past 4.9 software [82] and their distribution along the gradient of each environmental variable was shown with box-plot graphs created in SPSS software ( Figure 1). Furthermore, the descriptive statistic was calculated for all environmental variables for the Pannonian and Dinaric Ecoregion (Table A6, Appendix E), as well as for the Dinaric-Continental and Dinaric-Mediterranean subecoregion (Table A7, Appendix E).
To assess the relationship between the environmental variables and patterns in freshwater bryophyte species composition, a direct ordination method, canonical correspondence analysis (CCA), was used. After removing the outliers, vegetation and environmental data from 176 localities were included in the analysis. CCA was selected because the response data were compositional with a gradient longer than 4.2 SD units, meaning that analysis based on a unimodal, rather than the linear model, is preferred [83]. A step-forward selection procedure in CANOCO 5 [83,84] was used to identify the most-contributing subset of environmental predictors influencing the freshwater bryophytes. Eight variables with the highest conditional effect and with a 5% significance cut level (p < 0.05; Monte Carlo test, 499 permutations) were included. Prior to the analysis, species abundance values were square-rooted and rare species downweighted.
Generalized additive models (GAM) were employed to model the probability of occurrence of individual bryophyte species as a function of eight environmental variables highlighted in the CCA analysis. GAMs were selected as an efficient tool in ecology since they do not require an assumption about the shape of species response along the environmental gradient [83][84][85][86]. We used Poisson distribution with log link function and df = 2 in fitting the species response curves and Akaike Information Criterion (AIC) available in CANOCO 5 [83,84] to select the best model. AIC considers not only the goodness of fit but also selects the most parsimonious model. Species for which no candidate model had an AIC value lower than the null model were automatically detected and removed by this procedure. Furthermore, only statistically significant models (p < 0.05) were retained and shown in graphs (Table S3).

Conclusions
The present study revealed that freshwater bryophytes and their assemblages were segregated along the gradients of the water chemistry and the proportion of natural and urban area within the catchment. The two latter variables represent the degree of combined stress, because they are often related to the extent of eutrophication, and pollution in general as well as hydromorphological degradation, of the watercourses. Furthermore, the ecological responses of individual species were examined to determine their optima, degree of tolerance and bioindication potential regarding the studied variables. The results showed that most of the investigated species preferred natural, clean, well-oxygenated, oligotrophic watercourses, with low organic matter content and electrical conductivity. However, the widely distributed and most frequent species, such as Fontinalis antipyretica, Rhynchostegium riparioides, Cratoneuron filicinum, Fissidens crassipes, Cinclidotus fontinaloides and C. riprarius, showed wide ecological tolerance to studied water chemistry variables and are thus not reliable bioindicators concerning these variables. On the other hand, species such as Palustriella falcata, Eucladium verticillatum, Dichodontium flavescens and Jungermannia atrovirens had narrow ecological niches and were restricted to pristine watercourses. Riccia fluitans and Leptodyctium riparium were the most obviously separated from the rest of the species, being the most tolerant to poor water quality. These species had wide ecological ranges, but preferred neutral, hypereutrophic waters with high nutrient and organic content, as well as electrical conductivity. Furthermore, they were frequently associated with a higher share of intensive agriculture and a low share of natural land within the catchment.
The ecological responses of several species obtained from our study do not perfectly correlate with previous findings. This might be a result of different methodological approaches concerning data collection and analysis, differences in the gradient length encompassed within each study, or the existence of different ecotypes of particular freshwater bryophytes, which display different ecological behavior in different geographical areas. Nevertheless, our study covered a considerable geographic area, included many different types of watercourses, from ground-fed streams to eutrophic large rivers, and thus encompassed substantially long gradients of the environmental parameters investigated. This has provided new data on the ecology and bioindication potential of freshwater bryophytes, contributing to the existing body of knowledge on both subjects.

Supplementary Materials:
The following supporting information can be downloaded at: https: //www.mdpi.com/article/10.3390/plants11243451/s1, Table S1: Cover and abundance of the species at each surveyed locality with coordinates of localities in WGS84; Table S2: Yearly mean values of 18 environmental parameters for all surveyed localities, Table S3: Results of the generalized additive models (GAMs).

Data Availability Statement:
The data presented in this study are available on request from the corresponding author.

13.
Artifical waterbodies 10 <200-500 10-10,000 small, medium, large Table A2. Species frequencies across different types of watercourses. For abbreviations of species' names see Figure 2, and for types see Table A1 of Appendix A. Appendix C Table A4. Descriptive statistics of all environmental variables for 21 freshwater bryophytes. For abbreviations of species' names see Figure 2, and for environmental variables see Table 1.