How much is enough? Examining the sampling effort necessary to estimate mean eDNA concentrations in lentic systems

Abstract

The concentration of eDNA in an environment can provide important ecological information of relevance for management and conservation, but little research has explored optimizing sampling strategies to estimate mean eDNA concentrations in natural environments. Inter-replicate eDNA concentrations often exhibit right-skewed “clustered” or “clumped” distributions, likely due to the stochastic capture of large “aggregate” particles with high eDNA copy numbers. This has important potential implications for modeling the resulting sampling effort necessary to accurately quantify eDNA concentrations. In a previous study, 17–20 Brook Charr eDNA samples were collected from 28 lakes in Québec, Canada. We explored how variation in eDNA concentrations within a lake was affected by several habitat characteristics. We then conducted a power analysis to determine the sampling effort (“minimum n”) necessary to accurately quantify mean lake eDNA concentrations and, using simulations, explored how a bimodal distribution of eDNA particle copy count could affect inter-replicate variability. The median sample size such that 90% of sample mean estimates were within 20% of the “true” mean was 12.5; a sample size of 20 was sufficient to quantify mean concentrations in 21/28 lakes. We found no evidence that temperature or lake size impacted sample variability. We also found that variance among replicates was non-linearly related to mean lake eDNA concentration across years: variability was lowest at low and high concentrations and highest at intermediate concentrations. We hypothesize that this resulted from the stochastic capture of large “aggregate” particles at intermediate concentrations; at low concentrations, aggregates were likely rarely captured and at high concentrations may represent a consistent component of total eDNA. Simulations demonstrated that these patterns can emerge from some bimodal eDNA particle “size” distributions. Overall, we conclude that sampling efforts in many previous studies (notably including the authors' own) were potentially low, emphasizing the need to increase spatial replication in lentic surveys.

INTRODUCTION

Environmental DNA (eDNA) has revolutionized the capacity to detect aquatic organisms (Deiner et al., 2017) and offers significant advantages over conventional methods, including increased sensitivity and reduced impact on study organisms (Stewart, 2019). However, emergent research has also demonstrated that the concentration of eDNA present in an environment can provide important ecological information of potential relevance for management and conservation. A growing body of research, for example, has demonstrated a consistent positive correlation with the concentration of eDNA in an environment and organism abundance (Lacoursière-Roussel, Côté, et al., 2016; Rourke et al., 2021; Spear et al., 2021; Thomsen et al., 2012; Yates et al., 2019).

Very little research, however, has investigated strategies for optimally assessing the mean eDNA concentration in an environment. This issue is of particular importance for large aquatic systems, where obtaining spatially paired local estimates of absolute organism abundance can be challenging due to the movement of organisms within the system (Yates, Cristescu, & Derry, 2021). Estimating robust absolute organism abundance in lentic systems, for example, often necessitates capture–mark–recapture (CMR) techniques, and such abundance estimates are limited to the scale at which organism movement is largely restricted (Schwarz & Seber, 1999); in lentic systems, this often results in “whole-lake” estimates of organism abundance. When comparing eDNA concentrations to organism abundance, it is crucial to match the scale at which abundance estimates were obtained (Chambert et al., 2018). Such “whole-lake” estimates of organism abundance should therefore be compared only to “whole-lake” estimates of eDNA concentrations—that is, abundance (standardized for lake volume/surface area) should be compared to the mean eDNA concentration across replicates within a lake. Several studies have demonstrated strong correlations between mean lake eDNA concentrations and organism abundance (Klobucar et al., 2017; Spear et al., 2021; Yates, Glaser, et al., 2021). Quantifying the minimum effort necessary to accurately and reliably estimate the mean concentration of eDNA in a lentic system is thus critical for determining the effectiveness of eDNA sampling for such inferences.

When sampling eDNA, however, it is common to observe several samples with disproportionately high eDNA concentrations relative to other samples from the same environment (Lacoursière-Roussel, Côté, et al., 2016; Wilcox et al., 2015), manifesting in “skewed” or “clustered” patterns of inter-replicate variability (Chambert et al., 2018; Lacoursière-Roussel, Côté, et al., 2016; Wilcox et al., 2015). Many researchers have hypothesized that the capture of “clumps” or “aggregate” particles with high eDNA copy number (e.g., a particle of feces) can account for these observations (Lacoursière-Roussel, Côté, et al., 2016; Wilcox et al., 2015). Such “aggregates” may break apart rapidly into smaller particles with low eDNA copy count or settle quickly due to larger particle size (Turner et al., 2014; Wilcox et al., 2015). Consequently, the ratio of “large eDNA aggregates” to “small eDNA particles” in an environment may be very low but could potentially impact the skew and variability of eDNA concentrations observed across spatial replicates due to the stochastic capture of large particles.

The distribution of eDNA in an environment can thus often be described as “overdispersed,” that is, it exhibits a higher variance than would be theoretically expected for a randomly distributed “count” variable. As a result, the statistical distribution of eDNA samples can often best be modeled using non-Gaussian distributions, for example, the negative binomial distribution (Chambert et al., 2018; Furlan et al., 2016) or Weibull distribution (Barnes et al., 2021). The sampling effort necessary to accurately quantify a variable depends largely on the variance of the underlying distribution. Under the negative binomial distribution variance is dependent on the value of the mean (μ) and a dispersion parameter (θ). Environmental factors that affect the variance observed in eDNA concentrations could therefore represent useful a priori indicators of the potential sampling effort necessary to accurately quantify mean eDNA concentrations in a lentic system. Habitat size, for example, may be associated with increased habitat heterogeneity (Connor & McCoy, 1979). Increased habitat heterogeneity could result in increasingly clustered/patchy distributions of organism density as lentic systems increase in size, manifesting in increased variability in eDNA concentrations across spatial replicates and ultimately necessitating increased sampling effort to accurately quantify mean lake eDNA concentrations.

Herein, we utilized a previously published dataset on Brook Charr eDNA collected from 28 lakes between 5 and 91 ha (Gaudet-Boulay et al., 2023) to characterize spatial variation in eDNA concentrations observed within each lake and estimate the resulting sampling effort needed to accurately quantify mean eDNA concentrations in each system. We first explored the extent to which the distribution of sample eDNA concentrations within a lake deviated from normality. We then explored how eDNA concentrations varied across study lakes, predicting that variance in eDNA concentrations among our study lakes will positively correlate with lake size and negatively correlate with mean lake eDNA concentration. We also evaluated the effect of mean lake temperature on variance among spatial replicates within a lake. We then conducted a power analysis to estimate the distribution of “sampling effort” required to accurately (i.e., within pre-defined confidence thresholds) quantify the mean eDNA concentration of our study lakes and evaluated whether sampling effort necessary to accurately estimate mean eDNA concentrations among our study lakes was dependent on environmental conditions. We also evaluated the extent to which spatial autocorrelation was present among spatial replicates. We hypothesize that samples from smaller lakes (with closer eDNA sampling stations) will tend to exhibit greater spatial correlation than samples collected from large lakes; as a result, we predict that we will observe a negative correlation between lake size and eDNA sample spatial autocorrelation.

We also explored how non-linear patterns of sample variability could potentially emerge as a result of a bimodal distribution of eDNA particle copy count and total eDNA concentration. We generated simulated bimodal “pools” of eDNA with various ratios of small-to-large eDNA particle aggregates. We then repeatedly randomly sampled eDNA particles from these simulated “pools” and quantified how variation among “replicate samples” changed with the total number of particles sampled, thus simulating a gradient of samples with low-to-high “mean concentrations” of total eDNA copies.

MATERIALS AND METHODS Study species, eDNA sample collection, and quantification

As mentioned above, data were obtained from Gaudet-Boulay et al. (2023); for full methodological details, please refer to this manuscript. Environmental DNA samples were collected from 15 lakes in 2019 and 15 lakes in 2020 that, with the exception of a single lake also inhabited by Artic Charr, are inhabited solely by native Brook Charr (Gaudet-Boulay et al., 2023). Twenty 2 L samples were collected from each lake between 0.5 and 2 m in depth and 5 to 10 m from the shoreline. Samples were collected at systematic even intervals around the perimeter of each lake; equal distancing between samples was calculated by dividing the perimeter by the total sample number (20). Samples were collected between June 10–19 in 2019, but due to the COVID pandemic sampling was delayed to July 6–16 in 2020. Lake temperature data was collected during sample collection using a YSI ProDSS, and surface area was estimated using Google Earth.

For full decontamination, filtering, extraction, and quantitative PCR (qPCR) protocols, see Gaudet-Boulay et al. (2023). Environmental DNA concentrations were estimated based on a standard curve derived from known copy numbers of synthetic DNA of the target sequence (100,000, 10,000, 1000, 100, and 10 copies per reaction). Data from two lakes in 2019 were excluded due to potential problems that resulted in the exclusion of a number of replicates from the dataset and insufficient spatial replication for the power analysis (see Gaudet-Boulay et al., 2023); this resulted in a final dataset composed of observations from 28 lakes.

Statistical analyses Quantifying the distribution of eDNA concentrations within and across lakes

We first quantified the distribution of eDNA concentrations (e.g., mean and variance) within a lake in order to evaluate the extent to which lake eDNA distributions exhibited deviations from normality using the Shapiro–Wilks test implemented by the shapiro.test function in R (Core R Team, 2019). The “skewness” of the data was also tested using the skewness function from the R package e1071 (Meyer et al., 2023). The test implemented with this function evaluates the degree of skewness in the data, with negative values indicating “left-skew” and positive values indicating “right-skew”; values higher than |0.5| are typically interpreted as exhibiting skewness, with values greater than |1| indicating a high degree of skew.

We then explored whether degree of “skewness,” year of sampling, log-transformed mean lake eDNA concentration, and mean lake temperature affected whether the distribution of eDNA in a lake exhibited non-normality using a generalized linear model, with “non-normality” as a binary variable (non-normal = 1, “normal” = 0). Due to a limited binary sample size (28 lakes) it was not possible to include multiple predictors in a single model; the significance of variables was therefore tested separately. The significance of each variable was tested using Likelihood Ratio Tests (LRTs) against an intercept-only model. Year was included as a categorical variable, and log-transformed mean eDNA concentration, degree of “skewness,” and mean temperature were included as continuous predictors.

We then examined two alternative (to Gaussian) distributions to model eDNA concentrations: the Poisson and negative binomial, using gene-copy numbers obtained from a synthetic gene of known concentration. Under the negative binomial distribution, sample variance (s²) is a function of both the mean (μ) and a dispersion parameter (θ):[Image Omitted. See PDF]

As the value of θ increases, the observed distribution approaches the Poisson; small values of θ indicate high levels of overdispersion (i.e., larger variance than expected under a Poisson distribution).

To quantify trends in the mean and variance of eDNA concentrations observed for each lake and evaluate whether the negative binomial distribution improved models relative to the Poisson, we fit lake-specific Poisson and negative binomial distributions using the glm (base R) and glm.nb function from the MASS package (Venables & Ripley, 2002) in R. The presence of overdispersion in the Poisson models was tested using the dispersiontest function in the package AER (Kleiber & Zeileis, 2008), and the relative fit of the negative binomial model for each lake was compared to a similar Poisson model using a LRT.

Quantifying trends in the variance of eDNA concentrations within a lake

We quantified how the variance in eDNA concentrations observed within a lake varied with environmental correlates and the mean concentration of eDNA observed within a lake. To quantify variance in observed eDNA concentrations within a lake, we calculated the Geometric Coefficient of Variation (GCoV) for eDNA across spatial replicates:[Image Omitted. See PDF]Where s_ln is the standard deviation observed among samples collected within a lake after a natural-log transformation (Koopmans et al., 1964). The GCoV is a more robust and accurate standardized estimate of variation for log-normal/skewed data than the standard CoV.

To evaluate how variability in eDNA concentrations correlated with environmental variables and mean eDNA concentrations, we ran a linear regression with eDNA GCoV values for each lake as the dependent variable. However, due to unavoidable circumstances associated with delayed timing of sample collection during the COVID-19 pandemic (samples in 2020 were collected in July, but in June in 2019), several predictor variables exhibited strong collinearity (see Table S1 for VIFs). Log-transformed mean lake eDNA concentrations were significantly higher in 2019 relative to 2020 (detailed in Gaudet-Boulay et al., 2023). However, eDNA concentrations were strongly negatively correlated with lake temperature (r = 0.84) and temperatures were significantly warmer in 2020 vs. 2019; all lakes sampled in 2020 were at least 2°C warmer than the hottest lake sampled in 2019. Non-linear patterns between GCoV and predictor variables might emerge due to this experimental design, with a reduced ability to account for causal relationships among variables. To explore this issue, we ran two sets of models: (i) one set of models with a year-by-eDNA concentration term to investigate whether the linear relationship between eDNA concentration and predictor variables differed by sampling year, and (ii) another set of models with a quadratic term associated with mean eDNA concentration to explore whether, across years, variables exhibited a non-linear relationship with eDNA sample GCoV. Year (as a categorical predictor variable) exhibited a moderate level of multicollinearity with temperature (VIF >8). For the first set of models, we excluded temperature as a variable, and instead examined a year-by-eDNA concentration interaction term. For the second set of models, Year was dropped as a predictor variable and we instead included a continuous term for Temperature (including an additional quadratic term) and mean eDNA concentration. All continuous variables were mean-centered and scaled prior to analysis. For the second set of models with a temperature term, VIFs for temperature and log-transformed mean eDNA concentrations were problematic, but within tolerance levels (3 < VIFs <5, see Table S1). Collinearity inflates estimates of the standard error for variable effect sizes; at worst, significance tests for collinear variables will be overly conservative and prone to type-II error rates (Zuur et al., 2010). Both temperature, log-transformed mean eDNA concentration, and associated quadratic terms were therefore included in the analysis. Backward model selection using F-tests were employed to evaluate the significance of predictor terms, testing higher-order polynomial terms first; if a polynomial term was found to be significant, associated lower-order linear terms were automatically retained.

We also evaluated whether there was a significant relationship between the GCoV and the degree of skewness observed in the data using linear regression. We included the degree of “skewness” exhibited by each lake as a dependent variable and the GCoV as an independent variable; the significance of the GCoV term was evaluated using an F-test.

Power analysis – Determining optimal sampling sizes to quantify mean lake eDNA concentrations

To determine the number of samples necessary to estimate the mean concentration of eDNA within a pre-defined level of error, we conducted a power analysis on a set of 28 negative binomial distributions with “true” μ and θ values equal to the values fitted on our 28 study lakes in the previous analysis; as the largest published species-specific dataset on lentic eDNA systems, the μ and θ values estimated from the previous analysis represent our best estimate of the distribution of these values in naturally occurring waterbodies for Brook Charr eDNA. We determined sampling effort necessary to accurately estimate mean eDNA concentrations in these lakes by bootstrap sampling from simulated data. For each lake, we simulated 10,000 datapoints from a negative binomial distribution parameterized from the fitted μ and θ values from the previous analysis. We then estimated the sampling effort (“minimum n”) required to obtain a mean lake eDNA concentration estimate within the arbitrary thresholds of ±10% and ±20% of the “true” mean at an α-level of 0.1. For each simulated dataset of eDNA values, we bootstrapped 10,000 sample means for each value of sampling effort (n) ranging from 2 to 200. We then calculated the total percentage of the 10,000 sample means for each value of n that fell within ±10% and ±20% of the “true” mean for each lake and determined the “minimum n” for each lake in which 90% of simulated means fell within these thresholds (see Figure 1 for example).

FIGURE 1. Power analysis conducted on an example lake (Lake Boisvert). A total of 10,000 datapoints were simulated from a negative binomial distribution with μ and θ equal to the fitted values for Lake Boisvert presented as an example. We calculated 10,000 means from bootstrapped samples ranging in sample size (“n”) from each value ranging from 2 to 200 and determined the proportion of them that fell within ±20% of the mean of the original negative binomial distribution. Shown here are the distribution of mean values for n = 5, 10, 12, 14, 16, and 20. We would need to collect approximately 14 samples from a lake with a “true” mean and variance equivalent to the estimated mean and variance for Boisvert to be confident that 90% of our mean eDNA estimates fall within ±20% of the “true” mean. Plots were made using the R package WVPlots (Mount & Zumel 2023).

Identifying predictors of sampling effort required to estimate mean eDNA concentration

We then evaluated whether some environmental variables could be used to a priori estimate the sampling effort needed to accurately quantify mean eDNA concentrations within a lake. We used the “minimum n” to quantify eDNA concentrations within ±20% of the mean for each lake with an alpha-level error of 0.1, from the “power analysis” described above, as the dependent variable in our analysis. Lake surface area, lake perimeter divided by surface area (to evaluate the effect of lake “shape”), temperature, and log-transformed mean eDNA concentration were included as continuous predictors, and quadratic terms for temperature and log-transformed mean eDNA concentration were also included. All continuous variables were mean-centered and scaled prior to analyses. Unsurprisingly, “Year” was strongly collinear with temperature in this analysis, as well (VIF >9, see Table S2). As in previous analyses, we therefore ran two separate sets of models: (i) with a categorial “Year” predictor term and a Year-by-eDNA concentration term, and (ii) with temperature as a continuous predictor, including a quadratic term. For the second set of models, VIFs for temperature and log-transformed mean eDNA concentrations were again problematic, but within tolerance levels (3 < VIFs <5, Table S2).

“Minimum n” represents a discrete count variable and can thus be modeled using a Poisson or negative binomial distribution. For both candidate sets of models, we first evaluated whether a Poisson model exhibited significant overdispersion (see above), and the tested whether a negative binomial distribution improved the fit of the “full” model relative to the Poisson using a LRT. Backward model selection was then used to evaluate the significance of predictors using LRTs; higher-order polynomial terms were tested first and, if significant, associated lower-order terms were also retained.

Spatial autocorrelation in eDNA concentrations

The concentration of eDNA observed from environmental samples may not be spatially independent, either due to the horizontal transport of eDNA or underlying similarity in local habitat affecting patterns of organism abundance (and thus subsequent eDNA production). Environmental DNA samples were collected at relatively standardized intervals around the littoral zone of each lake (detailed in Gaudet-Boulay et al., 2023). As a result, samples in close proximity to one another may exhibit spatial autocorrelation; that is, closer eDNA samples may have more similar eDNA concentrations relative to samples farther away. We evaluated the degree of spatial autocorrelation observed among samples within a lake using Moran's I index values estimated from sample-station GPS coordinates using the ape package (Paradis & Schliep, 2019) in R. Z-scores can be calculated based on the Moran's I-value observed from the data minus the I-value expected under the null hypothesis, with significant positive values indicating positive spatial correlation (i.e., close samples exhibit similar values). To evaluate if lake size affected the degree of positive spatial autocorrelation observed, we ran a linear regression with z-scores for each lake included as the dependent variable. Lake surface area and temperature were included as continuous independent variables; “year” was not included due to its collinearity with temperature, as described above. The significance of model terms was evaluated using F-tests. The number of lakes that exhibited positive spatial autocorrelation was also compared to the proportion that exhibited positive spatial autocorrelation in a previous study examining Lake Trout (Salvelinus namaycush) eDNA concentrations (Lacoursière-Roussel, Côté, et al., 2016). Differences in study proportions were evaluated using a binomial generalized linear model with “study” as a fixed categorical variable; the significance of the “study” term was evaluated using a LRT.

Simulating “pools” of eDNA to evaluate the effect of hypothetical eDNA “state” on sample variance

Preliminary data analysis indicated non-linear patterns of sample variability with mean lake eDNA concentration across years. The capture of large aggregate particles containing high copy numbers of eDNA (i.e., “clumps” of eDNA) may be rare when total eDNA concentrations are low and, as a result, inter-replicate variability in eDNA concentrations may be relatively low. Inter-replicate variability in concentrations may similarly be low at high eDNA concentrations for opposite reasons: if multiple “aggregates” form a consistent component of total eDNA captured in a sample, any single large “aggregate” will have less of a disproportionate impact on the total eDNA concentration of a given sample. The stochastic capture of large “aggregates” may thus occur predominantly at intermediate eDNA concentrations; variability in eDNA concentrations, under some particle “size” distributions, may therefore be highest at intermediate pseudo-steady-state concentrations.

To explore whether non-linear relationships with sample variance could emerge as a function of bimodal eDNA particle copy count distributions combined with mean pseudo-steady-state eDNA concentrations we simulated a number of random pools of eDNA particles with varying distributions. We simulated eDNA particles from bimodal distributions based on varied small and large particle characteristics. First, we generated three pools of small, low copy number eDNA particles from a truncated Poisson distribution with means varying from 1, 3, and 5 copies of eDNA per particle. We then simulated a pool of large “aggregate” particles from a normal distribution with mean copy number from 1000, 3000, and 5000 eDNA copies per particle (maintaining copy number SD at 25% of the value of the mean) and varied the ratio of large-to-small particles from 1:500, 1:1000, 1:2000, and 1:4000.

As in our empirical study, we then drew 20 random samples from these “total” pools of particles to simulate collecting 20 spatial replicates. To simulate drawing samples from a gradient of environments with low-to-high mean eDNA concentrations, we varied the mean number of total particles collected per “sample” during each simulation iteration from 10 (simulating a low mean particle steady-state concentration) to 9000 (simulating a high mean particle steady-state concentration) and repeated this process 5000 times. To also simulate sampling from a “natural” environment, we also assumed that the concentration of particles among “sample sites” (i.e., for each of the 20 simulated samples per “site”) within an “environment” was heterogeneously distributed around the mean concentration of that “environment”; the number of particles collected per single “sample” (of 20) were thus drawn randomly from a Poisson distribution equal to the “mean” of the simulated environment. We then calculated the GCoV of samples for each mean “concentration” level and plotted mean sample GCoV against eDNA “concentration” to visualize whether patterns of sample variance were similar to what we observed in our empirical dataset. The smooth. spline function was used to reduce jitter in the final curves generated from the simulation.

RESULTS Distribution of lake eDNA concentrations, non-normality, and comparison of negative binomial vs. poisson models

Observed mean lake eDNA concentrations ranged from 184.11 to 137,958 copies/L (see Gaudet-Boulay et al., 2023). The distribution of eDNA in 17 of 28 lakes exhibited statistically significant (i.e., p ≤ 0.05) deviations from normality. All lakes exhibited positive (“right”) skew in the distribution of eDNA concentration data, with a small number of spatial replicates exhibiting eDNA concentrations substantially higher than most other replicates within a lake (Figure 2a,b). Notably, 27 out of 28 lakes exhibited skew values greater than 0.5. The extent to which a lake exhibited skewness was a strong and significant predictor of deviation from normality (p < 0.001). The log-transformed mean concentration of eDNA in a lake did not have a significant effect on the likelihood of a lake exhibiting non-normal eDNA concentrations (p = 0.102). However, lakes sampled in 2020 exhibited significantly higher rates of deviation from normality relative to lakes in 2019 (12/15 vs. 5/13, respectively, p = 0.023). Temperature was also positively correlated with deviation from normality in the distribution of eDNA within a lake (p = 0.019), although we notably cannot disentangle this effect from sampling year. The distribution of eDNA in 26 of the 28 lakes exhibited significant overdispersion relative to the Poisson distribution (p = 0.058 and 0.125 for the two lakes that did not) (see Table S1). The negative binomial model fitted to each lake exhibited, in all cases, a significantly better fit relative to a Poisson distribution (p < 0.001 for all lakes) (Table S1).

FIGURE 2. Distribution of eDNA values observed in each study lake. (a) depicts lakes sampled in 2019, and (b) depicts lakes sampled in 2020. See Figure 1 in Gaudet-Boulay et al. (2023) for a regional map of the study lakes.

Quantifying trends in the variance of eDNA concentrations within a lake

For the first set of candidate models with a Year and year-by-eDNA concentration term, the final model contained a year-by-eDNA concentration interaction (p < 0.001, Table 1, Figure 3); lower-order variables in this interaction were thus automatically retained. For the second set of models with a temperature term included, the final model contained only log-transformed mean eDNA concentration and its associated quadratic term, which was highly significant (Table 1, Figure 3). Lake eDNA GCoV was low at low concentrations, increased with mean eDNA concentration until it plateaued at intermediate concentrations, and declined at high mean concentrations. We found no evidence that lake surface area or lake temperature affected variability among the concentrations of eDNA observed for samples collected within each lake (Table 1), although this must be tempered by the fact that significance tests for the temperature term may have been over conservative due to moderate collinearity with log-transformed mean eDNA concentration. Inter-replicate GCoV was positively correlated with the degree of “skewness” present in eDNA concentrations from a lake (F_1,26 = 6.55, p = 0.025, Figure S1).

TABLE 1 Results of backward model selection (using F-tests) to assess significance of predictors of the geometric coefficient of variation for within-lake eDNA samples.

Model set 1				Model set 2
Term	F	df	p	Term	F	df	p
Year:Log Mean eDNA (copies/L)	16.37	1	<0.001	Temperature2	0.02	1.22	0.900
Year a	NA	NA	NA	Log Mean eDNA (copies/L) 2	9.21	1.23	0.006
Log Mean eDNA (copies/L) a	NA	NA	NA	Surface Area	0.01	1.23	0.921
Surface Area	0.48	1	0.495	Temperature	0.33	1.24	0.572
				Log Mean eDNA (copies/L) a	NA	NA	NA

Note: Statistically significant terms presented in bold.

^aLower-order term retained due to significance of associated higher-order term.

FIGURE 3. (a) Linear relationship between the geometric coefficient of variation (GCoV) among lake spatial replicates and the log-transformed mean lake eDNA concentration (copies/L) between years (2019 = red, 2020 = blue); (b) quadratic relationship between the GCoV among lake spatial replicates and the log-transformed mean lake eDNA concentration (copies/L).

Identifying predictors of sampling effort required to estimate mean eDNA concentration

The minimum sample sizes necessary to “accurately” quantify the mean eDNA concentration of the study lakes, such that 90% of sample means fall within ±10% and ±20% of the “true” mean, ranged from 5-173 (median = 50.5) and 2–44 (median = 12.5), respectively (Figure 4). Notably, however, these ranges were characterized by a small number of lakes that required large sample sizes to accurately estimate – only 5 of 28 of lakes required more than 25 samples, for example, to accurately estimate mean lake eDNA concentration within 20% of the “true” mean value, and 21 of 28 required less than 20.

FIGURE 4. Frequency distribution of sampling effort (“minimum n”) necessary to quantify the mean eDNA concentration of a lake within ±10% (a) and 20% (b) of the “true” mean. Sampling effort estimates were based on the negative binomial distribution, with μ and θ equal to the estimated values for each study lake.

Environmental predictors of sampling effort

For model sets 1 and 2, models with a Poisson distribution exhibited significant overdispersion (p = 0.005 and 0.004, respectively), and models with a negative binomial distribution provided a significantly better fit relative to the Poisson for the “full” models in each candidate set (LRT, p < 0.001 for model set 1 and 0.007 for model set 2).

For the first model set with a categorical “Year” predictor variable (and with “temperature” excluded), the final model contained a year-by-eDNA concentration interaction (p < 0.001, Table 2, Figure 5); lower-order variables in this interaction were thus automatically retained. For the second candidate model set, log-transformed mean eDNA concentration and its associated quadratic term had a strong and significant effect on sampling effort (Table 2). Similar to our analysis of inter-replicate GCoV, we found that the “minimum n” to accurately quantify mean eDNA concentrations was lowest at low and high eDNA concentrations, with required sampling effort peaking at intermediate concentrations (Figure 5).

TABLE 2 Results of backward model selection (using likelihood ratio tests) to assess significance of predictors of sampling effort needed to “accurately” quantify lake mean eDNA concentration (i.e., produce 90% of bootstrapped sample means within ±20% of the “true” mean of the distribution).

Model set 1				Model set 2
Term	χ 2	df	p	Term	χ 2	DF	p
Year:Log mean eDNA (copies/L)	15.19	1	<0.001	Temperature2	0.01	1	0.915
Year a	NA	NA	NA	Log mean eDNA (copies/L) 2	10.78	1	0.001
Log mean eDNA (copies/L) a	NA	NA	NA	Surface Area	0.41	1	0.523
Surface area	1.06	1	0.304	Temperature	1.69	1	0.193
Perimeter/Surface area	1.80	1	0.180	Log Mean eDNA (copies/L) a	NA	NA	NA
				Perimeter/Surface area	1.89	1	0.169

Note: Statistically significant terms presented in bold.

^aLower-order term retained due to significance of associated higher-order term.

FIGURE 5. Sampling effort (“minimum n”) necessary to quantify the mean eDNA concentration of a lake within ±20% of the “true” mean of a lake assuming (a) linear relationship among spatial replicates and the log-transformed mean lake eDNA concentration (copies/L) between years (2019 = red, 2020 = blue); (b) quadratic relationship among lake spatial replicates and the log-transformed mean lake eDNA concentration (copies/L).

We did not find evidence to support our hypothesis that the sampling effort necessary to quantify mean lake eDNA concentrations will increase with lake size or is affected by lake shape, as lake surface area and lake shape (i.e., perimeter/surface area) had no effect on “minimum n” for either set of models (Table 2). In the second set of models, we also found no statistically significant effect of temperature on sampling effort, although this result again must be interpreted cautiously due to its moderate collinearity with log-transformed eDNA concentration.

Spatial autocorrelation in eDNA concentrations

Only four of 28 lakes exhibited significant positive spatial autocorrelation. We found no evidence to support our hypothesis that spatial autocorrelation will be positively correlated with lake surface area (F_1,25 = 0.002, p = 0.966). Lake temperature (F_1,25 = 1.787, p = 0.193) also had no impact on spatial autocorrelation. The proportion of lakes exhibiting spatial autocorrelation was higher, but not statistically different, than a previous Lake Trout eDNA study (Lacoursière-Roussel, Côté, et al., 2016) in which one of 12 lakes exhibited significant spatial autocorrelation across spatial replicates (χ²₁ = 0.02, p = 0.590).

Evaluating the effect of eDNA particle copy count on sample GCoV from simulated eDNA pools

Our simulation studies demonstrated the emergence of a non-linear relationship between eDNA particle concentration and inter-replicate GCoV as a natural consequence of a bimodal distribution of eDNA particle copy count (Figure S2a–c). Inter-replicate GCoV typically was low at low and high concentrations and peaked at intermediate concentration levels. The shape and magnitude of the curve varied depending on the assumed properties of the eDNA particles. Higher ratios of large-to-small particles tended to increase sample GCoV, as did increasing the mean copy number of large particles. Conversely, increasing the mean copy number of “small” particles tended to decrease sample GCoV.

DISCUSSION Quantifying sampling effort necessary to estimate mean lake eDNA concentrations

Our study provides evidence that the distribution of eDNA concentrations in natural ecosystems can often deviate from normality. Almost two thirds of the lakes observed in our dataset exhibited statistically significant deviations from normality, and all lakes exhibited a “right-skew” in eDNA concentrations across spatial replicates. The degree of skewness in eDNA concentrations was also strongly correlated with statistically significant deviations from normality, indicating that the right-skew of the observed eDNA data was driving the extent to which the data exhibited non-normality. Alternative distributions were therefore necessary to model the mean and variance of eDNA concentrations among lakes. Poisson models fitted to each lake exhibited substantial overdispersion in the vast majority of cases; the negative binomial model, which has an additional term to model overdispersion, provided a much better fit due to the heterogenous variance in eDNA concentrations we observed. We thus add to a growing body of work (Chambert et al., 2018; Furlan et al., 2016) demonstrating that the negative binomial distribution may be particularly effective for modeling eDNA concentrations in natural ecosystems.

Our findings also have important implications for estimating accurate “lake-wide” mean eDNA concentrations; in seven of 28 of our study lakes, sampling effort >20 spatial replicates would be required to obtain estimates within ±20% of the “true” mean of theoretical negative binomial distributions with means and variances corresponding to our study lake estimates. The extent to which the distribution of eDNA within and across lakes observed in this study and species is generalizable to other species is unknown. Consequently, further studies on other species will be required to investigate how the spatial distribution of eDNA varies among taxa. However, the level of sampling observed in previous studies in lentic systems may not be adequate to accurately estimate mean lake eDNA concentrations for common/abundant species with relatively high pseudo-steady-state concentrations of eDNA. While the level of sampling performed in many previous studies (typically 8–10 samples, e.g., Lacoursière-Roussel, Côté, et al., 2016; Klobucar et al., 2017; Spear et al., 2021; Yates, Glaser, et al., 2021) may be adequate to characterize eDNA concentrations in seven of our 28 study lakes (“minimum n” <10), many lakes in our study required increased sampling effort to obtain similarly accurate estimates (median “minimum n” = 12.5), including the seven lakes that required more than 20 samples. This has important implications for a number of applications, but most importantly when attempting to study the correlation between organism abundance and mean lake eDNA concentrations.

The steady-state concentration of eDNA is largely driven by population production rates, of which organism abundance can be a dominant driver (Doi et al., 2017; Lacoursière-Roussel, Rosabal, & Bernatchez, 2016; Levi et al., 2019; Pochardt et al., 2020; Yates et al., 2019). Inaccurate estimates of mean lake eDNA concentrations will introduce additional unexplained variation in eDNA/abundance correlations, weakening the strength of observed relationships. We found that estimates of mean eDNA concentrations in an environment may be most accurate at extreme ranges of organism abundance (e.g., low and high organism densities), with the potential to observe significant heteroskedasticity in mean eDNA estimates at moderate levels of organism abundance unless more intensive sampling is conducted. Based on our observations, we suggest that sampling effort to quantify the average eDNA concentration in a natural environment may need to be higher when expected organism abundances and/or eDNA concentrations are low-to-moderate. We recommend that future studies should collect a minimum of 20–25 samples per lake. Collecting 20 samples, for example, would have ensured “accurate” estimates in 21/28 of our study lakes, and collecting 25 samples would have provided “accurate” estimates in 23 of 28 lakes. Alternatively, pilot studies to quantify the distribution and variance of eDNA concentrations in prospective study lakes may be warranted to optimize sampling effort.

There may also be important limitations to the extent to which mean lake eDNA concentrations can be accurately quantified from our survey methodology/design. Moving from estimates within ±20% of the “true” mean to ±10% of the “true” mean required realistically unattainable levels of sampling, increasing the median number of collected samples necessary from 12.5 to 50.5 per lake; this is a level of eDNA sampling that is not practical for most applications. Future research could explore the effect of different sampling designs on the accuracy of estimates; alternative methods could, for example, include samples generated from multiple spatial replicates that are pooled and/or homogenized (Chambert et al., 2018). Filtering much larger volume samples may also help reduce stochasticity in the capture of eDNA “aggregates” and thus minimize inter-replicate variability and right-skew across lake replicates and reducing error in eDNA concentration estimates (Pilliod et al., 2013).

Our assessment of the sampling requirements for lakes with disproportionately high variance, however, is further complicated by the fact that we “only” collected 20 samples per lake—although we would like to note that this is one of the largest datasets on the distribution of single-species eDNA yet published for lentic ecosystems. Estimates of the mean and variance of eDNA concentrations in seven of our 28 lakes were based on a sampling effort that is less than what would have been required to “accurately” estimate the mean of the “fitted” theoretical negative binomial distribution. As a result, our mean/variance estimates for these lakes may be inaccurate. At the very least, we can infer that the variance in eDNA concentrations in these lakes is likely large and would necessitate significant sampling effort to accurately quantify. Nevertheless, 20 samples were sufficient for most lakes to reach our pre-defined accuracy threshold. Notably, this level of sampling was sufficient to generate a significant correlation between eDNA concentration and fisheries data (Gaudet-Boulay et al., 2023) and may be sufficient in many management contexts for similar species/systems. It also important to note, however, that despite the stochastic occurrence of samples with high variability, variance in traditional capture methods (e.g., index netting) can be comparable to variance observed among eDNA samples (Lacoursière-Roussel, Côté, et al., 2016). Yet, sampling effort that is required using traditional capture methods (e.g., gillnets) to obtain reliable estimate of fish abundance (density or biomass) in lentic systems is rarely discussed or investigated.We also found no evidence to support our hypothesis that sampling effort necessary to accurately estimate mean eDNA concentrations will increase with habitat size. While we lack detailed habitat data to conclusively conclude that habitat heterogeneity is positively correlated with lake size in our study systems, at a minimum we can conclude that any increased heterogeneity in larger lakes is unlikely to impact the distribution of eDNA at a level that was detectable. Our study lakes were largely small-to-medium-sized lentic systems and varied in size from 5 to 91 ha; sampling effort may not need to vary in surveys of lentic systems at this scale for brook trout eDNA given that other correlates (e.g., eDNA concentration) explain substantially more variation in “minimum n.” The extent to which our findings can be generalized to larger systems inhabited by brook trout and other species is unknown. However, there is likely a spatial scale at which sampling effort will need to increase with lake size. Future research on larger lentic systems will be needed to determine optimal sampling effort for lentic systems larger than ~100 ha.

While we do not have corresponding local estimates of spatial abundance, we suspect that littoral Brook Charr abundance is driving trends in the concentrations of eDNA we observed across years (see Gaudet-Boulay et al. (2023) for discussion). In brief, sampling was primarily focused on the littoral zone since Brook Charr favor littoral habitats, occurring at higher densities in littoral zones relative to pelagic habitat (Tiberti et al., 2017). However, the species also prefer cool temperatures, and thus may retreat to deeper water during periods of elevated thermal conditions (Cote et al., 2020). Samples from 2020 were collected when water temperatures were substantially warmer than in 2019; it is likely that Brook Charr had, at that point in time, moved out of the littoral zone to favor deeper, colder refugia (Gaudet-Boulay et al., 2023). We therefore suspect that changes in the local density of Brook Charr in littoral/surface waters were predominantly driving the observed concentrations of eDNA and corresponding GCoV/“minimum” sampling effort values. Local activity of fish, for example, could be driving the likelihood of encountering “large” copy count particles, that is, large copy count particles may be more likely to be encountered if fish were recently swimming in the sampled area. High occupancy of the littoral zone could result in an increased probability of encountering such particles; conversely, if fish are inhabiting deeper areas this potential source of inter-replicate variability may be minimized. Our work thus highlights the importance of understanding the ecology of focal species when designing quantitative eDNA surveys, given that patterns in the sampling effort required to quantify mean lake “eDNA” were likely affected by the effect temperature on habitat usage, behavior, and physiology of our study species.

Spatial autocorrelation among sampling sites

We found that eDNA concentrations among sampling stations within four of our study lakes exhibited significant spatial autocorrelation. Spatial autocorrelation can arise from both “endogenous” (i.e., from processes associated with eDNA dynamics) or “exogenous” (i.e., spatially structured environmental gradients) processes (Miller, 2012). We cannot determine whether the spatial autocorrelation we observed in our study lakes arose from horizontal eDNA transport between sites (i.e., endogenous) or from underlying spatial autocorrelation in habitat preferred by Brook Charr (i.e., exogenous); further study is warranted. However, we did not observe a relationship between the degree of spatial autocorrelation and lake size (and thus the distance between sampling stations). The horizontal transport of eDNA in lentic systems can be very limited (Ghosal et al., 2018; Goldberg et al., 2018). We therefore speculate that the spatial autocorrelation in eDNA concentrations observed in this study resulted from underlying habitat distributions that affected the distribution of Brook Charr and the subsequent eDNA they produced. Additionally, while we observed spatial autocorrelation in a higher proportion of lakes than a previous study (4 of 28 in our study, 1 of 12 in (Lacoursière-Roussel, Côté, et al., 2016), this difference was not statistically significant. Future studies in lentic systems can likely expect to encounter a small proportion of lakes that exhibit some spatial autocorrelation in the distribution of eDNA. Further research is warranted to investigate whether spatial autocorrelation in eDNA concentration arises due to endogenous or exogenous processes, and to determine the spatial scale at which it is likely to encounter autocorrelation in eDNA samples from lentic systems.

A potential mechanism related to “eDNA state” to account for non-linear patterns of sample variance with mean lake eDNA concentrations

Phenomena that follow “clustered” or “clumped” distributions (i.e., distributions in which counts are right-skewed) are commonly observed in nature and are well-studied in population ecology (Krebs, 2014). Clustered eDNA distributions could arise for a number of reasons. In lentic systems, horizontal eDNA diffusion/transport is limited to a scale of tens of meters (Ghosal et al., 2018; Goldberg et al., 2018); as a result, eDNA concentrations in lentic systems can be positively correlated with local organism abundance (Ghosal et al., 2018; Littlefair et al., 2020). Variable inter-replicate eDNA concentrations may therefore reflect classic underlying clustering in the spatial distribution of organisms as a result of schooling behavior, habitat heterogeneity, etc. Changes in patterns of Brook Charr behavior/habitat usage associated with temperature could thus be driving the opposing trends in inter-replicate variability observed across the study years; we are unable to exclude this possible explanatory mechanism.

However, “clustering” in the pattern of eDNA concentrations may also arise due to the state in which eDNA is collected. Several experiments studying eDNA particle size distributions in natural or semi-natural conditions have found that eDNA tends to exhibit higher concentrations on filters with pores ranging from 1 to 10 μm, suggesting that eDNA may be collected from individual cells and/or mitochondria (Turner et al., 2014; Wilcox et al., 2015). However, both studies observed high inter-sample variability, a surprising finding if eDNA state is dominated by freely associated individual cells and mitochondria. Wilcox et al. (2015) proposed that fish eDNA may often occur in large, but easily broken-down, particle aggregations that break apart during the filtering process; this hypothesis is consistent with the observations that most eDNA in semi-natural conditions is captured on filters with pore sizes between 1 and 10 μm but spatial replicates can exhibit high inter-sample variability. Furthermore, Turner et al. (2014) hypothesized that large eDNA “clumps” break apart relatively quickly as they move from their source-of-origin into single cells or small particles with low copy numbers or, alternatively, rapidly “settle” due to their larger particle size, a hypothesis that has been supported by subsequent empirical experiments in lotic systems (Wood et al., 2020, 2021). A similar mechanism could be operating in our study system, because if brook trout were inhabiting deeper areas of the lakes eDNA would have had to travel to littoral zones, providing time for the particles to “break down.” Particles containing eDNA may therefore persist in two general “states” in nature with different particle size distributions/properties: a “background” pool of eDNA composed of small particles (potentially single cells/mitochondria or small aggregates) and a pool of large “clumps” or aggregates that occur at much lower relative densities but with high copy numbers of eDNA. Future studies on sample concentration variability in natural ecosystems could employ sequential filtering to evaluate how inter-replicate variability changes with filter size, as in Barnes et al. (2021).

Our (simplistic) simulations show that bimodal distributions of eDNA particle copy count have potential implications for inter-replicate variability. We demonstrated that, given a bimodal distribution of particle copy count described above, variability among replicate samples would, under some conditions, be expected to exhibit a non-linear pattern with eDNA concentration similar to the patterns observed across years. These findings support the hypothesis proposed by Jo & Yamanaka (2022) that the concentration of eDNA in smaller particles may more closely reflect organism abundance, presumably due to decreased inter-sample variability associated with larger particles (at least at small-to-moderate concentrations where such particles may not be regularly captured). We do note, however, that filter pore sizes varied slightly between years due to COVID pandemic-related supply issues (1.2 μM in 2019, 0.7 μM in 2020), and it is possible (but in our view unlikely, given that both pore sizes are adequate to capture eukaryotic cells) that this could have impacted concentrations/variability across years.

In our empirical data, we observed that the GCoV of eDNA concentrations exhibited a strong and significant non-linear (quadratic) relationship with total eDNA concentration in our study lakes. While we found broadly equivalent support for models with a year-by-concentration interaction term, we also found a statistically significant and positive relationship between the degree of right-skew and the CGoV among spatial replicates within a lake, providing further evidence that skew among replicates was driving inter-replicate variability within a lake. We therefore suspect that variability (and skew) among replicates was likely driven by the pseudo-stochastic capture of large eDNA “aggregates” with high eDNA copy numbers.

It is notable that other studies have conversely observed that variability in inter-replicate eDNA concentrations can increase at higher levels of relative organism abundance (Chambert et al., 2018). However, it is possible that these studies may not have achieved sufficient densities at which large “aggregates” of eDNA were regularly captured; this may particularly be the case in lotic systems, where constant transport of eDNA could result in lower overall particle numbers per L of water, limiting the capture of large “aggregates” to stochastic occurrences. Turner et al. (2014), conversely, observed higher inter-sample variability in carp eDNA particle size distribution and concentrations in a lake environment relative to a high-density pond environment.

However, we cannot rule out the alternative possibility that the inter-replicate variability we observed emerged as a result of spatial clustering patterns in our study organisms that were related to/dependent upon relative abundance patterns. Local spatial abundance exhibits a strong correlation with local eDNA concentrations in lentic systems (Ghosal et al., 2018; Goldberg et al., 2018; Littlefair et al., 2020). It is possible that in lakes with low littoral abundance levels (and thus low observed eDNA concentrations), Brook Charr may exhibit random or dispersed spatial distributions. At moderate littoral abundances, the species may cluster around particularly productive/suitable habitat patches, and at high littoral abundance levels it may again exhibit random or dispersed spatial distributions. Such a pattern of habitat occupancy/use, in relation to overall levels of organism abundance, could also account for the pattern of variability in the concentrations of eDNA we observed across lakes. Similarly, high inter-sample variability could be related to recent defecation events that produce temporary but highly localized high eDNA concentrations. At low density, organisms may be sparsely distributed enough that the likelihood of sampling in the immediate vicinity of a stochastic defecation event from an individual fish is low. At some hypothetical level of intermediate density, it may be more likely to stochastically sample near such an event and at high enough densities such events may be common enough that eDNA from recent defecation events may form a consistent component of eDNA. Future research with spatially paired abundance estimates may be necessary to further evaluate why eDNA concentrations exhibited increased variability at intermediate concentrations.

CONCLUSIONS

We found evidence that the distribution of Brook Charr eDNA concentrations across spatial replicates in lentic systems was non-normal, right-skewed, and could best be modeled using negative binomial distributions. We also found evidence that variance in eDNA concentrations among spatial replicates was highest at intermediate mean lake eDNA concentrations. This observation is parsimonious with the hypothesis that eDNA particles are captured in multiple “size” states, although we were not able to rule out the potential contrasting effects of temperature on eDNA concentrations observed between our two study years. Nevertheless, sampling effort to estimate a “mean lake eDNA concentration” may need to be higher if intermediate eDNA concentrations are anticipated. We also generally recommend expanding the typical number of spatial replicates collected during lentic surveys to 20–25 replicates; we note that most previous studies have typically collected fewer replicates, including some of the authors' previous work. We also recommend further studies on the effect of increasing the volume of water that is filtered and/or investigating the effectiveness of different survey designs (e.g., systematic vs. random sampling designs, homogenizing spatial replicates prior to filtering).

ACKNOWLEDGMENTS

We would like to thank the field technicians of Ministère de l'Environnement, de la Lutte contre les changements climatiques, d la Faune et des Parcs du Québec (MELCCFP) for their help in the eDNA and water sampling and all the lab members that helped in every stage of this project: A. Perreault-Payette, C. Babin, and everyone else. This project was supported by the MELCCFP, the Société des Établissements de Plein Air du Québec (SÉPAQ), and a Strategic Partnership Grants for Projects from the Natural Sciences and Engineering Research Council of Canada (NSERC) to LB. This study is part of the research program of Ressources Aquatiques Québec (RAQ).

CONFLICT OF INTEREST STATEMENT

None declared.

DATA AVAILABILITY STATEMENT

No new data were generated for this manuscript, analyses in this manuscript utilized data from Gaudet-Boulay et al. (2023). A curated dataset used for these analyses can be found at: https://doi.org/10.5061/dryad.5x69p8d8x.