Gut content metabarcoding of specialized feeders is not a replacement for environmental DNA assays of seawater in reef environments

Sample collection

Sampling was conducted over 7 days in April and May 2018 at 10 sites around Dongsha Atoll within the lagoon, outside the lagoon, and on the seagrass beds (Fig. 1). No more than two sampling sites were visited on any given day (Appendix S4). At sampling sites where seawater was collected, six 1 L replicates were obtained from a dive boat adjacent to the relevant habitat (at or within 2 m horizontal distance from the coral reef) and 30 cm below the surface using sterile Nalgene bottles (N = 36 in total). All water samples were immediately stored on ice in an enclosed cooler and filtered at the Dongsha Atoll Research Station (DARS). Each sample was individually filtered across Pall 0.2 µm Supor^® polyethersulfone membranes using a Pall Sentino^® Microbiology pump (Pall Corporation, Port Washington, NY, USA) within 6 h of collection. The filtration apparatus was cleaned by soaking in 10% bleach between samples for at least 15 min. We favoured a 1 L filtration volume as previous research has shown that this provides an acceptable compromise between the time to process samples (10 to 15 min per filter replicate) and the biodiversity uncovered in shallow coral reef environments (e.g., West et al., 2020; but also see Mathon et al., 2022). Filters were stored at −20 °C until processing in the Trace and Environmental DNA (TrEnD) Laboratory at Curtin University in Perth, Australia.

Adult specimens of a single species of butterflyfish (Chaetodon lunulatus), a known corallivore, were collected using pole spear on SCUBA taking care not to impale the stomach cavity, which would contaminate the gut contents (N = 40; Appendix S4). Fish collections were undertaken in accordance with the policies and procedures drafted by the Animal Care and Use Committee of National Sun Yat-Sen University. Most fish were dispatched by the action of the spear, but those that were not, were subject to immediate ike jime and placed in an ice slurry on board the vessel. Field permissions to undertake the research were obtained from the Ministry of Interior in Taiwan (Marine National Park Headquarters) under permit No. 1070001035. Following collection, fish were kept on ice in individual Ziploc^® bags until storing them at −20°C to mitigate further DNA degradation. To remove gut contents, the fish were partially thawed, the gut cavity was opened, and the entire semi-frozen stomach contents were placed into 2 ml sample tubes with sterile 80% ethanol taking care not to include the stomach lining; these were immediately stored at −20°C until processing in the Benthic Organisms and Molecular Ecology Lab at National Sun Yat-Sen University. Between each sample, all tools and trays were either replaced or rinsed in bleach (10%) and sterile ethanol (99%).

Due to the patchy distribution of Chaetodon lunulatus across Dongsha Atoll, fish were not sampled at all sites where seawater was initially collected. That said, the distributions of collected fishes from a single sampling site were unlikely to overlap with adjacent sampling sites given that the home range of this species is small (average = 117 m²; range = 2.3–976.3 m²; Cowlishaw, 2014).

DNA extraction

Following DiBattista et al. (2020), DNA bound to filter membranes was extracted using Qiagen DNeasy Blood and Tissue kits (Qiagen, Hilden, Germany) that was automated on a Qiacube (Qiagen, Hilden, Germany) to minimise human handling and cross-contamination in the TrEnD Laboratory with the following modifications to the standard protocol: 540 µl of ATL lysis buffer and 60 µl of Proteinase K. DNA extraction controls (blanks) were carried out for every set of 12 extractions. Half of each filter paper was manually shredded with sterilised dissecting scissors prior to immersion in the extraction buffer; the other half of the filter paper was re-frozen at −80 °C to serve as a permanent voucher.

Between 58 and 260 mg of gut contents were spooned into 2 or 7 ml Minilys homogenizing vials depending on the volume of the sample. Samples were homogenized manually by using a micropestle (Geneaid Biotech Ltd, Taiwan, China) for 30 s and then placed in a sterile 1.5 ml tube. Due to the co-purification of inhibitors in stomach samples, DNA extractions were performed using the QIAamp DNA stool kit (Qiagen, Hilden, Germany) according to the standard protocol in a laminar flow hood. DNA extraction controls, for which all steps remained the same except for the addition of stomach contents, were included.

Fusion-tag qPCR

Following DiBattista et al. (2020), a primer set targeting 18S (V1-3 hypervariable region; 18S_uni_1F: 5′–GCCAGTAGTCATATGCTTGTCT–3′; 18S_uni_400R: 5′–GCCTGCTGCCTTCCTT–3′; Pochon et al., 2013) was used to maximise detection of general eukaryotes, hereafter the 18Suni assay. We used a second assay targeting ITS2 (scl58SF: 5′–GARTCTTTGAACGCAAATGGC–3′; scl28SR: 5′–GCTTATTAATATGCTTAAATTCAGCG–3′; Brian, Davy & Wilkinson, 2019) to increase the taxonomic resolution of identified Anthozoa (Phylum Cnidaria) and Demospongiae (Phylum Porifera) within each environmental sample, hereafter the ITS2 assay (also referred to as “CP1 assay” in Alexander et al. (2020)).

As per DiBattista et al. (2020), quantitative PCR (qPCR) experiments for both seawater and gut content extracts were set up in a separate ultra-clean laboratory at Curtin University designed for trace DNA work using a QIAgility robotics platform (Qiagen, Hilden, Germany). For more details on contamination mitigation measures, please see DiBattista et al. (2019, 2020). In brief, fusion-tag qPCR was performed with each extract in triplicate on a StepOnePlus Real-Time PCR System (Applied Biosystems, Foster City, CA, USA) under the following conditions for 18Suni: initial denaturation at 95 °C for 5 min, followed by 45 cycles of 30 s at 95 °C, 30 s at 52 °C, and 45 s at 72 °C, with a final extension for 10 min at 72° C. PCR was performed under the following conditions for ITS2: initial denaturation at 95 °C for 5 min, followed by 45 cycles of 30 s at 95 °C, 30 s at 55 °C, and 45 s at 72 °C, with a final extension for 10 min at 72° C. These primers and PCR conditions have been optimized and applied elsewhere (Alexander et al., 2020; DiBattista et al., 2019, 2020, 2022; West et al., 2020). Moreover, the primer design incorporated indexes in both the forward and reverse primer, which allowed us to tag individual PCR replicates of individual samples by a unique combination of tags on the forward and reverse primers. To check for contamination, non-template control (labelled as NTC) PCR reactions were run alongside the template PCR reactions, which only contained master mix including the assay primers.

Libraries for sequencing were made by pooling amplicons into equimolar ratios based on PCR Ct values and the endpoint of amplification curves, which were then size selected using a Pippin Prep (Sage Science, Beverly, MA, USA; 100 to 600 bp) and purified using the Qiaquick PCR Purification Kit (Qiagen Inc, Venlo, The Netherlands). The volume of purified library added to the sequencing run was determined against DNA standards of known molarity on a LabChip GX Touch (PerkinElmer Health Sciences, Waltham, MA, USA). Final libraries were sequenced in a paired end approach using a 500 cycle MiSeq^® V2 Reagent Kit and standard flow cell on an Illumina MiSeq platform (Illumina, San Diego, CA, USA) located in the TrEnD Laboratory.

Bioinformatic filtering

Sequence data were quality filtered (QF) prior to ZOTU (Zero-radius Operational Taxonomic Unit) creation and taxonomic assignment. ZOTUs are equivalent to amplicon sequence variants (ASVs) and provide a more precise measurement of sequence variation (also see Edgar, 2018). Metabarcoding reads recovered by paired-end sequencing were first merged using the Illumina MiSeq analysis software under the default settings (i.e., minimum 10 bp overlap between read 1 and read 2). Sequences were then assigned to samples based on their unique index combinations and trimmed in Geneious^® Pro v 4.8.4 (Drummond et al., 2009). Only those sequences with 100% identity matches to Illumina adaptors, index barcodes, and template specific oligonucleotides were kept for downstream analyses. Sequences were further processed in Geneious by trimming for base quality (phred score > 30) and minimum length (18Suni, 200 base pairs; ITS2, 100 base pairs) using BBDuk v 36.92 (Bushnell, 2017). USEARCH v 11.0.667 (Edgar, 2010) was then used for dereplication into unique sequences, removing singletons, removing chimeric sequences, and generating a ZOTU (minimum number of sequences for ZOTU creation = 10) table based on the UNOISE algorithm, which performs denoising (error-correction) of the amplicon reads. USEARCH and the taxonomic assignment pipeline outlined below (BLAST, LULU, and LCA assignment) was facilitated by an in-house script developed by Mousavi‐Derazmahalleh et al. (2021) implemented on a high-performance computing platform at the Pawsey Supercomputing Centre in Western Australia.

Taxonomic assignment

ZOTUs were queried against the National Centre for Biotechnology Information’s (NCBI) GenBank nucleotide database (accessed in 2020/2021) using BLASTn v2.10.1 with the following conservative settings: percentage identity (99), query coverage of 100, best hit score edge of 0.05, best hit overhang of 0.25, and an E-value of 1e⁻³. LULU (Frøslev et al., 2017) was then run to curate the assignments and eliminate any remaining redundant sequences, prior to assigning the lowest common ancestor (LCA) with a custom Python script, with the default parameters: minimum_ratio_type = min, minimum_ratio = 1, minimum_match = 84, minimum_relative_cooccurence = 0.95. Based on the LCA process, some ZOTUs were assigned at the species level, some at the genus level, some at the family level, and the remainder only at higher taxonomic ranks, as noted. All additional information related to this bioinformatic workflow is in Mousavi‐Derazmahalleh et al. (2021) as well as freely available scripts at the following link (https://github.com/mahsa-mousavi/eDNAFlow). All ZOTUs that were detected from extraction controls were removed from further analyses (Appendix S1). This included ZOTUs assigned to Trichoderma, Malasseziomycetes, Chlorophyta, Oedogoniaceae, uncultured Chlorophyta, uncultured Chytridiomycota, Cultellothrix coemeterii, Eukaryota, Dictyoceratida, Ephydatia, Brassica rapa, Blastocystis sp., and uncultured eukaryotes for the 18Suni assay. This also included ZOTUs assigned to Montipora, Montipora digitata, Montipora flabellata, Poritidae, Porites, Porites lobata, Porites lutea, Eukaryota, and Dysideidae for the ITS2 assay. Note that ZOTU1841 was assigned to a primate sequence and thus also removed from the 18Suni analysis. Also note that “BLANK4” was the only extraction blank where the number of sequences passed the filtering thresholds for the ITS2 analysis. Order level dendrograms for the 18Suni and ITS2 assay taxonomic assignments were constructed using a phylogenetic tree generator based on NCBI taxonomy (phyloT v2022.3 and iTOL v5; Letunic & Bork, 2021).

Statistical analyses

Initial PERMANOVA tests PRIMER v7 (Clarke & Gorley, 2015) across all samples based on presence/absence Jaccard similarity matrices showed that the community assemblage recovered with each method (water and gut content ZOTUs) was significantly different (18Suni: df = 1, MS = 68,668, Pseudo-F = 43.426, p < 0.0001; ITS2: df = 1, MS = 33,147, Pseudo-F = 11.089, p < 0.0001), resulting in a near complete lack of overlap for the detected ZOTUs (Appendix S2). Subsequent analyses further investigating trends across all sites were performed on subsets of the data separated by sampling approach (Appendix S3 for PRIMER input). The site composition of ZOTUs were therefore analysed separately for 18Suni and ITS2 in presence/absence format using PRIMER v7.

To test whether eDNA detections from seawater and butterflyfish gut contents displayed similar biodiversity profiles, we only considered sites where both approaches were applied (DS4, DS6, DS13; also see Appendix S1–S4). For these analyses, Jaccard similarity matrices were calculated based on presence/absence transformed data. PERMANOVA (9,999 permutations) was subsequently used to quantify these differences and principal coordinate analysis (PCO) was used to visualise the data (Anderson & Willis, 2003).

Total biodiversity detected

A total of 3,206,449 amplicon reads pertaining to the 18S gene that passed QF were retrieved from 36 seawater samples (average per replicate = 89,068 ± 7,085 SEM) and 3,058,957 amplicon reads were retrieved from 38 butterflyfish gut samples (average per replicate = 80,498 ± 10,995 SEM; Appendix S4). Overall, the 18Suni metabarcoding data produced assignments pertaining to 12 eukaryotic phyla, 22 classes, 40 orders, 43 families, and 38 genera across the combined seawater and butterflyfish gut samples (Appendix S1; Fig. 2A). The average number of ZOTUs assigned per replicate of seawater was 115.94 ± 8.77 SEM and for gut content it was 9.16 ± 0.60 SEM. At sites where both seawater and gut contents were sampled, on average, 93.31 ± 1.60% SEM of the ZOTUs were unique to seawater, 3.39 ± 1.36% SEM of the ZOTUs were unique to gut contents, and 3.30 ± 0.37% SEM of the ZOTUs were shared between the two sampling approaches. Here, one family of stony corals was identified as unique to gut contents (Acroporidae), whereas 34 families were identified as unique to seawater, with the greatest representation included in the classes Bivalvia (marine and freshwater molluscs; five families), Demospongiae (sponges; five families), Dinophyceae (dinoflagellates; four families), Florideophyceae (red algae; three families), and Polychaeta (bristle worms; four families). Two families of stony coral (Agariciidae and Pocilloporidae), coral-associated symbiotic algae (Symbiodiniaceae), and a family of cyclopoid copepods (Xarifiidae) were shared between the two sampling approaches.

Using the ITS2 assay (designed to specifically target cnidarians and sponges), a total of 1,919,189 amplicon reads that passed QF were retrieved from 36 seawater samples (average per replicate = 54,644 ± 4,320 SEM or 56,206 ± 4,145 SEM excluding outlier sample DGX22) and 1,928,056 amplicon reads were retrieved from 39 butterflyfish gut samples (average per replicate = 50,335 ± 2,164 SEM; Appendix S4). In total, the ITS2 assay (combined gut contents and eDNA) produced eight eukaryotic phyla, 10 classes, 14 orders, 25 families, and 28 genera across the sample sites (Appendix S1; Fig. 2B). The average number of ZOTUs assigned per replicate of seawater was 21.28 ± 2.56 SEM and for gut content it was 27.05 ± 1.84 SEM. At sites where both seawater and gut contents were sampled, on average, 51.15 ± 17.69% SEM of the ZOTUs were unique to seawater, 27.43 ± 11.65% SEM of the ZOTUs were unique to gut contents, and 21.42 ± 8.85% SEM of the ZOTUs were shared between the two sampling approaches. Here, two families of stony corals were only found in gut contents (Agariciidae and Dendrophylliidae), whereas 16 families were unique to seawater, with the greatest representation included in the classes Anthozoa (sea anemones, stony corals and soft corals; six families) and Demospongiae (three families). Four families of stony coral were shared between the two sampling approaches, which included Acroporidae, Fungiidae, Merulinidae, and Poritidae. In many cases, taxonomic assignment at the level of species was incomplete or ambiguous, reflecting gaps or uncertainty in reference material and/or taxonomy (8.3% and 9.5% identified to species for 18Suni and ITS2, respectively; Fig. 3).

Seawater vs gut content metabarcoding

For 18Suni, we found that community assemblages identified in each sample based on ZOTUs was highly dependent on the approach (i.e., seawater vs gut content; PERMANOVA: Pseudo-F = 43.426, df = 1, p < 0.0001; Appendix S2). This effect of sampling approach was consistent for the ITS2 assay (PERMANOVA: Pseudo-F = 11.089 df = 1, p < 0.0001; Appendix S2).

Based on PERMANOVA tests when the two sampling approaches were analysed separately, we found significant differences between sites for 18Suni sequences generated from seawater samples (Pseudo-F = 6.6691, df = 5, p = 0.0001; Fig. 4B) and gut content samples (Pseudo-F = 2.5081, df = 6, p = 0.0001; Fig. 4A). Pairwise analyses however demonstrated that all sites were significantly different from each other for seawater samples (full results in Appendix S5), compared with only 33% of site comparisons being significantly different for gut content samples (full results in Appendix S5). Significant pairwise site differences for gut contents were found between 5 of the 7 sites: BG vs DS6, p = 0.0261; BG vs DS10, p = 0.0218; BG vs DS9, p = 0.014; BG vs DS13, p = 0.009; DS9 vs DS10, p = 0.0025; DS9 vs DS13, p = 0.0147; DS10 vs DS13, p = 0.0013. DS4 and DWC were not distinguishable from other sites or each other based strictly on gut contents (p > 0.05 in all cases) (Appendix S5). For 18Suni gut content samples, Pearson correlations (r > 0.7) indicated that site differences were driven by Acropora (ZOTU143, ZOTU553) and Symbiodinium (ZOTU696, ZOTU822) (Fig. 4A). In contrast, for 18Suni seawater samples, except for ZOTU168 (uncultured diatom) and ZOTU2487 (Thraustochytriaceae), all other ZOTUs driving site differences were identified as only Eukaryota (Fig. 4B).

Based on PERMANOVA tests, we found significant differences between sites for ITS2 sequences generated from seawater samples (Pseudo-F = 3.5597, df = 5, p = 0.0001; Fig. 4D) and gut content samples (Pseudo-F = 4.1188, df = 6, p = 0.0001; Fig. 4C). Pairwise analyses indicated that all but one between site comparison (DS4 vs NSB) was significantly different for seawater samples. ITS2 sequences in gut samples appeared to be more distinct among sites compared to 18S sequences in gut samples, with 71% of site comparisons significantly different for the ITS2 assay (full results in Appendix S5). The specific pairwise ITS2 site differences between all sites for gut contents are provided in Appendix S5. Moreover, for ITS2 gut content samples, Pearson correlations were driven by corals in Fungiidae (ZOTU58, ZOTU106, ZOTU344) and Montipora (ZOTU53, ZOTU201, ZOTU262, ZOTU783) (Fig. 4C). For ITS2 seawater samples, Pearson correlations were driven by corals in Poritidae (ZOTU18, ZOTU44, ZOTU574) and Montipora (ZOTU153, ZOTU227, ZOTU240, ZOTU253, ZOTU337, ZOTU403, ZOTU472, ZOTU783) (Fig. 4D).

To address the uneven sampling effort between seawater and gut content at individual sites, we repeated the above statistical comparisons, but this time including only those sample sites where both approaches were applied (DS4, DS6, DS13; Fig. 5). In this case, based on PERMANOVA tests with 18S sequence data, sampling approach (Pseudo-F = 23.986, df = 1, p = 0.0001; Fig. 5) and site (Pseudo-F = 2.6769, df = 2, p = 0.0011; Fig. 5) were significantly different. However, there was also a significant interaction between site and sampling approach (Pseudo-F = 2.889, df = 2, p = 0.0006; Fig. 5). For the ITS2 assay, both sampling approach (Pseudo-F = 6.7726, df = 1 p = 0.0001; Fig. 5) and site (Pseudo-F = 3.0838, df = 1, p = 0.0001; Fig. 5) were significantly different based on PERMANOVA, and again there was a significant interaction between site and sampling approach (Pseudo-F = 2.9675, df = 2, p = 0.0001; Fig. 5).

Differences in semi-quantitative metrics for both assays such as relative read abundance, the percentage of positive replicates, and the number of ZOTUs among sites (Fig. 6) were consistent with differences detected among sites in community assemblages (Fig. 5). Moreover, families at individual sites with a high abundance of reads tended to have a higher percentage of positive replicates and a greater number of ZOTUs, indicating that these three metrics infer similar patterns. The most apparent difference within each assay, at least for gut contents, was the prevalence of stony coral reads (18Suni = 98.61% or 2,718,582 reads; ITS2 = 86.77% or 565,347 reads) vs non stony coral reads (18Suni = 1.39% or 38,326 reads; ITS2 = 13.23% or 86,220 reads); a difference consistent across all sampling sites. This pattern was less apparent for seawater, noting that the proportion of stony coral reads (18Suni = 0.23% or 3,465 reads ITS2 = 64.33% or 226,792 reads) to non stony coral reads (18Suni = 99.77% or 1,482,882 reads; ITS2 = 35.67% or 125,766 reads) was dependent on assay.

Future directions

In this study, we explored the direct and indirect interactions of a corallivorous butterflyfish with their reef environment using two metabarcoding assays. An expanded treatment of gut content metabarcoding for fishes with different feeding preferences (e.g., planktivores; Leray et al., 2019) or other resident animals with less selective requirements (e.g., filter feeders; Mariani et al., 2019) could therefore be used to gather highly flexible datasets that address multiple ecological questions, particular in habitat that are difficult to survey using visual approaches. For example, highly productive kelp forests and seagrass meadows contain an abundance of organisms that are often small, epiphytic, symbiotic, or cryptic. If interested in a narrower range of taxa corresponding to the food spectrum, metabarcoding assays of the surrounding seawater (Port et al., 2016; Stat et al., 2019; Lamy et al., 2021; Momota, Hosokawa & Komuro, 2022) in tandem with gut contents of the animals that feed directly on them (Gajdzik et al., 2021; Díaz-Abad et al., 2022) can prove to be a more robust approach to validate the resulting detections. Moreover, given the sensitivities of DNA-based sampling approaches, we suspect that in the future, with further improvements to taxonomic assignment or a renewed focus on microbiomes, gut content metabarcoding may be able to estimate fish fitness as it relates to changing environmental conditions and/or habitat degradation (Clever et al., 2022).