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

We included sequencing depth (number of reads per sample per amplicon) as a continuous covariate (α1) at the level of probability of detection. Additionally, we included amplicon iden- tity as a categorical covariate (α2) at the level of probability of detection. We included water depth (m) as a continuous covariate (α3) at the level of occupancy. Covariates were included in the model as follows: logitðpijrkÞ ¼ lpk þ b1k � a1ijr þ β2k a2ijr logitðcikÞ ¼ lpsik þ β3k � a3i Continuous covariates were z-score standardized to have a mean of zero and a standard devia- tion of one to help with model convergence. We compared a null model (Model 1) with no covariates with seven models with different combinations of covariates (Table 2) using WAIC and CPO and using Brier and logarithmic scores for cross-validation. We assessed model fit using Bayesian p-values based on deviance residuals and by looking at diagnostic plots to examine model fit. We plotted deviance residuals for each species, site, and covariate. Species coefficients arise from additional community-level parameters: lpsik � Nðmlpsi; slpsiÞ lthetak � Nðmltheta; slthetaÞ lpk � Nðmlp; slpÞ β1k � Nðmβ1; sβ1Þ β2k a2 � Nðmβ2; sβ2Þ β3k � Nðmβ3; sβ3Þ Community-level parameters were described by weakly informative hyperpriors [31]. All mean values for the above prior distributions were selected from a normal distribution and all standard deviations were selected from a uniform distribution. m � Nð0; 10Þ s � Uniformð0; 5Þ Prior sensitivity was assessed by running the models with various prior parameterizations.

Posterior distributions were similar across all priors. All statistical analyses were conducted in R v3.5.1 [38]. MCMC sampling was achieved with JAGS [39], implemented using ‘jagsUI’ v1.5.0 [40]. The model was written for JAGS in the BUGS language (see S2 File for BUGS model structure). We fit models using known species richness to conduct our model comparisons, and assessed models and model fit to determine the best model. MCMC sampling was run in three chains, each with 50,000 iterations, a burn in of 10,000, and a thinning rate of 10. Convergence was verified using the Gelman-Rubin diagnostic [41] and by evaluating trace plots. For all models, we report parameter estimates as the mean of the posterior distribution with the 95% highest posterior density interval (HDI; [42]) calculated using ‘HDInterval’ v0.2.0 [43]. We conducted a data augmented model with unknown species richness for the best model at varying levels of augmentation to determine PLOS ONE Multi-species occupancy models for metabarcoding PLOS ONE | https://doi.org/10.1371/journal.pone.0224119 March 19, 2020 6 / 17