Hydrology and parasites: What divides the fish community of the lower Dniester and Dniester estuary into three?

Estuarine fish communities are the common object of ecological studies that aim to infer the environmental changes occurring in estuarine environments and connected catchments. However, the high variability of abiotic factors inherent in estuarine water areas often makes the effects of eutrophication hardly distinguishable from those generated by other natural and/or anthropogenic fluctuations of hydrological and hydrochemical parameters. In addition, estuarine environmental zones of primary marine, fresh water, and brackish water characteristics can support diverse species communities with different sensitivity to various environmental factors. In the present work we clarified the structure of fish fauna in the lower Dniester and Dniester estuary area (northwestern Black Sea coast) and assessed its changes as a response to various environmental factors. We applied hierarchical clustering with successive indicator species analysis, and Random Forest analysis on the data on fish community parameters (species richness, diversity, and evenness) to discover that the area under study contains three separate fish communities, each with a different response to fluctuations of environmental factors. Subsequent examination of parasitological data readily confirmed the conclusion about presence of three fish communities. There is no simple universal indicator, such as abundance of a single species in any of the three communities, which can be used to monitor the eutrophication impact in the area of research. Instead, further analyses revealed the species diversity as a perspective indicator of eutrophication with nitrates in this water area.


Introduction
The species community level of biological organisation has been suggested as the most suitable for the complex assessment of environmental factors (Adams et al., 2000;Attrill and Depledge, 1997;Levin, 1992) with particular advantages to be gained by using fish species in these types of studies (Warwick, 1993). Therefore, fish communities have been repeatedly used in estuarine aquatic ecosystems research (Jaureguizar et al., 2004;Lobry et al., 2003;Selleslagh et al., 2009).
Estuarine environment is naturally characterised by a high variability of abiotic factors, thus the effects of eutrophication in this water area are often difficult to distinguish from those of other natural stresses (Elliott and Quintino, 2007). Continuous observations suggest that estuarine systems, where one or several characteristics (physical, hydrochemical, hydromorphological, etc.) underwent sub-threshold deterioration, can still support native faunas, but become much more vulnerable to changes in other parameters, especially those caused by anthropogenic pressure, which rapidly shift them to alternative types of communities. For European waters this was clearly illustrated by comparative description of several estuarine and riparian ecosystems in Baltic Sea (Schiewer, 2008); similar comparative observations were also published for other regions (Malherbe et al., 2010). This may explain the apparent discrepancy in conclusions that arise from research of different European estuarine ichthyofaunas inhabiting areas of various stress level caused by natural and/or anthropogenic factors. In several studies authors have found that hydrochemical parameters such as temperature, salinity, nutrient and oxygen concentrations are the main predictors of fish abundance and species richness in a given estuary (Araújo et al., 1999;Marshall and Elliott, 1998;Snigirov et al., 2012;Thiel et al., 1995). However, other researchers have concluded that estuarine fish abundance cannot be predicted based only on physico-chemical variables (Maes et al., 1998;Power and Attrill, 2003; https://doi.org/10.1016/j.ecss.2018.11.022 Received 9 July 2018; Received in revised form 19 October 2018; Accepted 15 November 2018 Power et al., 2000a, b). Other studies have even associated fluctuations in estuarine fish assemblages with the phasing events of different periodicity that regulate spawning times and larval scattering rather than with variations in hydrological conditions (Potter et al., 1986(Potter et al., , 2001.
The estuarine environment includes zones that display primary marine, fresh water and brackish water characteristics and thus supports a variety of spatially distributed species communities (amphihaline, euryhaline and true estuarine), which play a vital role in the functioning of adjacent marine and river biocoenoses (Maes et al., 1998;Methven et al., 2001). On the other hand, due to the relatively small number of truly estuarine fish species, it is widely accepted that variability of the estuarine fish faunas is determined to a large extent by the recruitment of species from rivers and/or from the sea (Potter et al., 1986(Potter et al., , 1997. Therefore, a given estuarine habitat can potentially be divided between several communities with different modes of response to the local set of natural and anthropogenic factors.
Thus, maintaining a clear distinction between estuarine fish communities based on their response to environmental factors is vital for the environmental management of estuarine areas and for the prediction of trends in local ecosystems.
Biological diversity in fish parasites was repeatedly shown to be a reliable indicator of ecosystems wellbeing when used in combination with other quantitative methods (D'Amelio and Gerasi, 1997;Overstreet, 1997). Fish parasites commonly have complex life cycles depending on the presence of a wide variety of intermediate hosts such as members of zoobenthos, zooplankton, and vertebrates. In addition, life cycles of many parasite species possess free-living stages. Therefore, the dynamics of food web interactions have a significant impact on a fish parasite community, and such a community could potentially become a sensitive indicator of trophic structure, environmental stress, and biodiversity in a given area. With this knowledge, one should combine analyses of parasitological and environmental data in order to conduct a detailed research of fish community structure.
Eutrophication in the Dniester transboundary catchment, shared between Ukraine and Moldova, has been accepted as a main factor shaping the ecosystem of the lower Dniester and Dniester liman (Shevtsova, 2000). Further, the conditions of Dniester aquatic biological communities were shown to correlate with the amount of chemical pollution present (Lebedynets et al., 2005). Although the role of nutrients discharge in the productivity of estuarine and coastal fish communities has long been debated (Hanley, 1990;Philippart et al., 2007), eutrophication has repeatedly been shown to increase primary and secondary production (Beukema and Cadée, 1986;Colijn et al., 2002) and, as a consequence, fish biomass (Kuipers and van Noort, 2008;Teal et al., 2008). The combination of these factors makes the lower Dniester ecosystem an attractive monitoring point for developing a characterisation of the ecological situation in the whole Dniester basin.
The area studied in this work includes the Dniester liman and lower Dniester (see Fig. 2). Several decades of observation on the northwestern part of the Black Sea yielded a value of 2-4 cm as a mode of diurnal tide amplitudes, with a 18 cm prediction for the hypothetical 100-years maximum (Medvedev et al., 2016). These values are 1-2 orders of magnitude lower than characteristic amplitudes in micro-tidal estuaries (1-2 m) where tidal effects on ecosystem become significant (McLusky and Elliott, 2004). This makes negligible a tidal impact in an area under research. Therefore, the lower Dniester (upstream from the Dniester liman) is a fluvial area, whereas Dniester liman corresponds to the definition of true estuary: "Estuary is a semi-enclosed water body connected to the sea as the tidal limit or the salt intrusion limit and receiving freshwater runoff, recognising that the […] tidal influence may be negligible" (Wolanski and Elliott, 2015). However, lower Dniester and Dniester estuary (liman) share a number of fish species (see Table 1), which justifies their study as a single ichthyological complex.
Lower Dniester and Dniester Estuary (see Fig. 2), is one of the Table 1 Taxonomic composition of fish fauna of the Dniester estuary area. 1887: by Brauner 1887, 1949: by Berg 1949, 1968: by Zambriborshtch 1953, 19681992: by Sirenko et al. 1992, 2018: our data. "-" -Species not found; "+"species found; "?"species repeatedly reported by fishermen to be present in the area, but have not been detected in our catches. § -Introduced non-European species.  largest estuarine ecosystems on the Black Sea coastline with an area of ∼580 km 2 . Throughout the course of more than a hundred years starting from the end of the 19th century, this territory attracted the attention of many ichthyologists (Berg, 1949;Brauner, 1887;Sirenko et al., 1992;Zambriborshtch, 1953Zambriborshtch, , 1968. The effective management of the Dniester estuarine ecosystem is indeed critically dependent upon a deep understanding of ecological processes, a detailed characterisation of fish fauna, and the selection of adequate feedback indicators and operational approaches. However, to the best of our knowledge, there are still no quantitative studies of the fish community structure and its environmental determinants, not to mention the selection of key fish species (or species complexes) for environmental monitoring. Therefore, in our study we aimed to analyse the structure of fish fauna in the Dniester Estuary area and test a hypothesis about its heterogeneity, i.e. that it consists of significantly different species communities, rather than representing a unified complex. If the heterogeneity hypothesis is true, we sought to uncover the key indicator species for each community, as well as establishing which environmental factors shape the structural parameters of the fish fauna (species richness, diversity and evenness). Furthermore, we intended to identify which fish species, or parameters of the fish community structure, are most sensitive to eutrophication and, therefore, should function as indicators of this type of ecological impact. To accomplish these objectives, we set out to analyse fish community parameters in a combination with environmental and parasitological data.
Assessment of biological diversity in large and complex ecosystems with high-amplitude fluctuations of environmental parameters, as inherent to estuarine areas, ultimately requires simplification, such as clustering of species community into functional groups. These groups comprise species with similar response to environmental factors rather than with common evolutionary origin, and facilitate assessment of aquatic ecosystem diversity due to the simplified requirements for modelling (compared to modelling of the whole ecosystem). Such groupings can be used for the numerous types of tasks: to identify specialists (Dehling et al., 2016), to predict prey selection (Spitz et al., 2014), and/or analyse the use of particular habitat (Franco et al., 2008). However, such an approach (split of the fish community of given estuary into different groups or guilds) was applied to a plethora of different environments, ecosystem types and climate zones. This led to increasing overlap and confusion between different studies in a context of applicability of particular clustering techniques . Indeed, development of the universal method of fish community clustering seems to be unrealistic.
The classical approach to multivariate analysis of species community structure implies the following framework with pre-defined task hierarchy: (i) the display of community patterns through species clustering and ordination of samples; (ii) identification of species principally responsible for determining sample groupings; (iii) multivariate statistical tests for differences on parameters under research; and (iv) the linking of community differences to patterns in the physical, chemical and biological environment (Clarke, 1993). Such a hierarchy makes correct clustering a pivotal point for all subsequent analyses.
Therefore, in our work we have chosen the following approach to the stage (i) of data analysis. To assess clustering correctness, we applied two different methods of community clustering, which use mutually irreducible mathematical techniques: Hierarchical Cluster Analysis (HCA) with successive indicator species analysis, and Random Forest analysis (RF). The subsequently generated confusion matrix allowed quantification of differences between HCA and RF outputs, thus characterising correctness of community clustering. Next, proceeding within the same paradigm, we tested the correctness of initial HCA output by analysis of an alternative set of input data: fish macroparasite fauna rather than fish species abundance.

Sampling procedures
The Dniester Estuary area was tested at 29 sampling points, at depths from 0.5 to 2.5 m, during the warm period (May-October) for twelve years from 2006 to 2017. Each field sample consisted of fish catches, in situ measurements of water temperature, salinity, pH and   Data on the Dniester outflow volumes for the study period were obtained from the database of the hydrology station in Tighina, Moldova.
We have placed sampling points as to have as far as possibly full gradient coverage. Fish catches were performed with a moving seine (30 m long, 8 mm mesh size) and static gillnet (30 m long, 13 mm mesh size) at the same sampling points. Fishing gear was pulled for 30-50 m parallel to the coastline. Two sampling nets were used every time at every sampling. Therefore, we used as a unit of fishing effort sample collected by two fishing nets of different types at one sampling point. After identification and counting, most of the fishes were released, except 5-10% of individuals representing all collected species taken for the study of parasite infestation (if less than 10 individuals of the given species collected, 1 individual was taken); these fish were stored alive in water tanks for further dissection.
To analyse parasite community, the sampled fish were transported alive in aerated containers to the laboratory, then sacrificed and dissected within two days (Kvach et al., 2016). The fish were subjected to a full standardised parasitological dissection of fins, skin, gills, muscles and internal organs. Study of microparasites was avoided due to specific methods needed. Monogeneans were preserved in glycerine-ammonium-picrate as semi-permanent slides (Malmberg, 1957). Digeneans, cestodes and nematodes were preserved in hot 4% formaldehyde (Cribb and Bray, 2010). Then, digeneans and cestodes were stained with iron acetic carmine and mounted to slides in Canada balsam (Georgiev, 1986). Acanthocephalans were pressed between two slides with a 70%-ethanol drop, then preserved in 70% ethanol; glochidia and arthropods were preserved in 4% formaldehyde without specific fixation. The preserved acanthocephalans and nematodes were then mounted in glycerol as temporary slides. The acanthocephalans, nematodes, glochidia and arthropods were identified under light microscopy. All parasites were identified to species level or to the highest possible taxa. Data for each species are presented as abundance in accordance to (Bush et al., 1997).
In cases of doubtful initial identification of fish species, typical fish individuals were fixed with 4% formaldehyde and stored for further examination. Taxonomic identifications were made in accordance with the widely used field guides for European fish species (Froese andPauly, 2006-2018;Jennings, 1996;Maitland and Linsell, 2006;Miller and Loates, 1997;Nelson, 1994). The geographical coordinates of the sampling sites were controlled using a GPS navigator.

Data analysis
Since in our study area the abundance of species differed by twothree orders of magnitude, we fitted the species accumulation curve with a modified Clench's equation which was shown to be most adequate for this type of community (Soberon and Llorente, 1993): where N is the number of species detected, n is the number of fishing effort units, and a and b are fitting constants; wherefrom the total number of species in fish fauna was estimated as n To distinguish and characterise groups of sampling points with similar fish communities, we performed hierarchical cluster analysis and indicator species analysis in PC-ORD 5.2 (McCune and Grace, 2006). Groups of sampling points were first defined by hierarchical cluster analysis (HCA) using group-average linking based on Bray-Curtis similarities on quantitative catch data sets (Baldoa and Drake, 2002). Two-to five-group solutions were explored, however, higher than fivegroup solutions were not used further since they included groups with < 4 sampling points, which would be problematic in further analyses of relationships with ecosystem characteristics. To determine the significance level of clustering we performed 999 Monte-Carlo permutations. As a result, the three-group clustering had the highest level of significance.
Indicator species analysis was then performed for two-to five-group solutions by calculating the indicator value I v of each species for each group. A perfect (I v = 100) indicator species would be present in all of the samples in one group and in no other sample (McCune and Grace, 2006). The significance level of each indicator species was then determined by 999 Monte-Carlo permutations. Again, the three-group solution had the highest number of significant (α < 0.05) indicator species which therefore provided evidence of the greatest difference between groups. Thus both hierarchical cluster analysis and indicator species analysis suggest a three-group clustering as the most significant grouping of sampling points in the study area.
To reveal associations between hydrological characteristics and fish communities separated by the indicator species approach we used RF analysis (Breiman, 2001). This is a modification of classification and regression trees with predictions generated by creating multiple trees, each based on a bootstrapped subsample of the data and a random subset of predictor variables at each node of all trees. The advantage of this method is a reduced effect of potential spatial autocorrelation between ecosystem characteristics compared to other modelling approaches (Marmion et al., 2009). In our study the categorical response variable was a community type as defined by hierarchical cluster analysis, with ecosystem characteristics as predictor variables. RF was implemented with a Random Forest package in R (Liaw and Wiener, 2002). The default settings were used for m try , the number of predictor variables available for selection at each node (m try = 3; the square root of the total number of predictor variables rounded to the closest integer) and 9999 trees were generated. Classification accuracy was estimated by comparing the predicted fish community for each sample point to the actual community. The importance of the environmental variables was estimated by the mean decrease in accuracy of the model when each variable was randomized (Parravicini et al., 2012).
The Canonical Correspondence Analysis (CCA) was used to investigate the relative importance of environmental variables that determine the value of species abundance (Braak and Verdonschot, 1995;ter Braak, 1986). The method was chosen since it is less susceptible to nonlinearities in the data matrix than other ordination techniques, and because it explicitly assumes nonlinear (Gaussian) distributions of species quantity along the environmental gradients (Braak and Verdonschot, 1995). CCA was applied to the overall fish data matrix (dependent set) and the environmental data matrix (independent set) in order to clarify the relationships between biological assemblages of species and the environmental variables. The XLStat software provides the "Forward Selection" function with the built-in "Monte-Carlo permutation test" that ranks the environmental variables separately. The fit of a parameter is the eigenvalue of the CCA with this parameter as the only variable. In this procedure, the variable with the highest eigenvalue is separated, and at each further step the next variable is selected, which adds most to the explained variance of the species data. At each step, it is tested by 1000 random Monte-Carlo permutations to elucidate whether the parameter added is statistically significant. We set the limit at 1000 as this number of permutations yields an effective approximate test in that algorithm (Edgington, 1995). Only environmental characteristics identified by the Forward Selection algorithm as significant (α < 0.05) were included into CCA. CCA biplot was constructed to provide a visual assessment of the relationship between fish species and environmental characteristics. To assess significance of impact of environmental factors partial F-ratio (pF) for individual factors was calculated according to Verdonschot & Braak (Verdonschot and Ter Braak, 1994).
For graphical representation of the variations in parasite assemblage structure between different fish communities, cluster analysis was used. Initially, the Bray-Curtis dissimilarity index was calculated based on density data: specifically, the dissimilarity of species assemblage structure between communities. The species density values were obtained as a ratio of number of species per community to a number of sampling points related to this community (Gotelli and Colwell, 2001). Then, a dissimilarity matrix was generated and subjected to an average linkage clustering method in order to generate a dissimilarity dendrogram for communities (Horinouchi, 2009). Since in our study area the abundance of species varied by more than two orders of magnitude, for cluster analysis the abundance values were introduced as log (n p +1), where n p is the total number of detected individuals of a given species. Next, to quantify the difference between fish communities we built Bray-Curtis dissimilarity dendrograms based on the fish species density and parasite species density data pooled for three detected fish communities.
The fish community parameters: species richness, diversity and evenness were estimated as a summary of distribution and abundance of fish species. Species richness was considered to be the total number of species recorded in each community, and diversity was evaluated using the Simpson complement index, i.e. 1-D, where Here, S is the number of species, N is the total number of fishes caught, and n is the number of individuals of each species. We have chosen this index since it is less sensitive to sample size than other biodiversity indices (Lande et al., 2000;Smith and Wilson, 1996). To assess species evenness, we used the Simpson evenness index which has been shown to be independent of species richness (Smith and Wilson, 1996).
To examine the relationship between fish community parameters as dependent variables and environmental factors as independent variables, we used partial least squares regression analysis (PLSR) (Carrascal et al., 2009) implemented in XLStat-2014. For this analysis we used environmental factors which were found in CCA to cause significant impact on a fish community. This method allows treating, as a whole, multiple independent variables often themselves related, which maximizes the explained variability in dependent variables. PLSR is especially useful when the independent (predictor) variables are highly correlated i.e. strong collinearity is observed (Carrascal et al., 2009). In this analysis we used first two components which explained more than 5% of original variance of the dependent variable. From each simulation, we gathered the explanatory capacity (R 2 ) and the weight of each predictor within each component, which allowed us to understand the latent factors defined by each component. Therefore, the sum of the R 2 s of the two significant components gave the total explanatory capacity of the PLSR models. The addition of the squares of the statistical weights within each component sums to one, so the contribution of each predictor variable to the meaning of each component can be clearly estimated.
To find out if any fish species in the study area could be used as an indicator of eutrophication we applied similar PLSR analysis with each species abundance as a dependent variable.
Student's t-test and Pearson's correlation coefficient were calculated with Mathematica 10 software and applied as indicated.

Defining species communities
In our field work we found 74 species (Table 1). Fitting of the species accumulation curve yielded 87 as an overall species number for the study area with a coefficient of determination R 2 = 0.998. Values of fitting coefficients in the Clench's equation: a = 3.466, b = 0.0396; see Fig. 1. Taking into account that five more species were repeatedly reported by local fishermen to be present in the study area (Table 1), this looks as a plausible estimation. However, we did not include the aforementioned five species into the statistical analyses.
When we integrated our sampling data set for the full period of research, three groups of sampling points with similar fish communities were identified with HCA (Fig. 2) and described by significant indicator species (Table 2).
HCA based on group-average linking showed a highest level of significance (α = 0.031) when three-group solution was chosen. Similarly, the three-group solution had the highest number of significant (α < 0.05) indicator species which therefore provided evidence of the greatest difference between groups; thus the three-group solution was selected for further analyses.
The "Reedbeds" community was most widespread (found in 11 sampling points) and included Abramis brama, Rutilus rutilus, and Neogobius gymnotrachelus as significant indicator species. There were 52 species found in this community. The "Open-freshwater" community was also widespread (10 sampling points) and included Misgurnus fossilis and Cobitis rossomeridionalis as significant indicator species. This community, being inhabited by 61 species, was found to be most species-rich. The remaining "Open-brackishwater" community was smallest in terms of the number of sampling points (8) and total number of registered species (40), but had the highest number of indicator species: Atherina boyeri, Engraulis encrasicolus, Mesogobius batrachocephalus, and Syngnathus typhle.
Species which were not defined as significant indicators nevertheless provided additional insight into differences among the fish communities under study. The Open-brackishwater community was relatively species-poor, being mainly composed of marine and diadromous species; this community also had the highest number of species absent in the two other communities (Alosa maeotica, Gobius niger, Knipowitschia longecaudata, Liza aurata, Liza haematocheila, Neogobius eurycephalus, Percarina demidoffii, Platichthys flesus, Pungitius platygaster). The majority of species found in the two other communities were shared between them. However, species such as Abramis sapa, Blicca bjoerkna, Carassius carassius, Leuciscus idus, Sabanejewia aurata and Tinca tinca were detected only in the Reedbeds community. In turn, Pelecus cultratus and Romanogobio kesslerii were detected in the Openfreshwater community only. Several species (Carassius gibelio, Lepomis gibbosus, Syngnatus abaster) were found to be present in all three communities.

Characterisation of species communities
Upon being applied to the integrated data set for the full period of research, RF models correctly classified the fish community at 83% of sampling points (Table 3). The Reedbeds community had the smallest rate of misclassification, whereas Open-freshwater community spots were most commonly misclassified as having the Reedbeds community. Misclassifications often involved outlying sampling spots with communities which were different from the neighbourhood (e.g. Reedbeds community in Karagol bay, see Fig. 2). Variables associated with physical environmental characteristics (temperature, oxygen concentration) and eutrophication with NO 3 − ions had the greatest explanatory power (Fig. 3) with temperature as the most important factor. However,  . Dashed line denotes the start of this study. C: Fish communities by sampling point, as defined by hierarchical cluster analysis and described by indicator species analysis. D: Fish communities by sampling point, as predicted by random forest model. White circles: Reedbeds community; grey circles: Open-freshwater community; black circles: Open-brackishwater community; circles colour code and "5 km" scale bar apply to C and D. Area patterned in dots at panels C and D denotes reedbeds. the importance of other variables was much less clear. Therefore, to give more detailed characterisation of the interactions between ecosystem factors and the structure of fish fauna, we performed CCA on the data from the whole study area (Fig. 4A). Trends of the main environmental gradients were related simultaneously to first and second ordination axes. According to interset correlations, the most important environmental variables were temperature, oxygen concentration and salinity with smaller input of NO 3 − , PO 4 3and NH 4 + concentration. However, no clear distribution gradient was generated for the species under study. Therefore, we performed CCA separately for each of three species communities. According to interset correlations in both Reedbeds (Fig. 4B) and Open-freshwater (Fig. 4C) communities, the most important environmental variables were temperature and oxygen concentration. Together they formed the main sample ordination gradient from the origin of coordinates along their positive directions; however, a substantial number of species were found around the centre of the ordination biplot, implying relatively neutral relationships with ecosystem characteristics. In contrast, in the Openbrackishwater community salinity was clearly the main factor shaping the sample ordination gradient with a marginal effect from other environmental variables. Organic ions (NO 3 − , NH 4 + and PO 4 3-) were found to contribute significantly to the model in the Reedbeds community and Open-freshwater community (NO 3 − ) whereas in the Openbrackishwater community their effect on species distribution was not significant: see Table 4. To find one additional characteristic of the fish communities we analysed their tentative dissimilarity using metazoan parasites as biological markers. Altogether, in our catches we detected 31 parasite species (Table 5). Parasitofauna of the Dniester Estuary fishes revealed in our research was similar to that depicted in previous studies (Chernyshenko, 1960;Kvach, 2004); the only new species found in the fishes was a mite larva (Unionicola sp.) encysted in mesentery of Misgurnus fossilis. The Bray-Curtis dendrograms displayed similar structure with a convincingly higher value of dissimilarity index between the Open-brackishwater and the two other communities, than between the Reedbeds and Open-freshwater communities. Surprisingly, dissimilarity indices of the fish-species dendrogram were much lower than those of the parasite-species dendrogram (Fig. 5): 0.65 vs.0.84 for the difference between the Open-brackishwater and other communities, and 0.34 vs. 0.57 for the difference between the Reedbeds and Open-freshwater community, respectively.

Impact of environmental factors on community parameters
To clarify the influence of environmental factors on the community parameters (species richness, diversity and evenness) we performed a PLSR analysis separately for each of the three communities with factors that demonstrated significant impact on fish community in CCA. In all cases the two first significant components, being combined, explained more than 95% of the original variance in the response variable (Table 6). In the Reedbeds community temperature, oxygen concentration and concentration of NO 3 − ions had the highest weights shaping all three characteristics of fish community, whereas NH 4 + concentration had a relatively high impact on the species evenness. For the Open-freshwater community PLSR revealed temperature, oxygen, and concentration of all three organic ions included in the analysis as important factors that influence the community characteristics. In contrast, oxygen concentration had quite a little impact on the Openbrackishwater community characteristics.
To find out if any fish species in the study area can be used as a simple indicator of the eutrophication severity we repeated PLSR analysis separately for each species listed in CCA ordination plots, with abundance of this species as a dependent variable. However, no species demonstrated a significant correlation between its abundance and the concentrations of registered ions (P > 0.15 in all cases; data not shown).

Discussion
74 species detected in our field work and 87 as an estimation for the full number in study area, being compared with historical data (Table 1), suggest that the overall species richness was not lowered compared to that in XIX and XX centuries. This is to some extent counter-intuitive, taking into account a rapid increase of human population and, as a sequence, of anthropogenic pressure in study area during last 100-150 years. The plausible explanation, however, comes from the fact of arrival of several non-European species (see Table 1).
In this study we have found that fish fauna of the Dniester estuary and lower Dniester is divided between three communities along environmental gradients. Although some species were found to be ubiquitous throughout the study area (e.g. Carassius gibelio), the HCA supports the hypotheses regarding the presence of species communities that differ significantly amongst particular areas in the lower Dniester and Dniester estuary. The HCA-generated discrimination of fish fauna between the three communities was consonant with three other lines of evidence: RF algorithm, CCA, and analysis of parasite fauna. No species in the study area was a perfect indicator of a particular community, i.e. found in all sampling sites of one fish community and absent from all other sampling sites. This implies that the environmental characteristics work as a complex of filters affecting the likelihood of species presence rather than as strict determinants of species success in a given area.
An 83% overlap between the HCA and RF readouts (Table 3) Table 3 Confusion matrix displaying error estimates for RF analysis of the relationship between fish communities and environmental characteristics. Row headers represent the actual community according to HCA, and column headers represent the community predicted by RF. "Error" is the proportion of sampling points in each community which was misclassified. The "Community total" column lists the actual number of sampling points where the given community was detected. confirms correctness of our split of the fish fauna into three specific communities. This is, to some extent, an expectable result, since spatial distribution of the communities follows the apparent salinity gradient characteristic for estuarine systems (see Fig. 2C and D). More important, however, is a question about the nature of inconsistencies between HCA-and RF-generated community fragmentations; i.e. what factors may lead to discrepancies displayed in a confusion matrix. One of the discrepancy reasons may follow from the bootstrap re-sampling algorithm of RF, which is used for inputting of missing values. The bootstrap algorithm was shown to generate inconsistent results when estimated variances of randomized coefficients tend to zero (Andrews, 2000), or when two or more values are competing for the same rank in a regression tree (Hall et al., 1993). Both options become increasingly probable in a case of quantitative imbalance between primary sampling points, i.e. in inherent to field research situation when data for the some parts of the studied area are limited to few individual samples (McKenna Jr, 2003). This suggests prudence when RF or other bootstrap-based technique is used for analysis of field-collected datasets. On the other hand, comparison of RF-produced clustering with that generated with an alternative statistical approach (HCA in our work) provides a plausible accuracy estimation for other types of RF output, such as an explanatory power of environmental factors (see Fig. 3). The comparative importance of environmental variables as shaping factors of the local fish fauna indicated by RF was much less clear than that defined by CCA. This is consistent with previous observations regarding the high dependence of environmental variables' importance upon the Fig. 4. Ordination diagrams based on a CCA of species abundance throughout the study period. The significance of the effect of each environmental variable was verified by Monte-Carlo permutation test, P < 0.05 in all cases. Cumulative constrained variance is given on the graph axes. Length and direction of arrows represent the relative importance and direction of change in the environmental variables. Circles: fish species. Species codes correspond to those given in Table 1  tree-based ensemble analysis when fish community structure is studied (Knudby et al., 2010). In turn, such an observation leads to some uncertainty in interpretation of the importance of environmental variables when RF is in use. As exhibited by CCA (Fig. 4), temperature and oxygen concentration were the most important determinants in the Reedbeds and Openfreshwater communities (localized in the lower Dniester and the upper part of Dniester estuary, respectively), whereas the Open-brackishwater community (localized to south from the line between the towns of Bilgorod-Dnistrovskiy and Ovidiopol, see Fig. 2) was shaped mostly by the salinity gradient. Quite a low collinearity between the salinity and Dniester runoff, and comparatively low influence of the latter factor on species abundance revealed by the Open-brackishwater CCA biplot, probably indicates an important role of precipitation and wind-induced surges in the modulation of water salinity in the lower part of the estuary (Yurishinets and Romanenko, 2009). Lack of significant distribution gradient(s) at CCA biplot generated for the whole study area (Fig. 4A) suggests absence of common factors shaping all three species communities, and thus confirms indirectly the profound functional differences between them.
Surprisingly, variables representing mostly anthropogenic effects (eutrophication with organic ions) were not important predictors of fish community structure, though some were of secondary importance. The possible explanation for this observation is that anthropogenic pressure is a major determinant of fish community structure and sustainability in the areas where the mode of land use is dynamically changing (Leprieur et al., 2008;Meador and Goldstein, 2003). However, human-use intensity is of minor importance when its methods and means are similar at an interdecadal timescale due to the adaptation of water communities (Mehner et al., 2005). Therefore, the relatively low effect of anthropogenic eutrophication in our study area probably indicates their stable and relatively low affluent volumes during the last 15-20 years and/or low level of nutrients accumulation in Dniester estuary (Garkavaya et al., 2008).
In line with results generated by the RF and CCA algorithms, the analysis of fish parasite fauna of the Dniester Estuary area has confirmed the interrelation between three fish communities: the Openbrackishwater community was most unlike the other two. In turn, the Open-freshwater and Reedbeds were much less different from each other ( Table 5, Fig. 5). It is widely accepted that the distribution of fish parasites is determined to a large extent by hydrological characteristics of water masses (Esch and Fernández, 1993;Pietrock and Marcogliese, 2003). However, parasite communities in fish are also shaped by host behaviours, feeding habits, immunological responses etc. (Timi, 2007). Nevertheless, in the present study parasite fauna have shown quite a clear difference between three fish communities inhabiting water areas of different hydrological profiles. Moreover, the Bray-Curtis dissimilarity indices were higher in a dendrogram based on parasite species density data than in a dendrogram based on fish species abundance. There are two potential explanations for this observation. Firstly, this could be an effect of the high sensitivity of parasite fauna of the Dniester Estuary area to hydrological and hydrochemical factors associated with different water areas. Secondly, this may indicate a high level of separation between populations of the same fish species incorporated into different fish communities. To clarify this issue, further intensive parasite monitoring focused on fish species detected in all three communities (Lepomis gibbosus, Carassius gibelio) is needed.
In the present study, in all three species communities the stepwise regression models revealed the highest and most significant coefficients of correlation for the interaction between environmental variables and species diversity (Table 6). Again, this analysis displayed much higher similarity between the Reedbeds and Open-freshwater communities than between either of these and the Open-brackishwater community. The pattern of correlations suggests the use of the species diversity parameter of all three communities as a possible indicator of anthropogenic eutrophication by nitrates (NO 3 − ). This indicator can be derived from fish catchment statistics collected at different time intervals and can therefore serve as an assessment criterion for the temporal dynamics of anthropogenic pressure in the Dniester Estuary and in the whole catchment. Quantification of impact and clarification of action mechanisms of the factors that control species diversity in a given ecological community are of the key questions of ecological studies (Sommer, 2002). Thus a positive correlation between species diversity and NO 3 − concentration revealed by PLSR in our study, certainly requires interpretation.
Since ichthyofaunas of Dniester estuary undergo a substantial fishery pressure during more than a century, the possible explanation of such a correlation is a significant positive impact of nutrients' inflow on species diversity in aquatic communities which are under long-term anthropogenic pressure, as was repeatedly demonstrated by multivariate modelling for different ecosystems (Worm et al., 2002). However, low or moderate correlations observed in our study for the nutrients concentration and the community parameters of fish Table 5 Fish parasite species in estuarine communities of lower Dniester. "-" -parasite species not found; "+" -1-10 parasite individuals per fishing effort; "++" -11-100 parasite individuals per fishing effort; "+++" -> 100 parasite individuals per fishing effort.

Openbrackish
Acanthocephalus lucii - fauna are to some extent counter-intuitive, and thus rise a question about possible causes of this analyses outcome. The plausible explanation is that in temperate estuarine systems hydromorphological characteristics of study area and water discharge fluctuations in different periods of year cause a complex nonlinear modulatory impact of different type on integrated fluxes of different nutrients (Struyf et al., 2004). This impact may not be fully reflected by nutrients concentration in punctual samples, as collected in our study. If true, a massbalance model of water discharge and nutrients outflow, based on series of short-period hydrological data sets (Moore et al., 2006), is needed to reveal a detailed interrelation between fish community characteristics and amount of nutrients passing through the lower Dniester and Dniester estuary.
The key indicator of the current status and future perspectives of a given ecosystem is its stability, which has been shown to be critically dependent from species diversity (Ives and Carpenter, 2007). Therefore, the inflow of nutrients observed in our study may upregulate the species diversity in fish communities, thus being a vital factor supporting the overall ecosystem. However, large-scale studies on estuarine communities with different level of eutrophication but similar environmental conditions are needed to obtain a solid proof of this hypothesis. Nevertheless, sensitivity of local fish fauna to eutrophication, revealed in our study, suggests wariness when popular techniques of organic agriculture are installed in the catchment since these techniques imply a severe restriction on the use of fertilizers (Badgley et al., 2007).
The present analyses provide a first attempt to distinguish and characterise the environmental drivers of changes in the fish fauna of the Dniester Estuary area and to identify possible causal processes. Indeed, at the current stage of investigation the causes of these changes remain reasonably speculative. Division into several separate Table 6 Results of the PLSR carried out for fish community parameters (response variables) and predictor variables (environmental parameters) for three fish communities. CW 1 , CW 2 : weights of each variable in the first and second PLSR component, respectively. PLSR weights which squares are larger than 0.2 are shown in bold. communities, albeit sharing a number of species, alongside a direction from fluvial part of the river to the true maritime environment, is a common and most expectable output of ecological studies in estuarine areas (Araújo et al., 1999;Baldoa and Drake, 2002;Beukema and Cadée, 1986;Feyrer et al., 2015;James et al., 2013). Results of our fieldwork and subsequent analyses are indeed in line with this concept, but thus do not demonstrate unequivocally unique features inherent to local fish fauna in a context of response to environmental impacts.
Research of additional environmental parameters such as turbidity, current speed and sediments nature in the study area may reveal such unique features. On top of that, our research has clearly shown that no single environmental variable can explain the observed patterns; in turn, there is also no simple universal indicator, such as abundance of a single or a few species, which can be used to monitor the anthropogenic impact in the area of research.

Author contributions
S. Snigirov, I. Kvach and S. Sylantyev participated field work on fish sampling; A. Goncharov collected and interpreted hydrological data; I. Kvach collected, analysed and interpreted parasitological data; R. Sizo controlled localization of fish catchment points and prepared maps; S. Sylantyev performed numerical calculations and prepared the paper; S. Snigirov, I. Kvach, A. Goncharov and S. Sylantyev developed concept of the study and working hypotheses; S. Snigirov, I. Kvach and A. Goncharov provided critical review on the paper drafts.

Conflicts of interest
The authors declare that they have no conflict of interest.