Integrating multiple chemical tracers to elucidate the diet and habitat of Cookiecutter Sharks

The Cookiecutter shark (Isistius brasiliensis) is an ectoparasitic, mesopelagic shark that is known for removing plugs of tissue from larger prey, including teleosts, chondrichthyans, cephalopods, and marine mammals. Although this species is widely distributed throughout the world’s tropical and subtropical oceanic waters, like many deep-water species, it remains very poorly understood due to its mesopelagic distribution. We used a suite of biochemical tracers, including stable isotope analysis (SIA), fatty acid analysis (FAA), and environmental DNA (eDNA), to investigate the trophic ecology of this species in the Central Pacific around Hawaii. We found that large epipelagic prey constituted a relatively minor part of the overall diet. Surprisingly, small micronektonic and forage species (meso- and epipelagic) are the most important prey group for Cookiecutter sharks across the studied size range (17–43 cm total length), with larger mesopelagic species or species that exhibit diel vertical migration also being important prey. These results were consistent across all the tracer techniques employed. Our results indicate that Cookiecutter sharks play a unique role in pelagic food webs, feeding on prey ranging from the largest apex predators to small, low trophic level species, in particular those that overlap with the depth distribution of the sharks throughout the diel cycle. We also found evidence of a potential shift in diet and/or habitat with size and season. Environmental DNA metabarcoding revealed new prey items for Cookiecutter sharks while also demonstrating that eDNA can be used to identify recent prey in stomachs frozen for extended periods. Integrating across chemical tracers is a powerful tool for investigating the ecology of elusive and difficult to study species, such as meso- and bathypelagic chondrichthyans, and can increase the amount of information gained from small sample sizes. Better resolving the foraging ecology of these mesopelagic predators is critical for effective conservation and management of these taxa and ecosystems, which are intrinsically vulnerable to overfishing and exploitation.

www.nature.com/scientificreports/ Although it appears that sharks will consume small squid whole, many of the squid remains were estimated to be from relatively large squid, often as large as the sharks themselves. Given the lack of empirical data on the diet of Cookiecutter sharks, we only have the most basic understanding of their trophic ecology with little to no information on size, temporal, or spatial variability in diet. In addition, due to the high proportion of empty stomachs encountered in these studies, and the historical difficulty of identifying the origin of the various tissue plugs in shark stomachs, traditional approaches to assessing diet (e.g., stomach content analysis) have been of reduced utility. The goal of this study was to use multiple biochemical tracer approaches (SIA, FAA, and eDNA) to elucidate the trophic ecology of Cookiecutter sharks in the Central Pacific and demonstrate the utility of this integrated approach for investigating the trophic ecology of difficult to study species, where sample sizes are often limited. In particular, we were interested in using SIA and FAA to assess the relative importance of different prey groups (mesopelagic, epipelagic, and vertically migrating prey) to Cookiecutter sharks across a size range. We also investigated the utility of using eDNA approaches to identify prey from Cookiecutter shark stomachs.

Results
Stable isotope analysis: comparisons of size and sex. There was no significant relationship between length and δ 13 C or δ 15 N values ( Fig. 1) for liver (linear regression, δ 13 C: r 2 = 0.22, p > 0.09; δ 15 N: r 2 = 0.03, p > 0. 54) or muscle (linear regression, δ 13 C: r 2 = 0.17, p > 0.15; δ 15 N: r 2 = 0.04, p > 0.51). However, the smallest shark (17.1 cm TL) had one of the highest δ 15 N and lowest δ 13 C values of all sharks for both liver and muscle. There was no significant difference between sexes in muscle (δ 13 C: Mann-Whitney Rank Sum Test, p = 0.38; δ 15 N: T-test, p = 0.39). In liver, there was no significant difference in δ 15 N (t test, p = 0.12) though there was in δ 13 C (T test, p = 0.03), with females having higher δ 13 C than males. However, low samples sizes limit the power of these tests.
Stable isotope analysis: prey contributions. Potential Cookiecutter shark prey groups generally clustered based on their habitat type (see "Methods"). One group (hereafter "DVM group"), was primarily comprised of large species that exhibit diel vertical migration as well as marine mammals. Another group ("MESO group"), was entirely comprised of small mesopelagic and epipelagic micronektonic and forage species. The last group ("EPI group"), was primarily made of large epipelagic species.  (Table 1) and potential prey stable isotope values (not accounting for TDF) show that δ 13 C values of shark muscle were similar to the EPI values and generally higher than those of the prey other groups, whereas δ 15 N values of both tissues were higher than those of EPI, but similar to DVM values (Fig. 2, Table 1). When muscle was adjusted down one trophic level to account for TDF, adjusting the consumer to resemble potential prey, using TDFs of 1.7 ± 0.5 for δ 13 C and 3.7 ± 0.5 for δ 15 N 53 , δ 13 C values align with DVM and MESO groups, whereas δ 15 N values are intermediate to these two groups. Mixing model results for muscle suggest that Cookiecutter sharks primarily feed upon DVM and MESO prey (or prey with isotopic composition similar to these groups) (Fig. 3). The estimated median proportional contribution of MESO prey was 55% (40-68% 95% credible intervals), while DVM prey was 39% (22-53%) and EPI was 5% (0-2%).
Because there is a lack of TDFs for elasmobranch liver from controlled feeding studies, we used liver samples that had been archived from a previous study that characterized TDFs for different Leopard Shark (Triakis semifasciata) tissues 53 . We estimated TDFs for Leopard Shark liver fed on diets of squid (Doryteuthis opalescens) or Tilapia (Oreochromis aureus) (see Supplemental Information for full details on TDF estimates). Leopard Sharks fed squid had a Δ 13 C squid (TDF) of 2.7‰ ± 0.5 and a Δ 15 N squid of 3.5‰ ± 0.7 with error propagation. Sharks fed tilapia had a Δ 13 C Tilapia = 2.4‰ ± 0.9 and Δ 15 N Tilapia = 4.4‰ ± 0.4 with error propagation. Liver mixing model results indicated increased use of MESO prey relative to muscle using TDFs from both squid and Tilapia. The estimated median proportional contribution of MESO prey was 66% (47-82% credible intervals), while DVM prey was 23% (4-39%) and EPI was 9% (0-38%) using squid TDFs. The estimated median proportional contribution of MESO prey was 82% (65-95%), while DVM prey was 9% (1-24%) and EPI was 7% (0-27%) using Tilapia TDFs. The major difference between muscle and liver mixing model results was that liver, which has a shorter incorporation rate than muscle, showed higher contribution of MESO prey and reduced DVM relative to muscle for both squid and Tilapia TDFs. Table 1. Mean (SD) δ 13 C and δ 15 N for Cookiecutter shark tissues (bulk: lipid and urea extracted samples), and tissues adjusted down one trophic level to account for trophic discrimination factors (TDF). Liver samples were adjusted based on TDFs from Tilapia (T) and squid (S) diets (see "Methods").  Fig. S5).
Habitat type was also the significant driver of liver FA profiles (main test P(MC) = 0.001, all pairwise tests P(MC) values < 0.001, Fig. 4, Appendix S1: S3, Table S3, S5) but unlike the muscle (Table S2), liver FA profiles were most similar to deep sea demersal, followed by the pelagic and then the deep sea chondrichthyans (SIMPER dissimilarity scores 17.48, 20.63 and 23.23 respectively). Cookiecutter sharks were most like Little Gulper Sharks, Centrophorus uyato, and Birdbeak Dogfish, Deania calcea. As with the muscle profiles, Cookiecutter livers were extremely high in 18:1ω9, and relatively high levels of 16:0 and 16:1ω9, which drove the similarities between the Cookiecutter and Gulper Sharks and Birdbeak Dogfish (Fig. 4).
Fatty acid analysis: comparison to potential prey items. Each  Environmental DNA (eDNA). Of the 14 total sharks, 10 sharks had target DNA from gastric fluid successfully amplified and sequenced using eDNA metabarcoding. Visible tissue of prey was present in three out of these 10 stomachs, and these tissue samples were Sanger sequenced (Fig. S6). After merging paired-end reads,  respectively; these reads were annotated as human, chimeric sequences (i.e., artifact sequences formed during PCR when two or more biological sequences incorrectly join) or were unassigned based on no matches to our 12S database. The reads from the stomach samples and positive controls were rarefied to 100,000 per sample to account for unequal sequencing depths. All subsequent results presented below are based on the rarefied dataset. The three replicates of the Swordfish control produced a total of 299,990 out of 300,000 reads assigned to Swordfish, indicating negligible cross contamination across the samples. Sequencing of the mock community replicates identified 10 out of the 10 mock community taxa; there were no more than 5 reads assigned to taxa not present in the mock community. A total of 79,068 reads could not be taxonomically assigned to our 12S database or to a taxonomic rank within MEGAN, however these unassigned results represent just 2.8% of the rarefied reads. Of the 2,800,000 total reads in the rarefied dataset, 1,186,484 (42.4%) were assigned to Cookiecutter shark.  Table 2. SIMPER output of relative dissimilarity scores between Cookiecutter sharks (CCS) muscle FA profiles and prey group FA profiles. The top two FAs driving the dissimilarity are ordered by contribution, with the group containing the highest amount indicated in parenthesis. www.nature.com/scientificreports/ We identified three prey in the stomach samples of the 10 sharks that were processed: Ariomma sp. (genus of deep-water fish), Cololabis saira (Pacific saury) and the tribe Thunnini (tunas). These three taxa accounted for 93.2% of the rarefied reads when reads assigned to Cookiecutter shark are excluded. Within the 10 stomachs, Ariomma were detected in five, Cololabis were detected in two, and Thunnini were detected in four individuals (Table 3). Sanger sequencing of the tissue contents present in three stomachs confirmed the presence of Cololabis and Thunnini. The maximum number of prey detected in any stomach was two. Tunas were annotated to the tribe level of Thunnini because MEGAN was unable to identify the resolution at the genus level due to similarities between genera in the short fragment amplified. However, Thunnus and Auxis, two genera within Thunnini, were identified in the dataset and both these genera have geographic distributions within the Hawaiian Islands.

Discussion
Using an integrative approach that incorporated multiple biochemical tracers, we gained novel insights into the habitat and trophic ecology of the Cookiecutter shark. Our results indicate that small, micronektonic prey are an important resource for these sharks in the Central Pacific. Surprisingly, our results suggest that large epipelagic prey were less important to the sharks in this study, despite the common belief that they are dominant prey based on bite mark observations. This discrepancy suggests that the high observability of large epipelagic prey  Table 3. Heatmap of prey identified in the stomachs of Cookiecutter sharks (CC) via eDNA metabarcoding of gastric fluid contents. Data for a given stomach sample (e.g. CC1 is Cookiecutter shark one) is presented as the mean read count from the rarefied dataset. Note that while there was not sufficient DNA for metabarcoding of CC11, a tissue sample from this stomach was Sanger sequenced and identified as Thunnini. † Confirmed by Sanger sequencing of tissue contents. www.nature.com/scientificreports/ to humans may have skewed our understanding of their importance to Cookiecutter sharks. While biochemical tracers generally lack the resolution to identify particular prey species, they were able to identify different functional groups with a high degree of concordance between stable isotope, fatty acid, and eDNA results. Overall, our results indicate that Cookiecutter sharks have a broad diet and have a unique ecological role in pelagic ecosystems, feeding on consumers across all trophic levels, ranging from large, apex predators to small, low trophic level species, in particular species that overlap with the depth distribution of the sharks throughout the diel cycle. We also found evidence of a potential seasonal shift in diet as well as a potential size-based shift in diet and/or habitat, where sharks may start to vertically migrate once they attain a larger size (Fig. 6). This integrative biochemical tracer approach expanded our knowledge of Cookiecutter shark trophic ecology, despite the limited sample size. Only three sharks (20%) had any stomach contents that could be identified using traditional techniques, and those stomach contents were very limited. By using multiple tracers, we were able to use every animal to derive insight into their trophic ecology, increasing our available sample size by 400% in comparison to if we only used traditional stomach content analysis. In addition, because we used different tissues and tracers, we were able to obtain different temporal scales of data that allowed us to make inferences that would not be possible using traditional stomach content analysis without extensive sampling across time, which is generally prohibitively difficult in these ecosystems. eDNA metabarcoding of gastric fluids was able to detect the presence of prey DNA in shark stomachs, including prey with no identifiable or partially digested tissue material (i.e., Ariomma). It is likely that the eDNA results reflect a snapshot of the most recent diet. These results provide new insight into the diet of Cookiecutter sharks, in particular by identifying Ariomma and Cololabis as prey for this species. Sanger sequencing of partially digested material retrieved from inside two shark stomachs supported the eDNA findings. Our results also demonstrate that eDNA can be extracted and sequenced from stomachs frozen for extended periods of time to accurately identify prey consumed prior to death.

Stomach Sample
Importance of different prey groups to Cookiecutter sharks. The trophic ecology of Cookiecutter sharks will be influenced by the degree of spatial and temporal overlap between the sharks and potential prey. Cookiecutter sharks can potentially exploit mesopelagic and vertically migrating species throughout the diel cycle (depending upon the vertical migration of the sharks and the prey), but will only be able to exploit epipelagic species during the night when sharks are near the surface 41 .
Stable isotope, fatty acid, and eDNA results all suggested that Cookiecutter sharks feed more on small prey, whether mesopelagic or epipelagic, than previously thought. Over longer time scales, approximately half their www.nature.com/scientificreports/ diet is derived from small prey and FAA indicates mesopelagic and deep-sea demersal foraging (Meyer et al. 11 ). Large epipelagic species, though often seen with Cookiecutter shark bite marks in the Central Pacific and across the globe 44,45,54,55 , appear to comprise a relatively small proportion of the overall diet (< 10%) for individuals in this study. However, due to the limited number and geographical extent of samples used in this study, it remains unclear as to whether these patterns will hold in other populations of this globally distributed species, as dietary plasticity is increasingly being acknowledged as a common trait to consumers 56 .
Muscle FAA results indicated that sharks consume large mesopelagic and diel vertically migrating species, small meso-and epipelagic species, and marine mammals, which matches SIA results. These results were also concordant with eDNA results, as four sharks (30%) had Ariomma sp. (a small mesopelagic fish) DNA in their stomach, and one shark (8%) had Cololabis saira (a small epipelagic fish). Thunnus DNA was also found in two sharks (15%), but we were not able to resolve it to species level, so it's unclear whether it was a species that used mesopelagic (e.g., Bigeye, Albacore) or epipelagic (e.g., Yellowfin, Skipjack) resources. However, it is important to note that primers used in the eDNA study could not identify invertebrates. FAA results indicated the Cookiecutter sharks were most similar to Somniosus spp., based on high levels of 18:1ω9 and 16:1ω7, which suggest consumption of marine mammal blubber ( 57,58 ). Somniosus spp. are known to opportunistically feed upon marine mammals [59][60][61] , further suggesting that marine mammals are prey for larger Cookiecutter sharks. The exceptionally high 18:1 and 16:1 may also come from the skin and subdermal fat from teleost and elasmobranchs (Every et al. 33 ) as Cookiecutter sharks remove relatively high proportions of these tissues owing to their unique ectoparasitic hunting. While the FAA suggested that large epipelagic species and cephalopods were the most dissimilar to Cookiecutter sharks, which supports SIA results, the similarity between Cookiecutter and Mako sharks suggests squid, and a mix of epipelagic and mesopelagic teleosts may be important prey items (Preti et al. 62 ).
Our finding that large species that undergo diel vertical migration are an important prey resource for Cookiecutter sharks is also consistent with published data. Papastamatiou et al. 44  Mixing models are inherently sensitive to the selection of sources and TDFs 63 . We used the approach described by Smith et al. 64 to investigate how well the TDFs and three sources we used fit our consumer data. We found that our muscle data were well bounded by the mixing polygons, with only one animal falling outside the 95% contour of the mixing polygon, suggesting that our sources and TDFs were generally appropriate. For liver, four animals fell outside the 95% contour of the mixing polygon, suggesting a poorer fit, but one that was still generally suitable. It is possible that we have not accounted for all the potential sources given the diversity of potential prey items and lack of stable isotope data for many potential pelagic prey, in particular mesopelagic species. It appears that we may be missing a source with more depleted δ 13 C than any of the current sources. In addition, while the TDFs appear to be generally suitable, it is unclear how suitable muscle and liver TDFs from a small coastal Carcharhiniform shark will be for a mesopelagic Squaliform shark. However, until controlled feeding studies are done on a species similar to Cookiecutter sharks, we are restricted to using available TDFs. Furthermore, Hussey et al. 65 found that TDFs narrowed with increasing trophic level, yet how this might impact the Cookiecutter shark, which feeds across all trophic levels, is unclear. While there are important caveats associated with our analyses, the strong degree of concordance between the different approaches provides support for our model selection, analyses, and interpretations.
Evidence for seasonal dietary changes. Stable isotope mixing model results for liver suggest that recent diet of the sharks was dominated by MESO prey, a result that was consistent using either set of liver TDFs. Because liver represents the diet over the prior weeks to few months prior to death, it suggests small meso-and epi-pelagic prey were the primary prey for Cookiecutter sharks compared to the longer time frame integrated by muscle that includes substantial contributions from DVM prey. This interpretation is supported by FAA. The FA profile of Cookiecutter shark livers, which will reflect more recent feeding (< 3 and up to 12 weeks), is most similar Little Gulper Shark (Centrophorus uyato) and the Birdbeak Dogfish (Deania calcea), which are both mesopelagic and deep-water species that primarily feed upon mesopelagic micronekton. Interestingly, the low level of PUFAs (22:6ω3 and 20:5ω3) and high ω9 FAs (in this case 44.48 ± 3.64 18:1ω9) observed in shark liver have been used as indicator of essential FA deficiency in young bull sharks 32 , and may indicate poor nutrition. Liver δ 15 N and FAA values suggest sharks recently fed more on small mesopelagic species, which have a lower nutritional value than larger species such as marine mammals and other large diel vertically migrating and epipelagic species 66,67 . If true, it may suggest that this seasonal shift in diet is necessitated by a lack of availability of larger, more energetically rich prey, either as a result of seasonal movements of sharks or their prey. Indeed, Papastamatiou et al. 44 observed an apparent seasonal shift in frequency of fresh bites on fishes and suggested this may reflect a seasonal shift in the distribution of Cookiecutter sharks. However, this seasonal diet shift to less nutritional prey is speculative and more research is needed to address this hypothesis.
Shift in trophic ecology with length. We found no significant relationship between the shark length and stable isotope composition but recognize the size range of sharks in this study was limited, particularly for www.nature.com/scientificreports/ smaller individuals. The variability in δ 15 N values among larger sharks (> 30 cm) may reflect individual variability in diet or habitat use, with different individuals exhibiting differential reliance on mesopelagic vs epipelagic resources. However, the smallest shark had one of the highest δ 15 N and lowest δ 13 C values, which suggests increased use of mesopelagic resources. There is an increase in baseline ecosystem δ 15 N values with increasing depth so the elevated δ 15 N values could indicate increased reliance upon a mesopelagic suspended particulate based food web as opposed to prey directly connected to surface productivity [68][69][70] . It is possible that smaller individuals feed entirely in the midwater food web without diel vertical migration unlike larger sharks. This result is concordant with FAA that found high levels of 20:5ω3 in the smallest shark, which indicates feeding in deeper, colder waters in deep-sea food webs and on lower trophic level prey ( 71-73 ; Meyer et al. 11 ) including cephalopods ( 28,71 ). Estimated size at birth is 14-15 cm 74 , so it is possible that the small shark's stable isotope composition was influenced by maternal contributions 75,76 . Given that this shark was larger (17 cm) than the estimated size of birth and was isotopically distinct from the larger sharks, this interpretation seems less likely. While this result is interesting, it remains speculative due to the small sample size. While our sample size is limited, SIA and FAA results suggest that diel vertical migration in Cookiecutter sharks may develop as they grow, with small sharks remaining in and feeding in mesopelagic habitats throughout the diel cycle. Alternatively, they may vertically migrate but be restricted to consuming small mesopelagic prey that also vertically migrate. High levels of 18:1ω9, and 18:1ω7 also indicated that large sharks had the highest potential consumption of marine mammal blubber (Waugh et al. 58 ) and that sharks < 30.5 cm (n = 5) had significantly different FA profiles than sharks > 30.5 cm (n = 9). We hypothesize that these patterns may indicate that there is an ontogenetic shift in diet or habitat. This interpretation is consistent with Papastamatiou et al. 44 who found that Cookiecutter shark bites on fishes in the Honolulu Fish Market were almost entirely from sharks estimated to be > 25 cm TL based on measured bite sizes. Papastamatiou et al. 44 also reported that Swordfish had significantly higher occurrence of bites from larger sharks (> 40 cm TL), which they suggested may reflect a size or sex based shift in foraging ecology, as male sharks are rarely > 40 cm TL 41 . Unfortunately, we did not get eDNA data from sharks < 30 cm, so eDNA results could not provide any insight into this potential length-based shift.
Size-based shifts in diet or habitat are ubiquitous in chondrichthyan fishes, but poorly understood in mid-and deep-water species 77 . Sex and size-based segregation have been noted in other deep-water squaliform sharks, including species in the family Centrophoridae (e.g. Cetrophorus, Deania), Somniosidae (e.g. Scymnodon), and Etmopteridae (e.g. Etmopterus, Centroscyllium) [78][79][80][81][82] . In particular, Etmopterus princeps and Centroscyllium fabricii are noted to show an ontogenetic shift in habitat, with younger sharks having a deeper distribution than older ones 81 . In addition, some deep-water species may move to deeper waters and give birth 3 , suggesting that Cookiecutter sharks may be born at the deeper part of their depth distribution. While we provide some evidence that Cookiecutter sharks exhibit a length-based shift in habitat and/or diet, our biochemical tracer results do not allow us to further differentiate between these factors. Small sharks could vertically migrate but be constrained to consuming small prey due to gape limitations, though gape limitation seems unlikely given their dentition, or they may indeed remain at depth and feed primarily on mesopelagic micronekton. This size threshold may also reflect change in swimming ability, with larger sharks being able to swim more rapidly to feed upon larger, faster moving prey.
Understanding the trophic ecology of this poorly known predator sheds important insight into its ecological role in epi-and mesopelagic ecosystems. Understanding the ecological role of this globally distributed predator, and other similar mesopelagic predators, is fundamental to understanding food web dynamics in pelagic ecosystems and necessary for effective management. How relevant our findings are to other populations of Cookiecutter sharks around the world will remain unknown until additional studies are conducted. We found that a multi-tracer approach increases our ability to make ecological inferences, especially with small sample sizes, when studying these types of species and food webs.

Methods
Cookiecutter sharks (n = 15; 9 M, 6 F; mean TL = 33.5 cm ± 6.6 SD; min TL = 17.1 cm, max TL = 43.1 cm, mean mass = 213 g ± 113 SD) were collected off Hawaii in August and September of 2013 by the Monterey Bay Aquarium onboard the NOAA ship R/V Oscar Sette using a Isaacs-Kidd midwater trawl and Cobb trawl fishing primarily at depths between 60 and 150 m, but down to 600 m. While the size at maturity is not well characterized, it has been reported to be ~ 36 cm for males and ~ 39 cm for females 41 . Based on these criteria, three males (43%) and two females (33%) were considered to be mature. Sharks were held in an aquarium on board the vessel until they expired, generally after several hours, at which point they were stored in a -20 °C freezer until they were processed. Sharks were thawed and their stomachs removed for eDNA analysis, and liver samples and muscle samples from the dorsal musculature behind the gills were collected for SIA and FAA. Liver and muscle samples were stored frozen at − 20 °C until they were analyzed. All animals were collected under the auspices of the Monterey Bay Aquarium following the Association of Zoos and Aquariums Accreditation Standards and permits and ethical approvals were obtained from the Monterey Bay Aquarium's Animal Welfare and Research Oversight Committees. All methods are reported in accordance with ARRIVE guidelines.
Muscle and liver samples were collected because they have different isotopic incorporation rates and as a result integrate dietary information over different periods of time. Muscle has a slower incorporation rate and reflect a shark's diet integrated over long periods of time, whereas liver, being more metabolically active, has a more rapid incorporation rate and will integrate diet over shorter time frames 20,83,84 . There are no incorporation (turnover) rate studies for mesopelagic sharks or fishes, but we use the allometric relationship from Weidel et al. 85 , which was developed for teleosts but has been reported to be applicable to sharks 20 . The estimated half-life of δ 13 C based on the mass of sharks in this study is an average of ~ 78 days. If true, muscle will reflect the previous 3 -> 9 months (e.g., 50 to > 88% turnover), and liver would be a shorter time frame. Further, lipids are more Scientific Reports | (2021) 11:11809 | https://doi.org/10.1038/s41598-021-89903-z www.nature.com/scientificreports/ metabolically active than proteins, such that FAs reflect diet at shorter time frames than stable isotopes (weeks to months) (Beckmann et al. 30 ). Muscle FA samples begin to reflect a change in diet within 3 weeks and are fully transformed by 16 weeks (Beckmann et al. 30 ). As with isotopic incorporation, liver FAs have a quicker incorporation rate, reflecting diet from within 3-12 weeks (Beckmann et al. 30 ), though the rate likely varies across species. Because urea and lipid content influence the isotopic composition of elasmobranch tissue 10,86,87 , all shark tissues were both lipid and urea extracted following Kim and Koch 87 . Approximately 10 mg of liver was freeze dried, then lipid extracted three times. Each round of lipid extraction consisted of 20 mL 2:1 chloroform:methanol solution, 10 min of ultrasonication, and 24 h on a shaker table. After lipid extraction, liver samples were urea extracted with three rinses in deionized water with 10 min of ultrasonication. Samples were analyzed in triplicate at the Stable Isotope Lab of University of California Merced on a continuous flow isotope ratio mass spectrometer (ThermoScientific Delta V Plus) coupled to an Elemental Analyzer (Costech 4010) with a Conflo IV. Isotopic composition is expressed using standard δ notation, using Vienna Pee Dee Belemnite as the standard for δ 13 C values and atmospheric AIR nitrogen as the standard for δ 15 N values. The liver samples were run with three reference materials (Acetanilide [Costech] and two glutamic acids [USGS 40 and USGS 41a]) and analytical precision for δ 13 C and δ 15 N values was < 0.1‰ and 0.2‰, respectively. TDFs were estimated based on the average difference of captive feeding study shark liver and prey (i.e., δ 13 C liver -δ 13 C prey ) using the previously published isotopic composition of squid and tilapia (Kim et al. 20,53 ). The TDF standard deviation is based on the propagation of error from shark and prey δ 13 C and δ 15 N values.
Stable isotope analysis (SIA). Lipid was extracted using a 2:1 chloroform:methanol solution and urea was extracted using deionized water rinses 87 . Samples were then lyophilized and homogenized in a Spex/Cer-tiPrep 5100 mill. Approximately 500 μg of tissue were weighed into tin boats and analyzed at the Stable Isotope Laboratory at UCSC using an Elemental Analyzer (Carlo Erba) coupled to a continuous flow isotope ratio mass spectrometer (Thermo-Finnigan Delta XP). Isotopic composition is expressed using standard δ notation, using Vienna Pee Dee Belemnite as the standard for δ 13 C values and atmospheric AIR nitrogen as the standard for δ 15 N values. Analytical precision, based on an internal lab standard (Pugel), was < 0.2‰ for δ 15 N and < 0.1‰ for δ 13 C values.
The stable isotope composition of potential prey of Cookiecutter sharks around Hawaii was extracted from the literature 68,[88][89][90][91][92][93] or through sampling fishes (14 species, ~ 10 samples/species) at the Honolulu Fish Market (Appendix S1: Table S1, Fig. S1). We only included potential prey values from the literature that accounted for lipid content through either chemical extraction or mathematical adjustment based on tissue C:N, or did not require lipid extraction due to low lipid content of tissues (See Supplemental Materials for more detail). All potential prey species sampled as part of this study were lipid extracted using a 2:1 chloroform:methanol solution. When there were multiple stable isotope values for a particular prey group, we aggregated them into a single mean and standard deviation (SD) value by resampling 1000 values from each prey δ 13 C and δ 15 N distribution in the group and combining them into a cumulative distribution for that prey group. We grouped potential prey into three groups using k-means clustering. We limited the number of prey groups to three because we were only using two tracers (δ 13 C and δ 15 N), and stable isotope mixing models are most effective when the number of sources used is n + 1 (n = number of tracers) 63 . These groups generally clustered based on their habitat type. The DVM group (mean ± SD δ 13 C: − 17.9 ± 1.0, δ 15 N: 12.5 ± 0.5) was primarily comprised of large species that exhibit diel vertical migration as well as marine mammals, though there were a few large meso-and epipelagic species in this group as well. The MESO group (δ 13 C: − 17.9 ± 0.9, δ 15 N: 7.3 ± 0.7) was entirely comprised of small mesopelagic and epipelagic micronektonic and forage species. The EPI group (δ 13 C: − 17.1 ± 1.0, δ 15 N: 10.5 ± 0.5) was primarily made of large epipelagic species, though there were some mesopelagic and vertically migrating species in this group as well. We used the Bayesian stable isotope mixing models MixSIAR 94 , using the isotopic composition (δ 13 C, δ 15 N) of the sharks and their potential prey groups (i.e., DVM, MESO, EPI as described above), to estimate the relative importance of different prey groups to sharks over shorter time frames using liver and longer time frames using muscle. The model was run with 300,000 Markov chain Monte Carlo (MCMC) simulations (200,000 burn-in) and showed good convergence based on Gelman-Rubin (all variables < 1.01) and Geweke (no scores outside ± 1.96 in any chain) diagnostics. Uninformative priors were used.  Kartikasari et al. ( 96 ). FA results are expressed as a proportion of the total identified compounds and only those with means > 0.1% were included in the subsequent statistical analyses.
We investigated whether Cookiecutter shark muscle and liver FA profiles differed by sex and length, grouping animals into four size classes: extra small (< 29 cm TL, n = 1), small (29-30.5 cm TL, n = 4), medium (30.6-35 cm TL, n = 4), and large (> 35 cm TL, n = 5). We also compared Cookiecutter FA profiles to those from chondrichthyans occupying deep sea (meso-and baythpelagic > 300 m deep), deep demersal (> 300 m deep associated with the benthos), and epipelagic habitats (< 300 m deep) (as in Meyer et al. 11 98 were used to assess if and how FA profiles differed across biological factors (sex and size), habitat type, and prey groups. Permutational analysis of variance (PERMANOVA) tests with Monte Carlo simulations [denoted as p(MC)] were run on Bray-Curtis similarity matrices calculated from the square-root transformed profile data to determine if categorical groups were significantly distinct (p(MC) < 0.05). Where groups were significantly different (determined using PER-MANOVAs), similarity percentage analysis (SIMPER) was used on the square root transformed data, to quantify which groups within a factor were most dissimilar, and which FAs were driving pairwise dissimilarities. SIMPER analysis was also used to quantify the dissimilarity between the Cookiecutter sharks and individual shark species. Principal Coordinates Analysis (PCO) was used to visualize the differences between groups, with FA correlations > 0.7 overlaid as vectors on the PCO, indicating direction and magnitude of the correlation.
Environmental DNA (eDNA). For each shark (n = 14, one shark's stomach contents was not examined using eDNA), stomachs were cut open using sterile razor blades. Contents were then squeezed into a sterile weigh boat and homogenized using a razor blade. If prey tissue or bone was present, we sampled the material and extracted separately (n = 4 from 3 sharks). Gastric contents were transferred to a sterile 5 mL falcon tube and vortexed for 3-5 min. Between the treatment of each stomach sample, scissors and tweezers were sterilized with the following steps: wash with 0.05% bleach solution, submerge in 96% ethanol, rinse in RNAse Away and then flame.
DNA extraction and PCR amplification. DNA extractions were performed using the DNeasy Blood and Tissue Kit (Qiagen, USA) following the manufacturer's protocol with minor modifications (see Supplemental Methods). Of the 14 shark stomachs, 11 stomachs had sufficient amounts of starting material to proceed forward. Gastric fluid contents from the 11 stomachs were extracted in triplicate; four samples were taken from visible tissue present in three stomachs and were extracted in duplicate. Samples were randomized prior to extraction. DNA concentrations were determined using the Qubit dsDNA HS Assay (Invitrogen, CA, USA). In addition to the stomach samples, we also included negative controls (i.e., extraction blanks) and positive controls in our study (see Supplemental Methods).
We used a 106-bp fragment of the mitochondrial 12S rRNA gene 99 to assess the vertebrate diversity present in the samples. Primer sequences were F-5′ ACT GGG ATT AGA TAC CCC and R-5′ TAG AAC AGG CTC CTC TAG . This primer set has previously been validated in a seawater mesocosm study 100 and used in other eDNA metabarcoding studies [101][102][103] . These primers target vertebrates (e.g., fishes and marine mammals) but are less likely to amplify DNA from cartilaginous fishes (e.g., sharks, rays) or sea turtles due to base pair mismatches. The primer set does not amplify invertebrates. Unique tags were added to the 5' ends of the primers to allow for assignment of sequence reads to the correct sample during bioinformatic processing. Before amplification with the tagged primers, all samples were checked for inhibition (see Supplemental Methods). If deemed inhibited, the sample was removed from further processing (n = 1 shark). PCR reaction chemistry and cycling parameters were the same as Port et al. 102 (see Supplemental Methods). Each DNA extract was amplified in triplicate, along with no template controls (NTCs) using molecular-biology-grade water in lieu of DNA template for each sample.
Prior to our main study, we designed a blocking primer to inhibit amplification of Cookiecutter shark DNA. We tested the blocker with the 12S primers against DNA extracts from Cookiecutter shark and an in-house fish tissue repository as well as a Cookiecutter shark stomach sample and sequenced the amplicons. While the blocker was effective at inhibiting amplification of Cookiecutter shark DNA, there was differential amplification of some species (data not shown). Due to this primer bias, we did not employ the blocker primer going forward. Note that in these preliminary sequencing runs using Cookiecutter shark stomach samples, we detected Opah, Lampris guttatus, as a prey item when the blocker was both included and excluded (Opah was not detected in the main study).
Sanger sequencing of partially digested tissue samples. PCR products from tissue samples (n = 4 extracts from 3 tissue samples) were commercially Sanger-sequenced with an ABI 3730 × 1 sequencer (Elim Biopharmaceuticals, Inc., Hayward, CA, USA) in both the forward and reverse direction using the same primers as the initial PCR setup. The resulting sequences were trimmed, edited and aligned using Geneious Pro v. 7.1 (Biomatters, NZ) and then compared against the nonredundant GenBank sequence library using the nucleotide BLAST search engine provided by the National Center for Biotechnology Information (NCBI).
Next generation sequencing of eDNA in gastric fluid samples. Tagged PCR products for gastric fluid samples yielding sufficient target DNA (n = 30 across 10 sharks) were pooled in equimolar concentration (20 ng DNA per sample) along with controls (n = 8) to create a single library. Note that of the 14 total sharks, three sharks did not have enough starting material for DNA extraction and amplicons from one shark were not visible on a gel. The library was prepared for sequencing using the KAPA low-throughput library prep kit with real-time library amplification protocol (KAPA Biosystems, USA) and sequenced on an Illumina MiSeq platform (250 bp, paired-end) at the Stanford Functional Genomics facility using a 20% PhiX spike-in control (see Supplemental Methods for more information). www.nature.com/scientificreports/ Bioinformatics and data analysis. We processed the Illumina sequence reads with banzai, a custom Unix-based script (banzai, https:// github. com/ jimmy odonn ell/ banzai), following the approach described in Port et al. 102 .
Banzai calls third-party programs [104][105][106] to move from raw sequence data to a quality-controlled dataset of sequence counts from operational taxonomic units (OTUs). Sequences were demultiplexed and only retained if the tag added during amplification was present on both the forward and reverse reads to eliminate samples with tag jumping 107 . OTUs were compared to a local nucleotide database containing mitochondrial sequences from NCBI using BLAST+ 108 . This database-deposited in the Dryad Digital Repository-totaled 12,689 sequences and included the complete mitochondrial genomes as well partial 12S rRNA gene fragments of bony fishes (Actinopterygii), cartilaginous fishes (Chondrichthyes), true seals (Phocidae), sea lions (Otariidae), whales (Cetacea), marine dolphins (Delphinidae), sea otters (Enhydra) and birds (Aves) (sequences downloaded September 2014). We also Sanger sequenced Cookiecutter shark tissue using the 12S primers and deposited this 105 bp sequence in the local database. Default blast parameters were used with a few exceptions (see Supplemental Methods). We removed OTUs classified as non-marine vertebrates (e.g., Homo sapiens, Gallus gallus) as well as unassigned sequences and reads annotated as chimeras.
To account for uneven sequencing depths across samples, we rarefied each sample (stomachs and positive controls) to 100,000 reads using the "rrarefy" function in the R package vegan 109 and used the positive and negative controls to create conservative thresholds for filtering for false positives (see Supplemental Methods). Sequence counts were analyzed in terms of presence/absence as well as relative abundance. Given biases associated with PCR amplification (e.g. sequencing, rarefaction, etc.) our analysis was semi-quantitative. Rarefied read counts were binned on an order of magnitude scale (i.e., < 100, < 1000, < 10,000 and < 100,000).