Data was collected as part of a large citizen science project across Switzerland (Knüsel et al. 2024b). Specifically, drinking water providers were asked to sample groundwater for macroinvertebrates. The sampling was conducted in spring catchment boxes, where groundwater is captured passively from horizontal drainage pipes (Fig. 1). At the pipe inlets into the respective spring catchment box, filternets (mesh size 0.8 mm) were attached for one week and then all captured organisms were collected and stored in 80% ethanol. We used two datasets covering different temporal and spatial resolutions for this study.

Figure 1. 

Methodology for sampling groundwater macroinvertebrates at spring catchment boxes (groundwater extracted for drinking water usage). Open spring catchment box from outside (A) and from inside with filter nets attached to the inlets of the drainage pipes (B both modified from Knüsel et al. 2024b). Schematic representation (C modified from Couton et al. 2023a) of the passive groundwater collection (1), inlet into spring catchment box with filter net installed (2), and injection of the water into the local drinking water supply system (3).

Firstly, we used data from 17 inlets spread across Switzerland (Fig. 2A) that were sampled approximately monthly over one year between May 2021 and April 2022. The filtering was conducted during one week per month. The collected macroinvertebrates were morphologically identified to major taxa (Schneider et al. 2023) and then split into three groups. We used two groups covering the most commonly found, obligate groundwater macroinvertebrates (stygobites), namely groundwater amphipods (genera Niphargus and Crangonyx) and groundwater isopods (including genus Proasellus) (Schneider et al. 2023). All other collected macroinvertebrates were combined as a third group (hereafter named ‘other macroinvertebrates’). It consists predominantly of organisms that are occasionally (stygophiles) or accidentally (stygoxenes) entering groundwater ecosystems. For each sample, we calculated the filtering duration based on the given start and end date (using 7 days as default when the start date was not specified by the participant). If the outflow rate was available for each sample, we calculated the total volume of water filtered for the respective samples. In cases where only one measure of outflow rate was missing for a particular sampling date, but the previous and subsequent measurements were available, we used the mean between the two. We additionally extracted data on precipitation, as this affects groundwater recharge (Stoll et al. 2013) and thereby might have an effect on groundwater fauna wash out at spring catchment boxes (see also Di Lorenzo et al. 2018). Therefore we calculated each sampling sites’ local precipitation sum (mm, equivalent to liters per square meter) over two weeks prior to the end date of each filtering period from the RhiresD data provided by the Federal Office of Meteorology and Climatology MeteoSwiss (https://www.meteoschweiz.admin.ch/klima/klima-der-schweiz/raeumliche-klimaanalysen.html). For longer-term seasonal precipitation trends, we calculated the same 2-week precipitation sums from a measuring station in Bern (Federal Office of Meteorology and Climatology MeteoSwiss, https://data.geo.admin.ch/ch.meteoschweiz.klima/nbcn-tageswerte/nbcn-daily_BER_previous.csv) and then compiled the mean and standard deviation per day of year across the years 1990–2020.

Figure 2. 

Sampling sites across Switzerland. Pie charts depict sampling inlets from the monthly dataset with presence of groundwater amphipods, groundwater isopods and other macroinvertebrates marked in dark blue, turquoise, and orange, respectively (A). Absences of the corresponding organism group are marked as empty sectors in the pie charts. Sampling sites from the weekly dataset are shown enlarged (B), filled circles mark inlets with groundwater amphipod detection and empty circles without (remaining macroinvertebrate groups were not considered). The point size represents number of samples. Geodata from swisstopo (permission JA100119).

In the catchment area of the river Töss, 143 inlets were sampled 1–10 times each in subsequent weeks (mean 6.8 weekly samples per inlet) (Fig. 2B). The sampling was conducted in spring 2021 (Studer et al. 2022). All groundwater amphipods were sorted and identified to the species level by sequencing the Folmer fragment of the mtDNA COI gene (Folmer et al. 1994). A detailed protocol of the molecular identification can be found in Studer et al. (2022) and Knüsel et al. (2024b). Organisms other than groundwater amphipods were not further identified in this dataset.

All analyses were performed in R (ver. 4.2.2; R Core Team 2022b). We used a Generalized Additive Modeling (GAM) approach to assess effects of seasonality, outflow, and precipitation on the abundances of groundwater amphipods, groundwater isopods, and other macroinvertebrates in the monthly dataset. Both, outflow and precipitation were found to affect groundwater community assemblages over time in previous studies (Hutchins et al. 2021; Saccò et al. 2021). Hereafter, we refer to the pattern across the sampled year as “seasonality” for simplicity, yet we acknowledge that seasonal and event-based effects cannot be clearly distinguished from 12 monthly samples. For each organism group we included data from all inlets where the respective group was found at least once. Seasonality was modeled as day of year using cyclic cubic regression splines (to ensure matching ends for continuity, Wood 2023), while outflow [L/sec] and precipitation [mm] were modeled with thin plate regression splines (default, Wood 2023). We used splines varying by organism group for each of these covariates and the default 10 knots. Further, we added the inlet identifier as random effect, the organism group as a factor, and an offset to account for varying sampling effort. We ran two models using different log-transformed offsets, one using the sampling duration (in number of sampling days) and one using the total discharge volume (in megaliters). Models were fitted in the R package ‘mgcv’ (Wood 2023) using restricted maximum likelihood (REML) and negative binomial response distribution (Suppl. material 1: R-code). Abundances were then predicted for sampling efforts of one week (representing the commonly used filtering duration) and one megaliter, using varying levels of either seasonality, outflow, or precipitation while keeping the other parameters at their median value. Additionally to the GAM, we calculated temporal autocorrelation and corresponding 95% confidence intervals based on presence–absence of the three organism groups in consecutive samples grouped per inlet using the function acf() from R package ‘stats’ (R Core Team 2022a). We used the type ‘correlation’ and set lag.max to 11 (maximum number of sampling occasions per inlet is 12).

Using the weekly dataset, we calculated detection rates of groundwater amphipods based on their capture history per inlet:

$p_{i,j}=\frac{x_{i,j}}{n_j}$ (Eq. 1)

where p_i,j is the detection rate of species i at inlet j (given presence), x_i,j is the number of samples in which the species i was detected at inlet j, and n_j is the total number of sampling occasions at inlet j (see e.g., McArdle 1990; MacKenzie et al. 2002; MacKenzie 2018). We computed inlet-specific detection rates for four amphipod species, each of which was found at minimally 10 inlets (to ensure sufficiently large sample size for reliable estimates). Additionally, we calculated inlet-specific detection rates for all groundwater amphipods combined. We then calculated the mean detection and corresponding standard deviation across all inlets where the given species or combined amphipods were found at least once by:

${\overline{p}}_i=\frac{{\sum }_{j=1}^{N_i}{p}_{i,j}}{N_i}$ (Eq. 2)

where p̄_i is the mean detection of species i across inlets (given presence) and N_i the number of sampled inlets (given presence of species i). Higher sampling effort is expected to yield more precise detection rate estimates compared to lower sampling effort including few repeated samples per inlet. In a last step, we thus assessed how the detection rate estimates change under varying sampling effort and if they stabilize upon sufficient sampling. Therefore, we calculated the detection rates repeatedly based on the formula above, but using a subset of inlets that were sampled a certain minimum number of times by:

${\overline{p}}_{i,k}=\frac{{\sum }_{j=1}^{N_{i,k}}\left(p_{i,j}\right|n_j\ge k)}{N_{i,k}}$ (Eq. 3)

where k is the threshold of sampling occasions per inlet, ranging from 1 (all inlets included) up to 10 (inlets with 10 sampling occasions included only) and N_i,k corresponds to the number of sampled inlets that fulfill the given threshold criterion.

The sampled inlets showed pronounced differences in faunal composition with respect to the presence/absence of the two stygobite groups and other macroinvertebrates (Fig. 2A). Also, there was a large variety of measured outflow rates and their temporal patterns between inlets (Fig. 3A). Precipitation patterns were comparable across inlet locations, with a large peak in July and a smaller peak between December and February (Fig. 3B). Generally, precipitation in summer 2021 was higher than the long-term average for many areas, while in autumn 2021 it was lower than the long-term average (Fig. 3B, MeteoSchweiz 2022). The outflow rates of some inlets reflected the precipitation peak(s) with certain lags, whereas other inlets had rather consistent outflow rates throughout the sampling period, irrespective of precipitation fluctuations.

Figure 3. 

Environmental covariates along the sampling period (monthly dataset), plotted per inlet A outflow in L/s and B precipitation sum across two weeks prior to the sample collection. For comparison, long-term precipitation mean (black line) and standard deviation (grey shaded area) from Bern (1990 to 2020) are plotted in the background.

We found no significant effect of seasonality on groundwater amphipod and isopod abundances in the GAM (Table 1, Fig. 4). In contrast, seasonality significantly affected the other macroinvertebrates, with predicted abundances peaking in July (Fig. 4). This pattern was mostly driven by EPT (Ephemeroptera, Plecoptera, and Trichoptera) larvae (Suppl. material 1: fig. S1). The outflow rate did not have a significant effect when abundances where modeled per day, but was significant when abundances were modeled per megaliter (Table 1). In the latter case, abundances of all three organism groups per megaliter were predicted to be highest at low inlet outflow rates and to decrease with increasing outflow rates (Fig. 4). We did not identify any significant effect of precipitation on any of the organism groups’ abundances. However, there might be a confounding effect between precipitation and seasonality. In both models, the random inlet effect was highly significant. Still using the monthly dataset, the model with abundances per discharge volume had a higher deviance explained than the model with abundances per sampling duration (62.4% vs. 38.1%; Table 1).

Figure 4. 

Predicted abundances (mean and 95% CI) of groundwater amphipods (dark blue), groundwater isopods (turquoise), and other macroinvertebrates (orange), plotted along gradients of seasonality [day of year, labelled in months], groundwater outflow [liters per second], and precipitation sum over 14 days preceding the sample collection date [millimeters]. GAM predictions were computed based on abundances per discharge volume (megaliter, upper plots) and per sampling duration (week, lower plots). Significance levels based on Table 1 are indicated (ns for p > 0.05, * for p ≤ 0.05, ** for p ≤ 0.01, *** for p ≤ 0.001).

Temporal occurrence varied between organism groups (Fig. 5A) and their capture histories (detection–nondetection) indicated temporal autocorrelation patterns (Fig. 5B). For groundwater amphipods, temporal autocorrelation was significant up to the third lag (= three months). In contrast, the autocorrelation for groundwater isopods was significant up to the tenth lag (Fig. 5B), in line with either very high or very low detection rates across sampling occasions (Fig. 5A). The other macroinvertebrates showed a similar pattern as the groundwater amphipods, with significant and decreasing autocorrelation up to the second lag.

Figure 5. 

A capture histories across the sampling occasions for the three organism groups. Filled tiles mark presence and empty tiles mark absence of the corresponding group, while tiles marked in light grey depict occasions where no sample was taken. Sampling occasions (x-axis, approximated by month for comparability to other plots) consist of monthly one-week filtering periods B temporal autocorrelation of presence–absence of the organism groups. The lag is based on subsequent, monthly samples. The grey area (confined by red dashed line) marks 95% confidence band, autocorrelations larger than the band are significant.

Groundwater amphipods were detected in 55% of the 143 sampled inlets in the Töss catchment. In total, nine species were identified, of which five were only detected in 1–3 inlets each (Table 2). Detection rates varied substantially between species, both in terms of mean and variance across inlets (Fig. 6). Specifically, detection rates of Niphargus fontanus and N. tonywhitteni both indicated a very large heterogeneity across inlets. The same was found for groundwater amphipods combined; detection rates per inlet varied from 0.1 up to 1.0. In contrast, Crangonyx cf. subterraneus was found to have the lowest detection rate and the smallest variance across inlets.

Figure 6. 

Detection rates of amphipods at groundwater extraction sites. Data is shown for four species and for all groundwater amphipods combined. Each point marks the detection rate of a given species at a certain inlet. The filling indicates how many samples were available from the respective inlet to calculate the detection rate, with darker filling indicating more samples. Inlets with at least four sampling occasions were plotted. Boxes give the median and interquartile range (IQR, hinges at 25% and 75% quantiles), with whiskers extending from the box hinges to ±1.5 * IQR.

With increasing sampling effort, we found detection rate estimates to become more conservative (Fig. 7, Suppl. material 1: fig. S2). Because of the procedure how detection rate was estimated, detection rates at higher sampling effort are inevitably lower than estimations obtained at low sampling effort, but allegedly level off upon sufficient sampling. The detection rates of groundwater amphipods combined seemed to stabilize around a mean detection of 0.4, although standard deviations remained large (Suppl. material 1: fig. S3). The lowest mean detection of groundwater amphipods combined (i.e., 0.38, including inlets with at least 6 sampling occasions) would translate to a minimum of 6.3 samples necessary to infer the absence of groundwater amphipods with 95% probability (based on Kéry 2002). However, this calculation is based on the mean detection rate and does not account for the substantial variation in groundwater amphipod detection rates among different inlets. At the species level, we noticed considerable variation in whether detection rate estimates appeared to stabilize or not, given the sampling effort used in this study. Detection rates for N. fontanus and N. tonywhitteni stabilized around a detection of 0.37 and 0.3, respectively, at threshold levels set to 6 samples or higher. However, when setting a too restrictive threshold and thus including very few inlets per species, detection rates seemed to destabilize again. The detection rate for N. auerbachi did not stabilize and instead decreased almost linearly with increasing samples per inlet. As the maximum number of sampling occasions was 10 per inlet, we could not identify detection rates smaller than 0.1. The mean detection rate of C. cf. subterraneus converged towards this lower detection limit (Fig. 7).

Figure 7. 

Detection rates of groundwater amphipods depending on number of samples collected per inlet. Detection was calculated as mean over all inlets where the given species or group occurred and then repeated reducing the data to inlets with a given threshold of minimal sampling occasions. Standard deviations are reported in Suppl. material 1: fig. S3. Point sizes show how many inlets were included for calculating the mean. The area below a mean detection rate of 0.1 is shaded, because the maximum number of samples per inlet is 10, and thus smaller detection rates per species could not be captured.

Organisms from the surface are exposed to strong environmental fluctuations and many of them thus show seasonal patterns in detectability. We identified a peak in macroinvertebrates’ abundances (predominantly shaped by EPT larvae) around July, which may partly reflect the seasonal life cycle of these insects (see also Burgherr et al. 2002; Lubini et al. 2012). Their occurrence might indicate a close hydrological connection between the groundwater and nearby surface waters (Stanford and Ward 1988). In contrast, no seasonal effect was found for groundwater amphipod and isopod abundances, suggesting temporally more consistent detectability in these organisms, possibly reflecting the longer development time of stygobites (Hose et al. 2022). Precipitation did not have any significant effect on detectability of either of the studied organism groups, despite being identified as a key driver of community dynamics and energy flows in other groundwater systems (Saccò et al. 2020; Saccò et al. 2021; Saccò et al. 2022). However, a possible effect might be confounded by the inclusion of a seasonality covariate along with precipitation in the models, particularly as the main peak of EPT detection and (higher than the long-term average) precipitation coincide in July. Since our data is limited to only one year, some effects may be concealed. The non-significant precipitation result might also arise from varying time lags among aquifers, at which precipitation affects groundwater ecosystems (determined by factors such as geology). Generally, the models explained a moderate proportion of the deviance, suggesting that while the included covariates and smooth terms captured some of the variability in the organisms’ abundances, other factors not included may also play a role.

Out of nine groundwater amphipod species found in the Töss catchment, five species occurred at only 1–3 inlets each, a pattern that is characteristic for groundwater communities (Hahn and Fuchs 2009; Trontelj et al. 2009; Knüsel et al. 2024b). The other four species occurred at 10 or more sites each, allowing the calculation of detection rates with sufficient spatial replicas. The mean detection rates of N. fontanus and N. tonywhitteni stabilized when the effort consisted of six or more sampling occasions per inlet, indicating that detection might be overestimated when less repeated samples are available per inlet. In contrast, the detection rate of N. auerbachi did not stabilize and the one of C. cf. subterraneus showed an asymptotic shape towards the lower detection limit of 0.1. For both species, we expect true detection to be lower than 0.1. Results from Knüsel et al. (2024b) using an occupancy modeling approach based on a Swiss-wide dataset confirm the expected low detection probabilities of N. auerbachi and C. cf. subterraneus. Considering all groundwater amphipods combined, the lowest mean detection rate would translate to a minimum of 6.3 samples necessary to infer their absence with 95% probability. In comparison, findings from southwestern Germany indicated that 2–17 samples are necessary to collect 95% of species in vertical wells (boreholes) (Gutjahr et al. 2013). However, we found a high heterogeneity in detection rates between sampled inlets. Consequently, any inlet with an amphipod detection rate below the mean would require more than the computed 6.3 samples to infer the absence of groundwater amphipods. For example, an inlet with a detection rate of 0.1 would require 29 samples to infer the absence with 95% probability (Kéry 2002). However, it is not completely resolved yet to what extent the heterogeneity in detection rates reflects underlying environmental factors and community characteristics, or is an artefact of the sampling methodology (Knüsel et al. 2024b).

We acknowledge that we here only analyzed data from inlets where the respective species were present. Because the occurrence process was not modeled, neglecting false absences where the species is present but not detected likely results in overestimating true detection rates (see e.g., Kéry and Royle 2016; MacKenzie 2018; Knüsel et al. 2024b). Further, one must be careful to apply conclusions on detection rates and sampling effort from this region to other regions, as amphipod occurrence and likely also detection rates might be different, for example in mountainous areas (Knüsel et al. 2024a). The occurrence data analyzed here encloses a specific region in Switzerland, that exhibits a rather high amphipod occurrence and richness (Studer et al. 2022) and might have served as refugium during the Last Glacial Maximum (Knüsel et al. 2024a). If the necessary sampling effort based on detection rates is to be estimated, we advise to carefully consider the sampling area and species expected to occur, the heterogeneity in detection rates among sampling sites, and variability in environmental conditions among sites. Generally, when detection probabilities are low but spatially heterogeneous, monitoring programs should cover a relatively large time frame with sufficient spatial and temporal replicates and take into account possible heterogeneity in detection among species to ensure accurate assessments (Boulinier et al. 1998; MacKenzie et al. 2002).

Studies of groundwater communities in natural springs commonly standardize abundances by discharge volume (Di Lorenzo et al. 2005; Mori et al. 2015). Since the inlets sampled here were built by humans, the effects of outflow rate on organism abundances might not be directly comparable to that of natural springs. We here found a negative effect of outflow when abundances were standardized per discharge volume, which indicates a possible ‘dilution effect’ of the organism density. In other words, the density of organisms might be independent of an inlet’s discharge rate, and standardizing by volume could introduce a bias in the analysis. Hence, it may be sufficient to standardize sampling effort based on filtering duration, independent of the discharge volume of the given inlet, or to add the mean outflow rate as a random effect.

When occurrences are correlated across time, subsequent samples might be biased. Here, we detected temporal autocorrelation in groundwater amphipod occurrence up to three months apart, which is expected given the slow life history and limited mobility of stygobites (e.g., Hose et al. 2022). For standardized monitoring of groundwater amphipods collected at drinking water extraction facilities, given a limited number of sampling occasions per year, we recommend to place sampling occasions evenly throughout the year and to avoid clustering them in a short timeframe, as the latter could lead to correlated false absences. Generally, one could also consider extending the filtering periods from one up to multiple weeks, which would increase the chance of collecting groundwater organisms (e.g., Di Lorenzo et al. 2018). However, organisms will start to decompose at some point, which could hinder the identification of those organisms that remained in the filter nets for longer periods. When deciding on the passive sampling periods, it is important to consider whether the goal is to assess occurrence or abundance, as well as the desired level of taxonomic resolution and the identification method (as decomposition likely affects the accuracy of both, morphological and molecular identification).

Using citizen science data from drinking water providers has shown to be an effective approach to study groundwater macroinvertebrates (Alther et al. 2021; Knüsel et al. 2024b). However, the generally low abundances and low detection probabilities despite high sampling efforts, in combination with large spatial variability in detection, pose challenges. The question remains how monitoring strategies could be scaled up to cover regional to national scales, in order to assess groundwater ecosystem health. One approach might be to combine multiple sampling strategies and to focus on bioindicator species (Koch et al. 2024 and references therein). Apart from citizen science, one novel methodology that is being applied to study groundwater communities is eDNA metabarcoding (Deiner et al. 2017; Couton et al. 2023b). It is fast, less work-intensive, and provides information on the whole groundwater community. While eDNA samples generally integrate biodiversity of larger spatial and temporal scales than traditional samples (Deiner et al. 2016; Carraro et al. 2020), Couton et al. (2023a) found eDNA to be (yet) less effective in detecting groundwater amphipods than directly filtering the organisms. Nevertheless, combining multiple approaches might give complementary insights into the state and dynamics of groundwater ecosystems and could thus be one way forward.