© P. Lemmers et al., Published by EDP Sciences 2023

1 Introduction

Sculpins (Cottidae) are globally widespread benthic fish species occurring in both marine and freshwater environments. They are highly sensitive to water pollution and therefore thought to be useful as indicator to assess the ecological status of water bodies (Utzinger et al., 1998). In Europe, sculpins mostly occur in flowing and relatively cold freshwater systems where they have very limited recolonization possibilities due to natural and anthropogenic dispersal barriers. Isolated populations can be recorded a few kilometres apart in the same stream (Freyhof et al., 2005). A taxonomic revision of the genus Cottus revealed 15 freshwater Cottus species in Europe (Freyhof et al., 2005). Molecular studies confirmed that C. perifretum and C. rhenanus are different species (Czypionka et al., 2012; Fast et al., 2017). The revision prompted Crombaghs et al. (2007) to examine populations recorded as C. gobio in the Netherlands. They found that C. gobio in fact comprised two species, viz. the river bullhead (C. perifretum) and the Rhine sculpin (C. rhenanus). Cottus rhenanus is endemic to the River Rhine and Meuse catchment areas and it is native to Belgium, France, Germany, Luxembourg and the Netherlands (Freyhof et al., 2005; Kottelat and Freyhof 2007; Fig. 1). In the Netherlands, C. rhenanus is a relatively rare species and has a limited distribution. The largest currently known C. rhenanus population in the Netherlands occurs in the River Geul, a tributary of the River Meuse (Crombaghs et al., 2007). The River Geul with its tributaries represents the largest natural fast flowing river system with rheophilic fish assemblages of the Netherlands.

Analyses of museum collections revealed that C. rhenanus was the only sculpin species occurring in the River Geul and its tributaries during the period 1919–1961 (Aben, 2007). Populations of most fish species in this river and its tributaries strongly declined in the past century as a result of agricultural, communal, and industrial waste water pollution. Furthermore, habitat degradation, deterioration of the bed substrate (e.g., an increase of silt and colmation of gravel), and the construction of several dispersal barriers for migratory fish species attributed to the decline (Redeke, 1941; Steenvoorden, 1970; Crombaghs et al., 2000). Former lead and zinc mines upstream in the vicinity of Kelmis (Belgium) adjacent to the River Geul were substantial sources of heavy metal pollution (Leenaers and Rang, 1989; Kucha et al., 1996). The mining industry continued until the 1950s. Heavy metals were also deposited on the river banks and became exposed by erosion of the River Geul. The contaminated sediment was remobilized and transported downstream (Van der Werf, 2014). Land use in the river catchment area changed from extensive to intensive agriculture from the 1950s, which led to leaching of nutrients into surface waters (Provincie Limburg, 2017). It is unclear how this affected the biota. Cottus rhenanus was abundant and widely dispersed in the River Geul and its tributaries in the 1920s (Steenvoorden, 1970). Since the 1950s, however, local populations of the species started to disappear from this river system as a result of direct discharge of industrial effluents and manure into the streams. In 1996, it was reported that C. rhenanus still occurred in the River Geul catchment area, but less abundant than in the 1960s and early 1970s (De Nie, 1996). Crombaghs et al. (2000) described that only a few C. rhenanus specimens were caught in the River Geul in the period 1990–1999, despite an extensive sampling effort. Relatively large populations, however, were still present in two tributaries of this river.

Little was known about the occurrence, densities and habitat use of C. rhenanus in the River Geul and its tributaries since 2000. Only a few studies deal with the autecology of this species (Ovidio et al., 2009; Colleye et al., 2013). Furthermore, its habitat use in relation to its densities and co-occurring fish species has never been characterized. We monitored fish assemblages in its core distribution area during a period of 10 years as part of ecological status assessments of water systems in accordance with the EU Water Framework Directive (WFD; European Commission, 2000). Besides that, we carried out a habitat characterisation. Because of the sensitivity of Cottidae to water quality (Brown, 1989; Utzinger et al., 1997; Besser et al., 2007), we hypothesised that differences in C. rhenanus densities coincide with a changing water quality over the sampling years. Furthermore, based on previous observations (Gubbels, 1997; Crombaghs et al., 2000), we expected that C. rhenanus occurs in shallow parts of the stream only. It was unknown, however, to what extent depth and other habitat variables affect the occurrence of the species. This study integrates data of four monitoring investigations, which allowed to analyse the distribution, habitat, co-occurrence with other species and changes in the population status of C. rhenanus in the River Geul catchment area and the relation with water pollution.

Fig. 1 (A) The endemic range of Cottus rhenanus in green modified from Kottelat and Freyhof (2007), updated with validated distribution data from the Netherlands (NDFF, 2022) and Belgium (Waarnemingen.be, 2022). Country borders are indicated with a solid line. The location of the study area is indicated by a red frame. (B) The River Meuse (dashed dark blue line), River Geul and tributaries (solid light blue lines) on Dutch territory are depicted with only the names of rivers and tributaries that are included in this study. Dams and weirs are represented by a red hexagon.

2 Materials and methods

2.1 Study area

The River Geul in the Netherlands is located in the most southern part of the southern Province of Limburg (Fig. 1). The river originates in Belgium, is 58 km long and flows 21 km through Belgium with the last 37 km flowing through the Netherlands, before it discharges into the River Meuse. The mean width and depth of the River Geul are 10 m and 0.4 m respectively, and the mean velocity 0.3–0.6 m s⁻¹ (Tab. 1; Crombaghs et al., 2000). This river and its tributaries largely remain in their natural morphological state as their courses have not significantly changed (Kadaster, 2022). During the past 100 years, only local adaptations have been applied. The main channels are freely meandering and the channel morphologies are in a moderate to good state, in particular in the tributaries. The River Geul is a hillside chalk stream that has numerous smaller tributary streams that receive water from natural springs, ensuring a fairly constant discharge (Crombaghs et al., 2000). Hillside chalk streams in the Netherlands are relatively shallow and flow relatively fast (Makaske et al., 2020). They are characterised by a channel meandering strongly with short and sharp bends of low amplitude. The channel varies greatly in width and depth. The bed substrate is sandy and often gravelly, with silt occurring locally. Our study area includes four of these streams (i.e., the River Geul and its three tributaries Gulp, Selzerbeek and Mechelderbeek, Fig. 2). The Gulp originates in Belgium and is approximately 22 km long, of which 8 km flows on Dutch territory. The mean width and depth are 3.5 m and 0.4 m, respectively, and the mean velocity 0.3–0.6 m s⁻¹ (Tab. 1; Crombaghs et al., 2000). The Selzerbeek originates in Germany and after 0.5 km flowing in Germany, it forms the border between the Netherlands and Germany. It has a total length of 10.5 km, a mean width of 2.9 m, mean depth of 0.4 m and mean velocity of 0.3–0.4 m s⁻¹. The Mechelderbeek originates in the Netherlands and is approximately 3.75 km long with a mean width and depth of 1.3 m and 0.2 m, respectively. The Mechelderbeek has a mean velocity of 0.2–0.5 m s⁻¹. The streamside vegetation partly consists of trees such as common alder (Alnus glutinosa) and willow (Salix sp.), interspersed with grasses and annual herbs (Fig. 2).

Fig. 2 Photographic representations of the four streams that were investigated in the present study: (A) River Geul, (B) Gulp tributary, (C) Mechelderbeek tributary, and (D) Selzerbeek tributary. Photos were taken in January 2023.

2.2 Fish surveys

We conducted four surveys using handheld backpack electrofishing equipment (DEKA Lord 3000 and Bretschneider EFGI 650) to investigate fish assemblages. We also used fish data from an older survey (Crombaghs et al., 2000), but only for the analyses of the geographical distribution of C. rhenanus. The reason is that the older survey used a different sampling method (dip nets) and collected data on a different scale. Crombaghs et al. (2000) collected data between 1990 and 1999 at a resolution of 1 × 1 km but covered the Geul and tributary streams almost totally, which made the data suitable for geographical analysis. In September–October 2005, 2010, and 2015, we surveyed in total 110 transects in the River Geul and the tributaries Gulp, Selzerbeek and Mechelderbeek using a consistent monitoring methodology required for ecological status assessments of water systems in accordance with the WFD (European Commission, 2000). The fish assemblages were consistently surveyed using electrofishing following the method required by the WFD, and specified for Dutch inland waters by Bijkerk (2010). This description entailed that fishermen surveyed on average 300-metres transects in an upstream direction in one run using handheld backpack electrofishing equipment. A netter followed behind the fishermen to increase catch efficiency. Transects wider than eight metres were surveyed by two fishermen, each accompanied by a person with a dip net. After identification, counting and measuring total fish length (TL, from tip of snout to longer caudal fin lobe, accuracy 1 cm), fish were directly released. The sampling duration per running metre of each transect was equivalent. Additionally, in April–May 2014, we surveyed 73 transects in the same tributaries for fish assemblages and to characterize the physical habitat (see Sect. 2.3).

2.3 Habitat characterization

The density of C. rhenanus was related to six sediment types, the percentage area of large structures and riffles in a transect, and the mean transect depth. In order to do so, during April–May 2014 we sampled 73 transects varying in water depth and sediment types, using backpack electrofishing equipment. The transects were located outside the direct environment of dams and weirs (e.g., not in a slow-flowing impounded section of a stream) and concerned other transects than those sampled for the abovementioned ecological status assessments according to the WFD. After sampling, we estimated the relative abundance of each sediment type (bed substratum) per transect, which yielded the ordinal categories: 0%, <5%, 5–10%, 11–25%, 26–50%, 51–75%, and 76–100%. Distinguished sediment types were silt (mud/sludge), sand, fine granule (0.1–20 mm), pebble (21–65 mm), cobble (66–200 mm), and boulder (>200 mm). Furthermore, we estimated the percentage area of large structures (e.g., woody debris) and riffles (a relatively shallow, rough fast water with a substrate bed consisting of stony substrates such as fine granules, pebbles, and cobbles; Bisson et al., 2006) for each transect, as these large structures may provide shelter habitat for C. rhenanus. Riffles have been distinguished from other fast flowing habitats such as runs by the bed substrate and velocity (a run flows deeper, with lower flow velocity and has a bed substrate consisting of mainly sand with cobbles; Bisson et al., 2006). Finally, we measured stream width (m) and depth (cm). This was done by recording width with a tape measure and depth with a measuring rod at five (beginning, middle, end and two representatives in between) locations in a transect. The mean values were recorded as stream width and depth. The maximal depth that we were able to sample using handheld backpack electrofishing equipment was 150 cm. The mean values with standard deviations and sampled ranges of all tested parameters are included in Table S2.

2.4 Data analyses

Prior to our analyses, we corrected for transect length, width and sampling effort, in order to compare samples from various streams. Densities were calculated as number of individuals per unit of surface area, assuming an effective sampling width of three metres for handheld electrofishing equipment (e.g., transects wider than eight metres were sampled by two fishermen making the effective sampling width six metres). This yielded data on density of fish species per square metre, which was used in further calculations and statistical analyses. Generalized linear models (GLM) were performed to determine whether C. rhenanus densities in four streams varied between the years 2005, 2010, and 2015. Data collected in 2014 were only used for the habitat characterisation analysis and not included in this analysis as the sampling period (April–May) differed from those in the period 2005–2015 (September–October). Cottus rhenanus density (count data) was used as the response variable and sampling year as a fixed factor, separately for each stream. Due to overdispersed data, we used negative binomial regression models of the MASS package (Venables and Ripley, 2002). All assumptions and validations of negative binomial GLM regression models were met as described by Zuur et al. (2009). We used the habitat and C. rhenanus abundance data collected in 2014 for the habitat characterisation. Using GLM's, we tested all habitat variables (transformed from ordinal to continuous scale) in relation to C. rhenanus density. For computing all statistical tests, we used the program R version 4.2.2 (R Core Team, 2022). We created maps with the use of QGIS version 3.22.11 (QGIS Development Team, 2016).

To visualize whether the species assemblages changed across years, we performed non-metric multidimensional scaling (NMDS) based on a Bray-Curtis dissimilarity matrix, using the vegan package (Oksanen et al., 2019). This analysis was performed with abundance data. Fish species observed in less than five transects or represented in total by less than 25 individuals were considered rare and were not included in these analyses. In the dataset, columns represented species (16 in total). Each line represented one transect (110 in total) from one of the sampling years. The stress (a measure for the goodness of fit) for the two-dimensional solution stabilized at 0.174 after 20 iterations. Permutational multivariate analysis of variance (PERMANOVA) was performed to determine if the fish assemblage had significantly changed across years and per stream, based on abundances of all species. First, homogeneity of multivariate dispersion was tested and assumed equal. Subsequently, PERMANOVA was performed with 999 permutations and based on Bray-Curtis dissimilarity (Anderson, 2017).

Water quality data of 12 substances (5 heavy metals, 5 organic substances, 1 insecticide, and dissolved oxygen), was obtained from the water quality database of Water Board Limburg (Water Board Limburg, 2022). The Water Board Limburg collected these data between 1983 and 2022 for ecological status assessments of water systems in accordance with the EU WFD. Data were not available from every year for all of the parameters. We used the most downstream available sampling locations for the analyses. In the cases of the Mechelderbeek tributary, results of nearby sampling locations needed to be aggregated in order to have sufficient measurement points. The coordinates and a map of the sampling locations is presented in Table S3. Graphs of the data and trend lines were also included in electronic supplementary Table S3. The trend line (exponential, linear or logarithmic) with the highest explained variation was selected. As a reference value, the mean of the first five measurements was used. To determine whether concentrations of the 12 analysed substances decreased or increased, and to what extent, we divided the concentration by the reference value to calculate the percentage of concentration decrease or increase. We classified the percentages into classes of concentration decrease (colour scheme greens) and concentration increase (colour scheme reds): 0–25%, 26–50%, 51–75%, and 76–100% (Tab. S3).

3 Results

3.1 Cottus rhenanus in the fish assemblage

During fish sampling in 2005, 2010, and 2015, in total 1918; 4782 and 5369 Cottus rhenanus individuals were caught, respectively. During the additional sampling in April and May 2014, 7633 individuals were caught. The smallest and largest individual were 1 and 15 cm TL, respectively. During the sampling, C. rhenanus was often encountered together with stone loach (Barbatula barbatula) and young-of-the-year brown trout (Salmo trutta) (Tab. 2). Cottus rhenanus was also present in almost all transects where brook lamprey (Lampetra planeri) was caught. In total, 37 fish species from 11 families were caught (Tab. S1). Cottus rhenanus was found co-occurring in transects with 22 native fish species and nine alien species (Tab. 2). The three proportionally most abundant species in the River Geul were B. barbatula, C. rhenanus and Eurasian minnow (Phoxinus phoxinus) (Fig. 3). The three proportionally most abundant species in the three tributaries of the River Geul were B. barbatula, C. rhenanus and S. trutta.

From 2005, 2010, and 2015, the density of C. rhenanus increased in the River Geul and the tributaries Mechelderbeek and Selzerbeek but not in the Gulp. Species that displayed an overall decrease were European eel (Anguilla anguilla), barbel (Barbus barbus), gudgeon (Gobio gobio), roach (Rutilus rutilus) and chub (Squalius cephalus). A NMDS ordination analysis was conducted to visualise the change of the fish assemblages between 2005, 2010, and 2015. The NMDS shows that the fish assemblages have changed across years, as centroids have shifted from the year 2005 to 2015 towards the right (Fig. 4). This is supported by PERMANOVA, with a significant effect of sampling year (Pseudo F = 4.782, p < 0.001) and stronger significant effect of stream (Pseudo F = 7.297, p < 0.001) on species abundance, based on Bray-Curtis dissimilarity. Also, the interaction between year and stream was significant (Pseudo F = 2.443, p < 0.001). Thus, relative species abundance differed across years within each stream and within the entire river system.

Fig. 3 A representation of how species are proportionally part of the fish assemblage in each sampled stream and during each sampling year 2005–2010–2015, ordered by proportion. Species of which in total less than 200 individuals were caught, were accumulated as ‘Other’. Dashed lines demarcate proportions between species. The legend is ordered from top to bottom by highest to lowest cumulative proportion.

Fig. 4 Non-metric multidimensional scaling (NMDS) plot based on Bray–Curtis distance between relative species abundance (rare species not included) by year. Coloured dots represent each sampled transect of a year, black dots represent fish species. Ellipses are drawn around the centroid (the mass of points in multivariate space) of the fish assemblage per sampling year, which represent the 95% confidence interval of the centroid.

3.2 Recolonization and population densities

In the year 2000, individuals of C. rhenanus were found incidentally in the River Geul and large populations of this species were only found in the tributaries Gulp and Selzerbeek (Fig. 5; Crombaghs et al., 2000). During the WFD sampling in 2005, the species was found in the River Geul, and its tributaries Selzerbeek and Gulp, but now more widespread. In 2010, it was clear that the recolonization continued, and we found individuals in the confluence of the rivers Geul and Meuse. In 2015, C. rhenanus occurred widely spread through the River Geul and its tributaries (Fig. 5). In the River Geul, its densities increased significantly over the period 2005–2010 (GLM, z = 2.250, p < 0.05) (Fig. 6). Although the mean densities increased over the sampling years in the Mechelderbeek, too few transects were sampled for statistical comparison. In the Selzerbeek, densities have significantly increased over the period 2005–2010 (GLM, z = 2.426, p < 0.05) and 2005–2015 (GLM, z = 3.168, p < 0.01). In the Gulp, however, the density of C. rhenanus was significantly lower in 2015 than in 2005 (GLM, z = –2.809, p < 0.01) and 2010 (GLM, z = –3.717, p < 0.001).

Fig. 5 Distribution of Cottus rhenanus in the River Geul catchment area over the period 2000–2015. Data for the year 2000 were derived from Crombaghs et al. (2000). Data on the period 2005–2015 were collected in transects during standardized sampling for ecological status assessments according to the European Water Framework Directive and following the method of Bijkerk (2010). Sampled transects where the species was not found are depicted with open symbols in red.

Fig. 6 Mean density of Cottus rhenanus (number of individuals 300 m−2) in the River Geul catchment area per sampling year for each stream. Dots represent the density of each sampled transect. Error bars depict the standard error of the mean. Significant differences between sampling years tested with negative distributed regression models are indicated with letters. Sampling years of the Mechelderbeek could not be statistically compared due to too few sampled transects. *p < 0.05, **p < 0.01, ***p < 0.001.

3.3 Characterisation of the physical habitat

The negative binomial regression models revealed that the density of C. rhenanus was associated negatively to a significant degree with a relatively high abundance of boulders (>200 mm) (GLM, z = –2.642, p < 0.01; Fig. 7). In addition, the density of C. rhenanus was also associated negatively with a relatively high abundance of large structures (woody debris) in a transect (GLM, z = –2.349, p < 0.05; Fig. 7) and associated negatively with the mean depth of a transect (GLM, z = –2.834, p < 0.01; Fig. 7). No significant relations were observed between the C. rhenanus density and abundance of silt, sand, fine granule (0.1–20 mm), pebble (21–65 mm), cobble (66–200 mm), and percentage of riffles in a transect.

Fig. 7 Cottus rhenanus density (number of individuals 300 m−2) in relation to abundance of six sediment types (A-F) in transects, in relation to the percentage area of large structures and riffles in a transect (G-H), and the mean transect depth (I). p-values indicate significance of a relation, determined using negative binomial GLM's. Significant relations were determined between C. rhenanus density and (F) a relatively high abundance of boulders (> 200 mm), (H) large structures (woody debris), and (I) the mean depth of a transect. The solid lines represent the predicted values with corresponding 95% confidence intervals (shown with dashed lines), predicted by the models of a significant relationship.

3.4 Water quality

Table S3 shows that concentrations of cadmium, chrome, copper, and lead reduced by at least 51–71% over the period 1983–2022. Zinc reduced by at least 0–25%. Organic pollution has also declined. Ammonium and nitrogen-total decreased by 51–75% and 0–25% in all streams except the Selzerbeek, respectively. In the latter, ammonium and nitrogen values were relatively high in 2019. Nitrite and phosphate-total decreased by at least 26–50% and 51–75%, respectively. Nitrate decreased by at least 26–50%, with the exception of the Gulp where concentrations increased with regard to the reference value. Dissolved oxygen increased by 0–25% in all streams. The insecticide 4,4-DDT decreased by 76–100% in the River Geul. Only from the River Geul sufficient data on 4,4-DDT were available.

4 Discussion

4.1 Population recovery and recolonization

In this study, we analysed how fish assemblages and more particularly C. rhenanus populations changed between 2005 and 2015. We also related the species occurrence to its habitat. The results of the fish sampling for WFD ecological status assessments of water systems showed that fish assemblages have changed across time and streams. The NMDS shows that the centroids of the fish assemblages in 2005, 2010, and 2015, have shifted to the right. This shift is explained by the fact that a number of critical rheophilic species (including C. rhenanus) have increased during the three sampling campaigns, indicating that water quality has improved significantly (Tab. S3), in contrast to the habitat which largely remained intact (Kadaster, 2022).

It has been demonstrated that for at least a century, C. rhenanus has been present in the River Geul and its tributaries (Aben, 2007). From the second half of the 20th century its population and distribution has declined (Steenvoorden, 1970; De Nie, 1996; Crombaghs et al., 2000). Between 1990 and 1999, C. rhenanus was probably the rarest species in the catchment area (Crombaghs et al., 2000). Our results demonstrate that within a period of only 20 years, C. rhenanus recolonized the River Geul and three tributaries and showed increasing densities. Currently, this species is one of the most abundant fish species. This increase in density is also reflected in the diet of the cormorant (Phalacrocorax carbo) in the river catchment area (Van Rijn, 2013). Analysis of pellets, collected in 2012, showed that this bird preyed mostly on C. rhenanus (39.2% of total number of fishes consumed) and P. phoxinus (21.8%). This percentage was even higher during frost periods, when stagnant (or still) waterbodies were inaccessible for these piscivorous birds and they could only forage in streams. During such conditions C. rhenanus comprised 73.6% of their fish consumption (Van Rijn, 2013). Our fish assemblage surveys were carried out as part of the ecological status assessment required by the WFD. For this purpose, transects were sampled in one run (Bijkerk, 2010), which is not the optimal approach for determining density estimates of Cottus species. For instance, Abdoli et al. (2007) found that two to three consecutive runs are needed to obtain a correct density estimation of Cottidae. Thereby, the densities we determined are most likely underestimates of the actual density. We expect that fish from the tributaries and upstream part of the River Geul in Belgium, have recolonised the Dutch part of the River Geul. The recolonisation of C. rhenanus in the catchment most likely has occurred in a downstream direction, such as at high water discharge events, since obstructions (i.e., dams or weirs) higher than 18–20 cm are impassable for C. gobio s.l. upstream movement and thus act as dispersal barriers (Utzinger et al., 1998). Currently, there are at least eight weirs present in the River Geul (Fig. 1B) that inhibit upstream fish dispersal (Lemmers et al., 2020a). These weirs are all associated with historic water mills such as vertical water wheel mills.

Chemical pollution and habitat degradation in streams have been a major cause for the decline of rheophilic fish species in the Netherlands (Aarts et al., 2004). Since the morphological aspect of the studied streams has remained largely intact during the past 100 years (Kadaster, 2022), the C. rhenanus population increase and recolonization is most likely a long-term effect of water quality improvement measures since the late nineties. Cottidae are sensitive to water pollution such as heavy metals, organic pollution, oxygen deprivation and thermal pollution (Brown, 1989; Van den Brink et al., 1990; Admiraal et al., 1993; Utzinger et al., 1997; Besser et al., 2007). In the Emscher catchment area (North Rhine-Westphalia, Germany), C. rhenanus became locally extinct due to discharge of communal wastewater (Hempel et al., 2020). The organic load in the River Geul continued to increase until 1997, when sewage treatment plants were also constructed in the Belgian (upstream) part of the catchment area and as a result, the general water quality in the River Geul improved (Tolkamp, 2003). Analysis of 12 water quality parameters from the River Geul and the three tributaries between 1983 and 2022 show that pollutant values have overall decreased. All analysed heavy metals decreased in concentration. Organic pollution decreased structurally in all streams except nitrate in the Gulp. This may explain why the density of C. rhenanus decreased here between 2010 and 2015. The concentration of dissolved oxygen has increased in all streams. And the insecticide 4,4-DDT strongly decreased in the River Geul. Improvement of water quality in the River Geul is further supported by the reestablishment of a spirlin (Alburnoides bipunctatus) population in the River Geul since 2018 (Lemmers et al., 2020b). Alburnoides bipunctatus depends on clear, fast-flowing and oxygen-rich streams and is very sensitive to organic water pollution. Populations of other characteristic rheophilic (invertebrate) species are also expected to recover.

4.2 Species co-occurrence and habitat

The fish fauna in the River Geul largely consists of rheophilic species. The river's natural dynamic character, interspersed by anthropogenic historical structures (water mills), provides a wide range of continuously altering and complex microhabitat systems that supports a high diversity of habitat specialists. During samplings in 2005–2015, C. rhenanus was often encountered in transects with B. barbatula and young-of-the-year S. trutta, three-spined stickleback (Gasterosteus aculeatus), G. gobio, S. cephalus and P. phoxinus. These findings are equivalent to other research conducted in England (Prenda et al., 1997), France (Roussel and Bardonnet, 1997) and the Czech Republic (Vlach et al., 2005), where the occurrence of Cottus gobio s.l. has been associated with a similar fish species community. Prenda et al. (1997) determined that B. barbatula and G. aculeatus are associated with shallow, depositional areas of the stream. Salmo trutta, P. phoxinus, G. gobio and C. gobio s.l were classified as generalists. Cottus gobio s.l could not be associated to a certain system of microhabitats, however it was found more abundant in shallow sites (Prenda et al., 1997; Vlach et al., 2005), similar as with C. rhenanus in the present study. According to Gubbels (1997), C. rhenanus in the stream ‘Zieversbeek’ (a tributary of the River Geul in the Netherlands) mostly sheltered during the day and was never found exposed on sediment types such as solely bare sand or silt which is in agreement with our observations in the field in both Belgium, Germany, and the Netherlands. The negative relations between C. rhenanus density and abundance of boulders (>200 mm) and the abundance of large structures (woody debris) in a transect, indicate that these habitats are avoided (see section Predation risk). The presence of woody debris in streams often leads to the accumulation of fine sediments, such as sand and silt, which are evaded probably due to the degree of exposure (Gubbels, 1997) and by the fact that these habitats are preferred by S. trutta (Vlach et al., 2005). A combination of coarser gravel sediment types (viz. fine granule, 0.1–20 mm; pebble, 21–65 mm; cobble, 66–200 mm) may be used because C. rhenanus is a poor swimmer and shelters in coarser gravel sediment types on the bottom (Roussel and Bardonnet, 1997; Vlach et al., 2005). During the reproductive period, the ceilings of small cavities in these types of gravel sediment are used by females to deposit adhesive eggs, which are then guarded by the male (Kottelat and Freyhof, 2007). Furthermore, the species can detect vibrations of the bed substrate (Colleye et al., 2013). Moreover, in these type of streams the flow velocity is higher on coarser gravel sediments leading to higher dissolved oxygen levels (Sear et al., 2017). This can be beneficial for the oxygen requirements of C. rhenanus, but also for the presence of invertebrate species, such as gammarids, which it preys on. However, based on this study it cannot be completely ruled out that correlations between sediment types and C. rhenanus densities do exist, as the methodology we employed for this study may not be optimal for this purpose. More detailed local-scale and, for instance point sampling, studies should be carried out to link these sediments to the species densities.

The regression model showed that C. rhenanus density correlated negatively in a significant way with water depth. Unfortunately, these data might be slightly biased due to a lower catch per unit of sampling effort in deeper sections when using handheld backpack electrofishing equipment (only five out of 36 transects in the River Geul and out of 73 for all transects had a small section that was deeper than 150 cm). At these sites the fish densities might be underestimated. None of the sampled tributaries were too deep for reliable electrofishing. The presence of C. rhenanus in shallow parts of streams has also been demonstrated by other studies (Gubbels, 1997; Ovidio et al., 2009). Gubbels (1997) documented that 90% of the individuals in the stream ‘Zieversbeek’ (The Netherlands) were found between 10 and 15 cm where flow velocity on average was 0.22 m s⁻¹. The water depth varied between 1–130 cm and the maximum depth at which an individual was found concerned 70 cm. Ovidio et al. (2009) determined C. rhenanus home range in the stream ‘Falogne’ (Belgium) which had a mean depth of 15 cm. The established home ranges in periods of 25–27 days varied between 7 and 46 m. Moreover, between 2005 and 2015 we found the highest densities in tributary streams that were more shallow than the River Geul. It is known that sculpin (Cottus gobio) in the deep lakes Lucern and Thun (Switzerland) may occur at depths up to 150–200 m, where it has adapted to its preferred food sources (Lucek et al., 2018). However, as far as we know this data is lacking for C. rhenanus. Based on the presented results it cannot be ruled out that C. rhenanus does not also occur in deeper parts of larger streams elsewhere in its natural range, especially in meandering streams with riffle-pool habitats and during high discharges. We expect that the species will predominantly occur in the relatively shallow parts of larger and deeper streams possibly linked to the occurrence of its preferred food sources. Further studies will have to confirm this.

4.3 Future threats

Contrary to other rivers in the Netherlands, upstream parts of the River Geul do not yet harbour invasive Ponto-Caspian gobies that may strongly compete with Cottidae (Lemmers et al., 2022). The closely related C. perifretum has almost disappeared in less than a decade in the River Meuse main channel due to resource competition with invasive Ponto-Caspian gobies, primarily the round goby (Neogobius melanostomus) (Van Kessel et al., 2011, 2016; Cammaerts et al., 2012). Until 2016, three alien Ponto-Caspian goby species were found in the River Geul, but their occurrence was restricted to the confluence with the River Meuse. However, monitoring revealed that from 2017 until 2021, N. melanostomus started to colonise the River Geul at an approximate rate of 1.2 km year⁻¹ (Lemmers et al., 2022). There is a high risk that the establishment of N. melanostomus in the River Geul catchment area will negatively affect the C. rhenanus population via food and shelter competition (Van Kessel et al., 2011, 2016; Nagelkerke et al., 2018). A telemetric study found that a weir for a hydropower watermill in the River Geul, 5.5 kilometres upstream of the confluence with the River Meuse (Fig. 1B), acts as a solid barrier for upstream and downstream dispersal of fish (Lemmers et al., 2020a). It is assumed that the weir also acts as a barrier for invasive alien fish species, and therefore it was recommended to keep the weir (Gubbels, 2017; Lemmers et al., 2020a). The recent discovered invasion by the signal crayfish (Pacifastacus leniusculus) of upstream parts in the River Geul in Belgium might also negatively impact benthic fish species such as C. rhenanus by predation of eggs and juveniles (Lemmers et al., 2019).

Since Cottidae have low dispersal abilities (Utzinger et al., 1998), several migration barriers upstream of the above-mentioned weir for the hydropower watermill can lead to reduced genetic diversity of relatively isolated metapopulations (Junker et al., 2012). Removal of these fish migration barriers within the River Geul system is recommended to promote genetic exchange, with the exception of the weirs for water mills in the downstream part in order to prevent upstream expansion of invasive alien gobies from the River Meuse. Further water quality improvement and maintaining the habitat is recommended to strengthen populations within its natural range.

5 Conclusions

Analyses of historical and recent data reveal that C. rhenanus populations are resilient. Within two decades, the species has become one of the most abundant and widely spread species in this river catchment area. This event took place: (1) without major morphological restoration of the watercourses, (2) without (re)introduction of the species, and (3) in the presence of barriers for upstream dispersal and in the absence of invasive alien goby species. Therefore, the improvement of water quality remains the most probable explanation for population recovery. Cottus rhenanus is a benthic rheophilic fish species that requires flowing and shallow streams, and occurs on several coarse sediment types. The densities were negatively associated with a relatively high abundance of boulders (>200 mm) in the bed substrate and the presence of large structures such as woody debris. This knowledge can be used for management measures aiming at optimisation of habitat suitability. Cottus rhenanus can be used as an indicator species to assess the ecological status of small streams.

Data sharing

Source data for this study are available from: https://doi.org/10.5281/zenodo.7628451.

Supplementary Material

Table S1. Fish species and amounts that were caught during fish surveys in 2005, 2010, and 2015 in the sampled streams. ‡ indicates an alien fish species.

Table S2. Mean values with standard deviations and sampled ranges of characterised parameters for all sampled transects during the habitat survey in 2014.

Table S3. Analysis of water quality parameters of the River Geul and the tributaries Gulp, Mechelderbeek, Selzerbeek. Data has been collected between 1983 and 2022 (Water Board Limburg, 2022).

Acknowledgements

We sincerely thank two anonymous reviewers and the editor for providing valuable comments on an earlier version of this manuscript, which led to significant improvements. Furthermore, we would like to thank Dirk Heijkers, Paul van Hoof, Lars Huijnen, Jan Jeucken, Vincent de Jong, Nils van Kessel, Lambert Konings, Sascha Krysch, Didier Lemmens, Bart Niemeijer, Johan Pot, Douwe Schut, Guido Thewissen, and the members of the VBC Geul and Zijbeken for assisting with sampling. Maurice Gerits of Water Board Limburg provided us water quality data of the River Geul, for which we are grateful. We kindly thank dr. Jo Ebergen for revising the manuscript for correct use of the English language.

References

All Tables

All Figures

Fig. 1 (A) The endemic range of Cottus rhenanus in green modified from Kottelat and Freyhof (2007), updated with validated distribution data from the Netherlands (NDFF, 2022) and Belgium (Waarnemingen.be, 2022). Country borders are indicated with a solid line. The location of the study area is indicated by a red frame. (B) The River Meuse (dashed dark blue line), River Geul and tributaries (solid light blue lines) on Dutch territory are depicted with only the names of rivers and tributaries that are included in this study. Dams and weirs are represented by a red hexagon.

In the text

	Fig. 2 Photographic representations of the four streams that were investigated in the present study: (A) River Geul, (B) Gulp tributary, (C) Mechelderbeek tributary, and (D) Selzerbeek tributary. Photos were taken in January 2023.
In the text

Fig. 3 A representation of how species are proportionally part of the fish assemblage in each sampled stream and during each sampling year 2005–2010–2015, ordered by proportion. Species of which in total less than 200 individuals were caught, were accumulated as ‘Other’. Dashed lines demarcate proportions between species. The legend is ordered from top to bottom by highest to lowest cumulative proportion.

In the text

Fig. 4 Non-metric multidimensional scaling (NMDS) plot based on Bray–Curtis distance between relative species abundance (rare species not included) by year. Coloured dots represent each sampled transect of a year, black dots represent fish species. Ellipses are drawn around the centroid (the mass of points in multivariate space) of the fish assemblage per sampling year, which represent the 95% confidence interval of the centroid.

In the text

Fig. 5 Distribution of Cottus rhenanus in the River Geul catchment area over the period 2000–2015. Data for the year 2000 were derived from Crombaghs et al. (2000). Data on the period 2005–2015 were collected in transects during standardized sampling for ecological status assessments according to the European Water Framework Directive and following the method of Bijkerk (2010). Sampled transects where the species was not found are depicted with open symbols in red.

In the text

Fig. 6 Mean density of Cottus rhenanus (number of individuals 300 m−2) in the River Geul catchment area per sampling year for each stream. Dots represent the density of each sampled transect. Error bars depict the standard error of the mean. Significant differences between sampling years tested with negative distributed regression models are indicated with letters. Sampling years of the Mechelderbeek could not be statistically compared due to too few sampled transects. *p < 0.05, **p < 0.01, ***p < 0.001.

In the text

Fig. 7 Cottus rhenanus density (number of individuals 300 m−2) in relation to abundance of six sediment types (A-F) in transects, in relation to the percentage area of large structures and riffles in a transect (G-H), and the mean transect depth (I). p-values indicate significance of a relation, determined using negative binomial GLM's. Significant relations were determined between C. rhenanus density and (F) a relatively high abundance of boulders (> 200 mm), (H) large structures (woody debris), and (I) the mean depth of a transect. The solid lines represent the predicted values with corresponding 95% confidence intervals (shown with dashed lines), predicted by the models of a significant relationship.

In the text