Estimating species distribution and abundance in river networks using environmental DNA

• species distribution model
• ecohydrology
• proliferative kidney disease

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 byCj^=1Qj¯∑i∈γ(j)AS,i⁡exp−Lijvij¯τpi,[1]where Cj^ [NL−3] is the eDNA concentration at node j predicted by the model; Qj¯ [L3 T−1] is a characteristic (e.g., median) water discharge at node j; γ(j) indexes the set of nodes upstream of j connected by the river network; AS,i [L2] is the source area (or the spatial extent of the habitat) of the target species pertaining to node j; Lij [L] is the along-stream length of the path connecting i to j; τ [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); vij¯ [LT−1] is the average flow velocity along the path connecting i to j; and pi [NL−2 T−1] is the eDNA production at i per unit habitat area and unit time. We assume that pi is proportional to the target species density at node i (12, 24). In the case of sessile species dwelling in riverine habitats, AS,i can be considered equal to the riverbed area of the river stretch i; for terrestrial species, such as certain diatoms that colonize unchanneled areas and have been used to track surface runoff (25), AS,i=AL,i, where AL,i is the directly contributing area (26) to node i (that is, Aj=∑i∈γ(j)AL,i). Eq. 1 is based on the hypothesis that eDNA undergoes first-order exponential decay along the downstream path from the source i to the measurement node j (17). As an alternative characterization of decay, we may introduce the parameterization λ=vij¯τ, 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 vij¯, 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. 1 that take into account the structure of the network.

While not all of the above assumptions are equally valid in the general case (see SI Appendix 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 linkpi=p0⁡exp[βTX(i)],[2]where p0 is a baseline eDNA production value constant in space, X(i) is a vector of covariates evaluated at site i, 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 X(i) 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. 1C). 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.

• Download figure
• Open in new tab
• Download powerpoint

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 Materials and Methods. Note that we explicitly accounted for the possibility that samples with low eDNA concentration may be interpreted as zeros (Fig. 1C) by introducing the nondetection probability ϕj, a monotonically decreasing function of Cj^. The covariates included in vector X(i) 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 i. 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 3A) and frequency of nondetection of the target eDNA (compared with ϕj distributions in SI Appendix, Fig. S1). The median eDNA concentration of the different sampling sites (Fig. 3A) is well reproduced by the model (Pearson’s correlation coefficient r=0.84 for F. sultana, r=0.75 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.

• Download figure
• Open in new tab
• Download powerpoint

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 α=0.05 are characterized by a limited number of positive eDNA detections.

• Download figure
• Open in new tab
• Download powerpoint

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 r 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 SI Appendix, Fig. S2B; maps of covariate values are presented in SI Appendix, Fig. S3. (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 SI Appendix, Fig. S2A.

Predicted maps of eDNA production for F. sultana (Fig. 4A) 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 C and D). 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. 1C), because the estimated values of the decay time (Fig. 3D) 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. 4B) 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 C and D). The correlation between predicted densities of F. sultana and T. bryosalmonae is strong (Fig. 3B) 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. 3C) 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. 3D) 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⁻¹), in agreement with previous findings (19).

• Download figure
• Open in new tab
• Download powerpoint

Fig. 4.

Predicted species distributions. (A–D) Maps of predicted production p (A and B) and concentration C^ (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 SI Appendix. Another strength of our approach is the possibility of applying adequate parameterizations for p in an explicit manner (as in Eq. 2), 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.

• Copyright © 2018 the Author(s). Published by PNAS.