Quantifying wildlife conflicts with metabarcoding and traditional dietary analyses: applied to seabird predation by long-nosed fur seals

Introduction

New wildlife conservation and management scenarios are emerging during the Anthropocene as some species experience population increases through successful conservation efforts, while others continue to decline (Roman et al., 2015; Marshall et al., 2016; Cammen et al., 2019). Complex management scenarios arise when species recovery results in negative interactions with other species of value, whether that reflects a trophic role in the ecosystem, conservation status, community connection or economic opportunity (Marshall et al., 2016). For example between killer whales, sea otters and salmon (Estes et al., 1998; Williams et al., 2011); New Zealand sea lions and yellow-eyed penguins (Lalas et al., 2007); wolves and caribou (Hervieux et al., 2014); jaguar and marine turtles (Veríssimo et al., 2012); and, between Critically Endangered yellow-breasted capuchin monkeys and Endangered red-billed curassow (Canale and Bernardo, 2016). The interactions are natural but present a need for accurate information on natural predation levels and impacts to prey (Hervieux et al., 2014) to avoid unjustified persecution of the predator (Granquist et al., 2018), and for effective management of all species concerned (Marshall et al., 2016).

A key goal of investigating predator-prey interactions involves determining prey inter- and intra-specific diversity, dietary proportions, and abundances or biomass consumed by the predator (reviewed by Pompanon et al., 2012), while achieving reliable detection of predation (Deagle et al., 2005). Developments in eDNA extraction and metabarcoding techniques are demonstrating reliability for this level of environmental monitoring by: (i) identifying prey at high taxonomic resolution and when missed by other methods (Bowen and Iverson, 2013; Stat et al., 2019); (ii) estimating dietary proportions from relative genetic importance, and reconstructing approximate biomass consumed through genetic abundance data (Deagle et al., 2019); (iii) identifying species’ intraspecific genetic diversity within environmental samples for wildlife forensics and population estimation (Sigsgaard et al., 2016; Seersholm et al., 2018; Tsuji et al., 2020).

In southeastern Australia, the recovery of long-nosed fur seals (Arctocephalus forsteri) could cumulatively threaten populations of little penguins (Eudyptula minor) in addition to other known stressors (Reinhold et al., 2022). Both species are federally protected and garner significant cultural and conservation value (Introduced by Bill Heffernan, Ian Campbell., 1975 and Introduced by Bill Heffernan, Ian Campbell., 1999), albeit listed as ‘Least Concern’ by the IUCN Redlist (IUCN, 2020). Little penguins are a popular tourist attraction and valued species to communities across southern Australia (Tisdell and Wilson, 2012). Little penguin populations are estimated at 470,000 individuals (BirdLife International, 2021); yet 60% of sites have unknown population trends, 29% of colonies are declining, and most persist on offshore islands in southern Australia and are difficult to census (BirdLife International, 2021). Little penguins forage at sea but breed on land, making them susceptible to threats from both environments. Major contributors to declines in little penguin populations include terrestrial threats of: (i) changes in land-use and predators introduced by European settlers (Dann, 1991; Kirkwood et al., 2014) and (ii) susceptibility to hyperthermia during more frequent terrestrial heat waves (Tworkowski et al. (2018); Marker, PhD Thesis, 2016; also seen in African penguins, Welman and Pichegru, 2023). At sea, ocean warming and large-scale changes to food webs caused by climate change as well as competition with fisheries are threatening processes of concern (Ropert-Coudert et al., 2019).

The fur seals were decimated by fur trade through the 1800’s and culling into the late 1900’s due to perceived competition for resources with fishers (Shaughnessy, 1999). Long-nosed fur seals have experienced range-wide increasing population trends, reported at 97,200 in the state of South Australia (2013–14 census; Shaughnessy et al., 2015) where an estimated 83% of their recorded pup production occurs. While the original population size is unknown, harvesting records suggest that the current population represents a small fraction of that prior to European colonization and exploitation (Ling, 2014). The recovery of the species in Australia continues to conflict with some communities and fisheries (Shaughnessy et al., 2003; Goldsworthy and Page, 2007; Cummings et al., 2019) due to their potential predation impacts to prey species. While it is realistic to expect high levels of predation to affect prey population size or behavior (Visser et al., 2008), recurring, popular discussion to cull the long-nosed fur seal population in South Australia have grown despite an absence of quantitative information (Goldsworthy et al., 2019).

Little penguins and other seabirds were identified in juvenile, sub-adult, and adult male long-nosed fur seal diets, at two locations in southern Australia and at low occurrence frequencies (Page et al., 2005; Hardy et al., 2017; Goldsworthy et al., 2019). Across several sites in Bass Strait and Eastern Victoria, long-nosed fur seals were more likely to predate little penguins originating from sites in closest proximity to them (Reinhold et al., 2022). However, beyond the presence of diagnostic remains in scat samples, we fundamentally lacked any further quantitative information on this predator-prey interaction, particularly for ‘unquantifiable remains’, such as feathers and when prey are not consumed whole. Page et al. (2005) proposed a single scat containing feathers was equivalent to a single bird consumed; however, mass-balanced trophodynamic modelling in the Great Australian Bight marine ecosystem suggested that this likely overestimated predation and that multiple scats containing feathers could also represent a single bird consumed (Goldsworthy et al., 2013). This issue affects all predation systems where prey are not consumed whole, or where diagnostic hard parts do not survive digestion, and where they are lacking in the first place for soft-bodied animals; yet traditional morphological analyses of stomach or gut contents remains the standard in predator diet analyses.

To assess predation incidence on little penguins and seabirds across southeastern Australia, we apply two complementary surveillance techniques simultaneously, morphometric (prey hard-part) and avian-specific DNA metabarcoding (genetic) assays, to long-nosed fur seal scats. We consider these techniques complementary rather than overlapping, in that they are designed to detect evidence of predation from different sources of tissues when applied simultaneously to bulk sample processing, namely soft vs. hard remains of prey which are known to be differentially affected by digestion (Casper et al., 2007; Tollit et al., 2009; Jeanniard-du-Dot et al., 2017). Of note, it is important to understand that when hard parts are present in samples there are known differences in the quantities of DNA that can be obtained from chitinous compared to soft tissues (McDonald and Griffith, 2011). Therefore, (i) we compare overall seabird and little penguin detection rates across groups of samples, using hard-part and genetic assays; (ii) we investigate the diversity and relative importance of seabirds consumed by long-nosed fur seals apparent from both assays; and (iii) with the output of the genetic assay, we explore a minimum estimate of penguin abundance consumed by analyzing mitochondrial haplotype diversity among little penguin DNA obtained.

Methods

Key terminology & definitions

The goal of this study is to provide fully leveraged information from diet analyses which practitioners can use to make conservation and management decisions. To this end, it is important to clarify several terms and metrics reported in this study, for ease of understanding of the methods and results and their recommendations.

Regarding assays used, we refer to prey “hard parts” when describing their morphological remains, and to “hard-part” or “morphological” analyses when describing the assay technique that examines their morphological remains (see Jeanniard-du-Dot et al., 2017). We refer to “DNA-based” or “genetic” analyses when describing the assay technique that examines prey DNA extracted and sequenced from samples.

While incidence, frequency and occurrence are synonyms in the English language, we use “predation incidence” when discussing the impacts of our findings, because seabird predation by long-nosed fur seals to date has been incidental and infrequent (Page et al., 2005; Hardy et al., 2017; Goldsworthy et al., 2019). We use the frequency at which seabirds and little penguins were detected from two assay techniques to calculate the “frequency of occurrence” of these taxa across the range of samples to inform our reported range of predation incidence. There is concern that to refer to predation “frequency” implies that it is “frequent”, and we emphasize that this distinction requires additional monitoring of the prevalence of the patterns identified in this study.

As commonly practiced in diet studies, we report the “detection” of seabirds and little penguin remains when comparing the results of the genetic and morphological assays, and “occurrence” in comparison to other covariates (i.e., location) and studies. We note that diet studies commonly calculate frequency of occurrence based on the number of samples in which prey taxa were detected, and all studies understand that this is simply a representation of the true frequency of occurrence or incidence of predation. Finally, in addition to occurrence information, we report the relative and total abundances of prey sequences (Supplementary Information). This is because experiments performed on seals have shown useful relationships between the biomass of prey consumed and the relative sequence abundance, as well as the level of importance of prey species (Deagle et al., 2005; Thomas et al., 2014, 2016; Deagle et al., 2019), and highlight the potential for this information to be useful in the future, as further experiments enable the development of correction factors.

Sample collection

Individual predator scat samples (n = 99) were collected across multiple time points from four long-nosed fur seal breeding colonies in Victorian Bass Strait and New South Wales (NSW), southeastern Australia (Figures 1, 2). Pup abundances are illustrated as a conventional proxy for relative seal population abundance (Figure 1; Appendix S1.1). Most samples were collected from the two larger colonies, Barunguba and Cape Bridgewater, in the Austral spring (September) 2016 and summer (January) 2017, with additional samples included from spring 2015 and summer 2016 at Cape Bridgewater. Samples from Gabo Island were collected from summer 2017. One sample was opportunistically collected from a lactating female at Deen Maar Island and included in assays. Sample sizes resulted from balancing adequate replication per site with availability of fresh samples.

Figure 1 (A) Long-nosed fur seal scat collection sites (n = total number of samples), including pup abundances for sampling locations (McIntosh et al., 2014) as an index of seal population relative importance of sites. Sampled sites were Cape Bridgewater (38.3013° S, 141.4062° E) and nearby Deen Maar Island (formerly Lady Julia Percy Island, 38.4161° S, 142.0038° E) from western Bass Strait, Victoria; Gabo Island in eastern Bass Strait, Victoria (37.5649° S, 149.9133° E); and Barunguba (formerly known as Montague Island 36.2510° S, 150.2270° E) at the northeastern breeding range in New South Wales (NSW). Species distributions shown for (B) long-nosed fur seals and (C) little penguins (data from ALA, 2019).

Figure 2 (A) A long-nosed fur seal, Arctocephalus forsteri, from Barunguba, NSW; (B) the little penguin, Eudyptula minor, often burrowing near fur seal colonies; (C, D) seabird remains are conspicuous at long-nosed fur seal haul-outs and colonies, among scats and regurgitates.

Whole and moist scats were sampled to minimize bias from differential DNA degradation or partial loss of material (similar to Deagle et al., 2009). Whole scats were thoroughly mixed with individual disposable spatulas at point of collection and a 2 mL subsample was taken from each field-homogenized scat for genetic analyses of prey tissues (Hardy et al., 2017). The remaining whole scats were collected for analyses of morphological prey remains, using individual, zip-lock bags. Samples were stored within hours of collection between -10° and -20°C in portable freezers (WAECO) for up to 7 d in the field and transferred to -20°C freezer facilities.

Identification of seabird morphological and genetic remains

In total, 99 fecal DNA sample extracts were processed using morphological and genetic assays. For the morphological assay, all prey items recovered from scat samples were identified from hard parts using the methods described by Page et al. (2005). Birds were identified using feathers and other remains such as feet, flippers, and heads (Figure 2, Appendix S1.2).

For the genetic assay, prey DNA extractions used 250 mg of fecal subsamples and MoBio PowerSoil^® DNA Isolation Kits (now QIAGEN: www.qiagen.com) effective for DNA extraction from highly inhibited and mixed samples of fecal origin (similar to Carroll et al., 2019), with modifications to the manufacturer’s instructions made to optimize DNA extraction. These included an overnight digestion phase in cell lysis buffer (C2 solution) at 4°C, and the incubation step in inhibitor removal solution was extended from 5 to 60 minutes at 4°C. No host inhibitor step was required because the assays used do not detect mammals. Target DNA was then eluted in 100µL of 10 mM Tris buffer, MoBio PowerSoil^® C6 solution, (www.qiagen.com) and stored at -20°C. DNA extract concentrations were measured and verified using a NanoDrop™ Spectrophotometer (www.thermofisher.com/). For use in positive controls and to test primer specificity, nuclear DNA was extracted from the center of the muscle tissue matrix (25 mg) of a domestic chicken (Gallus gallus domesticus) and a little penguin carcass obtained by Phillip Island Nature Parks, using Bioline Isolate II Genomic DNA Kits (www.bioline.com/us/) as per manufacturer instructions. A dedicated controlled eDNA laboratory was used at RMIT University, Bundoora, Victoria, with separate spaces and rooms designated for the physical separation of eDNA extraction, pre-PCR preparations and post-PCR procedures. Positive and negative controls (extraction and PCR) were used to identify potential contamination at each laboratory procedural step from DNA extraction to diagnostic PCR steps.

The 99 fecal DNA sample extracts were screened in duplicate and at two DNA concentrations (neat and 1:10 dilutions), alongside extraction blanks (n = 5), PCR blanks (n = 4), and positive controls (n = 2) by diagnostic endpoint PCR (dPCR) using the Bird12sa/h assay (forward 5’ CTGGGATTAGATACCCCACTAT to 3’, reverse 5’ CCTTGACCTGTCTTGTTAGC to 3’), a conservative primer ‘Bird12sa/h’ targeting a ~230 base pair (bp) fragment of the avian 12S ribosomal RNA (rRNA) gene (Cooper, 1994) (Table 1). PCRs were run on Bio Rad C1000 Touch thermal cycler using cycling steps outlined in Table S2, and using the AmpliTaq Gold^® 360 Master Mix using reagents and concentrations provided by the manufacturers. All duplicate dPCR products were run on 1.5% agarose gels to determine the presence/absence of amplified target bird DNA. We obtained optimal amplification and low inhibition from neat DNA concentrations.

Table 1 Primers used, target taxa and genes, as well as designer references.

Of note, we conducted an initial pilot study using the Bird12sa/h and a second previously tested assay also for mitochondrial DNA, the AWCF1/R6 primer pair targets a longer fragment (~850 bp) of the cytochrome c oxidase 1 gene (COI) (Patel et al., 2010) (Table-S1, S2, described in Appendix S1.3). Both primers were tested using both positive controls (little penguin and chicken DNA) and run on a gradient PCR from 52–58°C (with 1°C increments) to optimize the annealing temperature. We selected 57°C for the Bird12sa/h primer (Cooper, 1994), and 54°C for AWCF1/R6 (Patel et al., 2010) as the optimal temperature for PCR for these primers and to further test their amplification success rates on a subset of our fecal DNA samples.

The Bird12sa/h assay targeting a shorter gene fragment was more sensitive than the COI primer targeting a long fragment environmental DNA, which was expected given the degraded nature of DNA found in fecal samples. Additionally, the Bird12sa/h primer produced a similar proportion of positive results for seabird eDNA to the overall detection rates of diagnostic seabird hard-parts in samples (~30%), while detection rates using COI primer were very low (~5% of samples) in our samples. Consequently, due to time constraints and costs, we chose to continue with only the Bird 12sa/h assay and discuss the availability and use of multiple specific genetic assaying tools in avian eDNA. Finally, because both primers are avian specific, there was no need to inhibit the host DNA (Cooper, 1994; Berry et al., 2017). We share this information to assist future studies in primer selection and optimization.

Using the avian specific Bird12sa/h assay, a total of 32 samples (of 99) showed target amplicons in both or a single duplicate at neat DNA concentration, all extraction and PCR controls were negative. DNA extracts of the 32 samples that tested positive for birds, and two extraction blanks and one positive control (n = 35 samples for sequencing) were therefore sent for quantitative PCR (qPCR), cleanup, sample-based rarefaction and extrapolation sampling curves, appropriate sequencing depth (< 10,000 reads per sample) and next generation sequencing performed on Illumina Miseq by Ramaciotti Centre for Genomics (RCG), University of New South Wales. There, a single-step fusion tagging PCR procedure was used to attach and assign unique MID (Multiplex IDentifier) tag combinations, next generation sequencing (NGS) adaptors and the Bird12sa/h assay. Amplicons were purified and blended at RCG in equimolar concentrations to form a library, which was sequenced with a 150 bp paired-end sequencing kit (Illumina Miseq v2 Nano 150 bp). We used the single-step fusion PCR procedure over the two-step PCR approach to reduce the risk of ‘tag jumping’ during the second amplification step where MID tags are assigned (Schnell et al., 2015; Taberlet et al., 2018). This type of error is difficult to detect and risks cross-contamination of amplified DNA among samples between initial PCR products and terminal PCR products. Single-step fusion PCR procedures therefore provide us with the least risk of sample cross-contamination over other procedures. After sequencing, samples were ‘demultiplexed’ and assigned to the correct original sample by their individual MID tags.

We used Geneious R8.1.5 (Kearse et al., 2012) to merge the paired-end forward and reverse sequences (2x ~150 bp fragments, with overlap of 70 bp) and retain only those with exact flanking MID tags, primers, and adapter sequences. Once paired, the MID tags, NGS adaptor sequences and the Bird12sa/h forward and reverse primers, were subsequently trimmed, leaving the complete target sequences for each sample. Sequences were discarded if they did not contain exact matches to both the forward and reverse PCR primers, tags, and adaptor sequences, failed to pair, or were > 10% shorter than the primer product length (expected 220 bp, discarded below 200 bp) (as in Berry et al., 2017, and Hardy et al., 2017).

These sequences were quality filtered and clustered into molecular operational taxonomic units (OTUs) using the UPARSE algorithm in USEARCH (Edgar, 2010; Edgar and Flyvbjerg, 2015) and using a 97% similarity criterion (as in Berry et al., 2017). Illumina’s Miseq has been found to have an error rate of about 0.1% (Fox et al., 2014), we chose a conservative 1% cut off of aggregated 100% identical, unique nucleotide sequences (hereafter ‘unique sequences’) to further minimizes the risk of erroneous sequences and false positives. We then mapped total filtered sequences (hereafter retained reads or sequences) for these identified OTUs back to individual samples. This can result in some samples containing very low to trace amounts of retained target sequences of high quality that belong to a more abundant cluster from the pool across samples.

Consensus, representative sequences for each out were queried against the National Center for Biotechnology Information’s (NCBI) GenBank nucleotide database using the algorithm BLASTn (Basic Local Alignment Search Tool; Benson et al., 2005). The resulting queried sequences were assigned to taxa, following criteria and taxonomic reference databases outlined in Hardy et al. (2017) and Appendix S1.4 (Table S3). These criteria maximized confidence in making a taxonomic identification by remaining conservative in our assignments (i.e., selecting an identification at genus level) where multiple species were found to be closely related on the Bird12s gene, or where other likely and related prey lacked genetic reference material, which could lead to false assignment to a genetically related taxon with representative genetic reference material. All the identified seabirds occurred within the geographic ranges of the LNFS and are considered viable prey species for LNFS.

Haplotype polymorphism analysis

As a species of conservation concern and key avian prey species, we sought to identify a minimum number of individual little penguins by exploring mtDNA haplotypes from 12S rRNA sequences obtained for little penguins. While dependent on sequence fidelity, such approaches have been used to explore intraspecific diversity in genetic data (similar to Seersholm et al., 2018). We imported the retained sequences for all seabirds, produced just prior to OTU clustering, we matched these sequences in relation to sample identifiers and formed these into clusters of unique sequences in Geneious (Kearse et al., 2012).

From each of the penguin DNA-abundant samples (Table S4), we conservatively selected the top 1–3 most abundant sequences to form the basis a conservative set of haplotypes identified in each sample and each represented > 5% of the total sequence abundance within each sample, and each diverging by 2–8 bases. Only one haplotype diverged from another by a single base, due to sequence abundance this would not have occurred due to sequencing error, and thus represents a truly divergent sequence cluster. We chose this conservative and manual selection method of haplotype identification over an in-depth sequence denoising method that would have maximized the identification of true haplotypes (see Nearing et al., 2018), due to the small sample size of penguin-positive and abundant penguin-DNA containing samples. We present a conservative number of individual haplotypes, while it is possible that additional true haplotypes exist within the data and were missed using this approach.

We produced a minimum spanning haplotype network using the software PopART (Leigh and Bryant, 2015) from an alignment of these identified penguin sequences to visualize relationships between haplotypes consumed, and between samples from the different locations sampled. All samples containing little penguin DNA, including lower abundance samples, were then searched for the presence of these five identified haplotypes. We estimated the number of penguins likely consumed based on the number of haplotypes within each sample. As samples were collected across multiple days from each location and sampling time, we treated each sample as from distinct predators or predation events.

Statistical analyses

To account for different sampling times and locations, samples were assigned one of seven unique grouping factors that combined location and time (e.g., Barunguba, January 2017). The single scat sample collected from Deen Maar Island was not included in statistical tests, but results were reported for future comparisons. The presence of seabird and specifically little penguin remains across samples using different dietary assay techniques (hard-part vs. DNA; Table S6, Appendix S2) were examined using two generalized linear models (GLMs), constructed in the stats package in R version 4.0.3 (R Core Team, 2020). A binomial error distribution for presence-absence data was used and additive term included to account for long-nosed fur seal group (location and time). Model fit was assessed using deviance explained and variable significance (p < 0.05).

Results

Bioinformatics results

From initial bioinformatic processing steps (sequence demultiplexing, read pairing and trimming barcodes), we obtained 64,361 individual sequences and 7,370 unique sequences. Following sequence quality filtering to 1% (n < 74) of the total abundance of unique sequences (n = 7,370) and clustering to OTUs (97% similarity). This filtering and quality control procedure resulted in 35,498 retained sequences, and 5 OTUs (of the thousands of unique sequences that were combined to within 97% similarity), for 25 samples.

Comparing diagnostic hard-part and genetic assays

Overall, the detection rates of seabirds and specifically penguins were statistically similar for each method (GLM seabird detection ~ metric: p-value = 0.6483; GLM penguin detection ~ metric: p-value = 0.2003; Table S6). The assays tended to detect seabird remains in different samples, with an overlap in detection of the same taxon using both assays for a third of the seabird positive samples, which was to be expected when using assays that target different tissue types (soft vs. hard parts) and with different passage times through predation digestion. Notably, DNA metabarcoding offered additional information: (i) in detecting multiple prey taxa within a single scat that never occurred with the hard-part methods (Figure 3), (ii) total and relative abundance information for amounts of DNA recovered (Figures S1, 4, Table S4), and (iii) exploration of within-sample genetic diversity assessed here for little penguin 12S rRNA enabling the estimation of a conservative number of predated penguins (Figure 5, Table S5).

Figure 3 The diversity of seabirds identified in long-nosed fur seal samples: (A) using hard-part analyses (n = 29) and (B) using DNA metabarcoding (n = 21). GI, Gabo Island.

Figure 4 Relative genetic contribution within seabird-positive long-nosed fur seal samples (n = 21) of seabird taxa identified and illustrated as a proportion of DNA abundance recovered for each taxon within each sample, and for locations and sampling times.

Seabirds were detected in 29.3% (n = 29) of samples using diagnostic hard-parts, and 21.2% (n = 21) of samples using DNA metabarcoding (Figure 6A). Most of these detections were little penguins, detected in 25.3% (n = 25) of samples with hard-parts and 10.1% (n = 10) of samples with DNA (Figure 6B). These detection rates, when calculated across all samples amount to reporting frequency of occurrence, provide the upper estimate of predation incidence in long-nosed fur seal samples for seabirds and little penguins. All positive detections were obtained after quality control and filtering; despite this, the majority of DNA sequences for all seabird taxa were identified in 9 out of 21 seabird-positive samples (Figures 6, S1, Table S4). Additionally out of the 21 seabird-positive samples, 10 were penguin-positive, while 6 contained the majority of the penguin DNA (Figures 6, S1, Table S4). This provided a lower and more conservative estimate for predation of seabirds in 9.1% and little penguins in 6.1%, across 99 samples.

Figure 5 Little penguin genetic diversity detected from sequences from the avian specific Bird12sa/h assay (for ~230 bp 12S rRNA gene): (A) presented as a minimum spanning network of five distinct haplotypes, and (B) number of haplotypes contained within each penguin-positive sample, including haplotype sequence abundances within samples.

While detection rates were statistically similar for both methods for seabirds and penguins, there was greater variability across locations using hard-parts compared to DNA (Table S6, Figures S2, S3). Minor, albeit statistically significant, differences were observed across sampling groups for seabird occurrences (binomial GLM for seabird occurrence ~ location: p-value = 0.017; Figure S3A, Table S6), but not for penguins (binomial GLM for penguin occurrence ~ location: p-value = 0.2291; Figure S3B, Table S6). Occurrence variability was high across sampling locations, times and methods used (Figures S2, S3). Little penguins were detected by the hard-part assay from Gabo Island samples, confirming predation occurs at this location, however only a few sequences of little penguin DNA were detected at Gabo Island and Deen Maar Island, and these sequences did not pass DNA quality filtering.

Seabird diversity in long-nosed fur seal diets

DNA metabarcoding resulted in the detection of taxonomic mixtures in scat samples compared to hard-part analysis (Figure 3), by detecting 2 distinct seabird taxa in 5 of the samples. Using hard-part analysis, no sample was found to contain more than one identified bird taxon (Figure 3A). Little penguins were the main seabird prey species detected using both analyses (Figure 3), and this was reflected in all datasets obtained from DNA (Figure 3), with greater total abundance of sequences (Figure S1) and greater relative sequence abundance (Figure 4) than any other seabird prey identified.

Figure 6 Detections of (A) seabird and (B) little penguin diagnostic hard-parts (‘HP’) and DNA (‘DNA’), as a percentage of all long-nosed fur seal samples (n = 99). We report all genetic sequences obtained from standard sequence quality control and filtering ‘DNA (all)’, and for samples that contained large sequence quantities of sequences ‘DNA (abundant)’ (> 90% of filtered sequences); and the number of samples that contained both the morphological and genetic remains of the same seabird (‘same taxon’).

Hard-part analysis revealed two other taxa in addition to little penguins: a shearwater family group (Procellaridae spp., n = 2), and Australasian gannet (Morus serrator, n = 1) (Figure 3). DNA metabarcoding detected abundant DNA from two distinct shearwater taxa, also identified at family-level, Procellarid spp. 1 in 9.1% (n = 9) and Procellarid spp. 2 in 5% (n = 5) and of samples (Tables S3, S4, Appendix S2). Black-browed albatross (Thalassarche melanophris) and greater crested tern (Sterna bergii) were each detected once (Figure 3, Tables S3, S4). Importantly, parallel use of both genetic and hard-part analyses revealed a greater diversity of taxa than either method alone.

Towards quantifying little penguin consumption

From the six scats containing abundant penguin DNA (Table S4), a total of five mtDNA haplotypes were identified in samples from Cape Bridgewater (haplotypes 1–2 and 4–5) and Barunguba (haplotypes 1–3) (Figure 5A). From each of the penguin DNA-abundant samples, we selected the most abundant unique sequences to form the basis of the first two haplotypes (haplotypes 1 & 2). We then examined the abundance of the remaining sequences in each sample. Within the Sept16:5 sample there were two potential haplotypes (haplotypes 4 & 5), each had an abundance of 400 or more sequences. These two haplotypes were unlikely to be the result of error as they diverged from the most abundant sequences by 4 and 8 bases. The final haplotype selected (haplotype 3) was from the Jan17:6 sample. While this haplotype differed only by a single base from the most abundant sequence, the abundance of haplotype 3 (262) was almost half of that of the main haplotype in that sample; haplotype 1 (589). This sequence ratio indicates that haplotype 3 is not the result of sequencing error (Table S5).

Only haplotype 3 was unique to Barunguba, haplotypes 4 & 5 were unique to Cape Bridgewater and the remaining two haplotypes (1 & 2) were detected at both sites (Figure 5B). Haplotype 1 matched the 12s rRNA gene fragments of two independent penguin samples (MF370525, MK761006) from whole genomes from the same region (Sarker et al., 2017; Vianna et al., 2020). Across all ten penguin positive scat samples, six contained a single haplotype, three contained 2 haplotypes and one contained 4 haplotypes (Figure 5B). It was reasonable to assume that a fur seal could consume at least part of more than one penguin per feeding opportunity. Therefore, by taking two distinct genetic haplotypes present within a sample to represent at least two distinct individual birds consumed, we estimated at least 16 individual penguins were consumed across the 99 scat samples, from two sampling locations and multiple seasons.

Discussion

We provide an updated and nuanced predation incidence range for seabirds overall (9.1–29.3%) and little penguins (6.1–25.3%) in long-nosed fur seal diets based on detection rates of hard-part and genetic assays, as well as the detection rate of genetically abundant samples, expanding previous information on occurrence rates (2–13%) (Page et al., 2005; Hardy et al., 2017; Goldsworthy et al., 2019). We confirm that little penguins remain the most commonly consumed seabird by long-nosed fur seals in comparison to other avian taxa as detected by both assay types. Here, DNA metabarcoding also offered key advantages over morphological analysis – quantitative information on: (i) greater diversity discovery by detection of multiple prey taxa within a single scat sample, (ii) absolute and relative abundances of taxa recovered using DNA to provide more conservative and nuanced estimates of predation rates, and (iii) preliminary estimates of predation impact, likely involving at least 16 individual penguins, by exploring the spatial and temporal distribution of haplotypes.

Implications for long-nosed fur seal and little penguin interactions

Previous studies using either assay technique alone have identified little penguin remains at relatively low frequencies in relation to seals’ total diets (as a percentage of samples: 5.9% in Page et al., 2005, <2% in Hardy et al., 2017, ~13% in Goldsworthy et al., 2019). These past values are similar to the lower range of estimates observed here (9.1% of samples for seabirds, 6.1% for penguins), and which corresponded to the number of samples containing large quantities of prey DNA. However, the detection from hard–part remains (29.3% of samples for seabirds, and 25.3% for penguins) and the upper range of estimates from genetic remains (21.2% of samples for seabirds, and 10.1% for penguins) is possible and is concerning if further monitoring reveals increases in penguin consumption by long-nosed fur seals. We report this range of predation occurrence because, regardless of the assay used, the persistence of feathers or soft tissues in scats remains an area of research and development (Goldsworthy et al., 2019), as does the relationship between genetic information and quantity of prey consumed (Thomas et al., 2014; Deagle et al., 2019).

Improving confidence in the results of these assays requires longer-term replication of the work and ongoing monitoring of this wildlife conflict using multiple surveillance techniques, paired with experimentation (for example Reinhold et al., 2022) to ascertain the closest to reality estimate of consumption of a species of conservation concern. Notably, when using relative read abundance information, we provided a more conservative estimate of predation incidence based on samples that contained large quantities of DNA and thus likely to represent a fresher meal and more recent predation event. These estimates, of 6.1% of samples containing little penguins, were most similar to past results from hard-part analyses generated from much larger datasets (> 1000 samples, <5% of samples; Page et al., 2005).

DNA metabarcoding was useful in detecting taxonomic mixtures in scat samples, here containing up to 2 distinct seabird prey taxa compared to hard-part analysis that did not detect more than one prey taxon at a time, and for little penguins using haplotype polymorphism analyses – the potential number of individual penguins consumed. Multiple individual birds are impossible to distinguish accurately from feathers and the number of scats required to pass a single bird is equally hard to quantify in the field (Page et al., 2005; Goldsworthy et al., 2013). The genetic assay provided multiple types of quantitative information on penguins consumed with which to overcome this issue, including a relative read abundance metric that confirmed little penguins as the most common seabird prey and a conservative estimate of individual penguin intraspecific diversity.

Our results confirmed that most scats with penguin DNA (n = 6 out of 10) contained a single haplotype, however, we found that four scats contained multiple haplotypes, with up to 4 distinct genetic haplotypes, suggesting that long-nosed fur seals can consume not only multiple seabird taxa but also multiple penguins in a single foraging trip (or, given the time it takes for prey to traverse the gastrointestinal system, within about 48h of sampling). This result is more consistent with pinniped foraging ecology than a previous assumption that each scat containing feathers corresponded to a single bird (Page et al., 2005), because fur seals typically process large prey by shaking it and breaking it up at the surface and may predate multiple prey in a single foraging event (Mumma et al., 2016; Hocking et al., 2016). Such feeding behavior can also result in partial consumption of the prey item and incomplete consumption of diagnostic hard parts. Therefore using genetic tools to fully extract intraspecific genetic diversity within samples warrants further application (see Nearing et al., 2018). Specifically, using additional genetic assays that target long fragments of less conserved genes could enable greater detection of intra-specific genetic diversity than the conservative 12s ribosomal gene (Banks et al., 2002), the use of which in this study, may have produced an underestimate of the number of penguins that could have been consumed.

While not their preferred prey in comparison to cephalopods and fishes (Page et al., 2005; Goldsworthy et al., 2013; Hardy et al., 2017), it is entirely plausible for long-nosed fur seals which can weigh from 50–150 kg (McKenzie et al., 2007) to consume more than 4 individual little penguins (or parts of) weighing less than 1.5 kg (Williams, 1995). If little penguin predation becomes an important individual foraging strategy even for a fraction of long-nosed fur seal populations, this could have serious negative impacts for isolated penguin populations. In the context of recent crashes in little penguin colonies (Sutherland et al., 2022), these data signal a need for broader and increased monitoring of predation mortality to little penguins to inform conservation and management strategies. Predatory behaviors could be transmitted to other predator populations such as the Australian fur seal (A. pusillus pusillus) or Australian sea lions (Neophoca cinerea) which partially overlap in geographic ranges. Despite little (Deagle et al., 2009; McIntosh et al., 2006) to no (Berry et al., 2017) prior identification of little penguins in these species’ diets, prey switching could become a possibility in response to food web disruption under ocean warming and changes in prey availability (Kliska et al., 2022) and this could have cascading effects on penguin populations.

We note that two thirds of the LNFS samples did not contain seabird DNA or hard-parts, and likely consist of fishes and cephalopods, their more common prey (Hardy et al., 2017). Further research of their total diet is ongoing using additional primers to investigate the importance of the non-avian composition of long-nosed fur seal diets in southeastern Australia. Additionally, the impacts of long-nosed fur seal predation need to be considered and managed within the wider forum, ideally using up-to-date cumulative impact assessments for threats to little penguins. For little penguins and in other wildlife conflict situations, endemic predation is natural whereas habitat degradation, environmental change and invasive species are often more significant sources of impact to susceptible species (Hervieux et al., 2014; Marshall et al., 2016; Ropert-Coudert et al., 2019).

Important methodological considerations and improvements

We highlight that for all assays (i.e., morphological vs. genetic) and metrics (i.e., detections vs. relative abundance data), there are important biological and methodological constraints and biases to consider when attempting to estimate predation incidence (detection and frequency-based) and magnitude (impacted number of prey or biomass) in reality.

Deagle et al. (2019) simulated the impacts of biases from metrics reported on diet summaries from 20 independent studies, and this indicated that relative read abundance information often provided a more accurate and nuanced view of population-level diet than frequency of occurrence-based data, despite known recovery biases (i.e., such as from amplification bias) – thus highlighting the importance of reporting relative and total abundance data from the genetic assay. Experiments performed on seals have shown relationships between relative read abundance and the ranked importance of prey species, and correlations with biomass consumed such that it is possible to develop correction factors (Deagle et al., 2005; Thomas et al., 2014, 2016; Deagle et al., 2019). It may become relevant to develop genetic tissue correction factors through experimental diet studies (similar to Reinhold et al., 2022) for this wildlife conflict scenario. Therefore, publishing all sound information on predation of seabirds by long-nosed fur seals will be useful to future investigations of this species interaction scenario.

Notably, we report sequence abundances after they have passed stringent sequence quality filtering steps, including a stringent threshold for excluding sequences from the recommended 0.1% (Fox et al., 2014) to 1% of total sequence abundance. This also follows many conservative decisions for the entire pipeline, including the use of a custom single-step fusion PCR over two-step PCR which reduces the risk of sample cross-contamination at a risky point in the workflow, that of amplification (Schnell et al., 2015; Taberlet et al., 2018). DNA-based diet analysis studies crucially rely on and have reported quantitative and semi-quantitative results from genetic assays for decades (Deagle et al., 2005). However stringency also represents a significant loss of data on wildlife interactions (Deagle et al., 2013) and thus we likely excluded several true positives that may have occurred with low DNA-abundance.

Due to financial constraints, technical replicates for each sample were not used in this study, a common practice in wildlife research where budgets are limited. Reported diet diversity would likely be higher when using multiple replicates for samples, thus our detection rate of seabird species and little penguin occurrences would likely have been higher when covering a larger quantity of the scat or results in slightly different dietary proportions than what the predator consumed (Mata et al., 2019). Our final study limitation is the use of manual sample homogenization and partial sub-sampling, while very common in DNA metabarcoding studies in wildlife biology, it remains possible that not all seabird prey species were detected from these samples and that dietary proportion information would be improved from the use of multiple replicate sub-samples. An experimental improvement where costs remain limiting would be to pool multiple sub-sampler replicates prior to sequencing (Shirazi et al., 2021; Van den Bulcke et al., 2021).

As a species of conservation concern and public value, the distinct variation in little penguin metabarcodes obtained in this study, beyond that attributed to sequencing error, led us to explore how many individuals might be within the data (similar to Seersholm et al., 2018 and Tsuji et al., 2020). However, haplotype diversity on this locus was likely underestimated because the 12S ribosomal RNA is conserved within species (Banks et al., 2002), which we selected for its proven reliability in detecting interspecific variability across seabirds (Berry et al., 2017; Hardy et al., 2017). Targeting alternative markers (variable barcodes or microsatellites) could reveal greater genetic diversity, however we did not succeed with a COI barcode (Appendix S1.2) due to the scarcity of longer fragments of DNA in samples.

In this study, we identified five individual haplotypes representing intra-specific diversity within little penguin DNA obtained after stringent and conservative sequence filtering and error mitigation were performed, then selected only in the samples containing high sequence abundance for this taxon. One of these matched the 12s rRNA gene fragments of two independent penguin samples (MF370525, MK761006) from whole genomes from the same region (Sarker et al., 2017; Vianna et al., 2020). However, these were the only two whole genomes for little penguins published for southeastern Australia at the time of analyses. The selection of haplotypes was conservative, and false positives are unlikely, however, we acknowledge that several true haplotypes exist in the data and were excluded. Given the successful identification of realistic haplotypes in this study, we recommend that future studies of this wildlife conflict and others, further develop methods to target intra-specific diversity for species of concern. Notably, we recommend larger-scale sampling across the broader geographic range of this predator, including multiple seasons and years as sampling depth is an extremely important factor in inter- and intra-biodiversity (Shirazi et al., 2021). We recommend pairing a greater sampling depth with more in-depth analyses of intra-specific diversity using both multi-locus assay metabarcoding protocols and post-clustering “denoising” pipelines to determine biological sequences at single nucleotide resolution from sequencing error noise (Nearing et al., 2018).

While both metabarcoding and morphological methods provided statistically similar results in overall detection rates, in applying these assays together or separately, it is not realistic to assume these methods should detect target taxa in the same samples or even at the same rates and further research to better understand and potentially convert information between these assay types is needed (Thomas et al., 2014). Here, a greater diversity of prey were detected using both assays than by either assay alone and both assays cumulatively provide a range of estimates of predation and a range of information with which to assess predation incidence (detection and frequency-based) and magnitude (impacted number or biomass of prey) more rigorously than either method would alone, a finding shared by over 40 similar comparative studies (Hardy, PhD Thesis, 2018). Different detection rates reported between these methods and on a sample-by-sample basis in this study are also common across a broad range of predator and prey taxa (Tollit et al., 2009; Zarzoso-Lacoste et al., 2013), with studies reporting as little as 10% of positive samples containing both the hard-parts and DNA of the same species. This is likely due to differential gastrointestinal passage times of prey soft tissues (< 48 h) compared to their hard tissues (up to 7 days in another otariid study; Tollit et al., 2009). Soft tissues likely contain greater concentrations of DNA particularly for mitochondrial DNA, than do hard parts (Mumma et al., 2016; Granquist et al., 2018); therefore, even when present in the same samples, genetic information from soft tissues would be preferentially amplified because DNA from chitinous tissues (i.e., hard-parts such as feather, fur, or bone) is of poorer quality and requires specific processes for extraction (McDonald and Griffith, 2011; Rothe and Nagy, 2016).

We therefore emphasize that genetic and morphological analyses represent complementary assays of predation. For many conservation practitioners and ecologists aiming to investigate population-level predation across a range of wildlife interaction scenarios, treating these assays as semi-independent sampling methods that provide complementary information on predation will be fit-for-purpose. However, in monitoring interactions and consumption between rare species, and species of high conservation concern, additional assays and experimental design procedures would enable direct sample-by-sample comparison of consumed taxa, such as greater sub-sample replication and potential sample-based replicate pooling, as well as amplification and sequencing of all samples beyond those containing a positive in initial screening PCR. In our case, this was not financially feasible and we posit for this study system that additional funding could be spent on other value-add assays such as: (i) including screening for predator genetic diversity to identify individuals in a population contributing to predation of a sensitive or valuable species (Wegge et al., 2012), (ii) developing species-specific probes using older and cheaper technology (Fox et al., 2012) including cross-validation of eDNA data with quantitative PCR (Murray et al., 2011), and/or (iii) development of penguin-specific DNA-to-tissue-based correction factors which could provide consumed biomass information (Thomas et al., 2014).

Conclusions

We leveraged the growing utility of DNA metabarcoding alongside traditional diet analysis methods to investigate a wildlife conflict in southeastern Australia and with the goal to better inform decision making and long-term monitoring of this conflict. Our multi-assay results provide updated information on predation incidence by long-nosed fur seals on culturally valued little penguins in southeastern Australia at a critical time in the conservation management of both species. The genetic assay provided more nuanced estimates of possible predation rates, including additional semi-quantitative information useful for conservation practitioners in understanding the potential range in predation incidence. The genetic assay enabled exploration of genetic diversity within samples and usable estimates of the number of penguins consumed within the cross-section of sampled locations and times. This study demonstrates a need for research and development of techniques at the nexus of population genetics and environmental sampling. Additionally, predator impacts need to be considered and managed within up-to-date cumulative impact assessments for threats. We have delivered an important step towards this for little penguins in south-eastern Australia.

Data availability statement

The raw sequenced data files are publicly available at this location to be cited as: Hardy, Natasha et al. (2024). Sequencing data for seabird eDNA in long-nosed fur seal diets from southeastern Australia [Dataset]. Dryad. https://doi.org/10.5061/dryad.stqjq2cb3. Processed data, code and figures drafted for this publication are publicly available on Github: https://github.com/NatashaAHardy/pinp_stats.

Ethics statement

The animal study was approved by University of Sydney ethics permit (L04/9-2013/4/6056); Phillip Island Nature Parks Ethics Permit (2.2016). The study was conducted in accordance with the local legislation and institutional requirements.

Author contributions

NH: Conceptualization, Data curation, Formal analysis, Funding acquisition, Investigation, Methodology, Project administration, Resources, Software, Validation, Visualization, Writing – original draft, Writing – review & editing. TB: Conceptualization, Investigation, Methodology, Validation, Writing – original draft, Writing – review & editing. MB: Conceptualization, Investigation, Methodology, Resources, Supervision, Visualization, Writing – review & editing. NB: Methodology, Resources, Software, Supervision, Writing – review & editing. WF: Funding acquisition, Project administration, Resources, Supervision, Writing – review & editing. RM: Conceptualization, Funding acquisition, Investigation, Project administration, Resources, Supervision, Writing – original draft, Writing – review & editing.

Funding

The author(s) declare financial support was received for the research, authorship, and/or publication of this article. This project was funded by Phillip Island Nature Parks and the Australian Research Council Linkage Grant to Will Figueira (LP120100228) for support in NSW sample collections and processing.

Acknowledgments

We are grateful to the assistance of Cecilia Power and Lisa Lee Nen That in the Bott Laboratory at RMIT Bundoora. We thank Dáithí Murray for the provision and use of his bioinformatics pipeline. We are grateful to all field volunteers for sample collections, especially Tony Mitchell (Department of Energy, Environment and Climate Action), and Francois and Clover (Underwater Safaris Narooma) for their dedication ensuring sample collections. Samples from NSW were collected for another project and sub-samples were submitted to Phillip Island Nature Parks for screening of bird DNA. Samples from NSW were collected under University of Sydney ethics permit (L04/9-2013/4/6056); Australian Government permits to conduct research under the EPBC Act (AU-COM2013-224), and from the Office of Environment and Heritage NSW Scientific License (SL101244). Victorian research was performed under Phillip Island Nature Parks Ethics Permit (2.2016) and Department of Environment, Land, Water and Planning Research Permit (10007974). Harley Schinagl at Phillip Island Nature Parks produced Figure 1. We thank Karling Roberts, Taylor McLeod, Jerry Moxley and peer-reviewers for vastly improving this manuscript. We acknowledge the following Nations and Traditional Owners on whose unceded lands we conducted this research: Yuin (Barunguba), Bunurong (Millowl, Phillip Island), Gunditjmara (Cape Bridgewater and Deen Maar Island), Eastern Maar (Deen Maar Island), Kulin (RMIT Bundoora) and Eora (USYD).

Conflict of interest

NH is employed by SMRU Consulting.

The remaining authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.

Publisher’s note

All claims expressed in this article are solely those of the authors and do not necessarily represent those of their affiliated organizations, or those of the publisher, the editors and the reviewers. Any product that may be evaluated in this article, or claim that may be made by its manufacturer, is not guaranteed or endorsed by the publisher.

Supplementary material

The Supplementary Material for this article can be found online at: https://www.frontiersin.org/articles/10.3389/fmars.2024.1288769/full#supplementary-material

References