© R. Ćuk et al., Published by EDP Sciences 2019

1 Introduction

Macroinvertebrate fauna of large European rivers has changed dramatically during the last two centuries (e.g. Bij de Vaate et al., 2006; Danube: Liška et al., 2008; Sommerwerk et al., 2009; Rhine: Nienhuis et al., 2002) as a result of multiple pressures, such as organic and chemical pollution, land use change, hydromorphological alteration and biological invasions (Petts et al., 1993; Tockner et al., 2009; Strayer et al., 2014). The interconnections of major European river basins via canals, transport by shipping, as well as intentional and unintentional introductions have synergistically facilitated the range expansions of great number of non-indigenous macroinvertebrate species (NIMS), many of which originate from Ponto-Caspian basin (Jazdzewski, 1980; Bij de Vaate et al., 2002; Leuven et al., 2009). Some of established NIMS spread and cause impact, i.e. become invasive and their impact may be interpreted as a decline in ecological quality resulting from changes in biological, chemical and physical properties of aquatic ecosystems (Olenin et al., 2007). Continuous rise of NIMS introductions, both in the number of species and pathways of introduction (Nunes et al., 2015), will further increase pressure from approximately 300 NIMS already recorded in European freshwaters (Strayer, 2010; EEA, 2012).

Gradual replacement of native macroinvertebrate fauna by NIMS in large European rivers have occurred since the beginning of the twentieth century (e.g. Van den Brink et al., 1990; Bij de Vaate, 2003 and references within). More recently, replacement of formerly established NIMS by mass invaders was observed in the Rhine (Haas et al., 2002; Bernauer and Jansen, 2006) and other European rivers (e.g. Jazdzewski et al., 2004). Consequently, most parts of major European rivers show high levels of biocontamination, i.e. they are dominated by NIMS in terms of abundance and taxonomic composition (Arbačiauskas et al., 2008, 2011; Arndt et al., 2009). Adverse ecological and economic effects of introduced NIMS are perceived as “biological pollution” sensu Elliott (2003). However, since quantitative effects of “biological pollution” are difficult to estimate (and usually are unmeasured in aquatic ecosystems), a more practical approach is to measure “biological contamination” (i.e. biocontamination), in which it is assumed that NIMS effect is proportional to their taxonomic and abundance proportion in composition of invaded assemblages (Arbačiauskas et al., 2008). Biocontamination of the aquatic ecosystem needs to be estimated as it represents a significant biological pressure (Arbačiauskas et al., 2008, 2011) threatening to the fundamental objectives of the EU Water Framework Directive (WFD) (e.g. Panov et al., 2009; MacNeil and Briffa, 2009). Therefore, it should receive special attention within the ecological status assessment required by the WFD (Cardoso and Free, 2008). A simple method of biocontamination assessment (“biocontamination index”) proposed by Arbačiauskas et al. (2008) is easily applicable alongside current routine water quality monitoring techniques required by WFD.

Benthic macroinvertebrates are the most prevalent biological group used in aquatic bioassessment (Friberg et al., 2006; Birk et al., 2012). Since the introduction of WFD in 2000, they represent one of five biological quality elements used in the assessment of ecological status of European aquatic ecosystems. A Multimetric Index was proposed for each biological quality element to assess a site's overall condition by combining different types of metrics, with a purpose of classifying a water body into one of five Ecological Quality Status categories ranging from high to bad (Hering et al., 2004,2006). However, it is not clear how to deal with NIMS during ecological quality assessments in European rivers (Orendt et al., 2010). On the one hand, NIMS may distort an assessment of true ecological status as they tend to be more tolerant of organic pollution than the natives they replace (MacNeil et al., 2012). Also, they can reduce taxa richness and distort assemblage composition by eliminating sensitive macroinvertebrate species which are indicators of higher ecological quality (Arbačiauskas et al., 2011). Several studies have shown that some of the frequently used indices (e.g. BMWP, ASPT and Biotic Indices) could provide unreliable ecological quality estimates when NIMS are abundant (Arbačiauskas et al., 2008, 2011; MacNeil and Briffa, 2009; MacNeil et al., 2010). MacNeil et al. (2012) argue that in water bodies subject to invasion, it may be necessary to rely solely on chemical water quality and disregard biological water quality assessments until there is a major overhaul of established biotic indices. On the other hand, in a review of the role of NIMS in ecological quality assessment, Orendt et al. (2010) concluded that they should be included in ecological quality and human impact assessments. Studies in Belgium (Gabriels et al., 2005; Boets et al., 2011) also support inclusion of NIMS. While practice of inclusion/exclusion differs between countries (Orendt et al., 2010), it seems obvious that there is no universal rule how to deal with NIMS in biological assessment, due to local specificity of established NIMS, their relative abundance and resulting different changes in composition of macroinvertebrate assemblages as well as their reaction to different pressures. Therefore, more studies are needed to test the suitability of using macroinvertebrate assemblages in large rivers, now mostly dominated by NIMS, to achieve reliable ecological quality assessments.

Large rivers in Croatia, located in Pannonian ecoregion (ER11) (Illies, 1978), are particularly exposed to invasion since they belong to the Danube catchment area and are, thus, part of the so-called southern corridor (Danube-Main-Rhine) (Bij de Vaate et al., 2002), which now represents the most important corridor for spread of Ponto-Caspian NIMS in Europe (Leuven et al., 2009). Except for studies focusing on particular groups of invasive species (Mollusca: Lajtner et al., 2004; Lajtner and Crnčan, 2011; Crustacea: Žganec et al., 2009, 2015, 2018; Maguire et al., 2011, 2018), only few studies examined macroinvertebrate assemblages of Croatian large rivers (Matoničkin et al., 1975; Paunović et al., 2012; Lucić et al., 2015). However, assessment of the biocontamination levels of macroinvertebrate assemblages in all major large rivers was never done. Also, the question of how NIMS effect ecological quality assessment remains unexplored and needs further research.

In this study, we aimed to (1) assess the composition and biocontamination of macroinvertebrate assemblages of four major large rivers in Croatia, (2) explore the relationship between key environmental (physicochemical, substrate composition) and spatial parameters and non-indigenous macroinvertebrate assemblages and (3) establish the relationship between biocontamination and biological metrics used in ecological quality assessment. By addressing these objectives, we aimed to provide a rigorous test of the practicalities of integrating assessment of biocontamination within established national monitoring programs for ecological status assessment of large rivers. Furthermore, by identifying key explanation factors responsible for the spatial variation of macroinvertebrate assemblages in large rivers, with focus on NIMS and biocontamination, we hope to contribute to better understanding of key processes which facilitate the spread of NIMS in Croatia.

2 Material and methods

Research was conducted at 48 sampling stations on four major large rivers in Croatia (the catchment area >10,000 km²): the Sava (SA) (21 sampling stations), the Drava (DR) (20 sampling stations), the Mura (MU) (3 sampling stations) and the Danube (DA) (4 sampling stations) rivers, all of which are situated in Pannonian lowland ecoregion (ER11) (Illies, 1978) (Fig. 1). Two out of three reservoirs at the Drava River were included in the research, both in the littoral zone and in the middle of the reservoir. According to Croatian national river typology, the sites SA1-SA5 at the Sava River (“Middle Sava”), all sites at the Mura River (MU1–MU3 = “Lower Mura”) and sites DR1–DR12 at the Drava River (“Middle Drava”) belong to “Very large lowland rivers − lower course of the Mura River and middle course of the Sava and Drava Rivers” (HR-R_5B). Sites SA6-SA21 (“Lower Sava”) at the Sava River and sites DR13–DR18 (“Lower Drava”) belong to the “Very large lowland rivers − lower course of the Sava and Drava Rivers” (HR-R_5C), while the studied reach of the Danube River is classified as “Very large lowland rivers − the Danube” (HR-R_5D) (Regulation on Water Quality Standard).

We carried out the sampling of benthic macroinvertebrates once in 2015 at each sampling station during a period of stable and low water levels. Samples were collected with a hand net (25 × 25 cm², mesh size 500 µm) according to AQEM sampling procedure (AQEM, 2002) in a littoral zone of rivers (44 sites) and reservoirs (2), where 20 quantitative replicative samples were sampled (0.0625 m²) from all microhabitats that covered more than 5%. Three replicate samples from the middle of two reservoirs on the Drava River were collected by Ekman grab (3 × 0.0225 m²). The collected material was preserved with 96% ethanol in the field and finally stored in 70% ethanol in the laboratory. We isolated and identified the benthic macroinvertebrates in the laboratory using a binocular stereomicroscope (Olympus SZX10). We identified all taxa to the lowest possible taxonomic level. All samples are deposited in the Central Water Management Laboratory of Hrvatske vode.

Water samples were collected on four sampling dates (March, July, October, December 2015) at each site and the following 27 physicochemical parameters were analysed according to standard analytical methods for assessment of surface water quality (ISO norms): water temperature (°C), pH, conductivity (µS cm⁻¹), total suspended solids (mg L⁻¹), alkalinity m-value (mg CaCO₃ L⁻¹), total hardness (mg CaCO₃ L⁻¹), dissolved oxygen (mg O₂ L⁻¹), oxygen saturation (%), biological oxygen demand (BOD₅) (mg O₂ L⁻¹), chemical oxygen demand (COD-Mn) (mg O₂ L⁻¹), ammonia (NH₄ ⁺) (mg N L⁻¹), nitrite (NO₂ ⁻) (mg N L⁻¹), nitrate (NO₃ ⁻) (mg N L⁻¹), total nitrogen (mg N L⁻¹), inorganic nitrogen (mg N L⁻¹), organic nitrogen (mg N L⁻¹), orthophosphate (PO₄ ³⁻) (mg P L⁻¹), total phosphorus (mg P L⁻¹), total organic carbon (TOC) (mg L⁻¹), dissolved organic carbon (DOC) (mg L⁻¹), calcium (mg L⁻¹), magnesium (mg L⁻¹), sodium (mg L⁻¹), potassium (mg L⁻¹), chloride (mg L⁻¹), fluoride (mg L⁻¹) and sulfate (mg L⁻¹).

Since measurements were not done for the bottom of the reservoirs, but only for the surface water samples, sites DR4 and DR7 (deep section of the reservoirs) were excluded from all analyses involving physicochemical parameters. Additionally, substrate composition was assessed at 46 sampling sites (excluding two sites at deep sections of reservoirs) using AQEM protocol (AQEM, 2002) and for analyses three aggregated fractions were used: hard substrate (mega-, macro-, microlithal and akal), soft substrate (psammal, psammopelal and argylal) and phytal. Instead of discharge data which were not available for each sampling site, size of subbasin was used for each site (calculated using GIS tools). To examine the effect of geographic location, geographic coordinates (X, Y Gauss-Krueger, zone 5), distance to the Danube (with four Danube sites having value 0) and altitude were used.

We assessed the biocontamination according to Arbačiauskas et al. (2008), i.e. a site-specific biocontamination index (SBCI) was derived from two metrics: an abundance contamination index (ACI) and a richness contamination index (RCI) at family rank. Obtained SBCI classifies the sampling site into one of five classes ranging from 0 (no contamination) to 4 (severe contamination).

We analysed the correlation between biocontamination indices and several biological metrics included in Croatian national methodology for ecological water quality assessment: the number of EPT taxa (EPT-S), the proportion of EPT taxa (%EPT), the total number of taxa, the total number of families, Biological Monitoring Working Party Index (BMWP) (Armitage et al., 1983), Altered Indice Biotico Esteso (IBE AQEM) (originally proposed by Ghetti, 1997, implemented in Asterics Software) and Croatian Indicator Saprobic System (HR-SI), which is based on Pantle-Buck Saprobic Index (Pantle and Buck, 1955), but with indicator values according to HRIS (Official Gazette, 2013) and absolute instead of relative taxa abundances.

We also calculated the saprobity module, which consists of two metrics: HR-SI and IBE AQEM. We used Asterics Software Version 4.0.4. for calculation of all biological metrics both with and without NIMS, except for the Croatian Saprobic Index (HR-SI) and Saprobity module, which were calculated separately in excel.

We analysed the spatial patterns of macroinvertebrate assemblages using non-parametric multidimensional scaling (NMDS) with software package PRIMER Version 6.1.13. and PERMANOVA+ Version 1.0.3 (PRIMER-E Ltd 2009). The average density of all taxa was square root transformed to control the influence of dominant species, and the Bray–Curtis index of similarity was used to calculate the similarity matrix.

Environmental variables for four months (March, July, October and December) in 2015 were averaged for each site, log-transformed and normalized and then analysed using principal component analysis (PCA) to examine spatial differences between sites, rivers and river sections. Before this analysis, Pearson correlation coefficients and scatter plots between all pairs of environmental variables were examined to eliminate co-linear variables and reduce redundancy, resulting in a subset of 18 selected environmental variables (of 27 in total) used for PCA. To examine the relationship between non-indigenous macroinvertebrate assemblages and environmental factors, Canonical Correspondence Analysis (CCA) was run using CANOCO 5 (Ter Braak and Šmilauer, 2012). Since preliminary analysis, using Detrended Correspondence Analysis (DCA) in CANOCO 5 revealed unimodal rather than the linear response of the species data, CCA was used instead of constrained multivariate analyses like RDA or dbRDA that assume linear response (Lepš and Šmilauer, 2003). To include seven sampling sites without NIMS in CCA, three native amphipods that were abundant were also included in CCA analysis. All variables were transformed (log (x + 1)) prior to analyses. Although 18 environmental variables and 24 species were initially used for CCA, finally a subset of 12 environmental variables (with inflation >5) and 14 species (11 widespread NIMS and 3 native amphipods, while 10 localized NIMS were excluded) were used for samples + species + environmental factors triplot. The statistical significance of the first (CCA1) axis and all axes in CCA were tested using Monte Carlo test (9999 permutations) to verify the significance of the models (Lepš and Šmilauer, 2003). Spearman rank (rs) correlations were used to test relationships between biocontamination indices, environmental parameters and biological metrics.

Fig. 1 Site-specific biocontamination index (SBCI-Fam, family level) values presented in different colours for all 48 studied sites sampled in 2015 along the Mura (MU), the Drava (DR), the Sava (SA) and the Danube (DA) Rivers. Right inset figure shows the position of Croatia in Europe with country codes (HR: Croatia, AU: Austria, BH: Bosnia and Herzegovina, HU: Hungary, IT: Italy, MG: Montenegro, SI: Slovenia, SR: Serbia), while left inset figure shows the enlarged area of the Drava reservoirs.

3 Results

3.1 Macroinvertebrate assemblages and biocontamination

Over 300,000 (315292) specimens belonging to 236 taxa of benthic macroinvertebrates were collected (37 orders and 92 families) at 48 sites (Fig. 1). Total number of taxa per river ranged from 51 in the Danube (average per site: 25) to 176 in the Drava River (average per site: 32) (Tab. 1), while total number of taxa per site ranged from 16 (at DR4 and DR7) to 48 (at DR2).

Non-indigenous species were represented with 21 species, belonging to the seven taxonomic groups (Tab. 2). The most NIMS were of the Ponto-Caspian origin (16), while representatives from other origins were scarce. Three NIMS have been recorded for the first time in Croatia: the flatworm Dendrocoelum romanodanubiale (Codreanu, 1949) in the Drava River, the leach Caspiobdella fadejewi (Epshtein, 1961) in the Sava River and the oligochaete Potamothrix moldaviensis (Vejdovsky and Mrazek, 1903) in the Danube and the Drava Rivers (Tab. 2).

The highest number of NIMS was recorded in the Danube (16), followed by the Drava River (15). However, the average proportion of NIMS in a total number of macroinvertebrate taxa was much higher at four Danube sites (30.8%) than at 19 sites in the Drava River (8.5%) where NIMS were found. In the Sava River, 10 NIMS were recorded at 17 sites with an average proportion of 6.9%, while in the Mura River, a single NIMS was identified (1.7%) at site MU2 (Tab. 1).

Out of all recorded NIMS, nine species were relatively widespread and recorded at larger number of sites in two or three large rivers (Fig. 2b), with variable densities (1–7964 ind. m⁻², average 242 ind. m⁻², extreme 127762 ind. m⁻² for Potamopyrgus antipodarum at site DR3 in the reservoir “Čakovec”). The most widespread invasive species were clam Corbicula fluminea (CorbFlum − abbreviation in Fig. 3b), isopod Jaera istri (JaerIstr), amphipods Dikerogammarus villosus (DikrVill) and Chelicorophium curvispinum (ChelCurv) recorded in the Danube, Sava and Drava at 27, 24, 22 and 19 sites, respectively. Out of five remaining widespread species, Chelicorophium sowinskyi (ChelSown) and Dikerogammarus haemobaphes (DikrHaem) were characteristic for the Sava River, Dreissena polymorpha (DreiPolm) and Echinogammarus ischnus (EchnIsch) for the Drava and Danube Rivers and mud snail P. antipodarum (PotmAntp) was found only in the Drava reservoirs, part of the Drava course downstream of the last reservoir and at one site in the Mura River. Remaining 10 NIMS had been recorded at a smaller number of sites (2–6), mostly in lower abundance (1–207 ind. m⁻²). The most upstream distributed invasive clams and crustaceans (amphipods) were C. fluminea and D. haemobaphes in the Sava River, and D. polymorpha and D. villosus in the Drava River (Fig. 2b).

There was no correlation between the number of NIMS and total number of taxa (Fig. 4; Spearman correlation rs = −0.04, p > 0.05), while the same correlation between the number of NIMS and the number of native taxa was negative and significantly different from zero (rs = −0.42, p < 0.05).

Out of all studied sites, only seven sites (14.6%) were not inhabited by NIMS (Figs. 1 and 2). On the other hand, most sites had high (35.4%) or severe (31.3%) family level site-specific biocontamination index (SBCI = 3 or 4). At four Danube sites, high (2) or severe (2) SBCI were observed, with similar richness and abundance contamination at three of four sites. In the Sava and Drava Rivers, severe SBCI mainly results from the high value of ACI (>0.5), while the percentage of sites with severe SBCI was more than two times higher in the Drava River (50.0%) than in the Sava River (19.5%). This was mainly due to the presence of invasive clam D. polymorpha and mud snail P. antipodarum in the Drava River, which were not recorded in the Sava River. In the Mura River, only low biocontamination was observed at one investigated site (MU2), where only P. antipodarum was found and two remaining sites showed no biocontamination.

In non-parametric multidimensional scaling (NMDS) analysis of whole macroinvertebrate assemblages (Fig. 5), individual sites were grouped according to the level of biocontamination and position at particular river section or two Drava reservoirs. Sites without or with a low level of biocontamination clearly separated from other sites, except for site SA5, due to the absence or low abundance of NIMS and similar composition of dominant taxa (Gammarus fossarum, Baetis sp. and Chironomini Gen. sp.). Among sites with medium to severe biocontamination, sites of the Lower Sava clearly differed from sites of the Middle Drava with the Lower Drava and Danube sites in-between. While the Lower Sava sites were dominated by four NIMS (two more abundant C. sowinskyi and D. haemobaphes, and also C. fluminea and J. istri), sites along Middle Drava differed from those at the Lower Drava and Danube sites, mainly due to different combinations of dominant NIMS (Middle Drava: C. curvispinum, D. villosus, D. polymorpha, J. istri and P. antipodarum; Lower Drava and Danube: same as previous but C. fluminea instead of P. antipodarum). Four sites at the two Drava reservoirs (2nd “Čakovec” and 3rd “Dubrava” in Croatia) had medium to severe biocontamination and impoverished macroinvertebrate assemblages, with extreme domination of P. antipodarum at DR3, while other three sites (DR4, DR6 and DR7) were more similar to other downstream Drava sites due to presence or domination of D. polymorpha.

Fig. 2 (a) Biocontamination indices (ACI: Abundance Contamination Index, RCI-Fam: Richness Contamination Index at family level) and Site-specific Biocontamination Index at family level (SBCI-Fam) at all 48 sites (colours as in Fig. 1). (b) Relative proportion of nine most widespread NIMS in total abundance of all NIMS (sum of NIMS = 100%, rare species summed as other NIMS) at each site.

Fig. 3 (a) PCA for 18 environmental parameters and 46 sites (sites DR4 and DR7 in the middle of second and third reservoirs were excluded). (b) Canonical correspondence analysis (CCA) triplot for 46 sites with different symbols and colours according to their positions at seven rivers sections (including two reservoirs), with 14 species (species abbreviation in orange: 11 NIMS-ObesObes: Obesogammarus obesus, ChelRobs: Chelicorophium robustum), other NIMS names abbreviations in text of Section 3.1; 3 native amphipods − GammFoss: Gammarus fossarum, GammRoe: Gammarus roeselii, SynrAmbl: Synurella ambulans) and 12 environmental factors as vectors.

Fig. 4 Relationships between a number of non-indigenous macroinvertebrate species (NIMS) and number of total or native taxa at all 48 sites with lines representing linear regression (full line: NIMS-nu. of native species, dotted line: NIMS-total taxa nu.) and rs-Spearman correlations indices (bold-significant correlations, p < 0.05).

Fig. 5 Non-parametric multidimensional scaling (NMDS) of macroinvertebrate assemblages (square-root transformed densities) at all 48 sites with markers according to site position at specific river section and colours showing the level of biocontamination at each site (colours correspond to those in Fig. 1).

3.2 Non-indigenous species and environmental parameters

Results of principal component analysis (PCA) for 46 sites using 19 environmental factors (13 physicochemical, 3 geographical and 2 for substrate composition) (Fig. 3a) had eigenvalues of 0.262 and 0.196 for first two axes, which explained 45.8% of total variation. Sites along the Sava River clearly separated from other sites due to higher conductivity and temperature, a higher proportion of soft substrate and concentration of orthophosphate and nitrogen nutrients. Sites along the Middle Drava were characterized by a higher proportion of hard substrate, higher oxygen concentration and were located at higher altitudes. The site DR5 represents outlier due to pollution from the town of Varaždin and maximum concentrations of ammonia and nitrites. Some other polluted sites include SA5–SA7 with high concentration of orthophosphate, nitrite and nitrate due to pollution from the towns of Zagreb and Velika Gorica; sites SA10 and SA11 downstream of the Lonja River confluence had increased orthophosphate, while sites SA19–SA21 had increased concentrations of chloride ions due to inflow of pollution from the right tributary of the Bosna River (upstream of site SA19) and from the sugar factory in the town of Županja (located upstream of site SA20, which affects sites SA20 and SA21). The four Danube sites and three Mura sites were characterized by higher biological oxygen demand (BOD₅) and higher concentrations of chloride ions and total organic carbon (TOC).

Results of canonical correspondence analysis (CCA) for selected 14 species, 46 sites and 12 environmental factors had eigenvalues 0.692 and 0.538 for the first two axes and explained 44.0% of total and 66.7% of fitted variation of species–environment relationships (Monte Carlo test: CCA1 F = 10.8, p = 0.0001, all axes F = 5.3, p = 0.0001). The CCA triplot showed a clear distinction of the Middle Sava sites where non-native amphipods D. haemobaphes and C. sowinskyi were abundant and had high conductivity, orthophosphate concentrations and proportions of the soft substrate. At sites without NIMS (SA1–SA4, DR2, MU1 and MU3), which were the most distantly located from the Danube, different combinations of abundant native amphipods G. fossarum, G. roeselii and Synurella ambulans were found. Sites at the Middle Drava with higher proportions of hard substrate, BOD₅, oxygen and organic nitrogen concentration were inhabited by different combinations of NIMS species that group with sites along the Lower Drava and Danube River reaches. Two sites in the Drava reservoirs and sites at the Drava immediately downstream of these reservoirs were characterized by higher abundances of D. polymorpha and P. antipodarum.

3.3 Biological metrics and non-indigenous species

Significant negative correlations of biological metrics used for ecological quality assessment and RCI were obtained in most cases except for Croatian saprobic index HR-SI (Tab. 3). Only two significant negative correlations with ACI were obtained for EPT-[%] and EPT-Taxa. Saprobity module showed a significant negative correlation with RCI-fam and SBCI-fam (Tab. 3).

Several biological metrics widely used for assessing organic pollution (BMWP, IBE AQEM, HR-SI) were calculated with and without NIMS and correlated with six physicochemical parameters that indicate organic pollution (BOD₅, organic nitrogen, TOC, concentrations of ammonium, total phosphorus and chloride) (Tab. 4). Correlations of BMWP and IBE AQEM with physicochemical parameters did not differ significantly when NIMS were included or excluded. On the other hand, HR-SI had little better correlation with physicochemical parameters when NIMS were excluded. Furthermore, saprobity module for all study river-types indicated the best correlation with three physicochemical parameters when NIMS were included in IBE AQEM but excluded from HR-SI (combined option in Tab. 4).

4 Discussion

4.1 Macroinvertebrate assemblages, non-indigenous species and biocontamination

This study provides the first comprehensive survey of NIMS in Croatian large rivers, with a total of 21 recorded species, most of which were crustaceans and molluscs of Ponto-Caspian origin. Out of a total number of recorded taxa (236), NIMS represented 8.9% and were found at most surveyed sites (85%). Similar proportions of NIMS were found in the Upper and Middle Danube during third Joint Danube Survey (JDS), 11.2 and 10.7%, respectively (Liška et al., 2015) and in the Rhine River (11.3%, Leuven et al., 2009), while significantly lower proportion was recorded in the Lower Danube (3.9%) (Liška et al., 2015). However, in our study, less than a half of the recorded NIMS were widespread and had higher densities and proportion in total density of benthic macroinvertebrates. RCI more or less regularly increased in downstream sections of the Sava and the Drava Rivers, with maximum values recorded in the Danube River. This suggests that the Danube River serves as a source of the most NIMS spreading to the western tributaries. Most of the nine widespread and abundant NIMS were probably introduced into the Sava and Drava Rivers by passive dispersal by ships or boats.

When compared to methodologically similar studies across Europe, the number of freshwater NIMS in Croatian large rivers was larger than in most studied waterways of all three major invasion corridors surveyed by Arbačiauskas et al., 2008, Arbačiauskas et al. (2008, 2011a) and comparable to the French section of the Moselle River (Devin et al., 2005) or to the middle section or the whole course of the Danube River (Paunović et al., 2007; Liška et al., 2015). The close proximity of studied reaches of the Sava and the Drava Rivers to the Middle Danube, which was probably invaded by most Ponto-Caspian NIMS during the first half of twentieth century, in combination with the existence of reservoirs on the Drava River, could be responsible for a relatively higher number of NIMS in Croatian large rivers.

Out of 27 NIMS recorded in the Middle course of the Danube in JDS3 (Liška et al., 2015), 14 species were identified in the Croatian section, which is comparable to 16 NIMS obtained in this study. According to the detailed study of the Danube River conducted in 2001, 2007 and 2013, it can be observed that the number of NIMS significantly increased, counting 12, 20 and 34 NIMS, respectively (Liška et al., 2015). On the other hand, the level of biocontamination calculated at the species level in 2013 decreased in every section of the Danube (Upper, Middle and Lower) in comparison to 2007 due to the reduced participation of non-native taxa in total abundance of macroinvertebrate assemblages (Liška et al., 2015). With this study as a baseline, future studies could examine the temporal and spatial dynamics of the relative proportion of NIMS in macroinvertebrate assemblages in Croatian large rivers.

Danube tributaries, the Sava and Drava Rivers in Croatia, contain a subset of Danube NIMS, mostly of Ponto-Caspian and Asian origin. Interestingly, the Drava River had a higher number of NIMS than the Sava River, although vice versa could be expected since the Sava is navigable up to the town of Sisak (rkm 594) (Komatina and Grošelj, 2015), while the Drava River only to the town of Osijek (rkm 22). Comparing the Sava and the Drava Rivers, out of nine more abundant and widespread species (Tab. 2), three NIMS, D. polymorpha, P. antipodarum and E. ischnus, were found only in the Drava River. Zebra mussel and mud snail presence in high density in the Drava River is a consequence of the existence of reservoirs along the upper part of the Middle section of the Drava River in Croatia. Distribution and dispersal of D. polymorpha, considered to be one of the most aggressive invasive species globally (Lowe et al., 2004), is heavily dependent on spatial positioning of standing waters such as reservoirs since it reproduces by planktonic larvae which are swept downstream and cannot swim against current in flowing waters (Havel et al., 2005, Johnson et al., 2008). Similarly, in habitats modified by human activities like reservoirs, P. antipodarum performs as a successful early colonizer (Quinn et al., 1998), reaching high abundance in low-speed waters with high densities in macrophytes (Richards et al., 2001). Similar findings are reported in this study: P. antipodarum was found in the Drava reservoirs with extremely high density (127762 ind. m⁻²) in macrophytes at site DR3 (littoral zone of the reservoir “Čakovec”). Further, the density of both zebra mussel and mud snail was increased immediately downstream from the last reservoir Dubrava, with a steady decrease of their densities in more downstream reaches of the Drava River.

Besides a different number of NIMS, there was a significant difference in species composition between the Sava and the Drava Rivers, with assemblages in the Drava River being more similar to those in the Danube (Figs. 2b and 5). Proximity to the Danube as the major source of NIMS and different colonization history could be responsible for observed differences. However, due to limited extent of previous similar studies in Croatian large rivers (Žganec et al., 2009, 2018; Lajtner and Crnčan, 2011; Lucić et al., 2015), colonization history of the most abundant and widespread NIMS can only partly be reconstructed. Zebra mussel and mud snail were probably introduced to Drava reservoirs shortly after their construction between 1975 and 1989 either by natural (birds) or some human-mediated vectors. Since D. polymorpha and P. antipodarum find optimal conditions for massive reproduction in reservoirs, the Drava reservoirs represent constant sources of larvae or juveniles that colonize reaches of the Drava River downstream of the reservoirs. Interestingly, although zebra mussel was found in some artificial lakes close to the Sava River but not directly connected to it, this species is very rare and found only at few sites and in low abundance in the Sava River (Lajtner, personal communication). This clearly shows that reservoirs are major source water bodies for this species in Croatian large rivers. However, failure of zebra mussel to establish denser population in the Sava River (few specimens were found in 2011, Lucić et al., 2015) could also be a consequence of some other factors as well, such as dominance of psammopelal and mud covering hard surfaces on the river bottom and severe pollution of the Sava River in the past (Matoničkin et al., 1975). Other NIMS have probably been spreading upstream by ship traffic, while in the Drava River they are spread also by small fishermen's boats or some other vectors (attachment on fisherman's equipment, by birds, mammals, etc.). A recent study of Peracarida fauna along the whole course of the Sava River (Žganec et al., 2018) demonstrated that stochastic processes of the passive upstream spread of “the killer shrimp” D. villosus, probably by ship traffic, could explain its absence in the Sava River until 2015. It is possible that similarity of physicochemical conditions between the Drava and the Danube Rivers (Fig. 3a) and proximity to the Danube were main factors responsible for earlier and faster upstream spread of D. villosus, which caused the disappearance of earlier colonizer “the demon shrimp” D. haemobaphes in the Drava River (Žganec et al., 2009, unpublished). Other differences in NIMS composition between Sava and Drava, dominance of invasive amphipods C. sowinskyi and D. haemobaphes in Sava and C. curvispinum and D. villosus in the Drava could be a consequence of different substrate composition and different levels of water pollution during the past.

In both, the Sava and the Drava Rivers, the decrease in the number of native taxa and increase of several NIMS and RCI in downstream direction was observed. The lowest number of native taxa and the highest RCI values were observed in the Danube River. Considering all sites together, there was a significantly negative correlation between the number of NIMS and native taxa (Fig. 4). Higher organic and chemical pollution of downstream sections of the Sava and the Drava Rivers, which was especially pronounced in the past (Sava: Matoničkin et al., 1975; Drava: Meštrov et al., 1978,1989; Ćuk, unpublished data), probably caused the decline or complete disappearance of native taxa. That probably facilitated colonization of these rivers by more pollution tolerant NIMS which could fill empty niches. However, it is also possible that NIMS were also the drivers of the change in macroinvertebrate assemblages of the Sava and the Drava Rivers. For example, some species like D. villosus were shown to be strong predators of other macroinvertebrates (Dick et al., 2002; Krisp and Maier, 2005) and could, therefore, cause the decline or disappearance of many native taxa, especially those from EPT groups. It seems to be the most probable scenario of multiple causal agents of the decline occurring synergistically (Didham et al., 2005; MacDougall and Turkington, 2005). That is, the decline in native populations probably occurred prior to invasion due to habitat disturbance and is additionally intensified by the invasion of NIMS. However, recent improvements in water quality in Croatian large rivers (Ćuk, unpublished) could facilitate further spreading of NIMS, but as some studies recently demonstrated (Hellmann et al., 2017), they could also enable recovery of native taxa.

4.2 Non-indigenous species and environmental parameters

Physicochemical parameters of the Sava River clearly differed from those in the Drava and the Mura Rivers, as higher recorded values of conductivity, alkalinity, calcium and magnesium ions concentration in Sava were due to the dominance of carbonates in basins of most of Sava's right tributaries, unlike catchment area of the Drava and Mura Rivers, where silicates dominate in the geological composition (Sommerwerk et al., 2009). Further, sites along the Sava River had a higher proportion of soft substrate, temperature, total suspended solids and nutrients (especially orthophosphate, ammonia and nitrite), while the Drava and the Mura sites had a higher proportion of hard substrate, oxygen and organic nitrogen concentrations. The differences in substrate composition, water pollution in the past and other factors mentioned above, together with previously mentioned colonization history, may have affected the composition and structure of NIMS. Obviously, the distance to the Danube River was an important environmental factor since upstream sites of the Sava and Drava Rivers were either not yet colonized by NIMS or only few NIMS were present. However, further upstream and downstream spreading of nine most important molluscs and peracarids could be expected in the Sava and Drava Rivers. Pollution of the Sava from towns of Zagreb and Velika Gorica, as well as the Drava from the town of Varaždin, offer the opportunity to conduct more detailed studies of tolerance of specific NIMS and whole macroinvertebrate assemblages to pollution.

4.3 Biological metrics and non-indigenous species

Frequently applied biological metrics used in water quality assessment, such as BMWP, biotic indices and Pantle–Buck saprobic index, are generally designed on the basis on native fauna. Unfortunately, invasion of NIMS irreversibly changed native communities in European large rivers (Arbačiauskas et al., 2008). In certain sites where NIMS are extremely abundant and native taxa scarce or even completely absent, standard biological metrics could, therefore, provide unreliable results and consequently, the assessment of ecological status of water bodies could be compromised.

Croatian Indicator Saprobic System (HR-SI) has the longest tradition in Croatia and many NIMS have already been included in operational taxa list with similar indicator values assigned (e.g. D. villosus − 2.1; C. curvispinum − 2.1; P. antipodarum − 2.8 etc.). Since it was observed that the level of biocontamination increases with time (Arbačiauskas et al., 2008) due to higher abundance of already established NIMS and newly arrived NIMS, the value of HR-SI may also change significantly even if the water quality remains constant, indicating either higher or lower saprobic status in comparison to compatible physicochemical parameters. For example, probable scenario is that NIMS decrease the index value (indicating lower saprobic status) in relatively polluted river stretches due to their lower indicator values comparing to resident polysaprobic taxa (oligochaetes, chironomids) which usually dominate in polluted sites. Furthermore, the tendency of many NIMS to appear in high abundance can diminish or completely disregard the presence and influence of the natives within the calculation of HR-SI with the absolute abundance.

The possible problems in distorted values of other biological metrics have been discussed previously by several authors, e.g. BMWP (MacNeil et al., 2004, 2009) or Belgian Biotic Index (BBI) (Gabriels et al., 2005). As NIMS contribute to the taxa richness, they may theoretically increase the index value (BMWP, biotic indices), particularly in the early phase of invasion when NIMS mostly appear in addition to the native species. For example, the presence of a single specimen of a tolerant “killer-shrimp” D. villosus (family Gammaridae) in a polluted surrounding increased the index BMWP value for 6 units in case no other gammarid is present in the sample! However, results obtained in our study indicate relatively high negative correlation between biocontamination indices (particularly RCI) and both BMWP and IBE AQEM (Tab. 3), which is in accordance with several previous studies (e.g. Arbačiauskas et al., 2008; MacNeil, 2014) and could imply that NIMS seriously decrease or even eliminate sensitive taxa (Ephemeroptera, Plecoptera, Trichoptera) that contribute to the higher BMWP or IBE AQEM score. If water chemistry remains constant, such biological results indicate certain discrepancy from the actual pressure they are designed to assess, thus providing unreliable results. Therefore, it would be most useful if these indices undergo certain revision in the presence of NIMS.

Practice differs in EU countries about whether NIMS should be included in standard metrics in water quality assessment. For example, Gabriels et al. (2005) suggested inclusion of NIMS into Belgian Biotic Index (BBI) as being members of benthic macroinvertebrate assemblage, thus reflecting biodiversity, however, with no tolerance class assigned; in the Netherlands, NIMS are excluded from water quality assessment (Orendt et al., 2010), while in Germany and Austria, NIMS are partly included. To check the validity of inclusion or exclusion of NIMS, we calculated the above-mentioned biological metrics with and without NIMS and then correlated the obtained values with five physicochemical parameters that indicate organic pollution (Tab. 4). In most cases, no significant difference was observed when NIMS were included/excluded from the metrics. However, the difference observed in HR-SI with NIMS included/excluded correlated with physicochemical parameters was small, but, nevertheless, at severely biocontaminated sites it would probably be advisable to exclude NIMS from this metric due to their unclear saprobic tolerance at this moment and minimal or even negligible influence of natives at sites with high ACI. This also refers to the same metric within the saprobity module.

In conclusion, more research is needed to fully understand the validity of inclusion/exclusion of NIMS into the existing biological metrics designed on the basis of native fauna. The suggested research direction is to examine the relationship of NIMS with pollution and other pressures. Within the phase of invasion, NIMS do not show a clear saprobic tolerance and more time and research is necessary to detect their saprobic tolerance (Friedrich and Herbst, 2004). As the presence of NIMS distorts the interpretation of community response to a certain pressure within the watercourse, assessment methods according to WFD can, therefore, provide incorrect and even inverted results regarding module “organic pollution” due to the dominance of invasive species, as was shown by Arndt et al. (2009). Detecting and understanding the real causes of changes in macroinvertebrate assemblages should be defined as a high priority. This is particularly important in water management where a certain set of measures should be applied to achieve a good water status required by WFD.

Acknowledgements

The authors would like to thank the employees of Central Water Management Laboratory of Hrvatske vode for providing the results of physicochemical analysis.

References

All Tables

All Figures

Fig. 1 Site-specific biocontamination index (SBCI-Fam, family level) values presented in different colours for all 48 studied sites sampled in 2015 along the Mura (MU), the Drava (DR), the Sava (SA) and the Danube (DA) Rivers. Right inset figure shows the position of Croatia in Europe with country codes (HR: Croatia, AU: Austria, BH: Bosnia and Herzegovina, HU: Hungary, IT: Italy, MG: Montenegro, SI: Slovenia, SR: Serbia), while left inset figure shows the enlarged area of the Drava reservoirs.

In the text

	Fig. 2 (a) Biocontamination indices (ACI: Abundance Contamination Index, RCI-Fam: Richness Contamination Index at family level) and Site-specific Biocontamination Index at family level (SBCI-Fam) at all 48 sites (colours as in Fig. 1). (b) Relative proportion of nine most widespread NIMS in total abundance of all NIMS (sum of NIMS = 100%, rare species summed as other NIMS) at each site.
In the text

Fig. 3 (a) PCA for 18 environmental parameters and 46 sites (sites DR4 and DR7 in the middle of second and third reservoirs were excluded). (b) Canonical correspondence analysis (CCA) triplot for 46 sites with different symbols and colours according to their positions at seven rivers sections (including two reservoirs), with 14 species (species abbreviation in orange: 11 NIMS-ObesObes: Obesogammarus obesus, ChelRobs: Chelicorophium robustum), other NIMS names abbreviations in text of Section 3.1; 3 native amphipods − GammFoss: Gammarus fossarum, GammRoe: Gammarus roeselii, SynrAmbl: Synurella ambulans) and 12 environmental factors as vectors.

In the text

	Fig. 4 Relationships between a number of non-indigenous macroinvertebrate species (NIMS) and number of total or native taxa at all 48 sites with lines representing linear regression (full line: NIMS-nu. of native species, dotted line: NIMS-total taxa nu.) and rs-Spearman correlations indices (bold-significant correlations, p < 0.05).
In the text

	Fig. 5 Non-parametric multidimensional scaling (NMDS) of macroinvertebrate assemblages (square-root transformed densities) at all 48 sites with markers according to site position at specific river section and colours showing the level of biocontamination at each site (colours correspond to those in Fig. 1).
In the text