Estimating species distribution and abundance in river networks using environmental DNA.

All organisms leave traces of DNA in their environment. This environmental DNA (eDNA) is often used to track occurrence patterns of target species. Applications are especially promising in rivers, where eDNA can integrate information about populations upstream. The dispersion of eDNA in rivers is modulated by complex processes of transport and decay through the dendritic river network, and we currently lack a method to extract quantitative information about the location and density of populations contributing to the eDNA signal. Here, we present a general framework to reconstruct the upstream distribution and abundance of a target species across a river network, based on observed eDNA concentrations and hydro-geomorphological features of the network. The model captures well the catchment-wide spatial biomass distribution of two target species: a sessile invertebrate (the bryozoan Fredericella sultana) and its parasite (the myxozoan Tetracapsuloides bryosalmonae). Our method is designed to easily integrate general biological and hydrological data and to enable spatially explicit estimates of the distribution of sessile and mobile species in fluvial ecosystems based on eDNA sampling.

Environmental DNA (eDNA), present as loose fragments, as shed cells (1, 2), or in microscopic organisms (3, 4), can be extracted from matrices such as water or soil and used to track the presence of target species or the composition of entire communities (5, 6). Approaches using eDNA for qualitative species detection have proved their value in management and conservation, improving the measurement of biodiversity in a replicable and consistent manner (7) and facilitating the detection of rare, invasive, or parasitic species (8–13).

Environmental DNA in river water carries a record of the species present upstream, but the interpretation of this signal is a complex issue (14). Once released to the environment, eDNA undergoes selective decay. Nucleic acids incur progressive damage [e.g., due to biological activity, temperature, or pH (15, 16)] during hydrological advection, retention, and resuspension (17, 18). These processes result in alterations that affect eDNA detection in environmental samples. The magnitude of the decay is highly dependent on the nature of the flow regime and the substrate type (19). Furthermore, eDNA has polydisperse properties due to its origin from diverse organic sources (e.g., spores, cells, tissues, feces), which complicates the evaluation of decay rates (20). The eDNA sampled at any point within a dendritic network of sources is the outcome of diffuse eDNA release from points upstream, modified by decay processes during transport that are governed by network connectivity, in which each path to the observation point may be described by different hydro-morphological conditions. As a result, while it is straightforward to link a positive PCR test with the presence of the target species at some (unknown) distance upstream, quantification of species densities and the location of populations is currently impossible because, besides a number of potentially confounding factors affecting eDNA shedding [e.g., animal behavior, movement, physiology, and size (21, 22)], it requires consideration of the effects of the dynamics of eDNA transport along river branches and the deconvolution of the hierarchical aggregation of the various network branches. Here, we establish a generally applicable framework to interpret quantitative eDNA point measurements in rivers and relate them to the spatial distribution of the DNA sources, jointly with estimates of the density distribution of the target species throughout the river basin.

The proposed framework stems from fundamental mass balance relationships. It is intended for use in river networks discretized into “nodes,” i.e., river stretches of suitable length within which hydrological conditions, as well as the target species density and hence its eDNA production, can be considered homogeneous (23). Within such nodes, the basis for the spatially explicit model contrasting measured eDNA concentrations is given bywhere [] is the eDNA concentration at node predicted by the model; [] is a characteristic (e.g., median) water discharge at node ; indexes the set of nodes upstream of connected by the river network; [] is the source area (or the spatial extent of the habitat) of the target species pertaining to node ; [L] is the along-stream length of the path connecting to ; [T] is the inverse of the decay rate (i.e., a characteristic decay time) subsuming the damage rates to genetic material experienced during hydrologic transport (18, 19), possibly constant for all eDNA fragments irrespective of the different hydrological and environmental conditions along all paths to the sampling site (24); [] is the average flow velocity along the path connecting to ; and [] is the eDNA production at per unit habitat area and unit time. We assume that is proportional to the target species density at node (12, 24). In the case of sessile species dwelling in riverine habitats, can be considered equal to the riverbed area of the river stretch ; for terrestrial species, such as certain diatoms that colonize unchanneled areas and have been used to track surface runoff (25), , where is the directly contributing area (26) to node (that is, ). Eq. is based on the hypothesis that eDNA undergoes first-order exponential decay along the downstream path from the source to the measurement node (17). As an alternative characterization of decay, we may introduce the parameterization , where is a decay length that is assumed to be constant irrespective of hydrological regimes and heterogeneities in morphological conditions across the watershed (23). This alternative formulation, allowing one to avoid the effective calculation of , is justified by the observation that water velocities in catchments for a large range of flow regimes generally show modest longitudinal gradients, as high velocities potentially prompted by steeper slopes in the upper reaches are limited by increased flow resistance (27). In previous studies, estimated decay lengths range from the order of magnitude of a few meters for experimental flumes (18) to that of kilometers observed in real catchments (19). Such scale dependence seems to indicate a role for the correlation scale of heterogeneous features, relevant to eDNA transport, likely to obey scaling arguments in river networks (28). Suffice here to note that all estimates of decay times (converted to travel lengths as noted above) from field measurements (19, 23) and mesocosm experiments (17, 24) point to values much greater than the mean distance between significant confluences in real catchments, at least in runoff-producing areas (26). This implies that a plurality of sources could be contributing to detectable eDNA concentrations at each measurement site, thus prompting the need to resort to approaches like Eq. that take into account the structure of the network.

While not all of the above assumptions are equally valid in the general case (see for a digression on caveats on model assumptions), the framework represents a flexible general theory of spatially explicit eDNA source tracking. Some of these assumptions may be easily relaxed. For instance, the purely convective treatment of the decay of genetic material outlined above may be the subject of more refined formulations of the transport problem. In particular, travel times from source to measurement site can be made explicitly dependent on the hydrodynamic and geomorphological dispersion induced by the hierarchical nature of the network (27).

To derive estimates of biomass distribution, we coupled the general source area model described above with a species distribution model. This combined approach provides a versatile method to capture the influence of the catchment-wide ecological, hydro-morphological, or geological drivers promoting eDNA production (and reflecting species density) within the defined river stretches. Local eDNA production can be expressed by means of the exponential linkwhere is a baseline eDNA production value constant in space, is a vector of covariates evaluated at site , and is a vector of parameters requiring calibration that identify the effect of such covariates on eDNA production and thereby on the distribution of the target species. Covariates included in the vector will depend on the particular system and data available. Their extent can be either local (e.g., pertaining only to the area directly contributing to a given stretch) or nonlocal, i.e., related to the whole catchment area upstream of the stretch.

We tested the framework with joint field measurements of eDNA concentrations of the myxozoan parasite Tetracapsuloides bryosalmonae and its primary host, the freshwater bryozoan Fredericella sultana, across various locations within the Wigger watershed (Switzerland) (Fig. 1). T. bryosalmonae is the causative agent of proliferative kidney disease (PKD), a high-mortality disease affecting salmonid fish populations. PKD is recognized as one of the leading causes of declines in brown trout populations in Europe, it affects diverse salmonid populations in North America, and it is a major aquaculture disease (29, 30). The water samples used for eDNA detection of both F. sultana and T. bryosalmonae were collected at roughly monthly intervals at 15 sites during 12 mo (one 500-mL sample per sampling occasion and site) (Fig. 1). T. bryosalmonae eDNA is likely to be largely derived from spores shed into the environment. Parasite spores, released into water by infected bryozoans, infect brown trout through skin and gills and proliferate in the kidney. To complete the life cycle, spores infective to bryozoans are excreted in the urine of infected fish. These two types of spores are genetically indistinguishable yet differentiated in terms of function, which poses further challenges for modeling. The T. bryosalmonae eDNA concentration may thus be a product of the genomic contents of the two types of spore originating from very different transport sources, i.e., from an immobile source (bryozoans) coupled with a mobile source (fish). In this particular case, a comparative analysis of field-measured eDNA for both F. sultana (sessile source of eDNA) and T. bryosalmonae (eDNA that could jointly originate from sessile and mobile hosts) will prove particularly instructive as a demonstration of the potential of the proposed framework.

Fig. 1.

Environmental DNA sampling of F. sultana and T. bryosalmonae. (A) Map and color-coded digital elevation map of the study region showing the extracted river network and location of the eDNA sampling sites. (B) Location of the study region within Switzerland. (C) Measured eDNA concentrations of F. sultana and T. bryosalmonae at the 15 sampling sites during the period May 2014–May 2015 (LOQ, limit of quantification).

We implemented and calibrated the model as described in . Note that we explicitly accounted for the possibility that samples with low eDNA concentration may be interpreted as zeros (Fig. 1) by introducing the nondetection probability , a monotonically decreasing function of . The covariates included in vector were local elevation, contributing area, and the fractions of contributing area covered by moraine, peat, or superficial water (e.g., lakes, ponds, or wetlands) upstream of site . These five covariates were chosen as representative of morphological and geological features of the catchment (23). Overall, our model proved to be considerably good at reproducing the observed data, with regard to both observed eDNA concentrations (Figs. 2 and 3) and frequency of nondetection of the target eDNA (compared with distributions in ). The median eDNA concentration of the different sampling sites (Fig. 3) is well reproduced by the model (Pearson’s correlation coefficient for F. sultana, for T. bryosalmonae). According to two-sample Kolmogorov–Smirnov tests (31) (Fig. 2), the null hypothesis that observed positive eDNA concentrations are drawn from the predicted distributions cannot be rejected at the 5% confidence level for 33% (F. sultana) and 71% (T. bryosalmonae) of the sampling sites. At the 1% confidence level, these percentages rise to 73% and 100%, respectively.

Fig. 2.

Comparison between the observed cumulative frequencies of measured eDNA concentrations at the 15 sampling sites for F. sultana (Fs; purple) and T. bryosalmonae (Tb; green) and the cumulative distribution function obtained by the model. values (color coded to match the solid lines) indicate the confidence level at which, according to a two-sample Kolmogorov–Smirnov test, the null hypothesis that the two samples (modeled and observed) come from the same distribution cannot be rejected. Higher values of indicate a better fit. Tested values for were 0.05, 0.01, 0.005, 0.001. N/A, not applicable (i.e., at site 15 no positive values of eDNA concentration for Tb were detected). Note that some of the cases where the two distributions are not different at are characterized by a limited number of positive eDNA detections.

Fig. 3.

Correlation between eDNA concentrations, predicted species density, and covariates. (A) Correlation between observed and modeled eDNA concentrations for F. sultana (Fs; purple) and T. bryosalmonae (Tb; green). Each circle corresponds to a sampling site. x and y values, respectively, are taken from the medians of the observed and modeled cumulative distribution functions displayed in Fig. 2. (B) Correlation between medians of predicted F. sultana and T. bryosalmonae eDNA production. Each circle represents one river stretch. The trend line is displayed in red. Pearson’s correlation coefficients are also reported. (C) Correlations between medians of predicted eDNA production for the two target species and the covariates representing the contributing catchment area and fraction of the upstream catchment covered by moraine. Correlations with other covariates are reported in ; maps of covariate values are presented in . (D) Posterior distributions for the decay time and values of the coefficients associated with the covariates reported on the axis labels. The posterior distributions for the other covariates are shown in .

Predicted maps of eDNA production for F. sultana (Fig. 4) identify the southeastern portion of the watershed as a hotspot for bryozoans. This is mainly due to the positive correlation between the presence of moraines and the production of bryozoan eDNA (Fig. 3 and ). This correlation was uncovered in previous work (23), although it is not clear yet what the underlying causes of this positive association are. The upstream bryozoan reservoir explains the eDNA concentration patterns observed at the downstream sites (Fig. 1), because the estimated values of the decay time (Fig. 3) allow for the detection of eDNA material at distances comparable with the maximum length from source to outlet of the river system. The predicted distribution of T. bryosalmonae (Fig. 4) mirrors that of F. sultana in headwaters and upper sites, while higher values of production are estimated toward the outlet. This is shown by a positive shift in the posterior distribution of the parameter expressing the effect of contributing area (Fig. 3 and ). The correlation between predicted densities of F. sultana and T. bryosalmonae is strong (Fig. 3) and suggests that bryozoans release disproportionately more spores compared with fish hosts. Thus, we suggest that a full description of the spatial distribution of PKD-infected fish (23) might be unnecessary to understand the bulk of the distribution of the parasite sources when the eDNA signal is dominated by locally abundant colonies of infected bryozoans like in the case at hand. We argue that the much stronger correlation observed between predicted densities of T. bryosalmonae and contributing area in comparison with that for F. sultana (Fig. 3) posits that the density of overtly infected (i.e., spore producing) bryozoans tends to increase along downstream directions. In fact, the positive correlation between PKD prevalence and total contributing area in river networks has been shown to be a byproduct of network connectivity (32). The fact that T. bryosalmonae is mostly shed by bryozoans rather than fish seems plausible in this specific host–parasite system, as parasite maturation in fish kidney tubuli is observed relatively rarely, compared with the prolific spore production within large parasite sacs inside the bryozoan host (33). We note further that the spatial match of the bryozoan and fish populations is unlikely to drive this relationship as the biomass of fish is expected to be higher in deeper, more downstream sections (34). Finally, estimated median values of the decay time (Fig. 3) were 4.0 h for T. bryosalmonae (with a 25–75% range of the posterior distribution of 2.7–7.0 h) and 6.9 h (25–75% range: 5.0–11.1 h) for F. sultana, corresponding to decay lengths of 14 km (25–75% range: 10–25 km) and 25 km (25–75% range: 18–40 km), respectively (obtained by assuming an average flow velocity of 1 ms−1), in agreement with previous findings (19).

Fig. 4.

Predicted species distributions. (A–D) Maps of predicted production (A and B) and concentration (C and D) of eDNA for F. sultana (Fs) and T. bryosalmonae (Tb). Estimates were obtained as the medians of the distributions predicted by the model.

Our framework proved capable of interpreting both eDNA data incidentally shed from benthic populations (F. sultana, with likely sources of eDNA from fecal pellets and sloughed cells) and eDNA from spores released into the water (T. bryosalmonae). Although different forms of eDNA may be differently impacted by environmental factors such as temperature and pH, the choice of formulation involving a single parameter expressing the decay time for both species appeared satisfactory for capturing the integrated eDNA transport dynamics at the catchment scale. Further considerations on this aspect are presented in . Another strength of our approach is the possibility of applying adequate parameterizations for in an explicit manner (as in Eq. ), to accommodate the nature of the link between the target species density and its biological and environmental filters along hydrologic pathways, such as the environmental conditions or the density of species with which it interacts. Such parameterization provides a simple and versatile means to assimilate field data and to integrate population or species distribution models.

Accurate field validations of the current assumptions are needed to generalize this framework, and this could be achieved by relatively simple experimental designs. For instance, the displacement and decay of genetic material from nonnative known biomasses placed in well-differentiated positions (say, within a catchment where hydrologic and geomorphologic drivers are known) could be key. Subsequent sampling at downstream sites, where eDNA would be contributed by sources at known distances, could then be used to assess the strengths and weaknesses of each assumption underlying the proposed approach.

Tracking the source area and the local biomass density of target species via downstream eDNA measurement is possible, provided that a suitable spatially explicit framework is used to interpret the field data. Key is accounting for the filtering produced by the progressive damage occurring during hydrological transport and harnessing it to recover spatial information on species distributions. The integration of quantitative eDNA measurements, hydro-geomorphological scaling, and ecological models presented here reveals another direction in ecohydrological studies by unlocking the huge potential of remote monitoring using eDNA.

Materials and Methods

eDNA Data Collection.

Stream water samples were collected in 15 locations (Fig. 1) along the river network of the Wigger watershed, Switzerland. For each site, a total of 21 500-mL samples were taken at approximately biweekly (or monthly during December, January, and February) intervals (except site 5, a connected pond which was artificially drained after 12 samples were taken).

Presterilized (10% bleach followed by UV-B treatment) plastic bottles were used to collect water from the river by submerging the bottle with a gloved hand. The samples were transported to the laboratory on ice and filtered using gentle vacuum within the same day onto 5-cm diameter, 0.45-m pore size individually packaged sterile membrane filters (Merck Millipore). A vacuum pump with a borosilicate glass filtration setup was used and sterilized between each sample in 10% bleach followed by three clean water rinses. Negative controls were created by filtering MilliQ water through a sterile filter at the start and end of each filtration session, as well as once during the filtration (after sample 7). Filter papers were placed in 2-mL bead beating tubes (obtained from the kit described below) and frozen at −80 oC until extraction. Before extraction, filter papers were cut with sterilized scissors to break them up. eDNA was extracted from all filter papers, including controls, using a PowerSoil DNA kit (MO BIO Laboratories) in a dedicated clean laboratory (free of PCR products). The kit includes a bead beating step and a separate inhibitor removal step. The eDNA was eluted in 60 L of Solution C6 and subsequently preserved at −20 oC. Samples were removed from the freezer only for analysis and remained at room temperature for a maximum of 2 h.

Details on F. sultana and T. bryosalmonae eDNA measurement and characterization of the study area are reported in .

Choice of Covariates and Model Settings.

Details on covariates and model settings are reported in . Chosen covariates were checked for multicollinearity. All variance inflation factors for the five considered covariates were below the rule-of-thumb threshold (35) of 10. The three geological covariates were obtained from the vectorized geological map of Switzerland provided by the Swiss Federal Office of Topography (Swisstopo). Covariates were normalized (i.e., linearly transformed into vectors within the range [−1; 1]). shows the values of the covariates plotted on the maps of the river network.

The average velocity along the path between nodes and was calculated as , where indexes the path connecting to , while and are, respectively, the length and the average water velocity of stretch . Velocities were calculated by assuming nearly uniform flow conditions for each stretch and at all times via Manning’s equation. The following further assumptions were made: Water discharges measured by the Swiss Federal Office for the Environment in Zofingen (corresponding to site 3 in Fig. 1) were used to calculate discharges across all network stretches based on the assumed proportionality between discharge and contributing area; river cross-sections were assumed to be rectangular with width much larger than depth; river widths were estimated for all stretches based on aerial images; Manning’s roughness coefficient was taken as m−1/3 s and deemed representative of the flow resistance in the whole river network. Source areas were assumed equal to river-bed surfaces .

Calibration Algorithm.

As sampled eDNA concentrations for both F. sultana and T. bryosalmonae did not show a clear temporal pattern (Fig. 1), measured values were considered as realizations of a random variable whose distribution is constant in time. Measured eDNA concentrations above the limit of quantification (LOQ) observed at site and time were assumed to be log-normally distributed; i.e., , where , or alternatively . Samples with low eDNA concentration may fall below the LOQ or induce sampling errors and be interpreted as indicating absence, owing to the small sampling volumes and the lack of replicates for a sample taken at a given site and time. Here, to account for the number of samples where the target eDNA goes undetected, we make the minimalist assumption that the probability of not detecting eDNA from a sample collected at site is a monotonically decreasing function of the eDNA concentration predicted by the model: , where is a concentration scale (species specific but constant in space) that requires calibration. The likelihood of nondetection at is therefore equal to , while the likelihood of a positive observation reads , where is the standard normal probability density function. Thus, the overall likelihood readswhere indicates the full set of eDNA concentrations observed at any time and site; is the number of sampling sites; and are, respectively, the number of null and positive observations at site ; and spans all positive observations at site . The sampling of the likelihood was performed by means of a Metropolis-within-Gibbs algorithm (36). For all free parameters, prior distributions were chosen as flat. Further details are reported in .