Environmental DNA allows upscaling spatial patterns of biodiversity in freshwater ecosystems

Global freshwater biodiversity is declining rapidly, necessitating monitoring tools that provide accurate, fine-scale spatial assessments to inform conservation. Environmental DNA (eDNA) metabarcoding enables efficient, non-invasive detection of multiple taxa from water samples, but downstream transport and decay of DNA complicate local interpretation and spatial upscaling. The study addresses whether coupling riverine eDNA data with a hydrology-based transport and decay model can disentangle upstream contributions, reconstruct local taxon distributions, and generate high-resolution biodiversity maps across entire river networks. Focusing on EPT (Ephemeroptera, Plecoptera, Trichoptera) indicator taxa, the goal is to predict presence and relative abundance at ~1-km stream segments across a 740-km2 basin and validate against traditional kicknet sampling.

Prior work established eDNA as a powerful tool for biodiversity assessment in diverse ecosystems and demonstrated metabarcoding can identify many species from single samples. In rivers, eDNA is transported downstream, integrating information from upstream sources; advection and decay dynamics (influenced by abiotic and biotic factors) affect detection. Traditional biomonitoring often relies on localized point estimates, limiting upscaling. Hydrological principles (e.g., power-law scaling of discharge, width, depth, velocity with drainage area) are well-established and can inform transport modeling. Earlier studies proposed frameworks to estimate species distribution and abundance using eDNA in river networks and highlighted EPT taxa as sensitive water quality indicators with high spatial variability, making them suitable test groups.

Study area and sampling: In June 2016, 61 sites across the Thur catchment (northeastern Switzerland; 740 km2) were sampled at base-flow. At each site, three independent 250 mL water samples were collected for eDNA, filtered on GF/F (0.7 µm) filters, and processed separately as technical replicates. Concurrently, standardized 3-min kicknet sampling of three microhabitats (pooled) collected benthic macroinvertebrates at 60 sites (one lost sample). Laboratory and sequencing: DNA was extracted (Qiagen DNeasy) in a clean eDNA lab. A short COI barcode region was amplified using a dual-barcoded two-step PCR and sequenced on Illumina MiSeq (two runs). Bioinformatics included quality control (FastQC), trimming (usearch), merging (Flash), primer removal (cutadapt), quality filtering (prinseq), error-correction and ZOTU inference (UNOISE3), clustering at 99% identity, ORF checking, and taxonomic assignment using NCBI COI databases (taxize, rentrez, Sintax). Outcomes: 423,043 reads across 183 water samples; 50 EPT genera detected by eDNA; median reads per site 3,406; median reads per genus 1,637. Kicknet detected 47 EPT genera, 36 overlapping with eDNA. River network and hydrology: A 25 m DEM (Swisstopo) and TauDEM (D8) extracted the perennial network using a 0.5 km2 contributing area threshold, yielding 760 reaches (median length 0.78 km; 0.07–3.18 km, 2.5–97.5th percentiles). Reaches were treated as graph nodes with edges along flow. Hydrological variables were estimated using power-law relationships with drainage area A based on four FOEN stations: Q = 0.072 A^1.056 (m3 s−1), w = 1.586 A^0.526 (m), D = 0.073 A^0.463 (m; fitted on three stations), and v = Q/(Dw) = 0.623 A^0.067 (m s−1). Covariates: 35 covariates were used: 18 morphological, geological, and land-cover variables (local vs. upstream-averaged as appropriate) plus 17 geographical cluster indicators capturing spatially structured effects. Covariates were z-normalized; multicollinearity was acceptable (|r|<0.8; VIF<10). Modeling framework (eDITH): eDNA concentration at a sampling site j is modeled as the sum of contributions from all upstream nodes i, accounting for source area, discharge normalization, advection, and first-order decay characterized by timescale τ. Taxon-specific eDNA production rates p_i are proportional to local relative density via a log-link with covariates: p_i = p0 exp(β' X(i)). Expected read numbers N_ij are assumed proportional to modeled eDNA concentration (scaling parameter estimated per genus). Measurement model: Observed read counts per replicate are modeled with a geometric distribution with mean equal to the expected read number, capturing discreteness and overdispersion with a single parameter and accommodating zero-inflation. Calibration and inference: Parameters (β, p, τ) were estimated per genus using Adaptive Metropolis MCMC. Goodness-of-fit was assessed with a bootstrap-based test comparing observed triplet residuals to those from simulated geometric-distributed triplets. Decay times (τ) were inferred per genus. Detection probability and presence: Modeled expected reads yield detection probabilities P_D = N/(1+N); maps of detection probability were thresholded at P_D ≥ 2/3 to produce presence maps. Genus richness maps were obtained by summing presence across genera. Validation and cross-validation: Presence predictions were compared to kicknet data (accuracy = fraction of matching presence/absence at 60 sites). Cross-validation trained models on subsets of sites (AS1: 80%; AS2: 60%; AS3: 40%), preserving stream order proportions via the D'Hondt allocation method, and evaluated loss in goodness-of-fit (GOF) and accuracy relative to the complete model (CM).

• Spatial prediction: The model upscaled eDNA data to produce space-filling maps of relative density, detection probability, presence, and genus richness at ~1-km stream segments across a 740 km2 basin.
• Decay times: Estimated eDNA decay times across genera had a mean of 1.5 h (distribution shown in Fig. 3), corresponding to decay distances on the order of kilometers at 1 m s−1 velocity.
• Biodiversity patterns: Predicted EPT richness was higher in central headwaters; downstream reaches and the Glatt tributary showed lower richness, consistent with geographical covariates having significant negative effects and known pollution in Glatt. Model identified biodiversity hotspots not captured by eDNA or kicknet alone.
• Agreement with kicknet by stream order: In headwaters (Strahler order 1), modeled richness distribution matched kicknet (2-sample Kolmogorov–Smirnov test p=0.43). eDNA alone underestimated richness in low-order reaches (2KS p<0.001). At high-order (≥4) reaches, modeled richness was lower than kicknet (2KS p=0.002), while eDNA and kicknet matched (2KS p=0.98).
• Goodness-of-fit: For all genera, at >90% of sites the null hypothesis that observed reads are geometric with model-predicted mean could not be rejected, indicating high GOF.
• Presence prediction accuracy: Average accuracy of modeled presence vs. kicknet across genera was 82.4% (range 40–100%). Counting false positives as plausible (for elusive genera) increased average accuracy to 92.8% (range 56.7–100%).
• Cross-validation robustness: Training on only 40% of sites led to modest average losses relative to the complete model: GOF loss 6.55% and accuracy loss 4.07%. Validation subsets mainly drove GOF loss; calibration subsets showed slight GOF improvements (e.g., −1.04%).
• Covariate effects and faunistics: Predicted covariate effects and spatial distributions for example genera (Habroleptoides, Protonemura, Athripsodes) aligned with ecological knowledge (e.g., habitat preferences, sensitivity to pollutants, elevation and stream size associations).

Coupling eDNA with a mechanistic hydrology-based transport and decay model (eDITH) resolves the key challenge of downstream advection confounding local presence, enabling reconstruction of fine-scale spatial distributions and biodiversity patterns across entire river networks. The framework transforms sparse point eDNA samples into high-resolution, space-filling predictions of presence and richness, highlighting hotspots that traditional methods or raw eDNA may miss. Validation shows strong goodness-of-fit to read counts and high presence/absence accuracy against kicknet, particularly in headwaters where upstream-integrated signals inform local patterns. The approach is robust to reduced sampling density and leverages universal hydraulic scaling, making it applicable where detailed hydrological data are scarce. However, predictions at the most downstream reaches are less accurate for local presence because downstream samples integrate multiple upstream sources that the model may attribute upstream rather than local. The method offers practical benefits for biomonitoring and conservation, supporting targeted management of sensitive, endangered, or invasive taxa over large networks.

This study demonstrates that environmental DNA, when integrated with hydrology-based transport and decay modeling, enables upscaling of aquatic biodiversity assessments to network-wide, high-resolution maps using limited sampling. Applied to EPT insects in a 740 km2 basin, eDITH produced accurate presence maps and richness estimates at ~1-km stream segments, identified biodiversity hotspots, and matched independent kicknet observations, especially in headwaters. The framework facilitates efficient, non-invasive, large-scale biomonitoring and can be extended to other aquatic and potentially terrestrial taxa in densely drained landscapes. Future work should refine the read-to-concentration relationship via laboratory validation, optimize sampling designs across networks, assess model sensitivity to site placement (especially downstream), incorporate improved primers tailored to target groups, and further disentangle decay vs. deposition processes to enhance local predictions.

• Read count interpretation: The assumed geometric distribution linking expected reads to eDNA concentration captures stochasticity but lacks independent lab validation; primer bias and sequencing platform errors can distort abundance-read relationships.
• Hydrological simplifications: Use of power-law scaling and rectangular cross-sections approximates hydraulics; more detailed hydrodynamics were not modeled. However, authors argue added hydrological complexity may not improve predictions given other uncertainties.
• Deposition vs. decay: eDNA removal processes are aggregated into a single decay time parameter; deposition on substrates is not modeled separately.
• Downstream prediction limitations: Local presence predictions at the most downstream reaches are less reliable due to aggregated upstream signals; model may attribute locally produced eDNA to upstream sources.
• Primer specificity: COI primers targeted broad eukaryotic diversity, potentially under-detecting EPT (conservative bias leading to false absences).
• Sampling design: Only three water replicates per site; eDNA alone underestimated richness in low-order streams. Site placement influences predictive power; further work is needed to optimize sampling strategies.
• Assumptions for scaling: Linear Q~A and daily-mean assumptions hold for the catchment scale considered; deviations could affect transferability to different hydrological regimes.
• Generalization to terrestrial taxa requires additional assumptions about transport from land to streams and source area definitions.