Validating metabarcoding-based biodiversity assessments with multi-species occupancy models: A case study using coastal marine eDNA

Materials and methods Model formulation The multi-species, multi-scale occupancy model. We used a Bayesian modeling frame- work to develop a multi-species, hierarchical occupancy model with three stochastic levels: occupancy (ψ), probability of capture (θ), and probability of detection (p) (Fig 1). The occu- pancy process describes whether sampling sites are occupied or not by a given species’ DNA. For eDNA sampling, there are often two levels of sampling replication within each site (e.g. [20,27]): biological replicates are samples collected from a single site in the field and technical replicates are repeated samples taken from a single biological replicate in the lab. The probabil- ity of capture refers to the probability that a species’ DNA is collected in a sample, given that the species was present at the site. The probability of detection refers to the probability that a species was detected in a technical replicate, given that the species’ DNA was collected in the sample. This model assumes no false positives occur in the data. While false positives may be a possibility in metabarcoding data [15], we used strict bioinformatic filtering to reduce this pos- sibility (see Bioinformatics below). Further comments on false positives can be found in the Discussion. This model can be fit to a dataset, yijrk, which is a binary indicator of whether a species k (k = 1,2,. . .K) was detected (1) or not detected (0) in a technical replicate r (r = 1,2,. . .R) from a given sample j (j = 1,2,. . .J) at a given site i (i = 1,2,. . .I). The model consists of three coupled

Bernoulli trials to describe a four-dimensional array of data yijrk. yijrkjwijk � BernoulliðpijrkwijkÞ wijkjzik � BernoulliðyijkzikÞ zik � BernoulliðcikÞ We model our observations of detection (yijrk = 1) or non-detection (yijrk = 0) of species k in replicate r in sample j at site i as a random variable having parameter pijrk, which describes the probability of detection. The second random variable, wijk, describes the capture (wijk = 1) or non-capture (wijk = 0) of species k in sample j at site i having parameter θijk (probability of cap- ture). The third random variable, zik, describes the occurrence (zik = 1) or non-occurrence (zik = 0) of species k at site i having parameter ψik (occupancy probability). Model assessment and comparison. The goodness-of-fit of multi-species, multi-scale occupancy models can be assessed using Bayesian p-values of the deviance residuals. We adapted Bayesian p-value calculations from [28] for a multi-scale model (S1 File) to assess goodness-of-fit, where values close to 0.5 indicate a good fit and values >0.95 or <0.05 indi- cate a poor fit. We also adapted model selection and cross-validation calculations from [28] for multi- scale, multi-species occupancy models to determine the best model. We calculated the Wata- nabe-Akaike information criterion (WAIC; [29]) and the conditional predictive ordinate crite- rion (CPO; [30]), and then evaluated the results of k-fold cross validation using the Brier score and the logarithmic score. Complete calculations for all model assessment and comparison methods can be found in S1 File.

Unknown species richness. An additional level can be added to the model described above for communities with unknown species richness [10]. This model uses data augmenta- tion to estimate species richness for the sampling area through the inclusion of another PLOS ONE Multi-species occupancy models for metabarcoding PLOS ONE | https://doi.org/10.1371/journal.pone.0224119 March 19, 2020 3 / 17