Secondary predation constrains DNA-based diet reconstruction in two threatened shark species

Increasing fishing effort, including bycatch and discard practices, are impacting marine biodiversity, particularly among slow-to-reproduce taxa such as elasmobranchs, and specifically sharks. While some fisheries involving sharks are sustainably managed, collateral mortalities continue, contributing towards > 35% of species being threatened with extinction. To effectively manage shark stocks, life-history information, including resource use and feeding ecologies is pivotal, especially among those species with wide-ranging distributions. Two cosmopolitan sharks bycaught off eastern Australia are the common blacktip shark (Carcharhinus limbatus; globally classified as Near Threatened) and great hammerhead (Sphyrna mokarran; Critically Endangered). We opportunistically sampled the digestive tracts of these two species (and also any whole prey; termed the ‘Russian-doll’ approach), caught in bather-protection gillnets off northern New South Wales, to investigate the capacity for DNA metabarcoding to simultaneously determine predator and prey regional feeding ecologies. While sample sizes were small, S. mokkaran fed predominantly on stingrays and skates (Myliobatiformes and Rajiformes), but also teleosts, while C. limbatus mostly consumed teleosts. Metabarcoding assays showed extensive intermixing of taxa from the digestive tracts of predators and their whole prey, likely via the predator’s stomach chyme, negating the opportunity to distinguish between primary and secondary predation. This Russian-doll effect requires further investigation in DNA metabarcoding studies focussing on dietary preferences and implies that any outcomes will need to be interpreted concomitant with traditional visual approaches.

Global biodiversity is under threat, with accelerating losses of species and increasing concern about ecosystem changes 1 . In the oceans, escalating fishing effort and associated bycatches and discarding practices are impacting on marine biodiversity 2 . Levels of extinction threat vary regionally and by taxonomic group, but owing to their low reproductive rates and late age-at-maturity, elasmobranchs, especially sharks, are among those taxa highly susceptible to increasing anthropogenic pressures 3 . More specifically, among the total known species of 'ground sharks' or the Carcharhiniformes (n = 287), 36% are listed as threatened with extinction by the International Union for Conservation of Nature (IUCN), and a further 18% are classified as 'Data Deficient' 4 . To mitigate this extinction threat, there is a critical need for more information on the ecology of sharks, particularly those that are threatened or endangered.
Two of the more abundant, cosmopolitan Carcharhiniformes found in Australia and well represented in bycatches of bather-protection fishing gears 7,8 are the common blacktip shark (Carcharhinus limbatus) and great hammerhead (Sphyrna mokarran), which are globally classified as Near Threatened and Endangered, respectively 4 . Sphyrna mokarran is also registered in Appendix II of the Convention on International Trade in Sample collection. The study was done using seven S. mokarran (three females and four males) and four C. limbatus (two of each sex) that died in bather-protection gillnets deployed off Lennox Head, Ballina, and Evans Head, NSW, Australia (28.77° S, 153.60° E to 29.10° S, 153.44° E) between 7 February and 13 March, 2018 ( Fig. 1). Each gillnet measured 150 m long × 4 or 6 m deep and comprised 600-or 800-mm mesh made from either 1.8-or 2.1-mm diameter twisted polyethylene, or 2.5-mm diameter polyamide twine (see 8 for details of the fishing gear).
Whole sharks were removed from the gillnets and stored at − 20 °C (within 3 h) until necropsied in May 2018. During the necropsy process, specimens were defrosted for 12 h and measured for TL before the stomach cavity was opened in a sterile field laboratory with bleach sterilized tools. All digestive contents were removed, with any animal matter identified (where possible), preserved in 100% ethanol, and stored at − 20 °C until further analyses. Extraction controls for DNA (ultrapure water samples, Invitrogen, Waltham, USA) were collected in sterile 1.5 mL microcentrifuge tubes alongside the stomach-content samples and subjected to the same workflow described below.
Stomach-content analyses and DNA extraction. Whole prey removed from S. mokkaran and C. limbatus were identified by visual inspection, most to the species level. Stomach contents of some of these whole prey items, comprising seven individual Urolophus sp. and two A. rostrata rays, predated by three S. mokarran individuals, were included in the metabarcoding analysis ( Table 1). Some of the incomplete specimens, typically teleost jawbones or scales, were identified to genus only (Fig. 2). For each S. mokkaran and C. limbatus, and each whole prey item detailed above, the remainder of its unidentifiable stomach contents (both liquids and solids) were separately homogenised using a previously sterilised (washed with detergent, followed by a 10% bleach wash, and rinsed thoroughly in MilliQ water) commercial food blender.
For all samples, DNA was extracted from approximately 5 mL of the homogenate using a QIAamp PowerFecal DNA kit (Qiagen, Sydney, Australia) according to manufacturer's instructions. This kit effectively removes PCR inhibitors from fecal and stomach content samples. The DNA extraction was carried out in a pre-PCR laboratory to minimise contamination, and clean-room protocols were followed with extensive bleaching and UV treatment of the area and equipment for all laboratory steps. Filter pipette tips were used in all instances and gloves were frequently changed, particularly between handling specimens or plasticware.
PCR amplification and Illumina sequencing. The DNA extracts from each stomach sample were amplified, tagged separately, and then pooled for sequencing. Two group-specific mini-barcode primers were selected for the amplification of teleost and crustacean DNA, targeting 12S (MiFish 40 ) and 16S (Crust16S short 41 ) mitochondrial DNA genes, respectively. We also used a universal 18S primer set 42 targeting the hypervariable V4 region of the nuclear small subunit ribosomal DNA to amplify templates from a broader fraction of marine metazoans. Polymerase chain reaction (PCR) was performed using the AmpliTaq Gold 360 protocol and thermocycling conditions recommended in 30 . , and Uni18S primer sets, respectively, and products were run on a 1% agarose gel to confirm amplification of the correct target size (MiFish = ± 170 bp; Crust16S = ± 170 bp, Uni18S = ± 220 bp). A second round of PCR was undertaken with the cleaned PCR products using unique dual-indexed primers for each sample, which included the Illumina-specific sequencing adaptors. PCR products were sent to the Ramaciotti Centre for Genomics at the University of NSW for cleaning, normalising, and pooling prior to paired-end sequencing, which was performed using a 500 cycle MiSeq V3 Reagent Kit on an Illumina MiSeq platform (Illumina, San Diego, CA, USA). Sample demultiplexing based on the incorporated indexes was conducted by the sequencing centre.
The ZOTUs were queried against the National Centre for Biotechnology Information's (NCBI) GenBank nucleotide database (accessed in 2020) using BLASTn with the following settings: percentage identity = 97, query coverage = 100, best hit score edge of 0.05, best hit overhang of 0.25, and an E-value of 1e-3. An additional, lessstringent data set was generated for each assay using the following BLASTn settings: percentage identity = 90,    44 was then run to curate the assignments assessing sequence similarity and their co-occurrence patterns with the default parameters: minimum_ratio_type = min, minimum_ratio = 1, minimum_match = 84, minimum_relative_cooccurence = 0.95. This entire process was completed on the Zeus SGI cluster based at the Pawsey Supercomputing Centre in Kensington, Western Australia using an abridged version of the fully automated 'eDNAFlow' pipeline 45 . All raw sequencing data needed to replicate the study are available from Dryad Digital Repository: https:// doi. org/ 10. 5061/ dryad. kd51c 5b4s.

Statistical analysis of the Russian-doll effect.
To test for the Russian-doll effect, and more specifically, to confirm no significant differences between metabarcoding reads from predators' stomachs (S. mokarran) and those of their whole prey (Urolophus sp.; A. rostrata), Kruskal's non-metric multidimensional scaling (NMDS) was applied to each metabarcoding assay (12S, 16S, 18S). These analyses were done using the isoMDS function in the Mass package 46 in R (R Core Team, 2017) 47 . For comparison, NMDS with stable solution from random starts was implemented, using the metaMDS function in the Vegan R package 48 . Permutational multivariate analysis of variance (PERMANOVA) was used to assess if the metabarcoding reads from the predator and its whole prey were equivalent, using the ADONIS function in Vegan 48 . Where the null hypothesis was rejected, a similarity percentages breakdown (SIMPER) was applied (using Vegan 48 ) to determine the contribution of each dietary species to the dissimilarity between metabarcoding reads from predator and prey stomach contents. For all of these analyses, we tested for differences between predator and whole prey using the number of reads for each taxon identified per assay, and also presence/absence of that taxon per assay.

Results
Four of the S. mokarran had identifiable whole or partial prey in their stomachs (Table 1, Supplementary Table 1). These prey taxa included stingarees (Urolophus sp.) (Fig. 2), eastern shovelnose rays (Aptychotrema rostrata), thornback cowfish (Lactoria fornasini), eastern smooth boxfish (Anoplocapros inermis), threebar porcupinefish (Dicotylichthys punctulatus), reef ocean perch (Helicolenus percoides) and a whiting (Sillago sp.). Notably, one S. mokarran stomach contained nine whole Urolophus sp., with a total stomach weight of 6.26 kg, and a second S. mokarran stomach contained three whole Urolophus sp. and three whole A. rostrata (6.15 kg) ( Fig. 2; Table 1).   Table 1). The PCR negative controls showed low levels of human contamination for the 18S assay only (Table 1), whereas the control water sample collected alongside the C. limbatus stomach samples demonstrated contamination from Carcharhinidae DNA (most likely from C. limbatus). Carcharhinidae and/or Chondrichthyes reads were ubiquitous across all three genetic assays, including the Urolophus sp. and A. rostrata rays removed from the S. mokarran stomachs, and most likely represent the host predator's DNA. Taxa identified in the S. mokarran metabarcoding assays included crustaceans, cartilaginous fish and teleosts. Crustaceans and teleosts were identified in the C. limbatus stomach contents (Table 1). In most cases, stomachcontent analyses on the Urolophus sp. and A. rostrata removed from the S. mokarran stomachs displayed high taxon similarity to that of the apex predators' stomachs. For example, a comparison of 4H (S. mokarran) and 4SR (Urolophus sp.) revealed Gnathophis spp., Lepidotrigla spp., Pseudorhombus spp., and Anoplocapros spp. in common ( Table 1).
Testing the Russian-doll effect. Both of the NMDS tests and PERMANOVA showed no significant differences between stomach contents from predator and prey for the 12S and 16S metabarcoding assays (p > 0.05) for number of reads, and for presence/absence comparisons, but did show a significant difference for number of reads for the 18S assay only (p < 0.05, Fig. 3, Table 2). The SIMPER analysis showed that the dissimilarities between number of reads for 18S predator and prey comparisons were driven predominantly by Chondrichthyes (38%) and Carcharhinidae (30%) (Fig. 4), which were elevated in the predator, and Actinopteri (4%), which was elevated in the prey. These results indicate that taxonomic composition, and even the number of metabarcoding reads from stomach contents of predator and their whole prey, were not statistically different from one another, at least for 12S and 16S. Differences in the number of reads, but not taxonomic composition, for the 18S assay likely reflects increased host DNA sequenced for this marker (i.e. Chondrichthyes, Carcharhinidae; but not Actinopteri) (Fig. 4). One caveat of these analyses is that the sample size was small, with three replicates for S. mokarran.

Discussion
The data collected here represent the first efforts at metabarcoding the stomach contents of S. mokarran and C. limbatus off eastern Australia. Our results support the general trends in the literature describing the diet of S. mokarran as being dominated by rays [14][15][16] , but also including teleosts and other sharks 12,13,17,18 , and reports of C. limbatus predominantly feeding on teleosts [20][21][22] . By 'Russian-dolling' 39 the stomach contents of these marine predators and their prey, we have also quantified some of the issues encountered when attempting to reconstruct trophic interactions from metabarcoding stomach contents.
Of the two species, S. mokarran dietary preferences were the most discernible, and clearly some had very recently fed in the general area (considering the regional abundances of some whole prey items). Based on whole prey items, S. mokarran consumed stingrays and skates (Myliobatiformes and Rajiformes), but also teleosts including both slow-(Anoplocapros inermis) and fast-moving species (Sillago sp.). Observations of prey handling by S. mokarran are limited, however the species has been observed to use its laterally expanded head (cephalofoil) to immobilise prey on the ocean floor during both the pursuit of bottom-dwelling rays (southern stingray Dasyatis americana) 16 , and the post-capture manipulation of pelagic rays (spotted eagle ray Aetobatus narinari) 14 . Most of the prey items of the S. mokarran individuals identified here were bottom dwellers (e.g. rays: Urolophus sp., A. rostrata and perch: H. percoides), which supports the above mechanisms of prey-capture and -manipulation (e.g. [49][50][51]. Unlike for S. mokarran, we found no whole fish in the digestive-tracts of the four C. limbatus. Nevertheless, metabarcoding indicated the presence of large species including Platycephalus spp. Several platycephalids occur nearshore off northern NSW, especially eastern bluespotted flathead (P. caeruleopunctatus up to 0.6 m TL) and, to a lesser extent, dusky flathead (P. fuscus up to 1.2 m TL). Large individuals of these species were unlikely to be prey of other animals in the stomachs of C. limbatus, but in the absence of whole specimens any conclusions based on the exact platycephalid species are speculative.
Other primary target species identified by metabarcoding in both S. mokarran and C. limbatus included Dexillus spp., probably tufted sole Dexillus muelleri; the only known species from this genus, occurring off more tropical areas of Australia. This might imply these sampled sharks (5H, 6H, 1BT, and 4BT) moved southwards from a northern more, tropical foraging area 13 , or less likely that an as yet unidentifed species of Dexillus occurs in northern NSW waters. Five of the S. mokarran and two of the C. limbatus also had Anoplocapros spp. metabarcoding reads sequenced from their stomach contents, and we observed A. inermis scales in the stomach of one S. mokarran, suggesting some dietary overlap with C. limbatus. Indeed, A. inermis is the only species in this genus (of three total species, including A. amygdaloides and A. lenticularis) found off eastern Australia.
While some larger primary prey items of S. mokarran and to a lesser extent C. limbatus could be identified, it was difficult to ascribe predation sources to the smaller taxa, such as the sandy sprat Hyperlophus vittatus or biddies Gerres sp. (likely silver biddy G. subfasciatus, although this species overlaps with at least two other congeners in northern NSW) detected in the 12S assay. These species could either be primary prey of the sharks or secondary prey consumed by the platycephalids. The Russian-doll effect 39 is exemplified well by comparing predator-prey stomach content metabarcoding reads. For those S. mokarran that contained whole prey items (e.g. Myliobatiformes), DNA metabarcoding assays yielded close matches between predator and prey stomach   www.nature.com/scientificreports/ contents. It is therefore difficult to conclude whether the S. mokkaran (less likely) or Urolophus sp. or Aptychotrema rostrata (most likely) were feeding on smaller prey items, such as decapod crustaceans and fish species. We suggest that the chyme of both predator and (whole) prey intermixed with the other. Similarly, our metabarcoding approach could not determine whether smaller sharks had been consumed by S. mokarran (as expected from the literature), due to the fact that all S. mokarran and their whole prey stomach samples included an abundance of 12S Carcharhinidae reads. However, the 16S assay did suggest these reads likely reflected the host predator's DNA, considering S. mokarran reads were identified both from S. mokarran and Urolophus sp. stomach contents, but at very low abundances. Blocking primers designed to exclude the host predator's DNA would be a worthwhile inclusion in future metabarcoding dietary studies (e.g. 36,52,53 ), while shark-specific primers could also prove more informative on this topic 30 . Alternatively, increasing the depth of sequencing may facilitate detecting the low template (prey) fraction of the metabarcoding chyme despite the presence of host DNA that would otherwise swamp the PCR amplification process.
In this study, metabarcoding assays provided some insights into the dietary preferences of S. mokarran and C. limbatus off eastern Australia. These species appear to have some dietary overlap, but with consistency among prey species identified in studies of their feeding ecology from across their broader ranges [12][13][14][15][16][17][18][20][21][22] . Sphyrna mokkaran fed predominantly on Myliobatiformes and Rajiformes, but also teleosts, whereas C. limbatus fed predominantly on teleosts, which is consistent with its smaller body size and lack of cephalofoil that allows specialized feeding of benthic prey.
The Russian-doll effect made the reconstruction of trophic interactions from stomach metabarcoding data problematic. Extensive intermixing in situ between predator-prey digestive tracts was evident, which also limited the ability to discriminate between primary and secondary predators of smaller teleosts and crustacean taxa. The literature on the Russian-doll effect 39 in metabarcoding dietary studies is limited (reviewed in 27,38 ) but certainly deserves further attention (e.g. 54 ). We suggest that while the approach offers some utility in identifying unseen taxa, confirmation via visual identification and across sufficient replication (ideally from various sampling methods) is still required to comprehensively understand primary dietary preferences.

Data availability
The datasets generated and analysed during the current study have been archived in Dryad: https:// doi. org/ 10. 5061/ dryad. kd51c 5b4s.  , showing a statistically significant increase (marked as ' > ' , p < 0.05) of the number of reads from the predators' stomach contents compared to that of their whole prey, likely reflecting elevated sequencing of host DNA for this assay (Chondricthyes and Actinopteri not presented).