Introduction

The EU Habitats Directive was adopted in 1992 (Council Directive 92/43/EEC) to ensure the conservation of a wide range of rare, threatened or endemic animal and plant species and habitats (European Commission 2023). The Directive’s annexes contain a list of protected species, one of them being the water beetle Graphoderus bilineatus. This species belongs to the family Dytiscidae, a family of predaceous aquatic beetles. It has a western Palearctic distribution (Löbl & Löbl 2017) and during the twentieth century it has considerably declined in Western and Central Europa (Kolar and Boukal 2020). It has also declined in the Netherlands (Cuppen et al. 2006) and today is mostly confined to large contiguous peat areas in the western and northern part of the country (Fig. 1). Here, the species mainly occurs in stagnant ditches and channels, with a depth of between 50 and 150 cm, and clear, fresh water with an often fairly sparse and varied vegetation, consisting of floating and submerged aquatic plants, a low cover of duckweed species (Lemna, Spirodella and Azolla species) and little or no algae (Cuppen et al. 2006). Eutrophication which has led to an increase in algae and duckweed cover in many water systems (Janse and Van Puijenbroek 1998), may have contributed to the long-term decline in the Netherlands.

Fig. 1
figure 1

Current known distribution of Graphoderus bilineatus in the Netherlands (grey), represented as 5 km × 5 km squares, and squares with monitoring locations (black)

Full size image

G. bilineatus is listed in both Annex II and Annex IV of the Habitats Directive. Species listed in Annex IV need to be protected in each EU Member State where they occur. For species listed in Annex II, core areas of their habitat need to be designated as Natura2000 sites which must be protected and managed in accordance with their ecological needs. In the Netherlands, seven Natura2000 sites have been designated for G. bilineatus.

Article 17 of the Habitats Directive prescribes that Member States report the distribution, the distribution trend, as well as the population status and trend of species listed in the annexes once every six years (DG Environment 2017). For Annex II species, the population trend must not only be reported at the national level, but also for the joint designated Natura2000 sites. In earlier Article 17 assessments for the Netherlands, the reported trends in G. bilineatus were based on expert judgement. However, the European Commission expects Member States to monitor protected species and here we report for the first time the trend derived from a monitoring scheme for G. bilineatus. This beetle is too scarce to be monitored by general surveillance schemes such as the EU Water Framework Directive (Council Directive 2000/60/EG). Therefore, a monitoring network tailored to this species has been developed to collect appropriate data.

G. bilineatus is a middle-sized beetle (14–16 mm). The larvae are free swimming predators, predominantly feeding on zooplankton such as small Crustaceae. The pupation predominantly takes place in July and most adults are found soon after. After mid-October their activity diminishes due to decreasing water temperature and the animals start to retreat for hibernation (Koese and Tienstra 2010). The adults are scavengers, so one way to monitor this species is to collect adult beetles using baited traps (Turić et al. 2017). Koese and Cuppen (2006), however, concluded that G. bilineatus could be captured with hand netting equally well. Experienced observers collect more beetles in less time with nets than with baited traps, but this method may be susceptible to observer differences. This variation can be reduced by using only trained professional observers. Yet, any systematic differences between observers can still impact the estimation of a population trend, especially if the number of observers is small. To allow statistical corrections for observer effects (see methods section), we applied double-visit sampling using hand netting, i.e., each sample location has been visited twice, by two independent observers. With hundreds of sample locations to survey, it was not feasible to survey all locations by professional observers every year. It was decided to visit only part of the locations each autumn, so that every six years all selected sample locations were examined once.

The scheme started in 2011 and since then two survey rounds of 6 years have been completed. We assessed trends in occupancy and population size and compared the results obtained for G. bilineatus with those of three other species (Cybister lateralimarginalis, Graphoderus cinereus and Hydrophilus piceus) which were often collected together with G. bilineatus. Here, we address the following questions:

  1. (i)

    Has occupancy and population size changed between the two survey rounds?

  2. (ii)

    Are the changes similar for locations within Natura2000 sites and locations outside Natura2000 sites?

  3. (iii)

    Are there any observer effects?

Material and methods

Selection of sample squares and locations

Prior to the start of the monitoring scheme, the knowledge on the distribution of G. bilineatus in the Netherlands was updated. First, locations of known observations were re-examined and data on aquatic plants and water quality were collected (Cuppen et al. 2006). Subsequently, a Species Distribution Model was constructed using these data which predicted the 1 km × 1 km squares in the Netherlands where the species likely occurred (Sierdsema and Cuppen 2006), and these squares were subsequently surveyed. The updated distribution map contained 85 squares of 1 km2 each with known presence of G. bilineatus during 10 years of active sampling in 2003–2013. The monitoring scheme was initiated 2011 but the final sampling area for the scheme was defined only in 2013. Therefore, 68 out of 85 km2 squares were chosen, covering 80% of the known distribution at the time, including all Natura2000 sites (41 squares). The remaining 27 squares were chosen at random from the pool of 44 squares outside Natura2000 sites, plus 12 squares with a high probability score according to Sierdsema and Cuppen (2006). Most selected squares were in the western and northern peat areas, with a few squares in the higher sandy soil region where the species still occurred (Fig. 1).

The 1 km2 squares have been used to evenly distribute the sample locations across the entire range of the species. Within each 1 km2 square, the first observer selected five sampling locations which were visually judged to be potentially suitable for G. bilineatus. This species mainly occurs in dead ends of waterways, 90° angle bends or inlets in the shoreline, particularly between floating organic matter (Koese and Cuppen 2006). The preference for these places is probably due to the accumulation of animal remnants on which the adults feed and the availability of shelter while resupplying oxygen. Furthermore, vegetation composition, such as the presence of Stratiotes aloides, Hydrocharis morsus-ranae and Lemna trisulca, presence of open water, water clarity, water depth, and natural banks for pupation, determine the suitability of a location for G. bilineatus. After sampling, coordinates and a description of the sampling location were shared with the second observer who surveyed the same location 5–8 days after the first visit. Six years later, the same locations were re-examined twice again.

Each year, 12 squares, representing 60 sampling locations, have been examined. Currently, data are available for 68 squares with 337 locations, each surveyed twice (2011–2016 vs 2017–2022), resulting in a total of 1336 field surveys (Table 1). Most of the surveys were done by Observer1 and Observer2 (Table 1). Observer3 joined during the second round, and a few others, taken together as Observer4, participated only occasionally. Observers alternated as first or second visitor to enable disentangling of observer effects from disturbance effects caused by the first visit. Preliminary analyses, however, showed that the observations between the first and second visit were highly similar, so there does not seem to be any disturbance effect caused by the first visit.

Table 1 Number of field visits to locations per observer. Observer4 represents a group of three persons. Round1: 2011–2016; round2: 2017–2022
Full size table

Field method

Each location was sampled for 15 min using a macrofauna-collecting net with a frame of 30 × 50 cm and mesh size of 3 × 3 mm. The contents of the net were either directly sorted by hand or emptied onto a white blanket and subsequently inspected. After identification, all specimens were released in the water. Sampling was done between mid-August and mid-October to collect adult individuals of all species.

Occupancy models

The ability to detect a species may vary between observers and this potential source of error occurs in most ecological research, even if standardized field methods are being applied (MacKenzie et al. 2018; Schmidt et al. 2023). Recent statistical models can take into account such variation in detection probability in a sophisticated manner, provided data from replicated visits are available (Kéry and Schaub 2012; MacKenzie et al. 2018). Occupancy models use presence/absence data to estimate the number of occupied locations. The replicated visits need to be done in a period of closure, which for occupancy models imply that the state of all locations (occupied or not occupied) has not changed during the period between visits. We assume that in our case the closure period assumption is met because of the short interval between the visits and sampling took place in a time of the year with low population dynamics.

Occupancy models consist of two coupled Generalized Linear Models (GLMs), one to describe the true state of occurrence of locations and one to describe the detection, given that the species does occur (MacKenzie et al. 2018). The GLMs can be extended with covariates to describe variation in occupancy and detection. Our state submodel is specified as follows:

$$z_{it} \sim \,{\text{Bernoulli}}(\psi_{it} )$$
(1)
$$\log it(\psi_{it} )\, = \,alpha.occ\, + \,beta.N2000 \, *N2000_i \, + \,r \, *\left( {t - 1} \right)$$
(2)

where zit is the estimated occupancy (0 or 1), ψit represents the occupancy probability of the species in location i and period t (t = 1 or t = 2), alpha.occ is the intercept in the GLM, beta.N2000 the effect of Natura2000 sites and occupancy trend rate r the effect of period 2 relative to period 1. N2000i indicates whether a location is situated in a Natura2000 site (N2000i = 1) or not (N2000i = 0). We followed Kéry et al. (2009) in the specification of a more extended model in which the overall trend is split into a trend in locations within and outside Natura2000 sites:

$$\log it(\psi_{it} )\, = \,alpha.occ\, + \,beta.N2000 \, *N2000_i \, + \,r_{N2000} *N2000_i * \, (t - 1)\, + \,r_{Non - N2000} *(1 - N2000_i )* \, (t - 1)$$
(3)

where rN2000 and rNon-N2000 represent the independent occupancy trends within and outside Natura2000 sites. By summing the posterior zi values over location cells, we estimated for each species the proportion of occupied locations, i.e., the occupancy estimate, as well as the total number of occupied locations as derived parameters.

The detection submodel estimates detection per visit and is specified as:

$$y_{ijt} \, \sim \,{\text{Bernoulli}}(p_{ijt} \, \times \,z_{it} )$$
(4)
$$\begin{aligned} \log it(p_{ijt} )\, = & \,alpha.p_t \, + \,beta.obs2* \, Observer2_{ijt} \, \\ + \,beta.obs3* \, Observer3_{ijt} \, + \,beta.obs4* \, Observer4_{ij} \\ \end{aligned}$$
(5)

The observations yijt indicate that the species was detected (yijt = 1) or not (yijt = 0) at location i during visit j (n = 1–2 visits) in period t. The observations are considered as the product of occupancy zit and the probability pijt of detecting the species in location i during visit j in period t. The GLM in eq. (5) contains the intercept alpha.pt and beta.obs2, beta.obs3 and beta.obs4 are the fixed effects of the categorical variable Observer2, Observer3 and Observer4 on detection probability, relative to Observer 1. The information on observer identity is represented in eq. (5) by Observer2 (Observer2 = 1 or 0), Observer3 (Observer3 = 1 or 0) and Observer4 (Observer4 = 1 or 0).

Population size models

In principle, double-visit sampling allows the application of binomial mixture models (Royle 2004). These models use abundance data to estimate the population size. The closure assumption here is more rigorous than for occupancy models, because for binomial mixture models the assumption is that local abundance has not changed between visits. Like occupancy models, a binomial mixture model consists of two coupled GLMs. The first GLM describes the local abundance, while the second one describes detection. This model makes it possible to estimate the true number of a species at a location and may produce similar results as obtained from mark-recapture studies (Ficetola et al. 2018). However, recent studies suggest being careful with the application of binomial mixture models as these prove sensitive to model assumptions and may produce unrealistic results (Duarte et al. 2018; Link et al. 2018). Indeed, our estimates of abundance and detection were frequently unrealistic and trend estimates did not correspond to results derived from cross-validation (see below). To obtain more reliable outcomes, we proceeded using a Poisson GLMM, as recommended by Link et al. (2018). In this model we included the same covariates as in the occupancy model, but all in a single GLM:

$$N_{ijt} \sim Pois\left( {\lambda_{ijt} } \right)$$
(6)
$$Log\left( {\lambda_{ijt} } \right)\, = \,\log lamb0\, + \,beta.N2000 \, *N2000_i \, + \,r \, *\left( {t - 1} \right)\, + \,beta.obs2* \, Observer2_{ijt} \, + \,beta.obs3* \, Observer3_{ijt} \, + \,beta.obs4* \, Observer4_{ijt} \, + \,lamb.sigma_{ijt}$$
(7)

Local abundance Nijt is described by a Poisson distribution with mean abundance λijt.

The logarithm of λ is the sum of intercept loglamb0, the effect of Natura2000 sites, population growth rate r and observer effects, all parameters being equivalent to the ones in the occupancy model described above. The Poisson GLMM model also contained a normally distributed random effect lamb.sigmaijt that differed between location, visit and period. This effect allows for larger unexplained variance in abundance than prescribed by the Poisson distribution (overdispersion; see Kéry & Schaub 2012). In our data the overdispersion was caused by many zero values. E.g., G. bilineatus was observed during only one-third of all field surveys. In contrast, high values (> 5 individuals) were found only in 41 surveys (3%). For the other three species, catch rates were equal or slightly greater.

The extended model in which growth rate is split into Natura2000 and other sites reads as:

$$Log\left( {\lambda_{ijt} } \right)\, = \,\log lamb0\, + \,beta.N2000 \, *N2000_i \, + \,r_{N2000} *N2000_i * \, (t - 1)\, + \,r_{Non - N2000} *(1 - N2000_i )* \, (t - 1)\, + \,beta.obs2* \, Observer2_{ijt} \, + \,beta.obs3* \, Observer3_{ijt} \, + \,beta.obs4* \, Observer4_{ijt} \, + \,lamb.sigma_{ijt}$$
(8)

Model runs

We report on the results of the extended occupancy model as described in (3) and the extended Poisson GLMM as described in (8). These models are formulated in a Bayesian framework using the JAGS language (Plummer 2009) and the R-package R2jags in R v. 4.1.3 (R Core Team 2022). We used vague priors for all parameters and ran the models with three parallel Markov chains. All parameters converged according to the Gelman-Rubin Rhat statistic (Rhat < 1.1) after 30,000 iterations, a burn-in of 15,000 and a thinning rate of 150 for the occupancy models. These values were 100,000, 50,000 and 500, respectively, for the Poisson GLMM. For the latter method, we also computed the Bayesian p-value using posterior predictive checks (Kéry & Schaub 2012). We found p-values between 0.41–0.46 for the four species, which suggested reasonable goodness-of-fit of the population trend models. We did not compute a goodness-of-fit measure for occupancy models because posterior predictive checks do not perform well for these models (Kéry & Schaub 2012). A parameter was considered significant when the 95% Bayesian credible interval (CRI) of the posterior sample of a parameter did not include 0.

The results were cross-validated by filtering the data to obtain subsets which can be directly compared without confounding influences. A large subset of locations has been surveyed by the same observer in both periods and that subset allowed checking the trends. A second subset contained only locations surveyed by both Observer1 and Observer2 in the same year. This subset was used to examine the differences between observers.

Results

Occupancy and population trends

Locations within Natura2000 sites had a higher probability to be occupied by G. bilineatus (see beta.N2000 in Table 2). C. lateralimarginalis also occupied locations within Natura2000 sites more often. In contrast, H. piceus had a higher occupancy outside these sites. The total number of locations occupied by G. bilineatus has significantly declined (Fig. 2). The decline in occupancy is stronger at locations outside Natura2000 sites (see rN2000 and rNon-N2000 in Table 2). The other three species had no evidence of a trend in occupancy.

Table 2 Mean ± SD of parameter estimates based on the posterior distributions under the extended occupancy model fitted to four water beetle species
Full size table
Fig. 2
figure 2

Number of occupied survey locations (± SD) per species in two survey rounds derived from the extended occupancy model (see text). Total number of survey locations is 337

Full size image

Generally, the results of the Poisson GLMM echoed the findings of the occupancy model: numbers of G. bilineatus and C. lateralimarginalis were higher at locations in Natura2000 sites than at non-Natura2000 ones, while it is the other way around for H. piceus and G. cinereus. Trends in population size were only found for G. bilineatus and its decline was most severe outside Natura2000 sites (Table 3).

Table 3 Mean ± SD of parameter estimates based on the posterior distributions under the extended Poisson GLMM model fitted to four water beetle species
Full size table

Observer effects

The probability to detect the presence of G. bilineatus at a location was about 0.7 for Observer1 (Table 2). Detection differed between the observers. Observer2 found C. lateralimarginalis and H. piceus at fewer locations than Observer1 (Table 2). Observer3 found C. lateralimarginalis fewer times than Observer1, but more often G. cinereus, while G. bilineatus was found at fewer locations by Observer4 (Table 2). Similar results were obtained with the Poisson GLMM: Observer2 found fewer individuals of C. lateralimarginalis and H. piceus as compared to Observer1 and Observer4 found fewer specimens of G. bilineatus (Table 3). However, no differences between Observer1 and Observer3 were detected. No significant change in detection between the periods was found for any species (see p1 and p2 in Table 2).

Cross-validation of trend estimates and observer effects

In locations surveyed by the same observer in both periods, G. bilineatus was found at fewer locations (n = 327; Chi2 test P = 0.05) and at significantly lower abundance (dependent-samples sign test P < 0.01) in the second period. In this subset, G. bilineatus has declined both at locations within and outside Natura2000 sites, but the decline was stronger outside Natura2000 sites: the species was found at 89 locations within Natura2000 sites in period 1 and at 78 in period 2, while outside Natura2000 sites, these values were 35 and 16 respectively. The other three species did not significantly differ in the number of occupied locations or abundance over time (Chi2 test P > 0.05; dependent-samples sign test P > 0.05). This confirmed the outcome of the models with respect to occupancy and population trends (Tables 2 and 3).

At locations surveyed by both Observer1 and Observer2 in the same year (n = 449 surveys), there was no between-observer difference in the presence of individual species. Yet, the combined presence of all four species was significantly higher for Observer1 (Chi2, P = 0.02), confirming the generally higher detection probability of Observer1 as compared to Observer2 in the occupancy model results. The observer effects obtained by the Poisson GLMM were confirmed at the level of individual species: Observer1 found more C. lateralimarginalis and H. piceus than Observer2 (dependent-samples sign test P = 0.03 and P = 0.02 respectively). According to the cross-validation, Observer1 also found more G. bilineatus and G. cinereus than Observer2 (dependent-samples sign test P = 0.05).

Discussion

Occupancy and population trends

According to the two models applied, G. bilineatus was declining in both occupancy and in population size. These trends were confirmed in the cross-validation test. G. bilineatus had already declined in the Netherlands before 2010 (Cuppen et al. 2006), and this study showed that the decline has continued. The recent decline has been strongest outside Natura2000 sites, but the species continued to decline also in Natura2000 sites. Occupancy and populations numbers were higher within Natura2000 sites than outside (Table 2 & 3), which is expected, because these sites are designated on the base of the occurrence of this species. Between 2013 and 2018, the conservation status of G. bilineatus has been favourable according to the Article 17 assessment for the Netherlands (Valainis 2022), but this status is now outdated. The decline is unfortunate because it means that the protection of this species has so far been unsuccessful.

The other three species are not declining in occupancy nor in population numbers. This indicates that the cause of the decline of G. bilineatus is species-specific. Factors that influence all species equally are less obvious as causes. This applies, for example, to the recent emergence of invasive crayfish species, such as Procambarus clarkii. These crayfish may impact aquatic beetles (Bameul 2013), either through their effect on submerged macrophytes or through resuspension of organic matter. But this is not a likely cause of the decline of G. bilineatus, as the other three species would have reacted similarly. A more likely driver of the decline is ongoing eutrophication, either due to run-off from farmland or through the intake of Rhinewater during years with summer drought. Compared to G. bilineatus, the other three species are considered less sensitive to eutrophication-related hypertrophic conditions (Drost et al. 1992). Also, climate change may have contributed to the decline; the other three species seem less sensitive to higher temperatures than G. bilineatus (Foster et al. 2016; Lenders 2017). The larvae of all four species are carnivorous and typically develop within a few weeks. During this time, larvae rely on seasonally high abundances of prey, such as blooms of Cladocera or larvae of insects with airborne adults (Chironomidae, Ephemoptera). Climate change could alter such synchronized trophic relationships.

Observer effects

Observer1 found higher numbers for all four species than Observer2. Observer1 deposited the samples on a white blanket and Observer2 in a white tray for closer examination and this may have contributed to the different results. This minor methodological difference may have resulted in a significant difference in detection. Both the occupancy model and the Poisson GLMM detected a difference between Observer1 and Observer2 but only in two of the four species, whereas the cross-validation detected differences in all four species.

Observer effects played a bigger role than we had anticipated. Observer effects can play a major role in trend assessment, especially when field methods are non-standardized (Johnston et al. 2018; Van Strien et al. 2022). This study confirms that standardization and working with experienced observers does not necessarily eliminate observer effects (Schmidt et al. 2023), and this may bias trend estimation. This confirms the necessity to statistically adjust for observer effects, as we did here, and the repeated visits to each location offered good opportunities to do so. In cases of single visits, observer and location effects are confounded and the Poisson GLMM would be unable to disentangle them. Even more unfortunate would be if a complete round is surveyed by different observers because that would make it difficult to statistically separate period and observer effects.

Design of monitoring scheme

The scheme has proved to be successful in the assessment of the trends in occupancy and population size. Yet, there are two reasons why we believe the scheme needs to be redesigned. First, at the start of the scheme we assumed that the distribution of the G. bilineatus has been adequately updated, so we used the known distribution of 85 1 km2 squares as the sample framework for the scheme. However, because of increased field efforts G. bilineatus has been found at 134 squares (Fig. 1). The distribution of these squares over Natura2000 vs. non-N2000 sites differs only slightly from the original 85 squares (69–31% vs 74–26%, respectively). We have no reason to assume that our current trend estimates are not representative, but it is better to survey the actual distribution of G. bilineatus.

A second reason has to do with the choice of fixed locations in the scheme. The selection of the five locations within each square is based on visual judgments of the physical characteristics and vegetation composition of waterways. Once selected, the same locations were revisited. However, the suitability of locations is a temporary characteristic which may decline while new suitable locations may emerge elsewhere. The current design of the scheme does not yet take this phenomenon into account.

The reporting obligation under the Habitats Directive has boosted the monitoring of animal species for which hardly any previous monitoring experience existed. Examples include not only the water beetle G. bilineatus, but the Cabrera vole Microtus cabrerae (Oliveira et al. 2023) and the Stag beetle Lucanus cervus (Bardiani et al. 2017). Not surprisingly, a lot of trial and error is needed to arrive at a proper monitoring scheme. For G. bilineatus, we found that observer effects played a bigger role than expected and that the originally planned statistical method (the binomial mixture model) did not work. Our experience also underlines that such schemes may need to be periodically modified to achieve reliable trend estimates.