Environmental DNA in a global biodiversity hotspot: Lessons from coral reef fish diversity across the Indonesian archipelago

Sampling Sites We collected eDNA samples across Indonesia, spanning a strong fish biodiversity gradient (Bellwood & Meyer, 2009; Roberts et al., 2002). Sampling focused on 7 reef systems within three regions (Fig. 1), including: 1) outside the Coral Triangle in Western Indonesia (Aceh and Batam-Bintan), 2) lower diversity regions of the Coral Triangle in Central Indonesia (Derawan and Wakatobi), and 3) high diversity regions of the Coral Triangle (Eastern Indonesia: Lembeh Strait, Ternate, and Raja Ampat) that have the world’s highest reef fish biodiversity (Allen & Werner, 2002; Bellwood & Meyer, 2009; Roberts et al., 2002). eDNA Sampling To assess marine fish diversity with eDNA, we employed a hierarchical sampling design (Table S1). To capture spatial variability in both habitat and eDNA signatures, within each of the 7 reef systems above, we sampled multiple sites separated by at least five kilometers (Andruszkiewicz et al., 2017; Miya et al., 2015), except for Wakatobi, where sites were 2.5 km distant. To maximize the capture of diversity at each sampling site and account for fine-scale spatial heterogeneity in local eDNA signatures, we collected three one-liter “biological replicates” at each site following standard sampling protocols used in temperate ecosystems (Miya et al., 2015); each replicate was no more than 3 meters from each other. To minimize variation in community composition associated with depth, each seawater sample was collected on SCUBA, 1 m off the reef, at depths between 11-15 m. To isolate eDNA from seawater samples, we filtered one liter of seawater through a 0.22-micron Sterivex™ filter (Millipore®, SIGMA MILLIPORE) following the methods of Miya et al. (2015) and Curd et al. (2018) with one key modification: we collected individual water samples in sterile one liter Kangaroo™ Gravity Feeding Bags (similar to intravenous drip bags) that allow for gravity filtration through the Sterivex columns, a method ideally suited to remote field locations. In addition, we filtered one liter of bottled drinking water at each location to serve as a “field blank” negative control. Filters were stored in a -20˚ C freezer until eDNA was extracted at the Indonesian Biodiversity Research Center (now “BIONESIA”). eDNA Extraction, Amplification, and Sequencing We extracted eDNA samples and blanks using the DNeasy Blood & Tissue Kit (QIAGEN, Hilden, Germany) following the modified extraction protocol of Spens et al. (2017), adding 720 µL of ATL buffer and 80 µL of proteinase K directly into the filter cartridge. Extracted DNA was then stored at -20˚C and transported to the University of California, Los Angeles for amplification and library preparation. We amplified extracted eDNA using the Multiplex PCR Kit (QIAGEN, Hilden, Germany), using the MiFish Teleost 12S rRNA mitochondrial DNA (mtDNA) primers (Miya et al., 2015). We conducted each individual eDNA sample via Polymerase Chain Reactions (PCR) in triplicate to account for potential PCR bias (Andruszkiewicz et al., 2017; Miya et al., 2015; Taberlet et al., 2012). Each PCR reaction consisted of 12.5 µL Qiagen 2x Master Mix, 2.5 µL (2 mM) of the primer, 6.5µL nuclease free water, and one µL the DNA extract. Thermocycling parameters utilized a touchdown protocol, beginning with a 15-minute pre-denaturation step at a 95 °C, followed by a touchdown thermocycling profile consisting of 30 seconds denaturing at 94 °C, 30 seconds annealing at 69.5 °C, and 30 seconds extension at 72 °C, with the annealing temperature dropping by 1.5 °C per cycle until 50 °C. Following this initial touchdown phase, the main cycle consisted of 25 cycles of 94 °C for 30 seconds for denaturation, 50 °C for 30 seconds for annealing and 72 °C for 45 seconds for extension, concluding with a 10-minute final extension at 72 °C. We confirmed successful PCR reactions through electrophoresis of 5µL products for 30 minutes at 150 volts on 2% agarose gels stained with 6x SYBR™ Green (ThermoFisher Scientific, Waltham, MA, USA) for visualization. To create the sequencing libraries, we pooled the triplicate PCR products from each individual eDNA sample into a single tube and purified these pooled PCR products using the Serapure bead protocol (Faircloth et al., 2014). Next, we quantified the DNA concentration (ng/µL) of each pooled PCR sample using the Qubit Broad Range dsDNA assay (ThermoFisher, Waltham, MA, USA) following the manufacturer protocol, and then then subsequently normalized pooled PCRs to equimolar concentration. We then used the Nextera DNA Library Preparation Kit (Illumina, San Diego, CA, USA) to index each pooled sample using a unique combination of Nextera i5 and i7 primers in a second PCR reaction, following the manufacturer protocol (Curd et al., 2018). The indexing PCR reaction consisted of 12.5 µL Kapa High Fidelity Master Mix (Roche, Basel, Switzerland) 0.625 µL of 1 µM i5 Illumina Nextera indices, 0.625 µL of 1 µM i7 Illumina Nextera indices, and 11.25 µL of PCR product for a total of 10ng of DNA. The indexing PCR thermocycling parameters began with an initial denaturation of 95 ˚C for 5 minutes, followed by 8 cycles of: 98 ˚C denaturation for 30 seconds, 56 ˚C annealing for 30 seconds, and 72 ˚C extension for 3 minutes, ending with a 72 ˚C extension for five minutes. To ensure the indexing PCR was successful, we electrophoresed indexed PCR products at 120 V for 45 minutes on a 2% agarose gel. We then created the final sequencing library by combining the cleaned and quantified indexed PCR products in equal concentration of 10 ng/µl per sample. Libraries were pooled by barcode resulting in two pooled libraries. The final libraries were sequenced at the University of California Berkeley sequencing core on an Illumina MiSeq platform utilizing 300 base pair paired end sequencing. Bioinformatics and Data Analysis To determine fish community composition, we used the Anacapa Toolkit (version: 1) to conduct quality control, amplicon sequence variant (ASV) parsing, and taxonomic assignment using user-generated custom reference databases (Curd et al., 2018). The Anacapa Toolkit sequence QC and ASV parsing module relies on cutadapt (version: 1.16) (Dobin et al., 2013), FastX-toolkit (version: 0.0.13) (Gordon & Hannon, 2010), and DADA2 (version 1.6) (Callahan et al., 2016) as dependencies, and the Anacapa classifier modules relies on Bowtie2 (version 2.3.5) (Langmead & Salzberg, 2012) and a modified version of BLCA (Gao, et al., 2017) as dependencies. We processed sequences using the default parameters and assigned taxonomy using a CRUX-generated reference database supplemented with additional 12S reference sequences generated by the Smithsonian (GenBank accessions MZ597881-MZ598481, Table S2). We note that CRUX relies on ecoPCR (version: 1.0.1) (Boyer et al., 2016), blastn (version: 2.6.0) (Altschul et al., 1990; Camacho et al., 2009), and Entrez-qiime (version: 2.0) (Baker, 2016) as dependencies. The raw ASV community table was decontaminated following Kelly et al. (2018) and McKnight et al. (2019) (See Supplemental Methods). However, we did not conduct site occupancy modeling as we were specifically interested in exploring under sampling of total ASV diversity. Because of noted issues with incomplete reference databases (Juhel et al., 2020), diversity analyses focused on ASVs rather than ASVs successfully assigned to species. In addition, because we focused on fish and the MiFish primers amplify vertebrates more broadly, we filtered out any ASVs that were not Actinopterygii and Chondrichthyes ASVs. Similarly, because of sequence similarities in 12S, ASVs can be erroneously assigned to taxa that are not native to local ecosystems (Gold et al., 2021). As such, we filtered out any ASVs assigned to fishes not known from the Coral Triangle. After applying these filters, we then transformed all read counts into an eDNA index for beta-diversity statistics (Kelly et al., 2019). Patterns of Biodiversity Across Indonesia We conducted analyses of fish biodiversity in a hierarchical fashion. First, we examined diversity at the level of individual sample replicate, where fish diversity was represented by eDNA sequences amplified from each of the three one-liter water samples at an individual sampling site. Second, we examined diversity at the site level, where fish diversity was represented by combining replicate eDNA sequences amplified from each site (Table S1). Lastly, we examined diversity in each of the seven sampled regions by combining diversity across all sites sampled. To better understand the efficacy of standard eDNA sampling protocols in high diversity marine ecosystems, we compared ASV recovery from our high diversity Indonesian reefs to Scorpion Bay, Santa Cruz Island, California, a lower diversity temperate marine ecosystem (Gold et al., 2021) as samples were collected, processed, an analyzed using the same methods outlined above. All analyses were conducted in in R (R Core Team, 2020). We examined patterns of alpha diversity by exploring ASV richness at each hierarchical level. We then compared sample coverage estimates by conducting ASV accumulation curves using the iNext package (version: 2.0.20) in R (Hsieh et al., 2016). Specifically, we explored the accumulation of ASVs 1) across all replicate samples within each site, 2) within all replicate samples within each region, and 3) within all replicate’s sites within each region. We then explored betadiversity patterns across sample replicates, sites, and regions by calculating Jaccard-binary dissimilarity indices and running PERMANOVA and companion homogeneity of dispersion test (adonis and betadisper functions from vegan (version: 2.5-6) in R). Results of the PERMANOVA were visualized using a nonmetric multidimensional scaling (NMDS) through the phyloseq (version: 1.32.0) and vegan packages in R (R Core Team, 2020). To explore the scale of species turnover across the biogeographic regions, we conducted zeta-diversity analyses using the zetadiv (version 1.2.0) package (Latombe et al., 2017). Zeta diversity measures the degree of overlap between any set of observed communities (Hui & McGeoch, 2014), in contrast to betadiversity which only compares pairwise overlap (Hui et al., 2018; Latombe et al., 2017; McGeoch et al., 2019). Importantly, the zeta diversity framework allows for the calculation of diversity metrics across multiple samples and sampling strata to better understand patterns of biodiversity (McGeoch et al., 2019), particularly the scale and shape of community turnover (Hui et al., 2018; McGeoch et al., 2019; Simon et al., 2019). From this framework there are two important zeta diversity metrics: zeta diversity decay and species retention rates (McGeoch et al., 2019). Zeta diversity decay measures the decreasing number of shared taxa across greater observed communities (McGeoch et al., 2019). For small numbers of samples (n<5) this can be visually represented as the center of the Venn Diagram between multiple observed communities (i.e., individual eDNA samples). We compared species retention rates and zeta diversity decay (i.e., the decreasing number of shared taxa across more observed communities) across 1) all replicate samples within each site, 2) within all replicate samples within each region, and 3) within all replicate sites within each region. In particular, we compare the rate of zeta diversity decay and species retention rates across the region to identify sites with lower community turnover across sampling strata. We further explore zeta diversity decay over distance to explore the scale of community turnover across sample replicates, sites, and regions through eDNA data. eDNA Survey and Visual Census Data Many eDNA studies do contemporaneous paired comparisons of eDNA to visual census data to compare the efficacy of each approach for biodiversity monitoring (e.g. Gold et al., 2021; Valdivia-Carrillo et al., 2019; Fernández et al., 2020, Stoeckle et al., 2021, Stat et al., 2019). However, comprehensive fish biodiversity studies in the Coral Triangle involve high-intensity visual census surveys involving multiple divers spending 6 hours per day over 10-21 days across a wide range of habitats to capture as much fish diversity as possible (G. Allen, M. Erdmann, pers comm.). Because we could not conduct contemporaneous paired studies, we explored the differences in eDNA and high-intensity visual surveys in capturing local fish diversity by sampling eDNA at three reef sites within Raja Ampat that G. Allen and M. Erdmann comprehensively sampled previously with visual surveys during multiple dives over multiple days. We then compared taxa recovered by eDNA to the visual census data, not to determine which was superior, but instead to better understand the relative strengths of weaknesses of each approach.

采样位点：本研究在印度尼西亚全境采集环境DNA（environmental DNA, eDNA）样本，覆盖了显著的鱼类生物多样性梯度（Bellwood & Meyer, 2009; Roberts et al., 2002）。采样聚焦于三个区域内的7个珊瑚礁系统（图1），具体包括：1）印度尼西亚西部珊瑚三角区（Coral Triangle）外海域（亚齐与巴淡-民丹）；2）印度尼西亚中部珊瑚三角区的低多样性海域（德拉旺与瓦卡托比）；3）珊瑚三角区的高多样性海域（印度尼西亚东部：勒姆贝海峡、特尔纳特与拉贾安帕特），该区域拥有全球最高的礁栖鱼类生物多样性（Allen & Werner, 2002; Bellwood & Meyer, 2009; Roberts et al., 2002）。

eDNA采样：为利用eDNA评估海洋鱼类多样性，本研究采用分层采样设计（补充表S1）。为捕捉栖息地与eDNA信号的空间异质性，在上述7个珊瑚礁系统中，除瓦卡托比海域的采样点间距为2.5 km外，其余每个系统内均设置至少间隔5 km的多个采样位点（Andruszkiewicz et al., 2017; Miya et al., 2015）。为最大化单一位点的多样性捕获量并抵消局部eDNA信号的精细空间异质性，本研究参照温带生态系统的标准采样规程（Miya et al., 2015），在每个采样位点收集3份1升的“生物学重复样本”，各重复样本间的间距不超过3 m。为最小化与水深相关的群落组成差异，所有海水样本均通过水肺潜水（SCUBA）在距礁体1 m、水深11~15 m处采集。

为从海水样本中分离eDNA，本研究参照Miya等（2015）与Curd等（2018）的方法，使用0.22 μm Sterivex™过滤器（密理博®，SIGMA MILLIPORE）过滤1升海水，仅做一处关键改进：采用无菌1升Kangaroo™重力喂养袋（类似静脉滴注袋）收集单份海水样本，可通过重力驱动过滤通过Sterivex柱，该方法极适用于偏远野外场景。此外，本研究在每个采样位点过滤1升瓶装饮用水，作为“野外空白”阴性对照。过滤后的滤膜保存于-20℃冰箱，直至在印度尼西亚生物多样性研究中心（现更名为“BIONESIA”）完成eDNA提取。

eDNA提取、扩增与测序：本研究使用DNeasy Blood & Tissue Kit（凯杰（QIAGEN），德国希尔登），参照Spens等（2017）的改良提取规程，直接向过滤柱中加入720 μL ATL缓冲液与80 μL蛋白酶K完成eDNA提取。提取得到的DNA保存于-20℃，随后转运至加州大学洛杉矶分校进行扩增与文库制备。

本研究采用Multiplex PCR Kit（凯杰（QIAGEN），德国希尔登），使用MiFish硬骨鱼12S rRNA线粒体DNA（mtDNA）引物（Miya et al., 2015）对提取的eDNA进行扩增。为抵消聚合酶链式反应（Polymerase Chain Reaction, PCR）潜在的偏差，每份eDNA样本均设置3次重复PCR（Andruszkiewicz et al., 2017; Miya et al., 2015; Taberlet et al., 2012）。单次PCR反应体系为：12.5 μL Qiagen 2×预混酶、2.5 μL（2 mM）引物、6.5 μL无核酸酶水与1 μL DNA提取物。热循环程序采用降落PCR方案：首先在95℃预变性15分钟；随后进入降落热循环阶段，包含94℃变性30秒、69.5℃退火30秒、72℃延伸30秒，每轮循环退火温度降低1.5℃直至降至50℃；初始降落阶段结束后，主循环包含25轮94℃变性30秒、50℃退火30秒、72℃延伸45秒，最终在72℃延伸10分钟终止反应。

本研究通过琼脂糖凝胶电泳验证PCR扩增成功与否：取5 μL PCR产物在2%琼脂糖凝胶中以150 V电泳30分钟，使用6×SYBR™ Green（赛默飞世尔科技，美国马萨诸塞州沃尔瑟姆）染色后成像。为构建测序文库，本研究将每份eDNA样本的3次重复PCR产物合并至同一管中，使用Serapure磁珠纯化方案（Faircloth et al., 2014）纯化合并后的PCR产物。

随后参照制造商规程，使用Qubit宽范围双链DNA定量试剂盒（赛默飞世尔科技，美国马萨诸塞州沃尔瑟姆）测定每份合并PCR产物的DNA浓度（ng/μL），并将所有合并PCR产物归一化至等摩尔浓度。本研究随后使用Nextera DNA文库制备试剂盒（伊鲁米纳（Illumina），美国加利福尼亚州圣地亚哥），通过第二轮PCR为每份合并样本添加唯一的Nextera i5与i7引物组合作为索引，具体规程参照Curd等（2018）。索引PCR反应体系为：12.5 μL Kapa高保真预混酶（罗氏（Roche），瑞士巴塞尔）、0.625 μL 1 μM i5伊鲁米纳Nextera索引引物、0.625 μL 1 μM i7伊鲁米纳Nextera索引引物与11.25 μL PCR产物，总DNA投入量为10 ng。

索引PCR热循环程序为：95℃初始变性5分钟；随后8轮循环：98℃变性30秒、56℃退火30秒、72℃延伸3分钟；最终72℃延伸5分钟终止反应。本研究通过琼脂糖凝胶电泳验证索引PCR成功与否：将索引PCR产物在2%琼脂糖凝胶中以120 V电泳45分钟。随后将纯化并定量后的索引PCR产物以10 ng/μL的等浓度混合，构建最终测序文库。本研究按条形码将文库合并为2个合并文库，最终文库在加州大学伯克利分校测序中心使用伊鲁米纳MiSeq平台完成300 bp双端测序。

生物信息学与数据分析：为确定鱼类群落组成，本研究使用Anacapa Toolkit（版本1）进行质量控制、扩增子序列变异（amplicon sequence variant, ASV）解析与分类学注释，所用参考数据库为自定义用户数据库（Curd et al., 2018）。Anacapa Toolkit的序列质量控制与ASV解析模块依赖cutadapt（版本1.16）（Dobin et al., 2013）、FastX-toolkit（版本0.0.13）（Gordon & Hannon, 2010）与DADA2（版本1.6）（Callahan et al., 2016）；Anacapa分类器模块依赖Bowtie2（版本2.3.5）（Langmead & Salzberg, 2012）与改良版BLCA（Gao et al., 2017）作为依赖工具。

本研究使用默认参数处理序列，并使用CRUX生成的参考数据库补充史密森学会提交的额外12S参考序列（GenBank登录号：MZ597881~MZ598481，补充表S2）完成分类学注释。本研究注意到CRUX依赖ecoPCR（版本1.0.1）（Boyer et al., 2016）、blastn（版本2.6.0）（Altschul et al., 1990; Camacho et al., 2009）与Entrez-qiime（版本2.0）（Baker, 2016）作为依赖工具。原始ASV群落表参照Kelly等（2018）与McKnight等（2019）的方法完成去污染（详见补充方法）。但本研究未进行位点占用建模，因为本研究的核心目标为探究总ASV多样性的采样不足问题。

鉴于参考数据库不完整的已知问题（Juhel et al., 2020），本研究的多样性分析以ASV为单位，而非成功注释至物种水平的ASV。此外，由于本研究聚焦于鱼类，且MiFish引物可广谱扩增脊椎动物，因此本研究过滤掉所有不属于辐鳍鱼纲（Actinopterygii）与软骨鱼纲（Chondrichthyes）的ASV。同理，由于12S序列的相似性，部分ASV可能被错误注释至非本地生态系统的类群（Gold et al., 2021），因此本研究过滤掉所有被注释为珊瑚三角区未记录鱼类的ASV。经过上述过滤步骤后，本研究将所有读长计数转换为eDNA指数，用于β多样性统计分析（Kelly et al., 2019）。

印度尼西亚海域生物多样性格局：本研究以分层方式开展鱼类生物多样性分析。首先，在单个样本重复水平分析多样性，以单个采样位点的3份1升海水样本扩增得到的eDNA序列代表鱼类多样性；其次，在采样位点水平分析多样性，以单个采样位点所有重复样本的eDNA序列合并结果代表鱼类多样性（补充表S1）；最后，在7个采样区域水平分别分析多样性，以对应区域所有采样位点的多样性合并结果代表区域鱼类多样性。

为更好地理解标准eDNA采样规程在高多样性海洋生态系统中的有效性，本研究将本研究在印度尼西亚高多样性珊瑚礁海域获得的ASV恢复结果，与美国加利福尼亚州圣克鲁斯岛蝎子湾的低多样性温带海洋生态系统的结果进行对比（Gold et al., 2021），两处样本的采集、处理与分析方法均遵循本研究前述规程。所有分析均在R语言（R核心团队, 2020）环境中完成。

本研究通过探究各分层水平的ASV丰富度分析α多样性格局；随后使用R语言中的iNext包（版本2.0.20）绘制ASV累积曲线，比较样本覆盖度估计值（Hsieh et al., 2016），具体包括：1）单个采样位点内所有重复样本的ASV累积情况；2）单个区域内所有采样位点的重复样本的ASV累积情况；3）单个区域内所有采样位点的ASV累积情况。

本研究通过计算Jaccard二元相异指数、进行置换多元方差分析（PERMANOVA）及配套的离散同质性检验（使用R语言vegan包版本2.5-6中的adonis与betadisper函数），探究采样重复、采样位点与区域间的β多样性格局。PERMANOVA的结果通过非度量多维标度（nonmetric multidimensional scaling, NMDS）可视化，所用工具为R语言中的phyloseq（版本1.32.0）与vegan包。

为探究生物地理区域间的物种周转尺度，本研究使用zetadiv包（版本1.2.0）开展ζ多样性分析（Latombe et al., 2017）。ζ多样性衡量任意一组观测群落间的共有类群重叠程度（Hui & McGeoch, 2014），与仅比较成对群落重叠的β多样性不同（Hui et al., 2018; Latombe et al., 2017; McGeoch et al., 2019）。重要的是，ζ多样性框架可计算多份样本与多个采样分层下的多样性指标，从而更好地解析生物多样性格局（McGeoch et al., 2019），尤其是群落周转的尺度与模式（Hui et al., 2018; McGeoch et al., 2019; Simon et al., 2019）。该框架包含两项关键的ζ多样性指标：ζ多样性衰减与物种保留率（McGeoch et al., 2019）。ζ多样性衰减衡量随着观测群落数量增加，共有类群数量的下降趋势（McGeoch et al., 2019）。对于少量样本（n<5），该趋势可通过多个观测群落间的韦恩图中心直观展示。

本研究比较了三类场景下的物种保留率与ζ多样性衰减（即随着观测群落数量增加，共有类群数量的下降趋势）：1）单个采样位点内的所有重复样本；2）单个区域内所有采样位点的重复样本；3）单个区域内的所有采样位点。本研究特别比较了各区域间的ζ多样性衰减速率与物种保留率，以识别采样分层间群落周转较低的位点。本研究进一步通过eDNA数据探究ζ多样性衰减随距离的变化，以解析采样重复、采样位点与区域间的群落周转尺度。

eDNA调查与目视普查数据：许多eDNA研究均同步开展eDNA与目视普查数据的配对对比，以评估两种方法在生物多样性监测中的有效性（例如Gold et al., 2021; Valdivia-Carrillo et al., 2019; Fernández et al., 2020; Stoeckle et al., 2021; Stat et al., 2019）。然而，珊瑚三角区的综合性鱼类生物多样性研究需要开展高强度目视普查：多名潜水员每日工作6小时，连续10~21天覆盖多种生境，以尽可能捕获全部鱼类多样性（G. Allen, M. Erdmann, 私人通信）。

由于本研究无法开展同步配对研究，因此本研究通过在拉贾安帕特的3个珊瑚礁位点采集eDNA，与G. Allen与M. Erdmann此前通过多日多次潜水开展的全面目视普查结果进行对比。本研究将eDNA检测得到的类群与目视普查数据进行对比，目的并非评判两种方法的优劣，而是更好地理解两种方法各自的相对优势与局限。