Evaluating the use of stable isotope analysis to infer the feeding ecology of a growing US gray seal (Halichoerus grypus) population

Gray seals (Halichoerus grypus) have been rapidly recolonizing the Northeast US coast, eliciting concern from the fishing industry. However, the ecological effect of this recovery is still unknown and as such, research is needed to better understand how the diet composition of gray seals in US waters will contribute to the ecological impact. While previous research on seal diets has focused on the analysis of hard prey remains, stable isotope analysis presents an alternative method that can be used to describe marine mammal diets when direct observation is impossible. To address this issue, we used stable isotope analysis of gray seal pup vibrissae and lanugo from Monomoy Island, Cape Cod, MA during the 2015/2016 winter breeding season to estimate adult female diet composition during pregnancy. Stable isotope mixing models (SIMM) suggested adult female gray seals were consuming greater amounts of cephalopod prey and less sand lance than previously indicated from analysis of hard prey remains. However, using SIMMs to estimate the diet composition of gray seals remains difficult due to the large number of isotopically similar prey species and uncertainty in tissue-specific, stable isotope trophic enrichment factors. Even so, by combining prey sources into ecologically informative groups and integrating prior information into SIMMs it is possible to obtain additional insights into the diet of this generalist predator.


Introduction
The gray seal's (Halichoerus grypus) return from near extirpation to the Northeastern United States has been celebrated by conservationists but has reignited controversy with some in the fishing industry [1][2][3][4][5]. This is because there is the potential for these seals to be predators of many fish species in economically significant fisheries [6,7]. Further, their recolonization comes at a critical time in the ecosystem, as the northwest Atlantic is experiencing rapid warming [8]. The gray seals in the US differ in distribution and foraging region from the nearby a1111111111 a1111111111 a1111111111 a1111111111 a1111111111 Science Center. We collected vibrissae from 25 weaned pups from Monomoy Island (41.5500N , 70.0000˚W) during January 2016. A single vibrissae was clipped from the base of the muzzle on the right side of each sampled pup. Of those 25 pups, 19 also had lanugo samples collected. Prey samples were acquired through the 2016 NOAA spring bottom trawl off the coast of Massachusetts. [5,9]. Pups were randomly sampled from stage III and stage IV pups (i.e. weaned pups that had yet to fully shed their lanugo) present on the island. Two to eight individuals each of the ten most important prey species were acquired based on previous scat/ stomach diet analyses of Northwest Atlantic gray seal populations [9,10,12] (Table 1). Specifically, these ten prey species represented more than 90% of dietary biomass in analyses of gray seal scat and stomach samples in New England [9]. As gray seal scat data is likely indicative of inshore foraging and stomach contents of offshore foraging, the selected species represent the dominate prey available to seals in this region, and it is unlikely a major prey item is missing from our analysis [9].

Sample preparation
We used pup vibrissae and lanugo tissues developed in-utero as proxies of maternal diet. Vibrissae begin growing in gray seal fetuses between 2-3 months of active gestation and continue to develop over the next 6-7 months until birth [35]. For pups born in December-January, vibrissae therefore reflect summer/fall foraging seasons. Pups are born with vibrissae roughly 2-3 cm long and grow at a rate of .0025cm/day to about 5-7cm during the lactation period [35][36][37]. The maximum growth rate recorded in phocid pups was 0.087 cm/day in the bearded seal (Erignathus barbatus) [36]. Therefore, for a 4-6 week old weaned pup with 5-7 cm long vibrissae, a maximum of 3.6 cm of new vibrissae growth could have occurred since birth, leaving the distal 2-3 cm to reflect in utero growth [16,36]. Thus for all vibrissae samples, we saved a distal 2 cm section for analysis. In gray seals, lanugo begins growing around 3-4 months of gestation and grows up until birth [35,38,39]. Like vibrissae, lanugo reflects a summer/fall foraging season. Prior to isotopic analysis, we rinsed vibrissae and lanugo samples in a 2:1 chloroform:methanol solution for 24 h and air dried for a further 24 h. We freeze-dried, and pulverized roughly 1.5g samples of muscle tissue from each of the potential prey species. Individuals for all species were sampled two times except elasmobranch species, which were sampled four times. Following the protocol of Kim and Koch [40] we performed a lipid extraction (LE) using a 2:1 chloroform:methanol on one of the samples. Both lipid extracted and non-lipid extracted (NLE) muscle tissue for all prey species were processed for SI analysis. Elasmobranch species underwent an HCl and non-HCl treatment following the protocol of Kim and Koch [40] to remove cartilaginous tissue from the muscle that could interfere with the δ 13 C signature. For most prey species LE treatments were used for δ 13 C analysis and NLE treatments for δ 15 N analysis, except in elasmobranch muscle where HCl and LE treated muscle was used for δ 13 C analysis.

Stable isotope analysis
We flash-combusted (Costech ECS4010 elemental analyzer) 2 cm vibrissae sections (0.3-1.6 mg) and approximately 0.8 mg of lanugo and homogenized prey muscle tissue samples loaded into tin cups and analyzed for carbon and nitrogen isotopes (δ 13 C and δ 15 N) through interfaced Thermo Finnigan Delta Plus XL continuous-flow stable isotope ratio mass spectrometer.
Raw δ values were normalized on a two-point scale using glutamic acid reference materials with low and high values (i.e. USGS-40 and USGS-41). Sample precision based on repeated reference material was 0.2‰ and 0.1‰ for δ 13 C, and δ 15 N, respectively.
Stable isotope ratios are expressed in δ notation in per mil units (‰), according to the following equation: δX = [(R sample / R standard ) -1] Á 1000. Where X is 13 C or 15 N and R is the corresponding ratio 13 C / 12 C or 15 N / 14 N. The R standard values were based on the Vienna PeeDee Belemnite for δ 13 C and atmospheric N 2 for δ 15 N.

Statistical analysis
We used Pearson's correlation tests to examine the relationship between vibrissae and lanugo SI values and paired t-tests to test for differences in SI values between tissues. Unpaired t-tests were used to test for differences in tissue SI values between sexes. An ANOVA with Tukey's HSD post-hoc comparisons were used to examine differences in the SI values of prey species. All statistical analyses were performed in R (Version 3.3.1).

SIMM analysis
Mixing models cannot be fit blind; background knowledge of consumer diet is necessary to choose appropriate food sources to fit into the model [32,41]. Even with this knowledge, models where the number of contributing sources is greater than or equal to one plus the number of isotopes used are fundamentally undetermined with multiple existing solutions [32,41]. As such, Bayesian SIMMs have been developed to estimate the probability distributions of multiple prey source contributions to a consumer while accounting for the observed variability in prey and consumer isotopic values, TEFs, and elemental concentration [42]. Even so, the discriminatory power of SIMMs generally decreases with the number of prey sources, is strongly influenced by the isotopic separation of prey sources, and can be biased by missing prey sources [32,43]. When the above issues lead to reduced discriminatory power, dietary data from independent sources can be incorporated as informative priors to improve the accuracy and precision of SIMMs [44].
Given the issues raised above we used a priori source-partitioning methods to combine prey species into statistically and ecologically relevant prey groups prior to their incorporation into SIMMs based on their SI values. First, similar to Hindell et al [36], a hierarchical clustering analysis utilizing Ward's method, accounting for variation in both δ 13 C and δ 15 N, was applied to cluster the ten prey species into six prey groups as the discriminatory power of SIMM can decline markedly above six or seven sources [32]. Second, following Edwards et al [45] and Philips et al [41], we further reduced these six prey clusters into three distinct groups based on ecological similarities in habitat and functional group to determine if condensing prey groups increased the discriminatory power of SIMMs predictions, albeit at a lower taxonomic resolution.
Stable isotope mixing models can be highly sensitive to the TEFs employed [46]. As species and tissue-specific TEFs are not currently available for fetal gray seals, δ 13 C and δ 15 N TEFs were chosen from three potential sources: averaged ecosystem TEFs developed from Post [47] derived from a meta-analysis of laboratory and field observations of trophic fractionation across a wide range of ecosystems, a mean of all experimentally determined TEFs for phocid seals [48,26,27], and meta-analysis derived TEFs from determined using the package Stable Isotope Discrimination Estimation in R (SIDER) [49] ( Table 2). The SIDER tool uses a Bayesian linear modeling approach to predict a species TEF value based on their physiology, phylogenetic relationships and experimental variation from published literature [49]. Following the method of Smith et al [50] a Monte Carlo simulation of stable isotope mixing polygons was used to select the most appropriate TEFs [16,42]. Consumer vibrissae and lanugo isotope values were independently fit to mixing model polygons created using prey cluster δ 13 C and δ 15 N values adjusted for each TEF, incorporating both TEF and prey SI value uncertainty, over 1500 iterations. TEFs were chosen for incorporation into SIMMs by selecting the values that provided the best fit (i.e. fewest number of consumer signatures near the edge of the polygon's 95% confidence interval).
Stable isotope mixing models were run using the MixSIAR package [51]. Separate models were run for each tissue (vibrissae and lanugo) and for each prey cluster grouping (six prey sources and three prey sources). Additional variants of the three prey sources models were run using diet composition data from Ampela's [9] analysis of gray seal scats in New England as informative priors as opposed to equally weighted priors to assess the ability of prior information to increase the resolution of SIMM outputs. Priors were adjusted for normalized dietary contributions of 1.4% for squid species, 53.3% for sand lance, and 40.6% for a combination of all remaining 'demersal fishes' in the analysis [9]. Each model ran over three Markov Chain Monte Carlo chains of 1,000,000 iterations, thinned by 500, and with an initial discard of the

Prey stable isotope values
Muscle tissue δ 13 C and δ 15 N values differed among the ten prey species (δ 15 N: F (9, 39) = 50.32, p<0.05; δ 13 C: F (9, 39) = 25.13, p<0.05). Northern shortfin squid (Illex illecebrosus) and longfin squid (Doryteuthis pealeii), had δ 13 C values near the middle of the prey distribution but had the first and third lowest δ 15 N values respectively, isolating them from other prey ( Table 1). Sand lance (Ammodytes spp.) had the lowest δ 13 C values and second lowest δ 15 N values distinguishing it from other prey (Table 1). At the opposite end of the spectrum, red hake (Urophycis chuss) and silver hake (Merluccius bilinearis) both had higher δ 13 C and δ 15 N values than all other species and were not significantly different from one another (p>0.05). Apart from these three distinct groups, the middle of the distribution contained both winter skate (Leucoraja ocellata), with high δ 15 N and low δ 13 C values, and thorny skate (Amblyraja radiata) with conversely high δ 13 C and a low δ 15 N values (Table 1). Also within this group were white hake (Urophycis tenuis), Atlantic cod (Gadus morhua) and winter flounder (Psuedopleuronectes americanus), which had δ 13 C and δ 15 N values near the upper-middle end of the distribution and, except for winter flounder δ 15 N values, were not significantly different from one another (p>0.05).

Cluster analyses and TEF choice
Our hierarchical cluster analysis assigned the ten prey species into six groupings based on their isotopic values: 1) sand lance, 2) northern longfin squid and shortfin squid, 3) red hake and silver hake, 4) winter flounder and winter skate, 5) thorny skate, and 6) Atlantic cod and white hake. Our second grouping method assigned prey species into three groups: 1) sand lance, 2) northern longfin squid and shortfin squid, and 3) a 'demersal fishes' group that included red hake, silver hake, winter flounder, winter skate, thorny skate, Atlantic cod, and white hake. TEF values from Post [47] were chosen as the most appropriate following a quantitative assessment of multiple potential TEF values for both vibrissae and lanugo using stable isotope mixing model polygons (Fig 1). For both tissues, the TEF values from Post [47] had the greatest number of consumer points within the confidence intervals of the mixing polygon created by prey source groups.

SIMM analyses
In general there was broad overlap in the predicted mean diet composition and 95% credibility intervals in all model variants regardless of the consumer tissue of interest (vibrissae vs. lanugo; Table 4). Models that included six prey clusters were poorly resolved with 95% credibility intervals overlapping substantially across all prey sources. Six prey source SIMMs predicted mean diet contribution ranging from 11.5 to 23.5% with a slightly more uniform distribution across species predicted by lanugo relative to vibrissae (Table 4).
Better resolution (i.e precision) was provided by SIMMs that only included three prey clusters. Three prey source SIMMs predicted that approximately half of seal's diets were comprised of other 'demersal fish' relative to squid and sand lance, with little to no overlap in 95% credibility intervals (Table 4). While the relative predicted mean contribution of squid vs. sand lance was consistently higher in three source models using vibrissae, 95% credibility intervals overlapped between these two prey sources in all models. Using scat data as a Bayesian prior to inform the three source models consistently resulted in reductions in the predicted mean contribution of squid vs. sand lance, though 95% credibility intervals still overlapped in these two prey sources (Table 4).

Gray seal and prey stable isotope values
We found no differences in tissue stable isotope values between male and female pups. Evidence from the eastern Atlantic suggests that male pups are born larger than female pups and are more energetically costly to gestating females [52,53]. However, studies from a more closely related gray seal population on Sable Island have failed to find similar results of greater male size at birth [54,55]. As such, our results may support the notion of equal male/female fetal investment in gray seals. Even so, given the level of isotopic overlap in possible prey species, we cannot completely rule out the possibility that diet composition might vary to some degree depending on the sex of the fetus.
The positive correlation found between tissues suggests that vibrissae and lanugo can serve interchangeably as proxies of maternal dietary isotopic signatures during gestation. However, significant differences between tissues for both δ 15 N and δ 13 C values indicate that they may be capturing different time frames of maternal foraging. Lanugo δ 13 C values were 0.41‰ (±3.0) lower than vibrissae, and lanugo δ 15 N were 0.36‰ (±1.8) higher than vibrissae. These two keratinous tissues presumably share similar TEFs due to their chemical similarity [39,48,26,27]. Lanugo is only grown in-utero and therefore should maintain a maternal signature [56]. Both tissues develop slowly throughout gestation and have major overlapping periods of growth yet vibrissae begin development 1-2 months prior to lanugo [35,39]. The discrepancies in SI values may reflect these different time frames of growth. In addition, lanugo values were more centrally located within the mixing model polygons in bivariate isotopic space (Fig 2). This made the lanugo model more sensitive to prior information (i.e posterior distribution more closely reflected prior distribution) than the vibrissae model (Table 4). Aside from this example, lack of congruity between the tissues was small enough not to cause major differences between the SIMMs (Table 4). Thus, lanugo may be preferred over vibrissae when trying to reduce invasiveness and maximize sample mass due to its ease of sampling. The stable isotope values of potential prey species sourced from the Northwest Atlantic illustrated substantial isotopic overlap between species. Among demersal ground fish, we found no statistical isotopic difference between Atlantic cod and white hake and no statistical difference between red hake and silver hake suggesting these species share similar trophic niches. Northern shortfin squid and longfin squid, though unique from many other species were also statistically isotopically identical to one another for δ 13 C but not for δ 15 N. No species were isotopically unique in bivariate isotopic space, though sand lance, with only two significant relationships, had the fewest statistical overlaps. These results indicate a general lack of spacing between many of the potential gray seal prey species in the Northwest Atlantic, and in some cases large variation within individual prey species. While difficult to assess, it is possible that the small sample sizes analyzed for some prey groups impacted the ability to provide precise estimates of prey SI values and contributed to higher uncertainty of the diet analysis. However, our findings generally agree with the results of an analysis of stable isotope values of common species from Georges Bank, Northwest Atlantic, albeit with slight differences in δ 15 N and δ 13 C values for white hake and silver hake and a larger sample size for species analyzed in our study [57]. Stable isotope analysis of the foraging ecology of the US gray seal (Halichoerus grypus) population

Combining prey sources
Our research illustrates one of the major challenges in using SIMMs to infer the diet composition of gray seals: They are a generalist predator feeding on several isotopically similar prey sources. Similar to studies of other generalist consumers, the large number of possible prey species and the degree of overlap between these species it was necessary to combine sources prior to including them in SIMMs [16,36,45]. However, while combining prey sources can increase the predictive strength of SIMMs, the loss of taxonomic resolution can limit the model's interpretation and applicability [41].
Our study used cluster analyses to isolate the most statistically similar of the ten likely prey species into six groups and further reduced that into three possible prey clusters. Models that included six prey clusters provide little interpretative value due to the large overlap in 95% credibility intervals. However, combining sources into three prey clusters allowed us to draw more robust conclusions, though at lesser taxonomic resolution. Specifically, three source models provide better constrained estimates into the consumption of the 'demersal fishes' group relative to the other possible prey sources. This suggests to us that for the gray seal, with its plentiful dietary inputs, combining sources into ecologically informative groups may be the best way to derive useful information from SIMMs. For example, the approach used here provides insights into seal's foraging into non-commercial species such as sand lance and squids.

Incorporating informative priors
One unique feature of Bayesian SIMMs is the ability to incorporate information from previous dietary knowledge as informative priors which can add additional refinement to model outputs [32,51,58,59]. Previous dietary knowledge incorporated into SIMMs may come from, but is not limited to, sources such as stomach content, scat, or prey DNA [44,60,61]. However, while integrating complementary dietary techniques in SIMMs can provide better estimates of the actual diet of consumers [44], there are potential trade-offs. For example, Franco-Trecu et al [60] found that incorporating scat analysis data into SIMMs used to estimate the diet compositions of the South American fur seals (Arctocephalus australis) and sea lions (Otaria flavescens) improved the precision in the estimated diet composition at the risk of inducing biases inherent to the original scat analysis in the estimates.
Incorporating data from scat analysis into our gray seal SIMMs consistently decreased the importance of squid relative to sand lance in comparisons to uninformed models but surprisingly did not improve the relative precision of these estimates (Table 4). This may have been driven by the large amount of isotopic variation incorporated into the prey groupings [44]. However, it is important to note that our source of prior information (e.g. scat analysis) may be biased towards inshore prey, such as sand lance, relative to other prey species [9,12]. Previous studies have shown that gray seal scat analysis in New England contains a dietary record of only 24-48 hours of foraging, occurring within 80km of sampled haul-out sites, and is likely reflective of inshore feeding [9]. As such, similar to Franco-Trecu et al [60] we caution that to draw reliable results from SIMMs using informative priors, there must be a high confidence in the accuracy of the complementary source data used to avoid inducing additional biases into model results.

Selecting TEF values
Another challenge when using SIMMs to infer the diet composition of gray seals is selecting appropriate TEFs. In our study, TEF selection was determined by the 'best-fit' mixing model polygon following the procedure outlined by Smith et al [50]. As the prey species used in this analysis reflect the dominate prey species consumed by gray seals in New England (i.e. greater than 90-95% of dietary biomass) [9] it is unlikely there were major unrepresented prey species in our analysis. Regardless, the high degree of variability in some prey taxa likely led to the large 95% confidence interval observed for our mixing model polygons (Fig 1). Surprisingly, the averaged ecosystem TEFs of Post [47] were the best fit and outperformed TEFs determined experimentally through a combination of captive feeding trials on phocid seals and those produced via the SIDER tool [48,26,27,49] (Fig 1).
This suggest the generic TEFs of Post [47] may be the most accurately calibrated to the isotopic modification between gray seal pup vibrissae and lanugo tissues relative to their mother's diet. There are several possible reasons for this result. First, it is possible that captive experimental conditions under which phocid TEFs have been determined in previous studies [48,26,27] are not representative of the environmental conditions that occur in the wild. This can occur either due to less varied diets and subsequent isotopic routing or lower energetic demands of animals in captivity [44,62]. Also, averaged TEFs from Post [47] may have been a better fit because the exchange of isotopes between maternal gray seals and pup tissues inutero may not be identical to the same relationship between the mother and her own tissues [62]. Stricker et al [62] found that pup vibrissae grown in-utero were 0.8‰ higher for δ 15 N and 0.4‰ lower for δ 13 C than paired maternal vibrissae in the otariid Steller sea lion, Eumetopias jubatus. These results were similar to a paired mother/pup vibrissae study by Lowther and Goldsworthy [63] on Australian sea lions, Neophoca cinera, which found differences between in-utero and maternal vibrissae of 1.2‰ higher for δ 15 N and 0.2‰ lower for δ 13 C [63]. Accordingly, the Post [47] TEFs had higher δ 15 N and lower δ 13 C values relative to experimentally and SIDER derived TEFs (Table 2). Thus this result may be expected if in-utero isotopic fractionation of tissues in gray seals was similar to Australian sea lions (i.e higher δ 15 N and lower δ 13 C values relative to adult TEFs). Also noteworthy, the Post [47] TEF values were in close agreement with the averaged TEF values chosen by Hindell et al [36] in their use of inutero pup vibrissae on another phocid, the bearded seal, suggesting their may be phylogenetic similarity in in-utero TEFS. Similar to Post [47], the TEFs from Hindell et al [36] are derived from a standard environmental system value and not from experimental data on in-utero tissue development in phocid pups [64].

Gray seal diet predictions
While both six and three prey source models results often had substantial overlap in 95% credibility intervals between prey groups, some inferences in the diets of gestating adult female gray seals are still possible. For example, all SIMM predictions suggest that female gray seals have a larger proportion of cephalopod prey in their diets and a lower reliance on sand lance than had been previously reported for gray seal populations [9,12,65] (Fig 2; Table 4). In addition, these results also support the notion that US gray seals are generalist consumers, foraging on a mixture of commercial and non-commercial fish stocks. For example, SIMMs with six prey sources suggested that between 30-34% of diet had the potential to be derived from commercially valuable species such as Atlantic cod and winter flounder, though there was substantial uncertainty around these predictions. When considering the three-source models it appears unlikely that gestating adult female gray seals specialize on commercially important fish as all 'demersal fishes' account for only 55-65% the predicted diet (Table 4).
For pregnant females, these results indicate a diverse prey assemblage geared towards opportunistic foraging. They also point to a potential bias of previous scat analysis: that they may underrepresent soft-bodied prey such as cephalopods and cartilaginous fishes such as winter and thorny skates [9,10,12,66]. This discrepancy may be accounted for in two ways. First, identification from hard parts can be difficult for both cephalopods and skates. Harvey [20] found that in a captive feeding study of the harbor seal, Phoca vitulina richardii, there was only a 37% recovery rate of cephalopod beaks. However, a more recent study of the harbor seal determined a successful recovery rate of 89% for cephalopod beaks [67]. Size of the squid prey may play a role, with large beaks more likely to remain in the stomach or be regurgitated [67]. Cephalopod beaks rarely fully break down but even if they are found in the feces they can be difficult to identify to species [1,13,20]. For skates, denticles are often the only hard part to survive digestion. It can also be difficult to accurately determine species, number of individuals consumed, or prey size from denticles alone [9,68].
Past studies of gray seal diet based on scat analysis have estimated up to half of diet may constitute sand lance [9,12,69]. In contrast, the results of our analysis suggest the contribution of sand lance to the diet of pregnant female gray seals is closer to 10-15% and likely does not exceed one third. This difference may be accounted for both by a more evenly distributed generalist diet and a greater reliance on previously underestimated prey items such as cephalopods and skates. In addition, this discrepancy may occur because scat analysis likely is representative of in-shore foraging near the collection site [9]. Sand lance is primarily an in-shore species in New England, and it therefore is unsurprising that scat analysis shows the greatest reliance on this species [9,12]. Foraging time frames captured by the longitudinal vibrissae and lanugo samples include time spent both offshore and inshore during the summer and fall [35,39,70]. Longfin squid and Shortfin squid are typically found offshore during the spring and come inshore during the summer [71,72]. Thorny skate are typically found offshore in waters 50-100m deep [73]. If throughout the timeline of gestation, gray seals are feeding both inshore and offshore it may have led to the observed lower reliance on highly inshore species such as sand lance in our analysis.

Conclusions
Determining specific prey sources for this generalist predator remains difficult with stable isotope analysis. Our work has shown that the large number of isotopically similar species gray seals consume limits the ability of SIMMs in predicting robust diet estimates at high taxonomic resolution. Even so, combining prey sources into ecologically informative groups and integrating prior information into SIMMs it is possible to obtain more robust insights into the diet of gray seals though with a loss of taxonomic resolution. Furthermore, while adding a priori information from independent estimates of gray seal diets has the potential to further refine SIMMs, they may also induce additional biases into model results. Finally, we were are able to establish the use of vibrissae and lanugo interchangeably as viable proxies for maternal signature and caution that the TEFs for in-utero growth for this phocid species likely differ from those determined experimentally in adults. With the gray seal population in New England continuing to grow, it is clear that accurate estimates of their diet will be necessary to accurately predict their impact on both regional fisheries and ecology [74]. As such we recommend that future studies in this region should seek to combine isotopic approaches with other, independent measures of diets such as fatty acid, scat, stomach, and prey DNA analyses to better account for the differing challenges and inherent biases found across methods.