Meta‐analysis supports further refinement of eDNA for monitoring aquatic species‐specific abundance in nature

The use of eDNA to detect the presence/absence of rare or invasive species is well documented and its use in biodiversity monitoring is expanding. Preliminary labora‐ tory research has also shown a positive correlation between the concentration of species‐specific eDNA particles and the density/biomass of a species in a given envi‐ ronment. However, the extent to which these results can be extended to natural environments has yet to be formally quantified. We collated data from experiments that examined the correlation between eDNA and two metrics of abundance (bio‐ mass and density) and, using mixed‐effects meta‐analysis, quantified the strength of that correlation across laboratory and natural environments. We found that eDNA particle concentration was more strongly correlated with abundance in laboratory environments compared to natural environments, accounting for approximately 82% and 57% of the observed variation in abundance, respectively. We found some evi ‐ dence of potential publication bias that may have impacted the estimation of the correlation in natural environments; after smaller studies were removed from the dataset, eDNA particle concentration accounted for approximately 50% of the ob‐ served variation in abundance in natural environments. No effect of abundance met‐ ric was found on the strength of correlation between eDNA particle concentration and abundance. Despite a weaker general correlation in natural environments, eDNA concentration often still explained substantial variation in abundance. eDNA re‐ search is still an emergent field of study; with only moderate improvements in tech‐ nology or technique, it could represent a powerful new tool for quantifying abundance.

essentially impossible to obtain (Luikart et al., 2010;Ovenden et al., 2016;Yates, Bernos, & Fraser, 2017). Novel technological and methodological developments have generated interest in using genetic techniques as a cost-effective alternative to indirectly estimate or infer N c in natural populations. Some recent attempts to indirectly monitor N c from genetic data used effective population size (N e ) (or effective number of breeders, N b ) (e.g., Whiteley et al., 2015, Bernos & Fraser, 2016, Perrier, April, Cote, Bernatchez, & Dionne, 2016), but variance in N b /N c ratios across populations can limit the utility of N b to infer N c , at least in species with highly variable and/or plastic life histories (Yates et al., 2017). Furthermore, even within a population, changes in N b often do not track changes in N c (Bernos & Fraser, 2016). Alternatively, environmental DNA (eDNA) is emerging as a powerful tool with which indirect estimates of abundance might be obtainable.
Environmental DNA refers to particles of DNA that are present in a sampled environmental substrate and can be reliably found in many environments (soil, air, and water). The utility of eDNA for detecting the presence or absence of rare or invasive species and for monitoring biodiversity has been well established (Goldberg, Strickler, & Pilliod, 2015). Preliminary research has encouragingly also demonstrated a positive correlation between the quantity of species-specific DNA particles in an environment and the density or biomass for many macroorganisms (Lacoursière-Roussel, Côté, Leclerc, Bernatchez, & Cadotte, 2016;Takahara, Minamoto, Yamanaka, Doi, & Kawabata, 2012;Wilcox et al., 2016).
Unsurprisingly, many of the studies that have demonstrated strong correlations between eDNA quantity and species-specific density or biomass have been done in controlled laboratory conditions (e.g., Takahara et al., 2012, Doi et al., 2015, Lacoursière-Roussel, Rosabal, & Bernatchez, 2016 where it is possible to control for the myriad of environmental factors that affect the production or degradation of eDNA and/or halt its transport from the environment. Yet, the extent to which these strong relationships can be generally extrapolated to natural environments remains questionable. Water volumes relative to organism abundance are much greater in natural ecosystems compared to most laboratory environments, where eDNA particles are likely to be more concentrated. Additionally, natural environments exhibit significant environmental stochasticity for factors that are known to affect eDNA production and degradation (temperature, pH, sunlight, microbial activity, etc.) (Lacoursière-Roussel, Rosabal, et al., 2016;Pilliod, Goldberg, Arkle, & Waits, 2014;Strickler, Fremier, & Goldberg, 2015). Animal behavior patterns may also fluctuate in a manner that affects eDNA production (activity related to diurnal patterns, seasonal spawning behavior, etc.) (Doi et al., 2017;Pilliod, Goldberg, Arkle, & Waits, 2013). Such stochasticity could introduce significant variation in the relationship between eDNA and abundance in natural systems, potentially limiting the utility of eDNA as an abundance quantifying tool in natural populations of conservation concern.
Studies in both artificial and natural systems often quantify species-specific "abundance" through either one (or both) of two metrics: (i) biomass or (ii) number of individuals (e.g. density). Which abundance metric might generally exhibit a stronger correlation with eDNA particle concentration is unknown. While per-gram eDNA shedding rates often differ among age-classes (Maruyama, Nakamura, Yamanaka, Kondoh, & Minamoto, 2014), the generally positive relationship between excretion rate and body size (Vanni & McIntyre, 2016) could result in biomass exhibiting a stronger predictive correlation with eDNA concentration (Maruyama et al., 2014;Takahara et al., 2012). However, studies in which the relationship was quantified for both density and biomass on the same species/system often indicate a weaker correlation between eDNA and biomass compared to eDNA and density (e.g., Lacoursière-Roussel, Côté, et al., 2016;Baldigo, Sporn, George, & Ball, 2017). Assuming a constant rate of eDNA production per unit biomass, a spatially large environment containing biomass that is distributed across a large number of individuals is more likely to have a more homogeneous distribution of eDNA concentration relative to an environment in which biomass is concentrated in a few individuals (Lacoursière-Roussel, Côté, et al., 2016). Resulting increases in sampling variability might reduce the direct correlation between biomass and eDNA concentrations across environments. Additionally, excretion is thought to be a primary source of eDNA in aquatic organisms (Klymus, Richter, Chapman, & Paukert, 2015). Yet, per-gram eDNA excretion rates in aquatic organisms often do not increase linearly with body mass; larger organisms, for example, often excrete proportionally less when controlling for body size (Maruyama et al., 2014). As a result, eDNA concentrations might exhibit stronger correlations across environments with the number of individuals, rather than their biomass. Overall, determining which metric (density or biomass) produces more consistent and reliable results could improve the focus of future eDNA studies attempting to quantify population abundance.
The application of eDNA technology to estimate abundance is currently in its infancy. Yet, enough studies have been conducted in both artificial and natural systems to contrast how correlations between eDNA particle concentration and abundance vary across these two environmental contexts. We collated correlation coefficients across studies and conducted a formal meta-analysis with two objectives: (i) to assess how the strength of the correlation between eDNA and species-specific abundance changed across artificial and natural aquatic environments; and (ii) to determine if the strength of eDNA/abundance correlations vary for different abundance metrics (biomass vs. density). We then converted mean correlation coefficients across environments into corresponding R 2 values to assess the capacity of eDNA to model biomass/density across environments and metrics.

| Survey of the literature
To find articles examining the correlation between density/biomass of a target species and the concentration of eDNA of that species in aquatic environments, keyword searches were conducted on the academic search engine ISI Web of Science TM . A complete keyword search for "eDNA + abundance" and "environmental DNA + abundance" was performed. Only studies examining macroscopic aquatic organisms from the kingdom Animalia were retained. Relevant references within studies were also searched for usable data. Studies conducted either in controlled artificial settings or in natural environments (e.g., streams, lakes) were retained.
Only studies that estimated absolute eDNA concentrations were included; such studies typically used quantitative PCR (qPCR) or digital droplet PCR (ddPCR) to estimate eDNA particle concentration. Notably, this criterion excluded metabarcoding experiments, which indirectly infer relative eDNA particle concentration based on read counts (e.g., Evans et al., 2016). For a study to be retained, raw eDNA concentration and biomass/density data must have been adequately reported in a figure or table either in the manuscript or supplementary files. Retained studies had to have at least four estimated eDNA/biomass/density datapoints, the minimum required to estimate variance under a Fisher z-transformation. To retain a study conducted in natural environments, standardized biomass/density estimates must also have been obtained through a method likely to generate high-quality abundance estimates (e.g., telemetry data, mark-recapture estimates, depletion estimates, standardized trap/ haul data, etc.).

| Statistical analysis
When meta-analyzing correlations between two continuous variables, Pearson correlation coefficients (Pearson, 1895) are an appropriate effect size statistic (Lipsey & Wilson, 2001). When possible, Pearson correlation coefficients were extracted directly from manuscripts; if unavailable, correlation coefficients were manually calculated between eDNA concentration data and associated biomass/ density estimates reported in the manuscript or supplementary file using R version 3.4.1 (R Development Core Team, 2017). In many studies each eDNA sample obtained from an experimental replicate (e.g., tank, field site, stream transect, etc.) was analyzed multiple times using qPCR techniques to generate a more robust estimate of concentration; to standardize how correlation coefficients were calculated across studies and to avoid pseudo-replication studies were only retained if multiple qPCR analyses were averaged to generate a single value for each individual replicate. If such data was available in the manuscript or supplementary files, replicate averages were manually calculated from extracted data. When eDNA and density/ biomass estimates were contained in figures, the package "digitize" (Poisot, 2011) in R version 3.4.1 (R Development Core Team, 2017) was used to extract the data. In many studies, logarithmic transformations of eDNA and/or biomass/density often improved linear fit.
Whenever possible, the data transformations presented in the original studies were retained except when inspection of the raw data indicated that additional transformations (typically natural logarithmic transformations [ln]) could further linearized data or when other data distributions were modeled (e.g., negative binomial, log-link, etc.).
When zeroes were present in datasets that were ln-transformed, "1" was added to all datapoints prior to transformation. Several papers estimated eDNA and density/biomass relationships across several experiments conducted over multiple environments or using different groups of organisms (e.g., different species, life history stages, etc.); separate correlation coefficients were calculated for each independent experimental grouping. In two studies, multiple filters were used to process water taken from the same replicates; correlation coefficients across filter types were averaged and the n corresponding to the final collated correlation coefficient set to the number of original replicates. Six studies also contained both density and biomass estimates associated with the same eDNA samples and generated from the same experimental grouping (12 total groupings); in such cases, separate correlation coefficients were calculated for biomass and density (see below for treatment of correlated sampling errors).
The escalc function from the R package "metafor" (Viechtbauer, 2010) was then used to calculate effect size and variance estimates using Fisher's z-transformation (Fisher, 1921).
The function rma.mv from the "metafor" package was then used to conduct an inverse-variance weighted mixed-effects metaanalysis. Two fixed categorical effects were tested: (i) environment ("Laboratory" vs. "Nature") and (ii) metric of correlation ("Biomass" vs. "Density"). Species largely covaried with study so only a study-level random effect was fitted to all models. An "experiment" random effect was also nested within study to account for studies that included multiple experiments conducted on independent subject groupings; this random effect was conditioned on "metric of correlation," assuming an unstructured variance/covariance matrix. Data from the six studies that estimated correlation coefficients for both biomass/eDNA concentration and density/eDNA concentration from the same group of study organisms violated assumptions of independence. Correlation coefficients were therefore calculated, whenever possible, between biomass and density for all experimental groupings within the six studies (excepting two experimental units); covariance between biomass/eDNA and density/eDNA correlation coefficients within experimental groupings was then calculated as in Olkin and Finn (1990) using the custom rmat function written by the author of the metafor package (Viechtbauer, 2018). It was then possible to calculate complete variance/covariance matrices for the experimental groupings within those six studies, which were then passed to the rma.mv function to account for the correlated sampling error. Cluster-robust variance estimation (clustering by experiment) was also used to account for non-independent sampling error implemented through the cluster function in "metafor." Fixed effects were tested using backwards stepwise selection conducted with likelihood ratio tests (LRTs), removing nonsignificant fixed-effects terms (p > 0.05). All random effects were included in models regardless of significance level. Final model coefficient estimates and confidence intervals were back-transformed to Pearson correlation coefficients and corresponding R 2 values were calculated. Degrees of freedom for t tests comparing moderator variable levels were conservatively set to the number of studies in the meta-analysis minus the number of estimated model parameters (including random effects).

Publication bias was assessed by conducting a modified
Egger's regression (Egger et al., 1998); standard errors of the effect size estimates were included as a covariate in a meta-regression, with a significant positive coefficient indicating a correlation between effect size estimate magnitude and study precision.
Visual examination of funnel plots (effect size estimate standard errors plotted against model residuals) were also used to assess potential publication bias, although caution should be used when interpreting funnel plots from mixed-effects meta-analyses due to the tendency for non-independent datapoints to cluster. If evidence of publication bias (i.e., positive residual asymmetry) was found, a sample size threshold was chosen below which all TA B L E 1 Published studies examining relationships between species-specific eDNA concentration and density/biomass

| Effect of environment and abundance metric on eDNA and density/biomass correlations
The final model included only environment as a fixed effect (Table 2). The environmental context of a study significantly affected the strength of the correlation coefficient between eDNA and abundance, with experiments conducted in laboratory or artificial pond settings exhibiting a much stronger correlation than experiments conducted in natural environments (t 15 = 2.707, p = 0.015, Figure 1). The correlation between eDNA particle concentration and abundance in laboratory experiments was extremely strong, explaining up to 81% of the variation observed when converted to a corresponding R 2 value. Although significantly weaker, eDNA particle concentration still explained substantial amounts of variation in abundance (57%). No significant effect of metric (i.e., density vs. biomass) was found (Table 2).
Observed correlation coefficients also exhibited significant residual heterogeneity (Q E,27 = 92.2, p < 0.001), indicating more residual variance remained than expected by chance despite the inclusion of random effects and significant moderator variables.
Some evidence for potential publication bias in some small studies was found; Egger's regression demonstrated a significant effect of standard error on mean effect size (χ 2 1 = 6.660, p = 0.010) and visual inspection of a funnel plot found that studies with small sample size (i.e., high variance) tended to exhibit positive residual asymmetry ( Figure 2).
Funnel plot inspection indicated that residual asymmetry was observed predominantly in experiments with six or fewer samples; these datapoints were removed and all models reanalyzed. A subsequent Egger's regression found that this adequately addressed publication bias issues (χ 2 1 = 0.020, p = 0.887). Re-analysis of the subsetted dataset exhibited similar trends, with the final model retaining a significant effect of experimental environment and no effect of metric (Table S1). Natural environments still exhibited a weaker correlation between eDNA concentration and abundance than laboratory environments (t 13 = 3.154, p = 0.008). However, natural environments in the subsetted dataset exhibited a slightly weaker correlation between eDNA particle concentration and abundance when compared to the full dataset (R 2 = 0.50 vs. 0.57, respectively) ( Figure S1).

| D ISCUSS I ON
To evaluate the reliability of eDNA as a tool to measure species-specific abundance or biomass, it is first necessary to quantify how consistent the relationship is across studies. By examining the strength TA B L E 2 Results of model selection for the meta-analysis, with statistically significant (retained) terms presented in bold. Environment refers to natural versus laboratory environments, and metric refers to the metric used to quantify abundance (biomass vs. density). Model testing was conducted using likelihood ratio tests  of the correlation coefficients between eDNA and two metrics of abundance, it was possible to broadly contrast its predictive utility across artificial and natural environments. We found that speciesspecific eDNA particle concentration was positively correlated with abundance, but that this correlation was substantially and significantly stronger in artificial experimental settings relative to natural environments.
It is unsurprising that eDNA better predicts abundance in controlled experimental settings. Multiple factors can affect the quantity of eDNA detected, which is a function of production, degradation, and transport of eDNA particles .
Field and laboratory studies have demonstrated that a variety of environmental variables (diffusion rate, temperature, pH, etc.) and life history characteristics (e.g., size/biomass) can affect the amount of eDNA particles available for detection (Barnes et al., 2014;Jane et al., 2015;Klymus et al., 2015;Lacoursière-Roussel, Rosabal, et al., 2016;Strickler et al., 2015). In an experimental laboratory setting, it is possible to control for such variation; replicates are typically held to standardized temperatures, lighting, flow, pH, and so on. Experimental organisms also are often derived from the same age-class, and experimental duration is similarly held constant. More importantly, laboratory experiments can precisely control levels of abundance (e.g., biomass or density of study organisms), whereas studies in natural environments must estimate, often with substantial error, such metrics. Expectedly, under controlled laboratory conditions variation in available eDNA results almost exclusively from variation in the density or biomass of eDNA producing organisms.
In natural environmental conditions, such control is not possible. Temperature, for example, significantly affects eDNA degradation (Pilliod et al., 2014) and likely production (Lacoursière-Roussel, Rosabal, et al., 2016), and yet can fluctuate both across aquatic environments and diurnally/seasonally within environments. Changes in eDNA concentration can also lag behind changes in organism occupancy  and can also fluctuate with flow, which can significantly affect the "settling" probability and/or transport of DNA particles (Jane et al., 2015). Additionally, different size/age-classes are likely to produce eDNA at different rates on a per-gram basis (Maruyama et al., 2014); variations in size and ageclass distributions across and within environments could affect the relationship between eDNA and biomass and/or density, especially when some age-class densities are not quantified during abundance surveys. As an example, streams with high densities of early juveniles relative to adults may exhibit greater average concentrations of eDNA particles relative to streams with lower juvenile densities, obscuring the relationship between adult abundance and eDNA concentration. Collectively, such factors may introduce significant variability across species, environments, and even between samples.
In comparing the relationship of eDNA concentration with metrics of species-specific abundance, we found no evidence that biomass was more strongly correlated with eDNA particle concentration compared with density. Ultimately, it may not matter much what metric abundance metric researchers choose to quantify because eDNA may be an equally useful abundance monitoring tool in either context, although this conclusion must be tempered with an acknowledgement that it is based on a contrast of 28 versus 16 correlation coefficients. The accumulation of future research could demonstrate a small but biologically meaningful general difference between the two metrics. However, until more data accumulates researchers should focus efforts on quantifying the most cost-efficient or least effort-intensive metric unless further research demonstrates a significant difference or experimental design (e.g., study species, question, etc.) requires the collection of a specific metric.
Finally, significant heterogeneity across studies remained despite accounting for random effects and multiple significant or marginally significant moderator variables. Numerous potential factors may have affected the heterogeneity observed among effect sizes.
eDNA capture and extraction technology and/or methodology can significantly affect eDNA quantification (Hinlo, Gleeson, Lintermans, & Furlan, 2017;Lacoursière-Roussel, Rosabal, et al., 2016), but the use of eDNA to predict abundance is a field in its infancy and tech- niques have yet to be both optimized and standardized. Studies often differed in water filtration techniques (e.g., automatic pump vs. hand filter vs. direct extraction on water samples themselves), eDNA storage, filter medium and/or pore size, extraction technique (salt-based vs. kit-based [e.g., Quiagen TM ]), PCR technique (e.g., qPCR vs. ddPCR), and the materials used in the generation of standardized curves during qPCR (copy number from synthetic plasmids vs. total eDNA from tissue samples). Furthermore, most eDNA studies rely on mitochondria DNA sequences (e.g., Wilcox et al., 2013Wilcox et al., ,2016; the number of copies of mitochondrial DNA sequences per cell can vary significantly across species and tissues, with many tissues possessing a large number of mitochondria (Robin & Wong, 1988). Due to the extreme sensitivity of the technology used to estimate eDNA, the presence or absence of even a small number of mitochondria rich cells could dramatically change the quantity of eDNA detected (Pilliod et al., 2013). Although contamination is a potential issue (Pilliod et al., 2013), standardizing qPCR curves from plasmids can allow researchers to directly estimate copy count number present in eDNA, bypassing this potential issue. Alternatively, it may be F I G U R E 2 Funnel plot indicating residual asymmetry in small studies with high variance and low sample size beneficial for future studies to explore the use of nuclear DNA sequences, as these would be less sensitive to variation in cell type present in eDNA samples.

| Future directions
Although laboratory studies clearly demonstrate the potential utility of eDNA as a tool to estimate species-specific abundance, attempts to duplicate these efforts in natural environments have met with limited success. Refining techniques used to estimate abundance from eDNA in field settings by incorporating methodological advances is required to improve the application of this technology.
Moreover, the development of standardized "best-practice" techniques is imperative to increase across-study accuracy and comparability. New studies should take advantage of the increasing number of methodological papers being published on eDNA collection techniques, extraction techniques, and the statistical analysis of eDNA data (e.g., Lacoursière-Roussel, Rosabal, et al., 2016;Hinlo et al., 2017;Chambert, Takahara, Pilliod, Goldberg, & Doi, 2018) to optimize and standardize field, laboratory, and analytical practices.
As more experiments are continually published, future meta-analyses could identify particular methods that might lead to more robust correlations between abundance and eDNA particle concentration, Despite experimental design limitations inherent in natural systems, there is substantial cause for optimism in its application as a tool to monitor abundance; eDNA particle concentration, depending on inclusion of smaller studies, explained approximately 50%-57% of variation in abundance on average in natural environments.
Whether this, in itself, is enough to accurately monitor abundance in populations of conservation concern remains unknown; substantial variation remains, on average, unaccounted for. However, the field of eDNA research is relatively new and methodologies are consistently improving; with just minimal improvements in technology or technique, eDNA could represent an extremely powerful tool for predicting species-specific population abundance for conservation applications in nature.

DATA A R C H I V I N G S TAT E M E N T
Data for this study are available from the Dryad Digital Repository: https://doi.org/10.5061/dryad.pm66g75

ACK N OWLED G M ENTS
We