eDNAPlus: A unifying modelling framework for DNA-based biodiversity monitoring

Ecology is undergoing a technology revolution that is making it possible to rapidly generate species inventories via automated and high-throughput DNA sequencers and via electronic sensors, such as drones, satellites, camera traps, and acoustic recorders. These techniques can, if coupled with appropriate algorithms and databases, simultaneously identify large numbers of target species, including those that are cryptic, difficult-to-access, tiny, and low-abundance (Tosa et al.,, 2021; van Klink et al.,, 2022; Bush et al.,, 2017; Besson et al.,, 2022; Piper et al.,, 2019; Ley,, 2022). So far, the most efficient method for generating species-resolution inventories is DNA-based surveys, which rely on reading DNA barcodes: short, standardized sections of the genome that can be compared to a reference library to enable taxonomic identifications without the need to examine organism morphologies (Ratnasingham and Hebert,, 2007).

DNA barcoding refers to the identification of single species (Hebert et al.,, 2003), and DNA metabarcoding refers to the detection of large numbers of species from environmental DNA (eDNA), which is the collective name for DNA isolated from environmental samples (Bohmann et al.,, 2014; Pawlowski et al.,, 2020; Taberlet et al.,, 2018). These environmental samples include water (Thomsen and Willerslev,, 2015), soil (Frøslev et al.,, 2019), air (Clare et al.,, 2022; Lynggaard et al.,, 2022), and bulk tissue (i.e. mass-trapped organisms) (Ji et al.,, 2013). For instance, Thomsen and Sigsgaard, (2019) demonstrated that traces of eDNA on flower petals could be analysed to describe the diversity of arthropods that visit wildflowers, including pollinators, parasitoids, predators, and herbivores. Ji et al., (2022) used the trace amounts of residual vertebrate blood left in 30,468 blood-sucking leeches to map vertebrate wildlife across a 677 km² nature reserve in China. Finally, Abrego et al., (2021) sequenced $542$ mixed-species, bulk-tissue samples of arctic arthropods captured over $14$ years and showed that species richness in the study site had declined by 50% during a time period in which local mean temperature had increased by 2C.

The potential of DNA-based surveys for monitoring and managing biodiversity comes with a number of statistical challenges. Firstly, species-specific absolute abundances cannot be estimated using DNA data alone. Secondly, DNA-based surveys yield data that have been subjected to several types of error and noise (see Section 1.1), some of which are species-specific. The framework presented in this paper addresses these challenges by developing a novel model and corresponding efficient inferential tools. Using our framework, we model within-species change in DNA biomasses (described in Section 1.1), which under certain conditions can be considered as a proxy for change in abundance, hence addressing the first challenge. To address the second challenge, we propose a hierarchical crossed-effects model that expresses all major sources of variation, error and noise in the data collection and analysis pipeline, whilst accounting for correlation across species and across sites, and for covariate effects on DNA biomass. We also model frequently employed controls at the PCR stage and evaluate their effect on inference.

We assume that $M_{i}$ physical samples are collected from site $i$, $i=1,\ldots,n$, and $K_{im}$ PCR replicates are performed on the $m$-th sample from site $i$. We denote by $y^{s}_{imk}$ the number of DNA reads of the $s$-th species, $s=1,\ldots,S$ in the $k$-th PCR replicate of the $m$-th sample collected at the $i$-th site. We have $n_{z}$ site covariates and $X_{i}^{z}$ represents their value at site $i$ and $n_{w}$ sample covariates, represented as $X_{im}^{w}$ for the $m$ sample at the $i$-th site. In what follows, $i$ indexes sites, $m$ samples, $k$ PCR replicates, and $s$ species.

Our proposed model (see Figure 2) is hierarchical, with three levels. The first level models the amount of DNA biomass of each species at the surveyed sites, which is a function of environmental and landscape covariates as well as between-species and between-sites correlation (Biomass availability). The second level models the amount of biomass collected for each species in each physical sample from each site (Biomass collection). Lastly, the third level models the number of reads obtained for each species in each PCR from each physical sample (Biomass analysis). Data are observed only at the third level, as the result of Stage 2 of the survey, with levels one and two corresponding to latent states.

We model the reads conditional on $c^{s}_{imk}$ as follows. Conditional on $c^{s}_{imk}=1$, $y^{s}_{imk}\sim\text{NB}(\exp(\lambda_{s}+v^{s}_{im}+u_{imk}+o_{imk}),r_{s})$, where $\lambda_{s}$ models the Stage 2 species effect on the amplification rate, $u_{imk}$ is the Stage 2 pipeline effect, with $u_{imk}\sim\text{N}(0,\sigma_{u}^{2})$, $o_{imk}$ is an offset modeling the normalisation steps described in Section 1.1, and $r_{s}$ is a species-specific variance of the Stage 2 noise. If more than one normalisation steps are employed, then they can all be incorporated into the same offset i.e. they can be added up. Conditional on $c^{s}_{imk}=0$, $y^{s}_{imk}\sim\pi\delta_{0}+(1-\pi)(1+\text{NB}(\mu_{0},n_{0}))$, that is, there are zero reads with probability $\pi$, and non-zero but negligible reads otherwise. Finally, conditional on $c^{s}_{imk}=2$, $y^{s}_{imk}\sim\text{Pois}(\tilde{\mu}_{s})$. The negative binomial is parameterised in terms of the mean and the number of failures. A visual representation of the PCR process when $c^{s}_{imk}=1$ is shown in Figure $1$ of the Supplementary material.

Stage 2 negative control samples (which are known to not contain DNA of any species) can be easily accounted for in our model by having additional samples for which $\tilde{\delta}^{s}_{l}=\tilde{\gamma}^{s}_{l}=0$. Accounting for spike-ins corresponds to having $S^{\star}$ additional species for which ($v^{S+1}_{im},\dots,v^{S+S^{\star}}_{im})$ is known. Since the pipeline effect is shared across all species (including spike-ins), the known values of $v^{s}_{im}$ for the spike-ins help to better estimate $u_{imk}$. We further investigate this effect in Section 4.

The model is summarised in Figure 2 (a), the directed acyclic graph of the model is shown in Figure 2 (b), while a graphical representation of the latent variables introduced across both stages is shown in Figure 2 (c).

The model presented in Figure 2 is not identifiable in its general form unless certain constraints are applied, as we discuss below. For example, if we define $\tilde{v}^{s}_{im}:=v^{s}_{im}-\eta_{s}-l^{s}_{i}$ and $\tilde{l}^{s}_{i}:=l^{s}_{i}-\beta^{s}_{0}$ the model for $\theta^{s}_{im}$ and $y^{s}_{imk}$ conditional on $c^{s}_{imk}=1$ and all offsets $o_{imk}$ set to $0$ can be expressed as

The reason for this unidentifiability is that data are observed only in the third level of the model, and hence the following sets of species-specific parameters are confounded: the baseline amount of biomass across all sites ($\beta_{0}^{s}$) with the baseline collection rate ($\eta_{s}$) and the baseline amplification rate ($\lambda_{s}$), and the former again with the baseline detection rate $\phi^{s}_{0}$. However, by assuming that all these baseline rates are constant across sites, samples, and PCRs, we are able to infer species-specific changes in biomass across sites and therefore covariate effects and correlations between species and between sites.

For inferential purposes, we set $\phi^{s}_{0}\equiv 0$ and $\eta_{s}\equiv 0$. Using these constraints, the new baseline (log) amount of biomass, $(\beta^{s}_{0})^{\star}$, is equal to $\beta^{s}_{0}+\eta_{s}$, which means that we can only estimate the sum of the baseline amount of available biomass and corresponding baseline collection. Similarly, the new baseline (logit) collection probability $(\phi^{s}_{0})^{\star}$ is equal to $\phi^{s}_{0}-\phi^{s}_{1}\eta_{s}$, and therefore the baseline collection probability is also confounded with the baseline collection rate. As a result, we cannot infer the amount of available biomass separately from the collection rate, and hence the estimates of log biomass obtained, as mentioned above, are only meaningful for comparison within each species. For the same reason, comparisons of absolute amount of biomass across species are not meaningful.

We also note that depending on the survey design in terms of the number of samples collected per site and the number of PCR replicates per sample, additional sets of parameters can be confounded and not estimable. Specifically, the following pairs of parameters are confounded:

Finally, we note that if the offsets $o_{imk}$ introduced in the model due to the several normalizations occurring in the pipeline are not recorded, estimates from the model would be biased. However, a potential way to mitigate this bias is the introduction of spike-ins, which contribute to the estimation of the “overall” pipeline effects $\tilde{u}_{imk}=u_{imk}+o_{imk}$.

Samples can be drawn from the posterior distribution of the parameters using a Gibbs sampler. Posterior sampling is greatly helped by representing the negative binomial distribution as a Poisson-gamma mixture, which allows many parameters to be updated in closed form from their full conditional distribution.

For the parameters $\sigma_{s}$, $\mu_{s}$, $B$ and $B_{0}$, the full conditional distribution is available in closed form, and therefore posterior sampling is straightforward. We use simple random walk Metropolis Hastings steps for parameters $\pi$, $\mu_{0}$, $n_{0}$, and $r_{s}$ and Metropolis-Hastings steps with a Laplace approximation proposal for the parameters $l^{s}_{i}$, $\lambda_{s}$, $v^{s}_{im}$, $u_{imk}$ and $r_{s}$. However, on its own, this naive Gibbs sampler will mix slowly since we have a complex hierarchical model with cross effects and many latent variables. We address this by updating parameters in blocks using re-parameterisation and an adaptive updating scheme for the discrete latent variables.

To illustrate our approach to blocking and re-parameterisation, we consider the error-free version of our model

A naive Gibbs sampler updating each parameter from its full conditional leads to prohibitively slow mixing, due to the form of the likelihood where $\lambda_{s}$, $v^{s}_{im}$ and $u_{imk}$ appear as a sum. To address the slow mixing in the nested effects, $\lambda_{s}$ and $v^{s}_{im}$, the use of a centred parameterisation (Papaspiliopoulos et al.,, 2007) has been suggested, which corresponds to defining $\bar{v}^{s}_{im}:=\lambda_{s}+v^{s}_{im}$ and $\bar{l}^{s}_{i}:=\lambda_{s}+l^{s}_{i}$. However, issues of slow mixing still exist between $\bar{v}^{s}_{im}$ and $u_{imk}$ and, as noted by Zanella and Roberts, (2021), re-parameterisation does not improve mixing in the case of crossed-effects models. In a classic crossed-effect model of the form $y_{jkl}\sim\text{N}(\lambda+v_{j}+u_{k},\sigma^{2})$, Papaspiliopoulos et al., (2020) propose a collapsed Gibbs sampler by first jointly sampling $\lambda$ with $v_{j}$ and then $\lambda$ jointly with $u_{k}$. However, this approach does not scale well in our setup, since it would involve sampling all the $\lambda_{s}$ and $u_{imk}$ jointly, which have dimensions $S$ and $\sum_{i,m}K_{im}$ respectively. Zanella and Roberts, (2021) propose the use of identifiability constraints on the model, which in Equation ((4)) correspond to assuming $\sum_{s}v^{s}_{im}=\sum_{k}u_{imk}=0$. Since sampling conditionally on constraints can be challenging, we propose a simpler strategy to improve mixing that is more suited to our framework. We consider re-parameterising to the factor averages $\hat{v}_{im}=\frac{1}{S}\sum_{s=1}^{S}\bar{v}^{s}_{im}$ and $\hat{u}_{im}=\frac{1}{K}\sum_{k=1}^{K}u_{imk}$ and the factor increments $\tilde{v}^{s}_{im}=\bar{v}^{s}_{im}-\hat{v}_{im}$ and $\tilde{u}_{imk}=u_{imk}-\hat{u}_{im}$ and performing an update by first sampling jointly the factor means conditional on the increments, that is, from $(\hat{v}_{im},\hat{u}_{im}|\tilde{v}^{s}_{im},\tilde{u}_{imk})$ and next using the standard updates $(u_{imk}|v^{1}_{im},\dots,v^{S}_{im})$ and $(v^{j}_{im}|u_{im1},\dots,u_{imK})$. In our simulations, we have found that jointly updating the factor means considerably improves mixing. The sampling scheme for the complete model is presented in the Supplementary material.

The indicator variables $(\delta^{s}_{im},\gamma^{s}_{im},c^{s}_{imk})$ can be updated directly from their full conditional distributions but, since there are $nMS(\bar{K}+2$) (where $\bar{K}$ is the average number of PCR replicates) of these variables and often one value of $(\delta_{imk},\gamma_{imk},c_{imk})$ has probability very close to $1$, evaluating every full conditional distribution in every iteration can be very time-consuming and computationally wasteful. Therefore, we use a cheap approximation as a proposal in an MH step. Specifically, every $B$ iterations, we update the approximation $\hat{p}((\delta^{s}_{im},\gamma^{s}_{im},c^{s}_{imk})=(\epsilon_{1},\epsilon_{2},\epsilon_{3}))=\frac{1}{T}\sum_{t=1}^{T}\mbox{I}\left((\delta^{s}_{im})^{(t)},(\gamma^{s}_{im})^{(t)},(c^{s}_{imk})^{(t)})=(\epsilon_{1},\epsilon_{2},\epsilon_{3})\right)$, where $(\delta^{s}_{im})^{(t)},(\gamma^{s}_{im})^{(t)},(c^{s}_{imk})^{(t)}$ is the value of $(\delta^{s}_{im},\gamma^{s}_{im}c^{s}_{imk})$ at the t-$th$ iteration, $\mbox{I}(A)$ is the indicator function which takes the value 1 if $A$ is true and 0 otherwise, and $T$ is the number of current iterations. Using this update scheme, we only need to evaluate the likelihood if the state is proposed to change. If the probability of one state is close to one, the adaptive scheme often proposes the current state, which can be accepted without computation. The adaptive scheme does not affect convergence of the MCMC algorithm since the approximation clearly has diminishing adaptation and the state space of the indicator variables is discrete (see e.g. Roberts and Rosenthal,, 2009, for more discussion of conditions for convergence of adaptive MCMC schemes).

In this section, we investigate the properties of the model under different study designs in terms of the number of sites, samples per site, and PCRs per sample, as well as the number of spike-ins. In each section, we consider the estimates of the differences in log-biomass within species, when that is not a function of site-specific covariates (no covariate case), and the estimates of the covariate coefficients when log-biomass is a function of a single continuous covariate (regression case).

In Section 4.1 we present theoretical results using a continuous version of our model that does not account for error in either stage. In Section 4.2 we fit our model as presented in Section 2 under different scenarios for study design by varying the number of sites, number of samples per site, and number of PCRs per sample. Finally, in Section 4.3, we explore the effect of spike-ins for different levels of noise in each stage of the process and different study designs.

We apply our model to an unpublished dataset of mostly arthropod invertebrates collected using $121$ Malaise-trap samples from $89$ sample sites in the H. J. Andrews Experimental Forest (HJA), Oregon USA (225 km²) in July 2018. Each trap was left to collect for seven days, and samples were transferred to fresh 100% ethanol to store at room temperature until extraction. The management objective that motivated the collection of this dataset is to interpolate continuous species distributions among the $89$ sample points so that areas of higher and lower conservation value at the HJA can be identified.

For each sample, the collected invertebrate biomass was combined with a lysis buffer, in an amount proportional to the starting biomass, to digest the tissue, and a fixed aliquot was then taken from the overall mixture for DNA extraction and subsequent PCR. This normalization, as described in Section 2, was accounted for in the model by setting the offset $o_{imk}$ equal to the log ratio between the aliquot and the overall amount of liquid mixture in each case.

We included $50$ species in the study by selecting the species that have the most non-zero counts across all PCR replicates. Log-biomass is modelled as a function of two environmental covariates: log elevation and log distance-to-road.

Figure 6 presents the 95% posterior credible intervals (PCIs) for the species-specific coefficients of log elevation and log distance-to-road in the model for log-biomass. The effects of the covariates on species biomass are not consistent within each taxonomic order, which suggests low phylogenetic inertia at this rank for response to these landscape characteristics. Elevation is a stronger predictor for species biomass than distance-to-road for this ecosystem. This makes ecological sense, since distance-to-road is only expected to exert an effect over a ca. hundred metres, via canopy openness, whereas elevation exerts a pervasive effect via its effects on temperature, precipitation, and vegetation.

Figure 7 (a) presents the posterior mean of the between-species residual correlations. Species in the Diptera (flies, spp. 14-30) exhibit higher positive correlations with each other, as well as with several species in the Hymenoptera (ants, bees, and wasps) and Lepidoptera (butterflies and months). We conservatively interpret these positive residual correlations as indicative of unmeasured environmental covariates, such as canopy openness, rather than of biotic interactions. We also note that two species in the Lepidoptera, (spp. 41, 43), one in the Hymenoptera (sp. 33), and one in the Psocodea (barklice, sp. 50) are among the few species showing strong negative residual correlation with many of the other species, and again, we interpret these correlations as indicative of unmeasured environmental covariates. There is a strongly positive, pairwise correlation between two tabanid fly species Hybomitra liorhina and Hybomitra sp. indet (spp. 12, 13), which might indicate the oversplitting of one biological species into two OTUs during the bioinformatic pipeline. Finally, there is also a strongly positive, pairwise correlation between the moth species Ceratodelia gueneata (sp. 44) and the predatory fly (Scathophagidae, Microprosopa sp. indet, sp. 20), which might indeed indicate a specialised predator-prey relationship.

In Figure 7 (b) we show the biodiversity map for the area, which is useful for identifying areas of higher species richness and compositional distinctiveness, which together can be used to identify areas of higher conservation value (i.e. higher ‘site irreplaceability’ sensu Baisero et al.,, 2022). The predicted mean log-biomasses on a continuous map over the HJA for all individual species are presented in the Supplementary Material. These can be used to identify species with a wide spatial range, such as the click beetle (Megapenthes caprella), or with a restricted range, such as the leafhopper (Osbornellus borealis).

Finally, Figure 8 (a) suggests that generally, there is a similar amount of variation between sites and between samples for these species. As suggested by Figure 8 (b), the species considered have similar collection probabilities across the several sites, possibly due to the fact that the most frequently detected species across PCRs have been selected. Figure 8 (c) demonstrates, as expected, that the Stage 2 true positive probability is close to 1 for all species. Similarly, the figure also suggests that the probability of a Stage 2 false negative error is very close to 0 for all but three species. One of these three (sp. 14) is in the fly family Tachinidae, which are parasitoids of other insects and thus might have been collected not only as adults but also occasionally as eggs attached to the adults of other (insect) host species, with the latter case being classified as false positives in Stage 2, given that an egg would contribute very low amounts of starting DNA biomass.

Over the last decade, DNA-based biodiversity studies, primarily using metabarcoding, have rapidly increased in popularity, and multivariate statistical models are now starting to be deployed to analyse metabarcoding data (e.g. Lin et al.,, 2021; Pichler and Hartig,, 2020; Abrego et al.,, 2021; Fukaya et al.,, 2022; Ji et al.,, 2022). Our paper provides the first unifying modelling framework that considers and quantifies all main sources of variation, error and noise in metabarcoding surveys (Table 1). As a result, our modelling framework should allow more reliable and more powerful biodiversity monitoring and inference on species responses to landscape characteristics than has been possible before. We have employed, extended, and developed a number of inferential tools to deal with the complexity of the proposed hierarchical model, which involves two latent stages and a large number of latent variables. Finally, this is the first modelling approach that accounts for spike-ins and negative controls (empty tubes), which are widely used quality-control methods in DNA-based biodiversity surveys but rarely explicitly considered within a modelling framework. We explored the benefits of spike-ins on inference and, for the first time, explored issues of study design in metabarcoding data, providing both theoretical and simulation results.

Our new framework allows us to infer and map species biomass as well as biodiversity across surveyed sites (Figure 7 (b)), and to link these to landscape characteristics (Figure 6). The resulting maps can be used to identify areas of high conservation value, as well as areas where particular species or groups of species are more or less prevalent, and to detect species-specific shifts, expansions, and shrinkage. We are also able to study pairwise correlations across large numbers of species (Figure 7 (a)), which is considerably more scalable using metabarcoding data than using standard observational data. We have shown that using spike-ins can substantially increase inference accuracy for parameters of interest (Figure 5). Our results also demonstrate that the current practice of collecting a single sample from each surveyed site considerably reduces our ability to infer changes in species biomass and that replication at both stages as well as the use of spike-ins is the optimal approach to designing metabarcoding studies (Figure 3).

In metabarcoding data, the baseline biomass of each species is confounded with its species-specific collection and amplification rates. Hence, we cannot infer absolute values of species-specific biomass across sites using metabarcoding data alone. However, by assuming that baseline species-specific collection and amplification rates are the same across sites, samples, and PCR replicates, we can infer species-specific biomass change across sites, species-specific covariate effects, and pairwise species correlations. These assumptions are justifiable in the case study, where eDNA release rates cannot differ across sites because the samples are whole organisms, and there is also no reason to believe that Malaise traps vary in trapping efficiency across sites. In other scenarios, eDNA release or degradation rates could vary across the environment, due to large variations in, for instance, food availability and water pH, respectively. In these cases, site covariates could be used to account for this suspected variation.

Generally, modelling changes in (proxies of) abundance, such as changes in biomass, is a more powerful monitoring tool than modelling changes in species presence across surveys sites (Joseph et al.,, 2006). Metabarcoding studies yield count data without any consequence on associated cost, and hence overcome the time and cost implications associated with collecting count data for multiple species. Our model uses the raw count data, and does not rely on ad-hoc rules about what constitutes a practically zero count for converting them to binary data, which has been the standard practice thus far (Ovaskainen et al.,, 2017; Bush et al.,, 2020).

The model can easily be extended to account for multiple primers, which essentially introduce different levels in Stage 2, with each level corresponding to a set of primers. For instance, vertebrate eDNA surveys can profitably use two primers (e.g. the 16S mammal and 12S vertebrate primers used in Ji et al.,, 2022), since the two primers together detect the largest set of vertebrates, with overlap. Additionally, our approach can be extended and applied to metagenomic data. In metagenomics, no PCR is employed. Instead, all the DNA in each mixed-species sample is sequenced, and the DNA barcodes (or other taxonomically informative genetic sequences) are discovered bioinformatically.

Metabarcoding studies, particularly when applied to microbiomes and meiofauna (e.g. nematodes, micro-eukaryotes), can detect 10000s of species, which leads to large numbers of latent variables and coefficients in the model. There are several ways that the inferential tools presented here could be further extended to scale to these cases. Firstly, the posterior distribution conditional on the $u_{imk}$ is independent across species. If $u_{imk}$ could be estimated at a first stage then inference across species could be easily parallelized. Secondly, variational Bayes methods could be applied to avoid the use of sampling methods. The choice of variational distribution will be important and can exploit the conditional normality of much of the model. Alternatively, the model could be adapted by assuming that the coefficient matrices such as $\beta^{z}=(\beta^{z}_{1},\dots,\beta^{z}_{S})$, have a low-dimensional representation.

The metabarcoding process produces large numbers of different DNA-barcode sequences, which represent not only different species but also within-species genetic diversity and outright errors generated during PCR and sequencing. Currently, these sequences are clustered into species hypotheses as part of the bioinformatic process, following a number of heuristic rules and thresholds (hence, the word ‘operational’ in OTU), but future work could explore building the OTU table during the model-fitting process.

We are not modelling species presence/absence and instead we have focused on modelling biomass on a continuous scale. As a result, we cannot infer whether a species is absent from a particular study site, but instead only if its biomass at given site is practically zero for modelling purposes. We have assumed that a sample which already contains biomass of a species cannot be further contaminated by the DNA of the same species from another sample or site in Stage 1. This is a reasonable but also necessary assumption. Under best practice, contamination in Stage 1 is expected to be relatively rare, and therefore, there is not enough information in the data to partition the collected biomass between that which was truly collected from the site and that which was contamination from elsewhere. If this assumption is violated, then the amount of collected biomass is no longer proportional to the amount of available biomass, and the amount of available biomass at sites where there has been contamination can be overestimated.

eDNA metabarcoding has revolutionised the cost-effectiveness, precision and scale at which biodiversity assessment can be performed. Nevertheless, the multiple stages at which imperfect detection of biomass can occur during the workflow are not insignificant. By facilitating estimates of within-site changes in biomass and covariates while accounting for workflow uncertainties, our modelling framework provides a step-change improvement in the design and analysis of eDNA metabarcoding data.