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

level of detection probability) and assessed these models using model comparison and cross- validation methods adapted for this multi-scale approach (Table 2). Three of the model com- parison methods (WAIC, CPO and one cross-validation score) were in agreement that Model 7 (ψ(water depth) θ(.) p(sequencing depth)) was the best model, while the Brier score from cross-validation suggested Model 1 (the null model) and Model 3 (ψ(.) θ(.) p(sequencing depth)) were the best models. We considered Model 7 our best model moving forward, given that most selection methods indicated this was the best model. We assessed model fit using Bayesian p-values and diagnostic plots for all models but pres- ent the results for the best model only. We obtained a Bayesian p-value of 0.51, suggesting that Model 7 (ψ(water depth) θ(.) p(sequencing depth)) provided a good fit to our data overall; diag- nostic plots, however, revealed higher deviance at sites with lower water depth, suggesting a poorer model fit at shallower sites (see S3 File). The community-wide estimate for occupancy was 0.29 (HDI: 0.22–0.36). Water depth had a significant effect on the community mean occu- pancy (Fig 2A and 2B), and we detected considerably more species at the shallowest sites com- pared to the other sites (274 species at two shallow water sites combined compared to 109 species across all six deep water sites). The community-wide probability of capture was 0.96 (HDI: 0.92–0.99) and the community-wide probability of detection was 0.14 (HDI: 0.12–0.17).

Sequencing depth did not have a significant effect on the probability of detection for most spe- cies in this case study (Fig 2C and 2D). Species-specific estimates of occupancy, capture proba- bility, and detection probability were also obtained from the model (S2 Table). We estimated species richness for the survey area by running the best model (Model 7) with data augmentation. This model used the probabilities of capture and detection to estimate the number of species missed in sampling efforts. We detected 231 species overall, and the esti- mated species richness for the survey area was 286 (HDI: 264–306), indicating that 55 (HDI: 33–75) species were undetected during our surveys. In other words, our survey detected ~81% of the estimated species in our study area. While it was not selected as our best model, we also ran Model 9 (ψ(.) θ(.) p(amplicon)) with a group-level effect (phylum) to investigate how the probability of each amplicon varied by phylum. Amplicons displayed relatively similar probabilities of detection at the commu- nity-level (Fig 3), however the probabilities of detection for each amplicon varied considerably when comparing between phyla, where some amplicons clearly failed to detect certain taxo- nomic groups (Fig 4). Discussion We applied a multi-species, multi-scale occupancy model to a DNA metabarcoding dataset generated from marine water samples and explored how the inclusion of categorial and con- tinuous covariates at different levels improved model performance. The best model included

sequencing depth as a covariate at the level of detection and water depth as a covariate at the level of occupancy, where we observed a higher species richness at shallower sites. One of the shallow water collection sites was within 1 km of a sewage outflow, which may have contrib- uted to this result, although a high species richness was also observed at the second, shallow water site located >10 km from the sewage outflow. While sequencing depth was included in the best model, we did not observe a strong effect of sequencing depth. However, the samples were all sequenced on a NovaSeq instrument, which generates an unprecedented number of reads, yielding very high sequencing depths (mean number of filtered sequences per sample ± standard deviation: 8,519,055 ± 2,514,998) compared to many other barcoding stud- ies (e.g. [44,45]). In studies where the mean sequencing depth is lower, differences in sequenc- ing depth are likely to have greater effects [5,46]. In such cases, analyzing data using PLOS ONE Multi-species occupancy models for metabarcoding PLOS ONE | https://doi.org/10.1371/journal.pone.0224119 March 19, 2020 8 / 17