Introduction

Invasive species affect ecosystem functions and services (Vilà et al. 2010) and can have a high economic impact (Cook et al. 2007). Fish species are important invaders, with a high number of species introduced by human activities outside their original biogeographic region (Strayer 2010). Predation, habitat degradation, competition for resources, hybridization and disease transmission are considered key factors in association with invasive alien fish introductions (Gozlan et al. 2010a). Threats caused by alien pathogens, often associated with other invasive species that serve as reservoirs, are poorly understood and deserve more attention in risk assessment and management of biological invasions (Roy et al. 2017). Pathogen pollution caused by spill-over from invasive fish species goes largely undetected in ongoing monitoring schemes, clouding the causes of their detrimental effects on naïve native species (Andreou and Gozlan 2016).

The cyprinid topmouth gudgeon (Pseudorasbora parva) and its associated pathogen Sphaerothecum destruens are an example of such a cryptic, invasive duo. Their invasion initially started via aquaculture trade (Sana et al. 2017), with P. parva serving as an unaffected reservoir. Pseudorasbora parva is listed as an invasive alien species of European Union concern (European Commission 2016). There is a growing need for evidence-based knowledge on its impact and effective management measures.

Pseudorasbora parva originates from China and South-East Asia (Bănărescu 1999). This species has rapidly colonised Europe from countries around the Black Sea, where it was introduced unintentionally in the 1960s (Gavriloaie et al. 2014; Simon et al. 2015). The range expansion of P. parva is regarded as one of the most compelling examples of fish invasion in the world (Gozlan et al. 2010b). As a ubiquitous, sedentary, hardy, omnivorous fish species (Bănărescu 1999), it occurs in both lentic and lotic waters, but is most abundant in lentic conditions such as well vegetated ponds, small lakes, and side channels of rivers (Gozlan et al. 2010b). It adapts well to new man-made habitats in both its introduced range and its native range (Onikura and Nakajima 2012). This high plasticity in life history traits and broad and opportunistic habitat use explains its success as an invader (Gozlan et al. 2010b). Pseudorasbora parva is considered to be a major threat to native fish communities. It was demonstrated that only four years after introduction into a lake, the native fish species composition and community structure had dramatically shifted due to competition (Britton et al. 2007). In addition, Gozlan et al. (2005) found that P. parva is a healthy carrier of S. destruens, a fish pathogen which has the capacity to threaten native fish populations.

Sphaerothecum destruens, previously named rosette agent and first discovered in 1984, belongs to the Rhinosporidiidae, a family of aquatic pathogens that belong to the Dermocystida, an order of the class Mesomycetozoea (Ichthyosporea) (Arkush et al. 2003).This class is phylogenetically classified on the boundary of animals and fungi (Mendoza et al. 2002). All known Mesomycetozoean species live in symbiotic relationships with animals and most of them are aquatic (Glockling et al. 2013). Sphaerothecum destruens lives as intracellular spores within its host that are released into the water via the host’s urine, bile or gut epithelium. Here it forms active, uniflagellate zoospores that can infect other hosts (Arkush et al. 2003; Andreou 2010; Al-Shorbaji et al. 2015). Sphaerothecum destruens tolerates a wide range of temperatures (Andreou et al. 2009) and can survive up to 26 days as a zoospore at 4 °C.

Sphaerothecum destruens became known for causing mortality in salmonid fishes in aquaculture (Arkush et al. 2003). It was subsequently identified as the cause of mortality and spawning inhibition in populations of Leucaspius delineatus (Gozlan et al. 2005; Andreou et al. 2011; Paley et al. 2012). Preliminary examination showed susceptibility of other cyprinids to the parasite as well (Gozlan et al. 2005). An increasing number of studies showed that S. destruens negatively impacts local native fish populations in Europe (Gozlan et al. 2009; Andreou 2010; Peeler et al. 2011; Andreou et al. 2012; Ercan et al. 2015). Abramis brama, Rutilus rutilus and L. delineatus, all species native to Western Europe, are highly susceptible to infection by S. destruens (Al-Shorbaji et al. 2015). However, the impact of P. parva or S. destruens on abundance, body condition and reproduction of native fish in situ are poorly known.

Our study aims to increase the body of knowledge on the impact of P. parva and S. destruens on native fish assemblages. We are addressing the following research questions: (1) how do fish assemblages differ in water bodies with or without P. parva in terms of abundance and reproductive success of native species and community diversity? (2) does P. parva presence have an effect on body condition of native fish? and (3) what is the prevalence of S. destruens in native fish species? As there is a need for quick and reliable methods to detect non-native pathogens, we also address the following research question: (4) can environmental DNA be used as a molecular method to detect S. destruens in situ?

Methods

Data were collected over the period of 2014–2016 for the analyses of effects of P. parva abundance on the fish assemblage, effects of habitat conditions in sampled water bodies, effects of P. parva presence on fish condition, prevalence of S. destruens in native fish and P. parva and finally the presence of S. destruens in water bodies (Table 1).

Table 1 Use of (sub)datasets for each section in the results
Full size table

Fish assemblage and habitat conditions

We selected 54 water bodies (oxbow lakes, shallow lakes and ponds) in river floodplains of the IJssel, Meuse, Nederrijn and Waal River. These water bodies were selected using the following criteria: (a) Potential presence of P. parva according to the Dutch National Database Flora and Fauna (NDFF 2015), (b) No permanent hydrological connection with the main stream or a side channel, (c) Similarity in habitat characteristics (e.g., depth and surface area, for habitat characteristics per sampling site see Online Resource 1), (d) Suitability for sampling with a seine net.

These criteria were set to reduce variance in the fish species composition created by habitat variables, as our aim was to detect effects caused by P. parva. The areas of sampled water bodies ranged from 100 to 80,000 m2. In total 54 sites located in the floodplains were visited (Fig. 1) and sampled using a seine net (21 m long, 2.4 m high, mesh size 4 × 4 to 10 × 10 mm). Fifteen sites could not be sampled sufficiently with this gear type due to high vegetation cover and/or water depth (Fig. 1). Hence, 39 sites were included in the analyses of effect on fish assemblages and body condition (Table 1). Sampling of the fish populations was carried out from October to December 2015. The seine net was used while wading and provided adequate data on juvenile and small fishes in shallow habitats. The sampling area ranged from 0.04 to 82.35% of the surface area of water bodies and was used to calculate fish densities (number of fish m−2). All caught fishes were identified, weighed (accuracy 0.05 g) and their total lengths (TL, from tip of snout to longer caudal fin lobe, accuracy 1 mm) measured in the field. Young of the year (YOY) were distinguished, based on length. Each individual was assigned to being a YOY, based on known YOY thresholds (Kranenbarg et al. 2010). Small fishes (< 35 mm) were pooled for weighting. In case a large number (> 50) of the same species and length were encountered, a representative number was weighed and measured and the remaining individuals were only counted. Subsequently, the fishes were released. Using the obtained data, Simpson index (1-D), Shannon–Wiener index, species richness (total number of fish species) and evenness (a distribution of abundance among species within a community) were calculated (Wilsey and Potvin 2000). Habitat and soil parameters which included coverage percentage of aquatic vegetation, littoral vegetation, and substrate (mud, sand, gravel and rocks), and tree branches in the water were visually estimated. The water transparency (cm) was determined using a Secchi disk (measured vertically). Water temperature (°C), conductivity (μS cm−1) and salinity (PSU) were measured at the site with the use of a Model 30 m (YSI incorporated). A water sample was taken and at the same day pH and alkalinity (eq l−1) were measured in the laboratory. Water samples in polyethylene bottles were stored in the freezer at a maximum storage time of 75 days until analysis. Metal ions were analysed using an ICP analyser (Thermo Electron corporation IRIS Intrepid ΙΙ XDL). Concentrations of nitrate (NO3−), ammonium (NH4+), phosphate (PO43−), chloride (Cl−) and potassium (K+) were determined using an Auto Analyzer 3 system (Bran and Luebbe, Norderstedt Germany).

Fig. 1
figure 1

Sampling sites with Pseudorasbora parva and/or Sphaerothecum destruens and without these species at the floodplains of the rivers IJssel, Meuse, Nederrijn and Waal in the Netherlands

Full size image

Prevalence of Sphaerothecum destruens

To determine infection rates of S. destruens in P. parva and three native fish species, eighty fish were caught in the Teelebeek stream (WGS84 coordinates 51.718, 5.929) in December 2014. In total 20 P. parva, 20 Gasterosteus aculeatus, 20 Barbatula barbatula and 20 Rutilus rutilus individuals were caught by means of dip-net fishing and euthanized with benzocaine (100 mg l−1). The fishes were individually preserved in tubes with ethanol (96%) and stored at − 20 °C. Of each individual fish, kidneys, liver and gonads were dissected out of the body and pooled prior to DNA extraction (Miller et al. 1998). DNA was extracted using the Powersoil DNA kit (MO BIO Laboratories, Carlsbad, CA) according to the manufacture’s protocol. Polymerase Chain Reaction (PCR) was performed on the extracted DNA using primers targeting the ribosomal 18S rRNA gene of S. destruens (Gozlan et al. 2005): (R1) 5′-GAA-GTC-ACA-GGC-GAT-TCG-G-3′, (F2) 5′-ACA-GGG-CTT-TTT-AAG-TCTTGT-3′ and (R2) 5′-ATG-GAG-TCA-TAG-AAT-TAA-CAT-CC-3′. For the first PCR, the F2 and R1 primer combination were used, yielding a 909 base pair product. All PCR results were examined using 1% agarose gel electrophoresis. To improve specificity of the PCR and exclude false positives, a nested PCR was carried out on the previously obtained PCR product using the F2 and R2 primer combination which yielded a 600 base pair product. Both primer combinations were also used separately (not as nested primer pair) for PCR. As a positive control, a plasmid containing a fragment encoding for the 18S rRNA of S. destruens was used (Gozlan et al. 2005, 2009). The PCR for both the first and the nested PCR was done using an initial denaturation of 5 min at 94 °C. Subsequently, 35 cycles of denaturation for 1 min at 94 °C, primer annealing for 1 min using a gradient of 50–60 °C and elongation for 1.5 min at 72 °C were performed. Finally, an elongation step of 10 min at 72 °C was performed. The analysis was performed in 2015.

The nested PCR products that yielded a 600 base pair product were purified with the QIAquick PCR Purification Kit (QIAGEN, Hilden, Germany). PCR products were directly Sanger sequenced by BaseClear B.V. (Leiden, Netherlands). The sequences were analysed using Chromas Lite version 2.4.4 (Technelysium Pty Ltd, Australia) and a BLASTN analysis (http://www.ncbi.nlm.nih.gov) was performed. Sequence alignment was performed using ClustalW (Larkin et al. 2007). Phylogenetic analysis was performed, to show possible genetic diversity within the amplified 18S fragments of the different samples obtained with above described primers targeting the 18S rRNA gene of S. destruens, linked to the location of sampling. The phylogenetic analysis was performed also to validate the outcomes of our eDNA experiment, to record and assess site specificity of populations of S. destruens. The phylogenetic analysis was performed with the software package MEGA 7.0.14 (Tamura et al. 2011) using the Neighbour-joining method and the p-distance model with pairwise deletion and 1000 bootstrap replications. The obtained sequences from the examined fishes were deposited in GenBank under the following accession numbers: MF319869 to MF319872.

Environmental DNA sampling and analyses

Each eDNA sample was composed of 40 water samples from a water body taken with a sterile 50 ml tube, up to 50 cm deep and later mixed in a sterile plastic bag. From the two litre homogeneous mixed water sample, a 15 ml sample was conserved in a sterile 50 ml tube, with 3 M 5 ml sodium acetate (pH 5.2) and 30 ml ethanol (100%) and separately stored at − 20° C (Dejean et al. 2011), in a sealed sterile plastic bag to avoid cross contamination between tubes. During the sampling, for each spot, fresh gloves were used and discarded after sampling. Great care was taken to avoid cross-contamination of the samples during eDNA sampling. DNA extraction was performed per sample separately on a decontaminated work spot. For the PCR, preparations were performed in a PCR UV chamber (Plas-Labs model 825) and negative controls were used (sterile Milli-Q).

Environmental DNA, biological material and solid particles were pelleted by 30 min centrifugation of the 50 ml sample tube at 4000 g at 4 °C. DNA was extracted from the pellet using the Powersoil DNA kit (MO BIO Laboratories, Carlsbad, CA) according to the manufacturer’s protocol. After DNA isolation, a fragment of the 18S rRNA gene was amplified with PCR using F2 and R1 primers, followed by a nested PCR using F2 and R2 primers and Sanger sequenced as described above. Since the primers target the 18S rRNA gene of S. destruens, only a similarity of the amplified product with known S. destruens sequences could be established. Obtained PCR products of a similar size as the amplification product of the positive control of S. destruens were cloned using the pGEM®-T Easy Cloning Vector (Promega, Madison, USA) and transformed in E. coli XL1Blue competent cells. Clones were subtracted to blue/white screening, positive clones were randomly picked, and the plasmids were subsequently isolated using GeneJET Plasmid Miniprep Kit (Thermo Fisher Scientific) and Sanger sequenced by BaseClear B.V. (Leiden, The Netherlands). The analysis was performed in 2016. The obtained sequences from the eDNA samples were deposited (provided with latitude and longitude coordinates) in GenBank under the accession numbers: MF319861 to MF319868.

Statistical analyses

Effect of Sphaerothecum destruens and Pseudorasbora parva on co-occurring fish species

To test for the expected effect of P. parva on native fish populations, generalized linear models (GLM) were used. First, we tested the relation between the abundance of P. parva and the abundance of each of the other fish species that were found living together with P. parva. Therefore, the abundance of each fish species was separately included as the response variable in the GLMs and the abundance of P. parva was included in the model as an explanatory variable. These analyses were performed for six species that were found in ≥ 7 sites. Furthermore, the relation between the abundance of P. parva and each diversity index (Simpson index, Shannon–Wiener index, species richness and evenness) was tested. Since the response variables were in many cases integers, we used the Poisson or negative binomial distribution, depending on over-dispersion. When the response variable data were non-integers, GLMs with normal distributions were applied when criteria for these models were met. The models best fit was chosen based on Akaike’s Information Criterion (AIC). These analyses were conducted with the use of R version 3.6.0 (R Core Team 2019) using the MASS package (Venables and Ripley 2002).

Multivariate analyses of habitat conditions

Multivariate analyses using Principal Component Analysis (PCA) were performed to provide insight in the physico-chemical characteristics of the sampled sites. One PCA was performed on the obtained chemical habitat data, and one PCA was performed on the physical habitat data. The data were first standardized (Z-score) by subtracting the obtained values by the mean. The resulted value was divided by the standard deviation. Permutational multivariate analysis of variance (PERMANOVA) was performed to determine if the difference between the presence or absence of P. parva is explained by ordering the sites along the first principal component (PC1). The PERMANOVA was performed with 999 permutations and based on Euclidean distance. Multivariate analyses were performed with the use of R version 3.6.0 (R Core Team 2019) using the vegan package (Oksanen et al. 2016).

Effect of Pseudorasbora parva on body condition

Body conditions of several species in the presence and absence of P. parva were compared, to test whether presence of P. parva decreases condition of native species. In this analysis, we used data from 54 sampled sites (Fig. 1). For testing differences in the length–weight relationship (also referred to as condition) of individuals, the isometric/allometric growth was determined (Santos et al. 2002). The length–weight relationship of individuals per species is expressed by:

$$ W = a \cdot L^{b} $$
(1)

Here, W is the recorded weight (g), L is the total length (from tip of snout to longer caudal fin lobe, cm), a is a body form coefficient and b is the exponent which indicates relative growth (Froese 2006). The values a and b were obtained by fitting least-squares regression on the natural log length as a function of the natural log weight for each species. Value b is the slope of the regression line in logarithmic form, value a is calculated as the exponent of the intercept. Subsequently, the mean relative weight (Wrm which is an indication of the body condition) per individual was calculated by:

$$ W_{rm} = 100 \cdot \frac{W}{{a_{m} L^{{b_{m} }} }} $$
(2)

in accordance with Froese (2006). Welch Two Sample t tests were used for comparing the relative weight per fish species between the sites with and the sites without P. parva using R version 3.6.0 (R Core Team 2019). However, the sample sizes of these two groups largely differed for most of the fish species and were regarded too unequal (when the sample size of the larger group outnumbered the smaller group by more than half). Therefore, only Carassius gibelio, Gasterosteus aculeatus, Pungitius pungitius and Tinca tinca could be tested.

eDNA related analyses

It is expected that the occurrence of S. destruens coincides with the occurrence of P. parva. To test whether S. destruens is randomly distributed in sites, a Chi square test was performed. The expected values for this test were determined using the number of sites with P. parva (P +), the number of sites without P. parva (P−), and the number of sites with S. destruens, assuming the S. destruens can be found in these sites proportionally to P+ and P−. A Welch Two Sample t test on surface area was performed between sites harbouring S. destruens and sites where it was not detected to test if S. destruens detection was related to surface area. This was tested as eDNA detection limits are affected by P. parva densities (Sana 2016), assuming to be related to surface area of water bodies. This test was performed with the use of R version 3.6.0 (R Core Team 2019).

Results

Fish assemblage and habitat conditions

In total, 12,607 individuals representing 21 fish species were captured, measured and weighed at 39 different sites (Table 2). Gasterosteus aculeatus was the most abundant species with 7,933 individuals found at 27 sites (Table 2). Densities of P. parva ranged from 0.01 to 9.85 individuals m−2 (Table 2, Online Resource 2). Pseudorasbora parva was found to be present at 21 sites (Fig. 1) and co-occurred with 14 other fish species: Abrama brama, Carassius auratus, Carassius gibelio, Cobitis taenia, Cyprinus carpio, G. aculeatus, Gobio gobio, L. delineatus, Neogobius melanostomus, Proterorhinus semilunaris, Pungitius pungitius, Rhodeus amarus, Scardinius erythrophthalmus and Tinca tinca (Online Resource 2). Pseudorasbora parva was the most abundant species based on densities in 10 out of 21 water bodies, while in 8 out of these 21 water bodies, G. aculeatus was the most abundant species. Where present, P. parva made up on average 48% (± SE 6.6%) of the total number of individuals (Online Resource 2). Young of the year (YOY) (≤ 35 mm, Gozlan et al. 2010b) of P. parva were present in 17 out of 21 (81%) water bodies and were not recorded in water bodies with low densities of P. parva of < 0.30 individuals m−2. Principal Component Analysis and PERMANOVA show that the presence or absence of P. parva is not explained by chemical (F = 0.40; p = 0.55) or physical habitat (F = 0.08; p = 0.80) variables (Online Resource 3). Online Resource 1 summarizes all physico-chemical habitat characteristics of sampled sites.

Table 2 Overview of fish data sorted by family
Full size table

Effects on fish assemblage and body condition

Pseudorasbora parva abundance showed significantly negative relations with L. delineatus abundance (Z = − 2.051; p < 0.05) and P. pungitius (Z = − 2.033; p < 0.05) (Fig. 2). Three diversity indices showed significantly negative relations with P. parva abundance, when including P. parva in the calculation of these indices: Simpson’s diversity index of fish assemblages (t = − 3.506; p < 0.01), Shannon–Wiener index (t = − 3.2; p < 0.01) and evenness (t = − 3.104; p < 0.01) (Fig. 2). A positive significant correlation was found with G. aculeatus abundance (Z = 2.03; p < 0.05) (Fig. 2). No correlation was found between P. parva and three other co-occurring species (P. semilunaris, R. amarus and T. tinca). Also, no correlation was found between P. parva abundance and species richness. None of the diversity indices showed a significant relation with P. parva abundance, when P. parva was excluded from the index calculations.

Fig. 2
figure 2

Obtained values (dots) and fitted values (solid line) with 95% confidence intervals (dotted line) of significant correlations by GLMs between Pseudorasbora parva abundance and aLeucaspius delineatus abundance (log scale, Z = − 2.051; p < 0.05), bPungitius pungitius abundance (log scale, Z = − 2.033; p < 0.05), cGasterosteus aculeatus abundance (Z = 2.03; p < 0.05), d Simpson diversity index (t = − 3.506; p < 0.01), e Shannon–Wiener index (t = − 3.2; p < 0.01) and f evenness (t = − 3.104; p < 0.01)

Full size image

YOY co-occurring with P. parva were found for the following species: C. gibelio, G. aculeatus, G. gobio, L. delineatus, N. melanostomus, P. pungitius, R. amarus, S. erythrophthalmus and T. tinca (Table 2).

No differences in condition (mean relative body weight Wrm) in the presence or absence of P. parva of co-occurring fish species could be demonstrated. Only C. gibelio, G. aculeatus, P. pungitius and T. tinca could be tested because sample sizes of other fish species were regarded too unequal (sample size threshold: when the sample size of the larger group outnumbered the smaller group by more than half). Due to the limited number of sites where S. destruens was detected (n = 5), body condition of individuals of other fish species than P. parva by the mean relative weight Wrm was only compared in the presence or absence of P. parva, and not for the presence or absence of S. destruens.

Prevalence of Sphaerothecum destruens in native fish and Pseudorasbora parva

Kidneys, liver and gonads of eighty individuals of fish, comprising four species (each n = 20), were tested for the presence of S. destruens to establish infection rates. For 13 out of 20 samples isolated from G. aculeatus, PCR resulted in fragments matching the expected size of the PCR product for S. destruens, when using the F2 × R1 primer combination (n = 10), F2 × R2 primer combination (n = 2). One of the fishes showed fragments for both primer combinations (n = 1). From these 13 samples, eight were randomly selected for sequencing. BLAST search analysis using BLASTN gave a positive match with S. destruens for five of the G. aculeatus individuals. For four of these samples (22R2, 25R1, 52R1 and 53R1), showing 97–99% similarity with the 18S rRNA gene of S. destruens, phylogeny could be performed. One sample that showed mixed sequences was excluded from the phylogenetic analysis. Prevalence of S. destruens in G. aculeatus is at least 25% (5 out of 20 specimens that were analysed), but might be as high as 50%, as only 8 out of 13 PCR products were sequenced. For at least 18 out of 20 sampled specimens (90%) from P. parva, PCR resulted in fragments matching the expected size of the PCR product for amplification product for 18S rRNA gene of S. destruens, when using either primer combination F2 × R1 or F2 × R2. From these samples, six (showing the highest intensity on agarose gel) were selected for sequencing. Sequence analysis and phylogeny showed that five out of six samples best matched S. destruens using BLASTN analysis (similarity 97–99% and 0.0 E-values). Prevalence of S. destruens in P. parva in this study is at least 25% (5 out of 20), but might be as high as 75%, as only 6 out of 18 PCR products were sequenced. These results were excluded from the phylogenetic analysis as they showed mixed sequences. Sphaerothecum destruens was not found in any of the R. rutilus and B. barbatula individuals.

Presence of Sphaerothecum destruens in water bodies

For five eDNA samples out of 39 water bodies, a nested PCR fragment of the expected size could be perceived, positively matching the 18S rRNA gene of S. destruens, with high similarity (97–99%) and 0.0 E-values, after BLASTN analysis (Online Resource 4).

In all five sites where S. destruens was found, presence of P. parva was confirmed by data of the fish surveys (Fig. 1; Online Resource 2). There is no clear pattern in the presence or absence of P. parva and S. destruens in relation to physico-chemical parameters, as supported by the PERMANOVA (Online Resource 3). Sites harbouring S. destruens were significantly smaller (t = − 3.4359, df = 33.074, p < 0.01) in surface area (\( {\bar{\text{x}}} \) = 218.1 m2; ± SE 82.1 m2; n = 5) than sites where S. destruens was not detected (\( {\bar{\text{x}}} \) = 9964.4 m2; ± SE 2775.2 m2; n = 34). The results suggest that S. destruens is more likely to be detected in smaller water bodies, possibly due to relatively higher DNA concentrations due to lower water volume. Of the 21 sites where P. parva is present, presence of S. destruens could not be confirmed in 16 sites. In addition, there are 18 sites without P. parva or S. destruens. To test whether the occurrence of S. destruens coincides with P. parva, we used a Chi square test with 5 sites with P. parva and S. destruens and 0 sites without P. parva and with S. destruens (observed values). Based on the number of sites with P. parva (21) and without P. parva (18), we assumed a random distribution of S. destruens over the sites and hence the expected values for this test were 2.7 sites with P. parva and S. destruens and 2.3 sites without P. parva and with S. destruens. The test showed that the distribution of S. destruens significantly deviates from the expected values (p = 0.0384, Chi square = 4.29, df = 1). This test suggests that the occurrence of S. destruens coincides with P. parva.

Discussion

Effect of Pseudorasbora parva on native fish species

Pseudorasbora parva is rapidly colonizing new water bodies since its first record in the Netherlands in 1992 (Spikmans et al. 2013). This species is likely to persist in newly colonised areas and reproduces successful in most water bodies it occupies, making further expansion of its distribution highly likely. The number of individuals of P. parva across sites comprised almost half of all fishes that could be caught (Online Resource 2). Dominance in densities of P. parva was earlier shown in Dutch floodplain lakes (Pollux and Korosi 2006), an English lake (Britton et al. 2007) and a Polish river (Witkowski 2009). Densities of P. parva in floodplain lakes with lentic conditions ranged from 0.01 to 9.85 individuals m−2 (Table 2), and are relatively high compared to densities found in lotic conditions in Italy (Carosi et al. 2016). Our study shows a negative relation between the abundance of P. parva and that of L. delineatus and P. pungitius, and the Simpsons diversity index, Shannon–Wiener index and evenness of the fish assemblages in the water bodies studied (Fig. 2). In laboratory conditions, P. parva and S. destruens have been associated with high mortality of L. delineatus (Gozlan et al. 2005; Andreou et al. 2012). Gozlan et al. (2005), Ercan et al. (2015) and Carosi et al. (2016) have recorded severe declines of several freshwater fish species associated with P. parva and S. destruens in the wild. In contrast, a population of L. delineatus in a French lake showed no decline after the introduction of P. parva (Carpentier et al. 2007), but information on the presence of S. destruens was not available for this site. Gasterosteus aculeatus was the only species with abundances positively related to P. parva. Combined with the molecular results, this suggests G. aculeatus may be a healthy carrier (see section Sphaerothecum destruens prevalence in co-occurring fish species).

The non-random selection of our sampling sites was focussed on temporarily isolated, shallow and stagnant water bodies in floodplains along large rivers because this type of water was known to harbour high densities of P. parva (Pollux and Korosi 2006). Our study revealed that this type of water is only favoured by a limited number of limnophilic and eurytopic species, such as L. delineatus, and is less favoured by other fish species, such as G. gobio, L. idus and A. brama. Most co-occurring species were only encountered in a few water bodies and effects of P. parva or S. destruens on their abundance could therefore not be established. Eleven species were found to reproduce in the presence of P. parva (Table 2), including L. delineatus. However, each water body was sampled at only one moment in time. Therefore, local extinction of species after P. parva colonisation might have taken place unnoticed as data on the historic composition of the fish assemblages are lacking. Factors other than S. destruens might also affect diversity in invaded water bodies, such as trophic interactions between P. parva and co-occurring fish species (Britton et al. 2010; Didenko and Kruzhylina 2015), habitat competition and hydrological events, such as flooding and droughts.

No effect of the presence of P. parva or S. destruens on the body condition of four co-occurring species (C. gibelio, G. aculeatus, P. pungitius and T. tinca) was found. Carosi et al. (2016) demonstrated a lower body condition of three native fish species where these species co-occurred with P. parva and suggested that this was caused by trophic overlap between species. Trophic overlap was indeed shown between P. parva and the native R. rutilus (Britton et al. 2010). On the other hand, Rola et al. (2019) showed high trophic plasticity for P. parva and a shift to a lower trophic level in the presence of competitors, resulting in a lower condition of P. parva individuals. Therefore, prediction of the trophic impact of invasive P. parva is difficult.

Sphaerothecum destruens prevalence in co-occurring fish species

Sphaerothecum destruens was found present in both P. parva and G. aculeatus individuals. Prevalence was at least 25% in G. aculeatus, possibly higher (up to 50%). For P. parva, the prevalence rate is at least 25% but possibly higher (up to 75%), since only 6 out of 18 PCR fragments were sequenced and 90% of the tested fishes resulted in a nested PCR fragment of the correct size. Prevalence in P. parva in the same study area was previously assessed at 67 to 74% (Spikmans et al. 2013). The prevalence in G. aculeatus might be an underestimation, as we did not perform sequencing on all samples with PCR product of matching size.

Sphaerothecum destruens is a generalist and known to have a broad host species spectrum (Arkush et al. 2003; Gozlan et al. 2009; Andreou et al. 2012; Andreou and Gozlan 2016) and it is being found in an increasing number of fish species in wild populations: Oncorhynchus tshawytscha, Oncorhynchus kisutch (Arkush et al. 1998), P. parva (Gozlan et al. 2005), L. delineatus (Andreou et al. 2011), Dicentrarchus labrax, Lepomis gibbosus, Oxynoemacheilus sp., Petroleuciscus smyrnaeus, Squalius fellowesii (Ercan et al. 2015), Leuciscus leuciscus, Rutilus rutilus, Squalius cephalus and Salmo trutta (Sana et al. 2018). In experimental or hatchery conditions, several other fish species have proven to act as hosts: Oncorhynchus mykiss (Arkush et al. 1998), Salmo salar (Paley et al. 2012), A. brama and C. carpio (Andreou et al. 2012). Several species mentioned here co-occur with P. parva in the Netherlands.

In our study no S. destruens infection of B. barbatula or R. rutilus could be found, although these samples were taken in the presence of P. parva and G. aculeatus known to be carrying S. destruens. However, Andreou et al. (2012) found high mortality rates in R. rutilus due to S. destruens. Rutilus rutilus might either have been immune to initial infection, able to clear early stages of parasitism or developed a latent infection (carrier state) that is more difficult to detect with PCR (Andreou et al. 2012). This is possibly due to co-existence with P. parva that has inhabited this area for 20 years. Sana et al. (2018) found that B. barbatula, Cottus gobio and G. aculeatus all tested negative as a host for S. destruens, even though S. destruens presence in the sampled stream was proven by eDNA detection.

In our study S. destruens was demonstrated for the first time in a population of G. aculeatus. While L. delineatus and P. pungitius abundances were lower in the presence of P. parva, for G. aculeatus we demonstrated an increasing abundance (Fig. 2). Gasterosteus aculeatus commonly co-occurred with P. parva (71% of water bodies, n = 21) and was the numerically dominant species in 38% of these water bodies. It also frequently co-occurred with S. destruens (80% of water bodies, n = 5). Their body condition was unaffected by the presence of P. parva and/or S. destruens. Presence of YOY indicates G. aculeatus reproduces in sites where S. destruens was found. Scharsack et al. (2007) studied the habitat-specific immune response of G. aculeatus and found that it is able to rapidly adapt its immune system to local conditions to cope with non-trophically transmitted parasites. Species-dependent susceptibility to S. destruens was shown in an experiment with A. brama, R. rutilus and C. carpio, viz. 53%, 37% and 8% mortality, respectively (Al-Shorbadji et al. 2015; Andreou et al. 2012). Since the wide host species range of S. destruens as described above, G. aculeatus might function as an asymptomatic carrier of this parasite. Therefore, G. aculeatus may act as a reservoir for S. destruens and thereby promoting disease outbreaks. This has direct implications for control programmes aimed at S. destruens. Further histopathological screening is needed to confirm the role of G. aculeatus as a reservoir.

Sphaerothecum destruens eDNA detection

Using eDNA, DNA fragments showing 97–99% similarity with the 18S rRNA gene of S. destruens were found to be present in five out of 21 (24%) of the sites where P. parva was found (Online Resource 4). Phylogenetic analysis shows minor genetic differences between more spatially separated sites, site 39 in relation to site 6, 30, 33 and 38 (Fig. 1; Online Resource 4). For example, more samples are required to determine spatial genetic variation of S. destruens in more detail. Environmental DNA is a promising, non-invasive, sensitive and specific method to detect S. destruens (Sana et al. 2018) as the species has a free-living zoospore stage with a broad temperature tolerance (Andreou et al. 2009) and (nested) PCR-amplification is extremely sensitive in detecting this particular zoospore life stage (Mendonca and Arkush 2004; Bass et al. 2015). Sampling techniques can be of great influence on the outcomes of eDNA studies (Hinlo et al. 2017; Goldberg et al. 2016). In our study the water bodies were sampled in late autumn/early winter at low water temperature. Zoospore concentrations and zoosporulation duration of S. destruens are known to increase at temperatures between 4 and 15 °C (Andreou et al. 2009). Detection of S. destruens eDNA is limited in turbid conditions, due to absorption of the DNA to soil particles and detection is also more difficult with low P. parva densities (Sana et al. 2018). Furthermore, we found that sites where S. destruens was detected were significantly smaller in surface area than sites where S. destruens was not detected. This indicates that the chance of detection is higher in the case of sites with a small surface area, likely due to relatively smaller volumes and higher DNA concentrations. In our study S. destruens was strictly found in the presence of P. parva. However, in 16 sites where P. parva was present, no S. destruens was found. Since we amplified a relatively large DNA fragment for detecting S. destruens, due to DNA degradation this might have yielded negatives meaning that our results are conservative. In addition, Sana et al. (2018) used longer primers that were more specific to detect S. destruens and real-time PCR, which is therefore recommended for use in future eDNA studies. Optimizing the eDNA sampling and analysis protocol in accordance with Mendonca and Arkush (2004) and Sana et al. (2018) is important as it can be useful to better understand the scale of S. destruens distribution, providing information for risk analysis and setting up control programs. Furthermore, eDNA can be used to screen fish hatcheries to help prevent further spread of the parasite.

Risk management

Effects of both P. parva and S. destruens easily go undetected (Andreou and Gozlan 2016). Ongoing monitoring programmes for ecological status assessments of water systems according to the EU Water Framework Directive, conservation status of protected areas according to the EU Habitat Directive and the Dutch national long-term fish monitoring programme (De Graaf et al. 2016) are unable to detect declines of native fish associated with P. parva or S. destruens. These programmes do not assess fish populations at a yearly basis and are limited to specific target species and habitats (lotic and large lentic water bodies). It is recommended to set up a long-term fish population monitoring programme in lentic water bodies, combined with an epidemiological survey using molecular techniques, to better understand the role of both P. parva and S. destruens in population changes of native fish species.

Since August 2016, P. parva is listed as an invasive alien species of European Union concern, making introduction, possession, trade and transport illegal in compliance with the EU Regulation 1143/2014. All EU member states are obliged to eradicate or at least control and contain its populations to prevent further spread and to mitigate harmful effects. Eradication of established populations of P. parva is not likely to be successful at large spatial scales in hydrologically connected water bodies. Therefore, cost-effective strategies seem to be limited to control population densities and prevent further spread of P. parva. Eradication of P. parva as a host, however, is not expected to prevent establishment of S. destruens in adjacent fish communities due to spore transmission through water (Al-Shorbaji et al. 2015). Restoring and conserving pristine ecosystems to prevent further loss of biodiversity can decrease transmission of pathogens (Keesing et al. 2010). In the current study P. parva was found to be the dominant species in many waters where it occurred, implying that lowering P. parva densities can limit the transmission of parasites such as S. destruens (Al-Shorbaji et al. 2015). Lemmens et al. (2015) show that the abundance of P. parva can be suppressed by Esox lucius stocking, enhancing the biotic resistance of fish communities. We also suggest focussing attention on minimizing the human assisted dispersal of P. parva via (ornamental) trade and angling associated stockings.

Conclusions

Since its introduction to Europe, P. parva has spread rapidly. Pseudorasbora parva poses risks to biodiversity due to dietary competition with native species, and hosting and spreading of the parasite S. destruens. This study confirms negative effects of P. parva on abundances of L. delineatus, P. pungitius and diversity of fish assemblages (i.e. Simpson index, Shannon–Wiener index and evenness) in floodplain water bodies along the rivers IJssel, Meuse, Nederrijn and Waal. However, a positive correlation was found between the abundance of P. parva and G. aculeatus. This relation, combined with the demonstrated infection of G. aculeatus with S. destruens, suggests G. aculeatus to be an asymptomatic carrier, further increasing the risk of S. destruens spreading to other water bodies. eDNA is a promising method for early detection of S. destruens and facilitates rapid response (e.g. eradication or control of its reservoir host). The ongoing colonisation of P. parva and its function as carrier of the parasite S. destruens will pose a significant risk to the native fish communities. The wide spread Gasterosteus aculeatus might substitute P. parva as a viable S. destruens host, implying that eradication or population control of P. parva does not necessarily lead to minimizing the risk of disease emergence caused by S. destruens.