Preliminary analysis of New Zealand scampi (Metanephrops challengeri) diet using metabarcoding

Selection of individuals for metabarcoding

The PCR products for sequencing were selected from six individuals (70.2, 70.3, 70.9, Fro1, Fro2 and Fro3), as they provided consistent and strong amplifications for COI and/or 18S genes. The PCR products encompassed the range of methodological factors that were to be addressed (Table 2). Individual 70.9 was used for comparison of the PCRs to determine which DNA polymerase and buffer mix was better (i.e., samples 1 and 2 versus 3 and 4). Preliminary data suggested the Platinum Taq reactions were inferior, so the remainder of the comparisons focussed on the Bioline reactions. The taxa that were identified from the digesta of the foregut and the hindgut were compared for the individuals 70.2 and 70.9 (i.e., samples 1 and 5 versus 2 and 6) to identify if the hindgut digesta could provide useful additional information about the diet, i.e., successive meals may be represented in the hindgut contents. The ethanol and frozen preservation methods were compared using the 18S sequences for three individuals each (i.e., samples 1, 2, 5, 6 and 7 versus 8 and 9), to assess the effects (if any) of differential DNA degradation of the digesta. All individuals were used for the comparison of the database identifications for the COI and 18S genes and thus a preliminary analysis of the diet was evaluated. Potential diet false positives were identified by evaluating the DNA negative.

DNA purification and pooling of selected samples

NucleoSpin^® Gel and PCR Clean-up (Macherey-Nagel, Düren, Germany) was used according to manufacturer’s instructions to purify the PCR products and remove fragments of less than 164 bp in length for all selected individuals. The amplified DNA concentration of those products was determined using Qubit™ dsDNA HS Assay Kit (Invitrogen™, Thermo Fisher Scientific Inc.) following the manufacturer’s protocol.

The purified COI and 18S PCR products from each sample were pooled and brought to equal molarity where possible. Sequencing was done through New Zealand Genomics Ltd (Auckland, New Zealand) at Massey University (Palmerston North, New Zealand) where indexing occurred using the Nextera™ DNA library Prep Kit (Illumina, San Diego, CA, USA) before sequencing on an Illumina MiSeq™ System (2 × 250 pair-end protocol).

Metabarcoding protocol

The raw Illumina sequences were analysed through the New Zealand eScience Infrastructure (NeSI; Auckland, New Zealand) high performance computing (HPC) facility. The resulting Illumina metabarcoding sequenced reads were processed by firstly removing the primers (no mismatch tolerated) using fastq-multx (version 1.3.1) and by pairing the reads with SolexaQA++. Low quality 3′ end sequences (Phred scores below 3) were truncated and reads merged using the sequence analysis tool VSEARCH version 2.3.0 (Rognes et al., 2016), allowing a maximum of five non-matching nucleotides in the overlap region. Merged reads were quality filtered on the expected error value (<1) and dereplicated. Chimera checking and removal was performed on the dereplicated sequences with the QIIME package (Caporaso et al., 2010) based on the de novo and reference based methods of Usearch61 (Edgar, 2010). Reference based chimera detection of 18S was performed with the largest curated database available (SILVA; Quast et al., 2013) and for COI, using the Midori database (Machida et al., 2017). The cleaned reads (Table S4 ) were clustered into OTUs with the Swarm methodology (clustering threshold of d2; Mahé et al., 2014).

Sequence assignment and analysis

Representative OTU sequences (seed sequence of each OTU) were taxonomically assigned with the Ribosomal Database Project (RDP) classifier (minimum confidence level of 80%; Wang et al., 2007) using the curated databases Protist Ribosomal Reference (PR2; Guillou et al., 2013) and SILVA for 18S, and Midori database for COI. Representative sequences were also assigned to reference sequences from the National Center for Biotechnology Information (NCBI) Genbank database (Benson et al., 2013), using the megablast option of BLASTn (Morgulis et al., 2008) with an e-value threshold of 0.001. Only the best hit for each sequence was kept. NCBI Genbank is an all-inclusive database and thus was used for both the 18S gene sequences and COI gene sequences, however, sequences were only assigned to a taxon if the sequence could be assigned to a genus or species. NCBI Genbank is also not a curated database. The Midori database is specifically for metazoan mitochondrial DNA sequences. The PR2 database consists mainly of nuclear-encoded protistan sequences, and SILVA is an aligned database of small and large subunit rRNA genes. Taxonomically assigned OTUs were filtered for potential PCR or sequencing artefacts by selecting OTUs with five or more sequence reads (“hits”). The resulting filtered OTUs for COI were then further refined to keep only the OTUs that could be identified to either a genus or species level by Midori and the NCBI databases. The same was done for 18S with PR2, SILVA and NCBI. The resulting data for genus or species were used for assessing the methodological factors and preliminary diet which was analysed using Rstudio^® (R Core Team, 2017) and the collection of R packages in tidyverse (Wickham, 2017).

To determine a list of identified taxa that were most likely to form part of the true diet of M. challengeri, taxa identified as contamination were removed from the final data. Any sequence reads identified as lobster (Astacidea) were presumed to be host contamination and removed. Sequence reads were also removed from further analysis if they matched to any taxa identified from the sequenced DNA negative sample, or to any taxa that could not be part of the M. challengeri diet (e.g., terrestrial species). The final “hit counts” pertaining to digesta are referred to as “diet hit counts”.

Taxa were defined as “exclusive” if they were identified in only a single source, such as taxa identified in the foregut digesta but not in the hindgut digesta, or solely identified in one PCR but not another.

PCR reagents (DNA polymerase and buffer)

A comparison between the Platinum Taq and the Bioline PCR reagents was made for individual 70.9, for which both sets of reagents were used (Table 2). Metabarcode sequences from both foregut and hindgut digesta were pooled for each set of reagents (thus comparing samples 1 and 2 with 3 and 4). The total number of taxa identified and the total number of exclusive taxa identified (i.e., those found only in one category) was greater for both the COI and 18S genes when using the Bioline reaction (11 different taxa in total and six exclusive taxa), compared to the Platinum Taq reaction (six different taxa identified and one exclusive taxon; Table 3). The total diet hit count (number of sequence reads matching a potential dietary reference sequence in that database—i.e., excluding host matches) was also greater for the Bioline reaction.

Five of the taxa identified by Platinum Taq were a subset of those found in Bioline, which had a low summed hit count of 55 compared to Bioline’s 468 hit count for the same taxa. Moreover, Bioline identified a higher number of exclusive taxa which averaged a hit count of ∼32 hits per taxa whereas Platinum Taq only identified one exclusive taxon with a hit count of 13.

Foregut versus hindgut digesta

The taxa identified from the hindgut and foregut digesta of two individuals were compared (i.e., samples 1 and 5 versus 2 and 6). The hindgut digesta contained a higher number of both total taxa and exclusive taxa compared to the foregut digesta, when the sequences were matched to the COI gene databases (Table 4). Conversely, the foregut digesta contained a higher number of total taxa and exclusive taxa with the 18S gene databases (Table 4).

Sample preservation method

The effects of potential differential DNA degradation due to the method of preservation were more difficult to assess directly, as the same individuals could not be compared across the two conditions, and the comparison was somewhat confounded by being collected in different locations. Instead, broad comparisons were made between the pooled results from the three individuals preserved for each of the ethanol and frozen methods. Firstly, the potential differential effects on preservation of diet DNA was assessed by calculating the diet hit count as a percentage of the overall hit count (for samples 1, 2, 5, 6 and 7 versus 8 and 9; Table 5). No difference was seen between the ethanol and frozen individuals for the 18S results, as both had a diet hit count percentage of 1% (Fro3 was the only frozen individual to have any diet hit counts; Table S6).

PCR inhibition

PCR inhibition is a substantial factor in this study. The spectrophotometer ratios of purity provided initial indicators of the presence of inhibitors in the extracted DNA (Watts, 2014). An optimal volume of 2 µl 1% BSA successfully reduced PCR inhibition by increasing the amplification efficiency. An additional step effective in reducing PCR inhibition, was 1:10 dilution of the template genomic DNA, thus diluting the PCR inhibitors it presumably contained.

The presence of organic compounds in the foregut (i.e., bile salts in the gastric juices) is thought to be one of the main contributors to PCR inhibition which cannot always be alleviated by BSA (Lorenz, 2012; Schrader et al., 2012; Harms-Tuohy, Schizas & Appeldoorn, 2016). Bile salts often found in vertebrate gastric juices, are similar to an organic compound in crustaceans such as crabs, lobsters and crayfish, and suggests why BSA is not an alleviator in all DNA amplifications tested in this study (Borgström, 1974). Different diets would also affect the secretion of the gastric juices, which may play a role in determining which PCR inhibitors are present and in what concentration in different parts of the gut (Rotllant et al., 2014). This demonstrates the problematic, high levels of PCR inhibition present in many genomic DNA extracts, and its influence on successfully amplifying gut content DNA. Increased quantitative monitoring of the PCR success, such as by utilising quantitative PCR, may be required in future.

One solution to PCR inhibition may be to pool replicate PCR products together. Another is the addition of the T4 gene 32 protein (gp32), which can be effective against inhibitors for which BSA is not, such as sodium dodecyl sulfate (SDS), sodium chloride and bile salts (Schrader et al., 2012).

PCR reagents (DNA polymerase and buffer)

Different DNA polymerases and buffers are available and are two of the key reagents for DNA amplification, each being susceptible to various components in the biological environment in which they amplify, thus affecting their performance (Rådström et al., 2004; Wolffs et al., 2004). A substantial difference could be seen visually in the DNA amplification when using the Bioline reaction compared to the Platinum Taq reaction. The Platinum Taq reaction appeared to be more susceptible to the effects of PCR inhibitors, and the DNA amplification was not as effective as in the Bioline reaction.

Amplification of the loci

The 18S gene region (∼425 bp), in general, amplified more successfully compared to the COI gene region (∼313 bp) targeted here. It is not common for shorter gene fragments to be more difficult to amplify, but it is relatively common to see some difficulties in amplification of the COI gene in some species, likely due to the greater sequence variation in the primer annealing sites (Chen, Jiang & Qiao, 2012; Lv et al., 2014). However, when both genes amplified, their PCR product intensity appeared consistent between sample conditions, i.e., if the PCR product was more intense for the foregut digesta for the 18S gene the same would be seen for the COI gene in the same individual.

PCR reagents (DNA polymerase and buffer)

Visualization of the PCR amplification can only suggest if amplification is occurring and not necessarily what the quality of amplification is. So although it was visually determined that the Bioline reaction generally resulted in stronger PCR amplifications, it was possible that it was amplifying only a subset of the possible templates compared to that of the Platinum Taq reaction. Individual 70.9 amplified equally well for both the Platinum Taq and Bioline reactions, and also for both the hindgut and foregut digesta. The DNA sequencing revealed that the Bioline reaction amplified a higher number of taxa from the same individual’s gut content DNA.

This was likely due to a higher tolerance of PCR inhibitors from the DNA polymerase and buffer in Bioline. This is suggested by the shared taxa, as Bioline (hit count 468) amplified approximately nine times more diet fragments than Platinum Taq (hit count 55). The six exclusive taxa, presumed to have the lower DNA concentrations, identified by Bioline have ∼32 hit counts per exclusive taxon, whereas, Platinum Taq only identified one exclusive taxon with a hit count of 13, thus indicating that the Bioline reaction appears to be more efficient at amplifying DNA from the digesta of M. challengeri.

Foregut versus hindgut digesta

The examination of both the foregut and hindgut digesta provided a broader overview of the diet of M. challengeri. There was an increase in the number of taxa when the entire gut digesta was used, compared to using digesta from only the hindgut or the foregut. As ingested food takes a few hours to pass through the foregut and into the hindgut, the examination of the digesta from both may permit the assessment of successive meals over a broader time period before capture, or differences in gut transit times for different dietary items (Simon & Jeffs, 2008; Lee, Hartstein & Jeffs, 2015; Kamio et al., 2016). Also, as the chemical environment (i.e., pH) and degree of degradation of digesta differ between the gut sections, an assessment of both foregut and hindgut digesta appears to maximise the opportunity to detect the widest possible breadth of diet per individual.

It was initially expected that the hindgut contents may provide lower quality sequencing results, because it was more difficult to eliminate host tissue contamination during dissection, the more advanced state of degradation of the food, and the higher ratio of sediment to dietary items. However, a higher number of taxa of smaller marine organisms (i.e., dinoflagellates) were identified in the foregut digesta than in the hindgut digesta, and a higher number of taxa of larger marine organisms (i.e., mostly fish) were identified in the hindgut digesta than in the foregut digesta. Thus, to dismiss sampling from either the foregut or the hindgut digesta would have narrowed the number of taxa that were able to be identified.

Host contamination

In most metabarcoding studies of diet in predators, host contamination of the PCR amplicons is an issue. However, the overall impact of host contamination can be reduced by introducing several molecular techniques that will help decrease any biased amplification of the host’s DNA, and may broaden the taxonomic analysis (Polz & Cavanaugh, 1998; Green & Minz, 2005; Vestheim & Jarman, 2008; Pompanon et al., 2012; Harms-Tuohy, Schizas & Appeldoorn, 2016; Devloo-Delva et al., 2018). PCR “clamping”, using additional DNA or peptide nucleic acid (PNA) oligonucleotide primers to bind to, and inhibit PCR amplification of host DNA, is a favoured method for suppressing host DNA amplification in PCRs, and has been successful in diet studies (Egholm et al., 1993; Orum et al., 1993; Vestheim & Jarman, 2008; Chow et al., 2011; O’Rorke et al., 2012). In this study, PCR clamping was not used, as previous research in lobsters had shown that it also restricted amplification from related crustaceans, many of which could be potential prey of M. challengeri. Other diet studies have also shown that the resulting sequences are sufficient to describe the diet without the use of PCR clamping (Piñol et al., 2014) and that there would not necessarily be any significant increase in diversity if they were used (Devloo-Delva et al., 2018). The metabarcoding results showed that host DNA amplification did not completely outcompete the amplification of prey DNA, and there were sufficient sequence reads to identify a broad range of diet taxa.

Preservation methods

The quality of DNA extracted from the digesta is affected by the state of degradation due to the digestion rate in the gut (Rådström et al., 2004; Deagle & Tollit, 2007; Troedsson et al., 2009). As the two different preservation methods used would have different rates in which the digesta was preserved, it was thought that one method may provide better preservation of digesta than the other. Determining which preservation method, frozen or ethanol, is better was relatively problematic to assess, due to the relatively high level of inter-individual variation and the confounded nature of the sample collection. Due to logistic constraints on the field sampling, different preservation methods could only be undertaken on different individuals taken from different locations.

The diet hit count percentage from the 18S databases were both ∼1% for ethanol and frozen individuals. Although not a particularly robust test due to sample size, it may suggest that both the frozen and ethanol preservation methods are similarly dependable.

Amplification of the loci and database review

Although there were difficulties in gaining equal success in the DNA amplification from both genes, together they provided a much broader overview of M. challengeri diet at the highest resolution possible.

The 18S gene and the COI gene identified a variety of taxa (Figs. 1 and 2). The COI gene is more commonly used than the 18S gene in studies of molecular taxonomy of large eukaryotes and is one of the pillars of DNA barcoding, making it a sought-after gene region for the discrimination of closely related Metazoan species (Hebert, Ratnasingham & DeWaard, 2003; Derycke et al., 2010; Tang et al., 2012). The 18S gene is more frequently used in molecular taxonomic studies for microbial organisms (i.e., micro-plankton and parasites; Wu, Xiong & Yu, 2015).

Using a combination of databases to identify the OTUs, from a selected gene, to genus or species level proved invaluable; where the curated databases (Midori, PR2 and SILVA) could not provide a high resolution identification (i.e., only to order or family), NCBI was usually able to. Also, targeting two genes has allowed for a more inclusive capture of the different taxonomic groups because of the more complete coverage from distinct reference databases. For further analyses of the diet of M. challengeri both genes should be used, as this lobster, being relatively small in size, appears to take advantage of the availability of both large and small prey, both dead and alive.