Introduction

The topmouth gudgeon Pseudorasbora parva (Temminck & Schlegel, 1846) is a freshwater fish with a natural distribution range in Northeast Asia but is introduced throughout Central Asia, Europe and North Africa (Perdices & Doadrio 1992, Gozlan et al. 2002). In its non-native distributional range, the species is highly invasive and recognised as one of the worst aquatic invaders in Europe and was listed on the European List of Invasive Alien Species of Union concern since 2016 (European Commission 2016). The adverse impacts of P. parva are related to the transmission of parasites and pathogens (Ahne & Thomsen 1986, Gozlan et al. 2005, Spikmans et al. 2020), suppression of zooplankton populations due to its feeding (Musil et al. 2014), competition for food and habitat (Britton et al. 2010), predation on eggs and early stages of native species, and even facultative parasitic behaviour on other fish species (Gozlan et al. 2010). It can become a dominant species when introduced to new areas (Britton et al. 2007), and a negative relationship was observed between P. parva abundance and fish diversity indices (Spikmans et al. 2020).

The rapid dispersal of P. parva and its ability to establish itself in new areas is a result of its life-history traits, such as high morphological variability, which facilitates the utilisation of different types of habitat (Záhorská et al. 2009) and enables tolerance of extreme environmental conditions (Pollux & Korosi 2006, Carosi et al. 2016). In addition, the species has a short generation time, small body size, and high reproductive capacity, traits of species with high invasiveness risk (Gozlan et al. 2010, Radočaj et al. 2021). It is also noted that even restricted localised introductions of a few individuals of P. parva can lead to their successful establishment and formation of invasive populations (Simon et al. 2015).

Pseudorasbora parva was accidentally introduced to Europe in the 1960s with Chinese carp species' imports and then dispersed through the River Danube basin (Bănărescu, 1964 as cited in Gozlan et al. 2002, Piria et al. 2018). With the help of several human-mediated introductions (Caiola & De Sostoa 2002, Gozlan et al. 2002, Piria et al. 2018) and by the connected river systems, e.g. Danube and Rhine (Pollux & Korosi 2006), the species was able to achieve pan-continental colonisation of much of Europe (Gozlan et al. 2010). Once the species established self-sustaining populations in its non-native range, most often colonisation followed up by dispersal utilising ‘stepping stones’ (Gozlan et al. 2010, Simon et al. 2015, Carosi et al. 2016), a network of habitat patches providing food and shelter within unsuitable habitat (Saura et al. 2014, Rocha et al. 2021). In the early 1980s, P. parva was recorded in the River Danube in Serbia (Cakić et al. 2004), where the species likely gradually invaded Croatian watercourses via upstream dispersal (Piria et al. 2018). Since then, the species likely has been accidentally translocated through numerous fish farms in Croatia. Finally, P. parva was found in inland waters in 1985 in the River Sava, around Zagreb town, very far from the initial entry point (Habeković & Popović 1991). Today, P. parva has a wide distribution of established populations across Croatia (Boršić et al. 2018, Mihinjač et al. 2019, Piria et al. 2019).

Pseudorasbora parva occurs in various lentic and lotic habitats (Gozlan et al. 2010). However, the species likely favours well-vegetated lentic habitats and slow sections of lowland rivers (Britton et al. 2007, Gozlan et al. 2010, Kottelat & Freyhof 2007) rather than lotic habitats that are mainly utilised as dispersal corridors (Pollux & Korosi 2006) and where it may be found in coexistence with other non-native fish species (Carosi et al. 2016).

In a laboratory experiment, Britton (2012) determined that native cyprinids, as facultative predators, could reduce the chance of P. parva establishment. On the contrary, some authors found that in nature, this species has positive associations with native fish species (Beyer et al. 2007) or that the presence of P. parva coincides with an increase in species richness and diversity index, although it reduces the abundance of some native fish species (Rechulicz 2019).

Most of the records on P. parva in Croatia are known from regularly-monitored larger rivers and lakes (Boršić et al. 2018), but knowledge of its distribution in numerous lowland streams and canals is sparse and limited (Delić 1993, Mustafić et al. 2020), whereas environmental factors associated with the occurrence of this species remains unknown. Environmental factors that shape fish communities (including the occurrence of each species) differ between regions (Wang et al. 2006, Bierschenk et al. 2019, Czeglédi et al. 2020). Thus, we hypothesised that the occurrence of P. parva is associated with a set of specific environmental predictors in lowland waterbodies of the River Sava basin in Croatia.

This study aimed to explore the exact extent of the distribution of P. parva in lowland streams and canals in the River Sava basin and determine the environmental abiotic and biotic factors associated with its occurrence. An information theoretic approach was used to compare different plausible hypotheses about the association between environmental factors and the presence of P. parva. Because streams and canals, as an interconnected network of aquatic habitats, are a likely source of further continuous invasions of P. parva, the obtained results could be helpful for the management of this highly invasive species.

Material and Methods

Study area and sampling

The sampling sites in the River Sava basin were chosen across central and eastern Croatia. The selection was based on the scientific and popular literature and the results of our previous research to fill the existing gaps in distributional data of self-sustaining populations of non-native fish species, including P. parva. Sites where no fish were found, and sites where it was impossible to collect complete environmental data (intermittent or inaccessible water bodies) were excluded from further analysis.

Table 1.

Factors associated with Pseudorasbora parva presence/abundance listed in the literature.

Sampling was performed from mid-June 2019 to mid-November 2019, mid-March 2020 (only one site sampling) and from the end of May 2020 to the end of August 2020. Winter months (December, January, February) were excluded intentionally due to the low activity of fish and the absence of aquatic vegetation cover (Table S1). Sampling was not performed from the second half of March to almost the end of May 2020 because of the lockdown due to Covid-19 (see Jakovljević et al. 2021). Fish were collected by wading in a single pass on various substrates and mesohabitats (e.g. pool-riffle-run) alternating upstream using a 2.2 kW (Hans Grassl) electrofishing device in both natural streams and canals, without using block nets. Electrofishing parameters were set specifically for every location to minimise the negative impact on the fish. The mesh size of the dip net was 6 mm. Sampling depths ranged from 0.1 to 1 m in both types of water bodies. Natural streams were sampled along two 100 m transects with at least 500 m between them where possible. Canals were sampled along 300 m transects at one sampling site where possible. All transects on the same water bodies were included as separate locations because of the possible differences in environmental variables. Fish identification was performed immediately after sampling following Kottelat & Freyhof (2007), while the most recent scientific nomenclature was used according to Ćaleta et al. (2019). The standard length (SL, ±1 mm) was measured, after which the native fish were released, and non-native fish were removed from the habitat as required by national legislation. Only 1+ age fish, estimated by their SL, were considered in further analysis. The sampling of the fish complied with Croatian Nature Protection Act and Croatian Freshwater Fisheries Act as approved by Croatian Ministry of Economy and Sustainable Development (permit class UP/l-612-07/19-48/140, number 517-05-1-1-19-2) and Croatian Ministry of Agriculture (permit class UP/I-324-05119-03/03, number 525-1311855-19-2 and class UP/I-324-01/20-01/04, number 525-13/0733-20-2). Because P. parva was not universally present in sampled locations, species presence-absence data was chosen as a more informative approach in this study.

At each sampling site, essential physio-chemical water quality variables (temperature in °C, oxygen in mg/l, conductivity in µS) were measured with a multiparameter device (SI Analytics HandyLab 680). The altitude and coordinates of the sampling site were determined with a GPS device (Garmin GPSMAP 78S) with an accuracy of 3 to 5 m to a 95% certainty level. The water depth was measured while wading and expressed as the maximum depth in meters at the sampling site. Sediment type cover was estimated while wading and expressed as a percentage of each sediment type at the sampling location (silt, sand, gravel, rock). Aquatic vegetation cover (including submersed, emergent and floating plants) was estimated while wading and was divided into six categories: 0 – not present; 1 – less than 5%; 2 – 5-25%; 3 – 25-50%; 4 – 50-75%; 5 – more than 75% cover. Distance from the nearest fish farm was determined using satellite imaging and expressed in meters.

Fig. 1.

Sampling sites used in the analysis (black dots) in the River Sava basin (dark grey area) in Croatia (light grey area) (for details, see Figs. S1-S3, Table S1).

Data preparation

An information-theoretic (IT) model selection approach was used to predict P. parva occurrence. The IT approach uses prior knowledge to formulate a set of biologically plausible models and selects the best supported based on the data (Burnham & Anderson 2002). Results from published studies on the ecology of P. parva (Table 1) were combined with hypotheses based on the experience of the researchers in this study to formulate 11 a priori models (Table 2). Such an approach was selected to avoid “data fishing” (Burnham & Anderson 2002).

Since P. parva prefers lentic habitats, stream sections where gravel was absent and only a silt and/or sand substrate was present were used to indicate lentic-characteristic conditions; in low flow conditions, silt and sand particles covered the streambed; thus, the absence of a gravel substrate was used to distinguish this habitat type (Coulombe-Pontbriand & Lapointe 2004, Earle 2019). Based on this approach, a substrate variable was replaced with a simple presence-absence of gravel categorical variable, which also had a good balance among sites. The Shannon-Wiener diversity index (SWI) was calculated from the sampling data as the most commonly used fish diversity measure (Clarke & Warwick 2001) and was included as a diversity predictor in the analysis, along with species richness. Northern pike Esox lucius L., huchen Hucho hucho L., European perch Perca fluviatilis L., brown trout Salmo trutta L., zander Sander lucioperca L., European catfish Silurus glanis L., and chub Squalius cephalus L. were grouped as a predatory species category; black bullhead Ameiurus melas (Rafinesque, 1820), gibel carp Carassius gibelio (Bloch, 1782), pumpkinseed Lepomis gibbosus L., monkey goby Neogobius fluviatilis (Pallas, 1814) and round goby Neogobius melanostomus (Pallas, 1814) comprised a non-native species predictor. The natural, altered or artificial state of the waterbody, presence of predators, presence of non-native species and vegetation cover were used as categorical variables and all other predictors as continuous variables. All environmental variables included in the formulated models were tested for possible collinearities and outliers. There was high collinearity between SWI and species richness; SWI was retained because it had an even distribution with no outliers.

Table 2.

Formulated models for the IT model selection approach with an explanation for every included model (alt – altitude of the sampling location; gravel – presence/absence of gravel in the stream bed as an indicator of low flow conditions; pred – the presence of predatory fish species: E. lucius, H. hucho, P. fluviatilis, S. trutta, S. lucioperca, S. glanis, S. cephalus; distf – distance from the nearest fish farm pond; veget – level of aquatic vegetation; SWI – Shannon Wiener index, a measure of fish diversity; cond – conductivity of water, a proxy for water quality; NNS – the presence of other non-native fish species: A. melas, C. gibelio, L. gibbosus, N. fluviatilis, N. melanostomus; alter – natural or regulated state of the waterbody).

Results

Ultimately, 111 sites on 46 streams and 35 canals that met the criteria were included in the analysis (Fig. 1, Table S1). Of 111 sampled locations, P. parva's presence was confirmed at 33 sites in 12 streams and 18 canals, where its distribution was previously unknown. There was no particular pattern in the distribution, as it was continuously present in the entire basin, from east to west (Figs. S1-S3, Table S1). In addition to P. parva, 40 other fish species were found during sampling. The most abundant species were S. cephalus. Pseudorasbora parva co-occurred with 33 of them; among non-native species and predators, it most frequently co-occurred with C. gibelio and S. cephalus, respectively (Table S1).

Fig. 2.

Plotted probability of occurrence of Pseudorasbora parva according to the best-fitting model.

The two best-fitting models (M7 and M8) had an almost identical fit (one AIC unit apart), so a new model was constructed that comprised variables common to M7 and M8 (designated M13). Further investigation of this model showed that including a covariate for conductivity did not substantially improve model fit, so simpler model M14 was recognised as the final best-fitting model (Table 3). Finally, the model was validated by plotting model residuals against fitted values and model covariates.

The best fitting model to predict P. parva occurrence in these habitats included two fixed effects: the presence of gravel and biodiversity index (SWI) and was formulated as:

Where Pparva_i is the probability of P. parva occurrence at location i, which follows a binomial distribution that has an expected probability E of P. parva occurrence with mean n_i × π_i and variance ni × πi × (1 – πi). A logit link function enables the mapping of changes in the predictors to the probability scale of a Bernoulli distribution, with an interval between 0 and 1. The variable gravel_i is a categorical independent variable with two outcomes, the presence and absence of gravel at location i, and it is used as a proxy for flow conditions. Variable SWI_iis a continuous independent variable describing fish biodiversity (Shannon Wiener index) at location i.

The probability of P. parva occurrence was negatively associated with the presence of gravel and positively associated with SWI (Table 4, Fig. 2). There were equal numbers of sites with and without a gravel substrate where the species was absent (39 vs. 39), but P. parva was present at 25 sites without gravel and only eight sites with a gravel substrate. The SWI index of the sampled sites ranged from 0 to 2.17, and the sites where P. parva occurred ranged in SWI index from 0.52 to 2.17 (Table S1).

Discussion

The known distribution data on P. parva in continental Croatia is based chiefly on its presence in fish farm ponds, popular angling venues and major rivers (Boršić et al. 2018). This study confirmed its occurrence in several lowland streams and canals throughout the River Sava basin. Although it mainly inhabits lentic ecosystems (Beyer et al. 2007, Kapusta et al. 2008, Mihinjač et al. 2019), P. parva's presence in a wide range of lotic water bodies such as canals and streams confirms its ecological versatility, facilitating its high rate of its dispersal (Pollux & Korosi 2006). Nevertheless, even in lotic habitats, P. parva primarily inhabits stretches with environmental conditions characteristic of lentic habitats (Carosi et al. 2016). Indeed, in this study, the species mainly occurred in lower reaches with low flow rates containing sand and silt on the bottom.

Table 3.

Models ranked by relative goodness-of-fit from the best to worst (see Table 2 for abbreviations).

Table 4.

Numerical summary of Bernoulli GLM to predict the probability of Pseudorasbora parva occurrence in lowland streams and canals in the River Sava basin, Croatia (SWI – Shannon Wiener index).

Pseudorasbora parva was also significantly associated with high local fish diversity, measured by SWI. The top six best-fitting a priori models included both SWI and the presence of gravel predictors. Notably, the two models combined to produce a best-fitting model had SWI values included with different expectations. In the ‘Poor habitat’ model (M7), it was predicted that the occurrence of P. parva would coincide with low SWI values due to habitat degradation. The final result showed the opposite, as predicted by the ‘Lentic model’ (M8); the species occurred in low-flow conditions that support species-rich fish communities. Jackson et al. (2001) mention a pattern of continual addition of fish species along downstream gradients: headwater streams of low order support low diversity of specialists, while larger, higher-order streams situated downstream along the watershed gradient provide more stable conditions for numerous species. A similar effect is observed after flow alterations, causing a shift to lentic-type habitat, with fluvial specialists replaced by lentic generalists (Perkin & Bonner 2011). This pattern explains both predictors included in the model: the occurrence of P. parva is expected within downstream, low-flow sections of streams that generally support richer generalist fish species communities. This prediction is also supported by the fact that P. parva in this study mostly co-occurred with another lentic non-native species, C. gibelio. There are two opposite hypotheses regarding non-native species’ interactions with native communities after introduction. First is the biotic resistance hypothesis, presented by Elton (1958), which states that species-rich communities are more resistant to the successful invasion of non-native species. The basis for this hypothesis is that rich native communities efficiently utilise almost all habitat resources, leaving few ecological niches for non-native species to use. On the other hand, the acceptance hypothesis presumes a positive relationship exists between native and non-native species richness (Jeschke 2014). This relationship is usually not a direct effect of native on non-native species but an indirect effect of rich habitats, providing enough resources to sustain the establishment of new non-native species in an already rich native community (Stohlgren et al. 2003). The present study best supports the latter hypothesis and was supported by several previous studies of freshwater fish assemblages (Kennard et al. 2005, Muniz et al. 2021). Additionally, Kennard et al. (2005) state that non-native fish species’ presence could result from random introduction and dispersion, as they can persist in degraded and undisturbed habitats and coexist with rich native communities. This situation could be the case for the P. parva population investigated in this research: it has been established for a long time, it is effectively spreading in all suitable habitats, and is assimilating within respective fish communities. It should be noted that the presence of P. parva at some of these sites may not have been long enough to express its potential impact on native communities, as suggested by Rechulicz (2019).

The presence of P. parva was confirmed at some sites with only a single individual (Table S1), so it was impossible to confirm a viable population in these habitats, but again, this presence could be indicative of possible migration routes used by the species. Another drawback of the study is the need for more sampling seasonality and timing; sampling was done throughout the year (excluding winter), and each site was sampled only once, no matter the season and the time of the day. For this reason, temperature and oxygen levels could not be included in the models since they were not comparable, and the balance of the sampling months was not good enough to include seasonality as a random factor influencing the species' occurrence. Also, substrate composition was used as a proxy for flow conditions, presuming that it is directly correlated to average flow conditions of the site via sedimentation rates (Earle 2019) and flow velocity would not be an average representation of the entire site but just of the exact measuring point in the time of the sampling.

Positive associations of P. parva with macrophytes were observed in stagnant water bodies in introduced areas (Wolfram-Wais et al. 1999, Kapusta et al. 2008), but researchers studying this species in lentic environments in its native range (Ye et al. 2006) and in lotic environments in introduced range (Pollux & Korosi 2006) did not find any significant relationship. In this study, macrophyte cover was also not a significant predictor of the species' presence. Regarding these results, associations between macrophyte cover, P. parva presence and abundance are more complex and result from more complex interactions. Furthermore, the unbalanced seasonality of sampling in this study may have influenced the results, especially early and late in the year when macrophytes are likely not abundant. However, the samplings conducted in autumn resulted from mild weather and relatively high water temperatures, and aquatic macrophytes were present at sampling sites during the earliest sampling in March and the last sampling in November (Table S1).

Several studies mention fish farms as a primary source of P. parva invasions in new areas (Caiola & De Sostoa 2002, Aparicio et al. 2012, Rakauskas et al. 2021). Boršić et al. (2018) stated that P. parva occurs in high abundance near fish farms in the Danube River basin of Croatia. However, such a pattern is not confirmed in this study, and the reason for that may be a large number of cyprinid fish farms in this part of Croatia (Food and Agriculture Organization 2005); even if the farms are the primary source of the invasion, they are close to each other and present across the entire area so there is no pattern that a model could interpret as significant. It could be argued that the abundance of small-sized non-native species increases when the water from farming ponds is discharged into the adjacent canals after the farming season, but in Croatia, the farming period of warm-water cyprinids usually ends just before winter, in mid-November (Debeljak & Fašaić 1999), after the sampling period used in this study.

The impact of predators has previously been described as unfavourable to P. parva abundance (Csorbai et al. 2014, Verhelst et al. 2016) and even establishment (Lemmens et al. 2015), but the latter study was done in experimental ponds with high predator density. Although the presence of predators can reduce the abundance of P. parva, the impact on the occurrence of the species was not observed in this study.

Regarding abiotic factors, all of the investigated sites in this study were located below 200 meters of elevation, so they were all considered lowland streams and canals (Mihaljević et al. 2011). This narrow range of altitude was probably not enough to show any significant impact on the occurrence of P. parva, as the authors describing this factor as relevant were investigating habitats ranging from 42 to 416 meters above sea level (Carosi et al. 2016). Carosi et al. (2016) also mention water quality in the context of upper reaches of rivers with a salmonid character and low conductivity and nutrient load vs. lowland slow river sections with increased conductivity and nutrient levels. The authors suggested that the occurrence of P. parva may be determined by river hydraulic properties and its isolated and fragmented reaches rather than water quality.

The state of the waterbody (natural vs. regulated) was not recognised as a significant predictor, probably because of the similarity of conditions and communities between numerous natural streams with high sedimentation rates and altered similar streams and canals. Despite this, it has to be highlighted that there is a negative trend in human impact on different types of natural stream habitats in this area (Mrakovčić et al. 2006). Some of the natural streams in this study still had a diverse bed substrate and channel morphology, implying the natural hydrologic regime and supporting different microhabitats (Lau et al. 2006). However, it seems that the high degree of regulation, such as weir and dam construction, leads to a decrease in such natural habitats, resulting in stream bed siltation (Auerswald & Geist 2018) and causing habitat homogenisation (Lau et al. 2006), possibly making it more susceptible to P. parva invasion.

Conclusion

The prolonged presence of P. parva, its high phenotypic plasticity, tolerance and high dispersal rate allowed it to spread across different habitats throughout the lowlands of the Sava River basin. Previous studies confirm that it has become established in many locations in this basin, and this study provides new data on possible established populations and dispersal corridors of the species. As this species is on the European list of invasive alien species of Union concern (European Commission 2016), a management plan has been developed for Croatia, and the results of this study may help in its improvements. The effective containment or control of P. parva abundance is highly unlikely due to its wide distribution in Croatia. Habitats with firm substrate and low species diversity, such as undisturbed streams with natural water regimes and specialised species, seem to be where this species is least expected, so efforts should be directed towards protecting such habitats from further modifications. This study was conducted only in lowland streams and channels of the Sava River basin. It would be interesting to study the occurrence of this species in natural and modified streams at higher elevations and especially in the karst ecosystems of the Adriatic Sea basin, where habitat conditions are much more unstable.