An examination of seasonal variation in taxonomic richness and community composition using eDNA on a tropical coral reef

Small volumes of water containing environmental DNA (eDNA) are increasingly combined with metabarcoding to generate biodiversity data for specific fractions of marine flora and fauna. To date, however, few studies have utilized this technique to assess how well it captures seasonal patterns in coral reef communities or how environmental or methodological factors influence eDNA detections. In our study, we used three eDNA metabarcoding assays primarily targeting bony fish, elasmobranchs, as well as cnidarians and sponges (Cnidaria/Porifera) combined with monthly seawater sampling to 1) investigate temporal variation in taxonomic detections and 2) statistically test the potential effect of season, sea surface temperature, timing of spawning (using moon phase as a proxy), and sample preservation on taxon detection across a 12-month period in a model coral reef system (Big Vicki’s Reef, Lizard Island, Great Barrier Reef, Australia). Species-level fish and genus-level scleractinian coral detections from standardised visual surveys conducted at the same coral reef, in addition to a curated list of all known fishes recorded from the more expansive coral reef system, were used to validate eDNA detections. Our eDNA dataset indicated that the number of taxa detected were consistently highest in September for fish, and in February followed by September for Cnidaria/Porifera. Conversely, detections were lowest in June and July for all taxa. Some, but not all, of the environmental and methodological variables explained the observed temporal pattern in biological communities or systematic changes in the number of taxa, and in some cases, this effect was taxon dependent. Our study also highlights the significance of timing in eDNA biodiversity surveys conducted on tropical coral reefs in the Southern Hemisphere. To obtain the most meaningful estimates of site diversity, we recommend focusing sampling efforts between early spring to early autumn. Alternatively, allocating an entire year to sampling would better capture seasonal variation and provide more comprehensive insights into coral reef biodiversity.

Methods
Laboratory processing

DNA was extracted from half of each filter membrane (including filtration controls) within six months of collection using the DNeasy Blood and Tissue Kit (QIAGEN) in the TrEnD Laboratory at Curtin University in Western Australia with the following modifications: 540 µl of ATL lysis buffer, 60 µl of Proteinase K, and a 3-hr digest at 56°C. An extraction blank was processed in parallel with every set of eleven samples to detect any cross-contamination. Three previously published primer sets (16S Fish, CoralITS2, and CoralITS2_acro; Table 1) were employed in this study to amplify primarily teleost fish and scleractinian coral taxa, respectively. Quantitative PCR (qPCR) was carried out in 25 μl containing the following concentrations: 1X AmpliTaq Gold® PCR buffer (Life Technologies, Massachusetts, USA), 2 mM MgCl2 (Fisher Biotec, Australia), 0.4μM dNTPs, 0.1mg BSA (Fisher Biotec, Australia), 0.4 μM each of forward and reverse primers (Integrated DNA Technologies, Australia), 0.6 μl of 5X SYBR® Green (Life Technologies), 1U AmpliTaq Gold® DNA Polymerase (Life Technologies), 2 μl of eDNA template, and made to volume with Ultrapure™ Distilled Water (Life Technologies). All qPCRs were prepared in dedicated trace DNA (clean room) facilities at the TrEnD Laboratory, Curtin University, and amplified on a StepOnePlus Real-Time PCR System (Applied Biosystems, Massachusetts, USA) with the following PCR cycling conditions: initial denaturation at 95°C for 5 min, followed by 40 cycles of 95°C for 30 s, 52-55°C for 30 s (see respective annealing temperatures in Table 1), and 72°C for 45 s, with a final extension of 72°C for 10 min. Each sample was amplified in duplicate in a single-step process via the use of fusion-tagged primer architecture, whereby original primers are flanked by a unique 6–8 bp multiplex identifier tag (MID-tag) and an Illumina compatible sequencing adapter. 

MID-tagged PCR amplicons were then pooled into two libraries (the first containing 16S Fish amplicons, and the second CoralITS2 and CoralITS2_acro amplicons combined) at equimolar ratios based on qPCR ΔRn values. The two libraries were then size-selected using a Pippin-Prep (Sage Science, Beverly, USA) to remove any amplicons outside of the target range (150-450 bp, and 160-600 bp, respectively). Size-selected libraries were then purified using the QIAquick PCR Purification Kit (Qiagen, Venlo, Netherlands), quantified using a Qubit 4.0 Fluorometer (Invitrogen, Carlsbad, USA), and diluted to 2 nM for loading onto the sequencing platform. The 16S Fish library was loaded onto a 300 cycle Illumina MiSeq® V2 Standard Flow Cell for unidirectional sequencing. The CoralITS2 library was loaded onto a 500 cycle Illumina MiSeq® V2 Standard Flow Cells for paired-end sequencing. The sequencing of both libraries was conducted on an Illumina MiSeq platform (Illumina, San Diego, USA) housed in the TrEnD Laboratory.

Bioinformatics

Sequences were demultiplexed into their respective samples based on their MID-tags using the ngsfilter and obisplit commands in OBITools (v1.2.9; Boyer et al. 2016) for the single-end (16S Fish) sequencing reads, and the insect package (Wilkinson et al. 2018) for the paired-end (CoralITS2 and CoralITS2_acro) sequencing reads. Data from the two sequencing libraries were separately quality filtered (minimum sequence length = 100 bp, maximum expected errors = 2, no ambiguous nucleotides), denoised, merged (paired-end reads only, 20 bp overlap, no mismatches), filtered for chimeras and dereplicated (pool = TRUE) using the DADA2 package (Callahan et al. 2016) in R (v3.5.3; R Core Team 2015). To maximize biodiversity detections, we retained high-quality singletons (that passed strict quality filtering, denoising, and chimera removal) if they were observed in at least two independent samples (DADA2 Pool=TRUE function). If the singleton was only found in one sample, it was removed from the dataset. This generated an amplicon sequence variant (ASV) .fasta file and count table; the former was queried against NCBI’s GenBank nucleotide database (accessed September 2020; Benson et al. 2005) using BLASTn (minimum percentage identity = 90, maximum target sequences = 10, reward value = 1) with Zeus, an SGI cluster, based at the Pawsey Supercomputing Centre in Kensington, Western Australia. Taxonomic assignments of ASVs were curated and, when necessary, collapsed to the lowest common ancestor (LCA) using a custom LCA python script (Mousavi‐Derazmahalleh et al. 2021), with a threshold query coverage (qCov) of 100%, a minimum percentage identity of 90%, and a difference (Diff) of 1. ASVs detected in extraction blanks with a read count >10 were removed across the entire dataset, prior to ASVs being merged by taxonomy using the phyloseq v1.24.2 ‘tax_glom’ function (McMurdie and Holmes 2013) in R. Sample replicates were then pooled by their temporal sampling period and taxa accumulation curves were plotted by sequencing depth using the vegan v2.5-7 rarecurve function (Oksanen et al. 2019) in R. These rarefaction analyses (Fig. S2 and S3 in Online Resource 1) indicated that subsampling of sequencing depth was not required, given that a plateau was apparent in the taxa accumulation curves.