Large scale molecular barcoding of prey DNA reveals environmental predictors of intrapopulation feeding diversity in a marine predator

,


Introduction
Predator-prey relations are an integral force in many communities and are often key to understanding how communities function.Predators occupy a wide spectrum of foraging strategies, ranging from generalists to specialists.Individuals using a smaller subset of resources than the population as a whole are defined as individual specialists (Van Valen, 1965) whereas individuals consuming a wider range of resources than used on average by the population are defined as generalists (Hanski, Hansson, & Henttonen, 1991).In some cases, preying on one species may preclude consumption of a different species as there may be trade-offs in skills required to utilize different resources (Arthur et al. 2016;Wilson & Yoshimura 1994).Interindividual differences in resource use that are transient and the result of short-term choices in habitat or hunting strategies are best described as intrapopulation feeding diversity while permanent differences between individuals based on sex, size, or personality are better described as individual specialization (Van Valen 1965).Both types of inter-individual differences in resource use need to be examined to understand predator-prey interactions.
The level of individual specialization and/or intrapopulation feeding diversity can affect food web dynamics, responses to changes in prey availability, and the accuracy of predictive models (Bolnick et al., 2003).For example, in a population of bluegill sunfish Lepomis macrochirus prior experience foraging on a single prey type increases the likelihood of an individual using that resource, even when another resource becomes more profitable (Werner, Mittelbach, & Hall,1981).Further, theoretical models predict that changes in prey abundance in a system with highly specialized individuals (i.e.slow to switch preferred prey) are much more likely to have chaotic dynamics (Abrams & Matsuda, 2004).In general, dynamics can be sensitive to small variations in the speed of predators changing prey preferences (Abrams & Matsuda, 2004), which would be affected by the predator's level of specialization.More generally, diversification within a population can have significant impacts on ecosystems functions, such as prey community structure (Harmon et al., 2009).In addition, differences in a single species' population structure can have larger impacts on community composition than differences between species (Rudolf & Rasmussen, 2013).These findings imply that differences between individuals of the same species can be important drivers of ecosystem functions (Harmon et al., 2009, Rudolf & Rasmussen, 2013).Thus, including metrics of variation in foraging decisions between individuals of the same population in ecosystem studies provides a more clear and accurate description of the system and ignoring them can be an oversimplification of the ecological interactions in the community (Araújo, Bolnick, & Layman, 2011;Bolnick et al., 2003;Bolnick et al., 2011;Dall, Bell, Bolnick, & Ratnieks, 2012).Unfortunately, most foraging studies do not describe the level of intraspecific variation in the predator population.
Variation in foraging decisions within populations of predators is difficult to describe empirically because they require observing a large number of predation events in multiple individuals across many ecological contexts.The empirical problems are even greater when studying predators that forage in environments where it is difficult to directly observe predation events (e.g., marine environments) and that prey on a large diversity of taxonomically similar prey species that make it difficult to determine which species has been consumed.Here, we demonstrate how the application of an individual diet specialization metric to a large set of molecular prey barcoding data from scat can be used to describe intrapopulation feeding diversity by examining the short-term variation in individual foraging decisions in a marine predator.Our analysis allowed us to explore 1) correlations of individual diet diversity with the sex of the predator and time of year in which the predation occurred, and 2) to test whether short-term diet diversity was related to the consumption of particular prey species and their preferred habitat.
Due to the large diversity of prey species that harbor seal populations eat, the species has historically been considered a generalist predator (Teilmann & Galatius, 2018).However, prey composition and foraging dive behavior of harbor seals in the Salish Sea vary relative to habitat, sex, and time of year (Lance et al., 2012;Olesiuk et al., 1990;Wilson, Lance, Jeffries, & Acevedo-Gutiérrez, 2014).Harbor seals also eat different types of prey depending on the type of environment in which they forage.Scat samples from haul-outs located in estuaries have higher prey diversity than those coming from outside estuaries (Lance et al., 2012;Luxa & Acevedo-Gutiérrez, 2013).Further, males and females consume different prey (Bjorland et al., 2015;Schwarz et al., 2018) and have different foraging dive patterns (Wilson et al., 2014).Specifically, females frequently perform longer and deeper foraging dives than males, and more commonly consume benthic species (Schwarz et al., 2018;Wilson et al., 2014).These traits of high abundance and differences in diet and foraging patterns between males and females, suggest that harbor seals could display intrapopulation feeding diversity..As such, ecosystem dynamics with regard to the effect of harbor seals on prey species are likely more complex than described in current models of the system which assume consistent generalized behavior (e.g.Chasco et al. 2017, Howard et al. 2013).As harbor seals are the most abundant mammalian predator in the Salish Sea, and prey on many species of economic and conservation concern, an accurate understanding of their role in ecosystem dynamics is important, for which high quality diet data are required.While current bioenergetics based models are useful descriptors of harbor seal consumption (e.g.Chasco et al. 2017, Howard et al. 2013), they can be improved by including the effects of different foraging strategies across sexes and individuals.
Obtaining high quality diet data from large, mobile organisms, such as marine mammals, can be costly and time consuming (Rothstein, McLaughlin, Acevedo-Gutiérrez, & Schwarz, 2017).Analysis of prey contents in scat via metabarcoding is a relatively cheap, non-invasive, and time efficient way to obtain large sample sizes with species-level taxonomic resolution (Deagle et al., 2005;Deagle et al., 2019;Rothstein et al., 2017;Tollit et al., 2009).However, to our knowledge, these molecular techniques have not previously been used to quantify intrapopulation feeding diversity at large spatial and temporal scales.
Here, we use molecular barcoding of prey DNA from scat in a novel way to examine intrapopulation feeding diversity, in harbor seals by answering the following questions: 1.How do the factors Sex, Time of year, Location, and Year affect cross sectional estimation of a specialization metric? 2. What prey items correlate with high levels of specialization relative to sex and environment?To answer these questions, we collected and analyzed scat from wild harbor seals in the Salish Sea.Diet of harbor seals was determined from the scat using both molecular and traditional techniques.Sex of the depositor was determined using molecular techniques.Diet data were analyzed with a proportional similarity index (Bolnick, Yang, Fordyce, Davis, & Svanbäck, 2002) to describe the variation in individual foraging decisions in harbor seals.

Collection and processing of scat samples
Scat collections were conducted by multiple researcher groups at five harbor seal haul-outs in the Salish Sea over a period of four non-sequential years (Figure 1).Haul-outs varied in seal population size as well as by habitat type (Table S1).Not all sites were visited every year and the months during which each site was visited varied between years (Table S1).Collections at Belle Chain, Cowichan, Comox, and Fraser River were conducted by teams from University of British Columbia under Fisheries and Oceans Canada's Marine Mammal Research Licenses (MML 2011-10 andMML 2014-07) and University of British Columbia's Animal Care Permits (A11-0072 and A14-0068) awarded to University of British Columbia Marine Mammal Research Unit.Collections at Baby Island were conducted by a team from Western Washington University under Federal Permit 18002 from the United States Office of Protected Resources, National Marine Fisheries Service, and a Western Washington University's Animal Care and Use Committee exception awarded to Alejandro Acevedo-Gutiérrez.
Collection of scat followed the general procedure described in Thomas, Deagle, Eveson, Harsch, and Trites (2016).Briefly, upon arrival at a haul-out we searched the entire area for scat.When scat was found, the entire scat was collected into a 126 µl nylon strainer inside of a 500 ml sealable container using a wooden tongue dispenser and plastic spoon.The container was then stored in a cooler with ice until transfer to a -20 freezer later that day.At Baby Island and Cowichan Bay in 2014 the entire outside of the scat was swabbed before collection.Swabbing focused on any mucus material, as it likely contains higher proportions of seal DNA (Rothstein, 2015).The swab was then placed in a vial of ethanol and stored in a cooler with ice until transfer to a -20 freezer later that day.
A DNA slurry of homogenized scat in ethanol was prepared for each sample to obtain a representative set of DNA following the procedure described in Thomas et al. (2016).Briefly, the entire scat was thawed in ethanol and homogenized within the mesh bag.After homogenization, a representative sample of DNA slurry was allowed to pass through the bag.The mesh bag was then removed, zip-tied, and stored at -20 for later use in prey hard part analysis.We then let the DNA slurry settle in the containers on the bench top overnight.The next day we pipetted the settled slurry into 20 mL scintillation vials that were subsequently stored at -20 until further analysis.

Sex determination of harbor seals via scat
To obtain DNA for sex determination, DNA was extracted from the scat matrix-ethanol slurry for all locations, except Cowichan 2014 and Baby Island.For these last two sites, DNA was extracted from the swabs.To extract DNA from swabs, the excess ethanol from the vial was poured off and the swab was dried in a vacuum centrifuge at 39°C until all ethanol had evaporated, approximately one hour.We then used QIAGEN DNeasy Blood and Tissue Kit to extract DNA from the dried swabs.DNA was extracted from slurry matrixes using QIAGEN QIAamp DNA Stool Mini Kit.Extracted DNA, from either the ethanol slurry or swab, was used in Taqman quantitative polymerase chain reactions (qPCR) to determine the presence and absence of X and Y chromosomes.The procedure was modified from Matejusová et al. (2013) and is described in depth in Rothstein (2015).The two probes that we used targeted the paralogous zinc finger X (ZFX) and zinc finger Y (ZFY) genes.Both probes are described in Matejusavá et al. (2013).Two reactions were run for each sample with each probe (four reactions total per sample).Each reaction consisted of: 4.5µl of ABI Taqman gene expression master mix, 0.5µl of either the ZFX or ZFY probe, and 5µl of DNA template.Reactions were run on a quantitative thermocycler with the following protocol: one holding cycle (50°C for 2 min, 95°C for 10 min) followed by 60 cycles of denaturation and annealing/extension (95°C for 15 sec, 60°C for 1 min).Four positive (two reactions for each sex, one ZFX and one ZFY probe each) and four negative controls (two reactions for each ZFX and ZFY probe) were run with each set of reactions.Positive controls came from captive harbor seals of known sex at the Vancouver Aquarium in Vancouver, BC, and Point Defiance Zoo & Aquarium in Tacoma, WA.Negative controls consisted of PCR grade water in place of a DNA template.
If no amplification occurred in either ZFX reactions, the sample was excluded from further analysis.If no amplification occurred in either ZFY reaction, but amplification did occur in either or both reactions with the ZFX probe, the sample was tallied as deposited by a female.If amplification was observed in either or both ZFY reactions, as well as in either or both ZFX reactions, the sample was tallied as deposited by a male.The false negative rate for two failed ZFY reactions (and thereby incorrectly classifying a male as a female) was 1.35%.This value was calculated from the occurrence of only one of the two ZFY reactions having positive amplification within a sample that was classified as male.Although this false negative rate is low, we excluded any samples that amplified in ZFY reactions but failed to amplify in ZFX reactions.

Prey determination in harbor seal scat
The diet of harbor seals was determined by combining DNA and hard part data.The DNA prey identification and quantification were completed following the procedure outlined in Thomas et al. (2016).Briefly, for all locations the scat matrix DNA (obtained as described above) for each sample underwent a multiplex PCR using primers for a 16s mtDNA barcoding fragment (˜260 bp) described by Deagle, Chiaradia, McInnes, and Jarman, (2010).Amplicons were labelled using a combination of unique F and R primer tags, in addition to indexed, post-PCR ligated Illumina TruSeq adapter sequences (for details see Thomas et al. 2016).An Illumina MiSeq was then used to sequence the amplified DNA fragments.Lastly, a custom BLAST database comprised of publicly available reference sequences specific for known prey species was used to produce identifications to the lowest taxonomic level possible for each amplified sequence.
Extraction and preparation of prey hard parts were completed by Thomas et al. (2017) for Belle Chain, Comox, Cowichan Bay, and Fraser River 2012 and 2013 samples, by one of the authors (BAN) for Cowichan Bay 2014, and by the first author (MRV) for Baby Island samples.Each scat was placed in a set of nested sieves, and then rinsed and stirred until all that was left in the sieves were prey hard parts.All hard parts, except cephalopod beaks, were transferred to 20 ml scintillation vials with 70% ethanol.They were allowed to sit for a minimum of two weeks before the liquid was poured off and the hard parts were allowed to dry.The cephalopod beaks were transferred to separate 20 ml scintillation vials with ethanol.All diagnostic prey hard parts were identified to the lowest taxonomic level possible using reference sets of prey bones from Washington and British Columbia by Thomas et al. (2017) for Belle Chain, Comox, Cowichan Bay, and Fraser River, and by collaborators at Long Live The Kings for Baby Island samples.Published keys for both fish bones and cephalopod beaks were used as described in Thomas et al. (2017).Notably, this analysis allowed differentiation between the proportion of adult and juvenile Oncorhynchus spp.consumed.The percentage of juvenile versus adult salmon was determined using the method described in Thomas et al. (2017).

Quantification of specialization metrics
Cross-sectional sampling, which only requires data from a single time point, can be used to estimate intrapopulation feeding diversity at large spatial and temporal scales.For pinnipeds, a single time point can be examined via scat collection and analysis with a single scat being indicative of the last few foraging bouts (Bowen & Iverson, 2012).,For this case, due to various limitations, one cannot calculate absolute individual specialization using cross-sectional sampling (Araujo et al., 2011;Novak and Tinker 2015).However, one can compare relative specialization between samples collected in the same manner (i.e., within the dataset).Thus, by calculating specialization metrics based on prey proportions in harbor seal scat we were able to deepen our understanding of intrapopulation feeding diversity and uncover patterns at the level of sex, season, and location.
To this end, we quantified the level of specialization represented by each sample using the proportional similarity index (PS i ) function in the R package RInSp (Zaccarelli, Bolnick, & Mancinelli, 2013).PS i calculates the overlap between what an individual is eating and what the population is eating using the following formula: Where p ij represents the proportion of resource j used by the individual i and q j represents the proportion of resource j used by the population.PS i is bounded by a theoretical minimum, which is population dependent as described below, and one.The variable population dependent minimum indicates a complete specialist and a PS i of one indicates a generalist (Bolnick et al., 2002).Because PS i is bounded, we report the overall average value with 95% confidence intervals calculated using Monte-Carlo resampling in the R 3.3.1 package "resample".Traditionally, prey counts have been used for calculatingPS i , not proportions, as each count is assumed to represent an independent prey capture decision (Araujo et al., 2011;Bolnick et al., 2002).Proportions of prey metabarcoding reads are representations of the prey biomass proportions that were consumed by the predator, and similar proportions can result from consuming a few large or many small prey individuals.Correspondingly, calculatingPS i , using the proportions of prey metabarcoding reads will produce a metric of "biomass specialization" that does not necessarily reflect independent prey capture decisions.Nevertheless, it describes intrapopulation variations in the utilization of different prey species.In addition, we calculated PS i relative to groups of samples from a certain point in space and time.If individuals in that particular group are encountering the same size distribution of prey, then diet proportions may represent the same relative relationship of prey capture decisions as counts of individual prey items.Despite the potential limitations, there are several benefits to using this type of data.Coupled with scat collection, it allows for large samples sizes, is non-invasive, and gives high taxonomic resolution.
To define our groups for analysis, samples were separated by Location, Sex, Year, and Month of collection, yielding a total of 89 groups (Table S1).PS i values for each sample were then calculated for each one of these groups.Within each group, each sample was treated as coming from a different individual due to the low probability of resampling the same individual (Rothstein et al., 2017).
Because different groups for analysis can have different theoretical minima, there is potential bias when comparing specialization metrics across groups.Differences in theoretical minima occur due to differences in sample size (the number of scat in each group) and/or differences in minimum prey densities (the smallest occurring proportion of a prey species in a group's diet).Due to very low minimum prey densities in our data set, the theoretical minima are determined by sample size (Table S1).We addressed this potential bias in multiple ways.First, we excluded from analysis the smallest groups (those with < 5 samples) as they have the highest theoretical minimum and thus the most potential for bias.We also used Spearman's rank correlation to estimate how much variance was explained by differences in sample size.This correlation was accomplished by comparing sample size to the average PS i for each group we kept.We also calculated the theoretical minima for each group by dividing one by the number of samples in the group and then examined the range, average, and median of those minima.Additionally, sample sizes of each group were included in modeling of the data, which is described below.Lastly, the seasonal changes in sample size were visually compared with the seasonal patterns in PS i values.

Comparison of factors influencing relative specialization
We analyzed the relative influence of the factors Sex, Month, Location, Year, and Sample Size on PS i using generalized linear mixed models (GLMMs).We chose mixed models because they allowed us to include Sample Size, Location, and Year as random variables.Restricted maximum likelihood estimation was used because it considers the loss of degrees of freedom when estimating fixed effects and thus offers a more unbiased estimate than maximum likelihood methods (West, Welch, & Galecki, 2015).Before modeling the data, we performed a logit transformation (log( PSi 1−PSi )) on thePS i values to normalize them.This transformation was necessary because PS i is bounded by a theoretical minimum and one (Bolnick et al., 2002).When numbers are bounded, the variance distribution is shifted towards the mean (Sokal & Rohlf, 2012).A logit transformation is an excellent choice for addressing this because it extends the tails of the distribution more than other alternatives (Warton & Hui, 2011).
All models were tested in the R 3.3.1 package lme4 (Bates, Machler, Bolker, & Walker, 2015).This package provides basic measurements of goodness-of-fit including AIC and coefficients.The R 3.3.1 package MuMIn was used to determine the r2 values for mixed models.Subsequent calculations of AIC, and w i (positive Akaike weights or likelihood of being the best model (Anderson, 2008)) were completed using Excel.AIC was calculated as the difference between two AIC scores; w i was calculated following Burnham and Anderson (2010).
To more clearly understand the relationship between sex ratio of the population and PS i , sex ratios were produced for every paired group (groups of males and females from the same location, month, and year) by calculating the percent of scat identified as female.The average PS i for each paired group was then compared with this female percentage using a Spearman's rank correlation.Spearman's rank correlation was used to account for the heteroscedasticity of the dataset, and was completed using R 3.3.1.Additionally, the average proportion of female scat for each month and location are visualized in the supplemental material (Figure S1).

Correlations between prey items and relative specialization
For each scat, prey item proportions (number of sequences per prey species / total number of sequences generated) were lumped into orders and summed.We then performed correlations between the proportion of the diet that each order comprised in each sample and the PS i for that sample.To examine if similar patterns occurred in both sexes, the analysis above was completed for male and females separately.We used the Bonferroni method to correct alpha for multiple comparisons (alpha values are specified in table captions).Additionally, to determine if benthic species were associated with a more specialist diet, a correlation was run between the proportion benthic prey and specialization value for each scat.Due to the heteroscedasticity of the dataset, we used Spearman's rank method for all correlations (Sokal & Rohlf, 2012).All correlation analysis was conducted in R 3.3.1.Because smaller PS i values indicate higher levels of specialization, a negative correlation value suggests a positive relationship with specialization.

Quantification of specialization metrics
Over the course of four non-sequential years, at five different locations, we quantified the diet of 1,520 scat samples.We successfully determined the sex of the depositor for 1,145 of those scats (75% success rate) (Voelker et al. 2018).The number of scat with successful sex determination varied by location and month (Table 1).Samples with successful sex determination were then binned by the factors Sex, Location, Month, and Year to form unique groups for analysis (Table S1).After eliminating samples without sex determination and with small sample sizes (< 5 samples), we were left with 1,083 samples in 89 groups.Only these 1,083 samples were used in further analyses.PS i was calculated relative to the samples in a specific group, and the average PS i across all samples and groups was 0.399 (95% CI = 0.026, R = 100,000).The PS i values of the 1,083 samples were not normally distributed (kurtosis = 2.66, skewness = 0.65, Figure S2).Therefore, a logit transformation was used to adjust the variance distribution (kurtosis = 5.21, skewness = 1.01, Figure S2).These transformed PS i values were used to run the GLMMs.Additionally, the range of theoretical minima across the 89 groups was 0.027 -0.2 (average = 0.103, median = .091);there was also a correlation between average PS i and theoretical minimumPS i (rho = -0.231,p = 0.03).This potential bias is addressed in the discussion.

Comparison of factors influencing relative specialization
Based on AIC values, r 2 results, and model likelihood, the best fit GLMM was Month*Sex + (1|Sample Size) + (1|Location) + (1|Year) (Table 2).The r 2 value and residual plots indicate that this model fit the data well (Table 2, Figure S3).The random factors of Sample Size, Location, and Year explained 0.39, 0.36, and 0.002 of the variance (SD = 0.62, 0.597, 0.05), respectively.The r 2 value calculated with fixed and random effects was over four times that of the r 2 value calculated using just fixed effects.Removing Month from the model caused a larger decrease in goodness-of-fit measurements than removing Sex (Table 2).Removal of the interaction term also caused a decrease in goodness-of-fit measurements (Table 2).Further, the interaction terms for Sex and the Months of August and October were significant (t = 2.86, 2.68, p = 0.004, 0.007 respectively).However, correlation analysis between the percent female scat collected for each paired group (which acted as a proxy for the effect of sex ratio in the population) and the average PS i for that pairing revealed no significant trend (rho = -0.071,p = 0.655).
To further examine the interaction between Sex and Month, the factor of Month was grouped into three levels: spring (April and May), summer (June, July, August), and fall (September, October, and November).The specialization metric showed a distinct shift throughout the year in males but not females (Figure 2).In summer and fall, males had relatively lower levels of specialization than females (Figure 2).To address the potential bias introduced by sample size for this mode of data analysis, we plotted the sample size for each group by season.The pattern observed in PS i values was not reflected in sample size (Figure 3).
Visual inspection of the data by month suggested males had a decrease in relative specialization in July through October (Figure 4).Based on 95% confidence intervals of logit transformed PS i values,PS i during these months only overlapped with April (Figure 4).The same pattern was not apparent in females because the 95% confidence interval for logit transformed PS i of female groups overlapped for all months (Figure 4).This trend varied in intensity by location (Figure 5).The described pattern was reflected most strongly in Belle Chain, Comox, and Fraser River (Figure 5).However, because scat were not collected at Baby Island after July, no comparison could be made at that location (Figure 5, Table S1).

Correlations between prey items and relative specialization
Correlation analysis between diet proportions of prey orders andPS i revealed that eight orders out of twelve showed significant correlations (alpha = 0.0038, Bonferroni correction, Table 3).Adult Salmoniformes, a largely pelagic group, were correlated with a generalist diet (rho = 0.27, p < 0.001) while juvenile Salmoniformes showed no significant correlation.Clupeiformes, another largely pelagic group, correlated with a generalist diet as well (rho = 0.24, p < 0.001).Conversely, Pleuronectiformes, a demersal and benthic group, correlated with a specialist diet (rho = -0.38,p < 0.001).Further, Gadiformes, which has both pelagic and demersal representatives, showed no correlation (rho = -0.04,p = 0.38).
Correlations performed with just data from female scat showed similar patterns.All orders of prey showed the same relationship withPS i , or were no longer significant, such as Salmoniformes (Table 4).The only new order to show significance was Batrachoidiformes.Correlations performed with only male scat once again showed similar patterns.New orders to show significance were Chimaeriformes (with a specialist diet) and adult Salmoniformes (with a generalist diet) (Table 5).
The correlation between proportion of benthic species in each scat andPS i suggested a positive relationship between the amount of benthic species consumed and the level of relative specialization (rho = -0.384,p > 0.001).A similar relationship was observed when the female and male data sets were examined separately (rho = -0.407,p > 0.001; rho = -0.35,p > 0.001 respectively).

Discussion
We successfully assigned PS i values to 1,083 scat collected from five different locations over the course of four non-sequential years (Table S1, Figure 1) (Voelker et al. 2018).As measured by repeated cross-sectional sampling and a specialization metric (PS i ), the overall level of intrapopulation feeding diversity in the region was high (PS i = 0.399, 95% CI = 0.026, R = 100,000).Further, Month, Sex, and Location were all important factors influencing this feeding diversity.Interestingly, Month and Sex had a significant interaction.Habitat of an individual's primary prey also seemed to have an impact on relative specialization suggesting that seasonal and sex-specific patterns in the use of benthic vs pelagic were the underlying cause for the observed intrapopulation feeding diversity.These indications of intrapopulation feeding diversity suggest the feeding ecology of harbor seals in the Salish Sea is complex and that prey species of concern may be impacted differently by each sex.

Estimated level of intrapopulation feeding diversity
We used cross-sectional sampling to calculate PS i , a specialization metric, as an estimate of intrapopulation feeding diversity in Salish Sea harbor seals.Our data confirmed intrapopulation feeding diversity across the spatial (hundreds of km) and temporal (years) scales that the scat samples represented (averagePS i = 0.399, 95% CI = 0.026, R = 100,000).These data demonstrate intrapopulation feed diversity but leave room for two alternative hypotheses (which cannot be separated in this case), regarding the absolute level of individual specialization: the occurrence of congruent long-term generalists and short-term specialists, or the occurrence of long-term specialists.
4.2 Importance of Time of year, Sex, Location, and Year on relative specialization ( PS i ) Month was an important predictor of relative specialization because removing it from the model caused a large drop in goodness-of-fit measurements (Table 2).This pattern makes intuitive sense as the type of prey eaten by harbor seals (Lance et al., 2012;Olesiuk et al., 1990) as well as their dive foraging behavior (Wilson et al., 2014) vary throughout the year.Therefore, changes in foraging behavior (both prey choice and dive type) were likely mechanisms behind the observed change in relative specialization throughout the year.However, there were likely other factors influencing relative specialization in addition to month.
Sex also had an impact on relative specialization, yet smaller than that of Month (Table 1).Differences in the level of relative specialization between female and male harbor seals were likely due to females and males in the region eating different prey items and having different foraging strategies (Bjorland et al., 2015;Schwarz et al., 2018;Wilson et al., 2014).For instance, females more often perform deeper foraging dives, eat benthic prey more commonly, and have smaller home ranges than males ( (Peterson, Lance, Jeffries, & Acevedo-Gutierrez, 2012;Schwarz et al., 2018;Wilson et al., 2014).Therefore, we propose the following theoretical resource distribution: males have more overlap between individuals while the females have less overlap between individuals; variations in this pattern appear to be associated with prey type (which will be addressed in the following section) (Figure 6).
Including an interaction term between Month and Sex increased the goodness-of-fit of the model (Table 2).This result indicates that differences between male and female seals likely varied throughout the year.Specifically, there were clear decreases in relative specialization in male harbor seals during the summer and fall months that were not reflected in females (Figure 2), indicating that the behavior of both sexes was similar in the spring but diverged in the summer and fall.This behavior was likely due to changes in feeding patterns of females and males throughout the year (Lance et al., 2012;Wilson et al., 2014).A possible reason for the different feeding patterns in the summer months is female change in behavior due to pupping (Ternte, Bigg, & Wigg, 1991).While nursing, females spend most of their time on the haul-out and make short foraging trips (Boness, Bowen, & Oftedal, 1994;D'Agnese, 2015).A similar difference was seen between sexes during the fall; however, both sexes were relatively least specialized during the fall.During the fall, there is a large influx of returning adult Salmoniformes (Quinn, 2005) that are preyed upon by both female and male harbor seals (Schwarz et al., 2018).In the Salish Sea, Salmoniformes can compose >50% of the population diet in the summer and fall (Lance et al., 2012).This resource could be rich enough that it is beneficial for a majority of seals, both males and females, resulting in less need for specialization.This explanation is further supported by the correlation between feeding on adult Salmonifomes and a relatively less specialized diet (Table 3), indicating it was a widely used resource in the region.
Our data also suggest that location explained a large amount of variance in relative specialization.The random factors of Year and Location increased the r 2 by more than four times, indicating that both had a large influence on relative specialization.However, because Sample Size, Location, and Year explained 0.39, 0.36, and 0.002 of the variance (SD = 0.62, 0.597, 0.05), respectively, one can assume that Sample Size and Location were the random factors responsible for the increase in goodness of fit of the model, not Year.This result indicates that where the seals were foraging impacted the level of relative specialization in the population, without noticeable changes from year to year.Our results also indicate that there was likely some bias introduced by the number of samples in a group.For instance, there was a correlation between average PS i and theoretical minimum PS i (rho = -0.231,p = 0.03).However, this potential bias is unlikely to have had a substantial effect on the outcome of our study because we included sample size as a random variable in the model and variation in sample size does not appear to explain the seasonal changes in PS i (Figures 2, 3).

Correlation between relative specialization and prey species composition
Our data suggest that the higher proportion of benthic species consumed, the relatively more specialized the diet of the predator.This pattern was observed in the full dataset, as well as when female and male data were considered separately (Tables 2-5).This information ties to our knowledge of the foraging patterns of male versus female harbor seals in the region.Females more often perform deeper foraging dives (Wilson et al., 2014) and eat more benthic species than males, who eat more pelagic species (Schwarz et al., 2018).In Scotland, harbor seal scat samples represented either a largely pelagic foraging strategy or largely benthic foraging strategy (Tollit, Greenstreet, & Thompson, 1997), and males had larger range and duration in foraging trips (Thompson, Mackay, Tollit, Enderby, & Hammond, 1998), suggesting that the separation between the two foraging strategies is not just a regional phenomenon.
The consistency across sexes indicates this pattern is reflective of foraging strategies specific to the ecology of prey species, and not just indicative of different diet preferences between males and females.We hypothesize that this pattern was caused by higher variability in benthic environments (Lalli & Parsons, 1997).If prey have more variable life strategies, a single foraging strategy will not suffice to catch them all.Because an organism is likely limited in the number of foraging strategies at which it can be effective, an individual could be limited in the number of prey species it can exploit.
There is the possibility that the sex of the individual determines the level of specialization regardless of the prey consumed.However, the consistency of benthic prey being associated with a relatively specialist diet, and pelagic prey being associated with a relatively generalist diet, in the complete, only female, and only male datasets suggests that the ecology of the prey species and the sex of the seal was driving the observed pattern.This idea is supported in other literature as well.Individual male harbor seals in Nova Scotia use different behaviors when foraging for benthic versus pelagic prey (Bowen, Tully, Boness, Bulheier, & Marshall, 2002) and larger seals are more likely to forage in pelagic environments regardless of sex (Bjorkland et al., 2015).
If prey species ecology is driving specialization levels, it is especially interesting to consider harbor seal consumption of juvenile Salmoniformes.For example, juvenile sockeye (Oncorhynchus nerka , Walbaum 1792) correlated with a generalist diet (rho = 0.22, p = 0.004).This could indicate that seals were not seeking out juvenile Salmoniformes specifically but rather eating them as a byproduct of focusing on fish that match the image of forage fish (e.g.small and silver) while conducting pelagic foraging strategies.This is just one example of how understanding the level of specialization could deepen our scope of knowledge regarding harbor seal impacts on prey species of concern.

Study limitations
There are a few notable limitations to this study.First, there was the potential for variation in sample size to introduce bias.However, there were no discernable patterns between sample size and average relative specialization by season (Figures 2, 3).We also included sample size as a random factor in the model to account for any bias introduced there.Hence, any bias introduced by sample size was likely minimal.Second, because scat were collected from the same haul-out multiple times there is a chance that some scat collected came from the same individual.However, this chance is low as Rothstein et al. (2017) estimated the sampling scheme to track five individuals at Cowichan Bay (i.e. a single haul-out) as 440 samples over 22 sampling bouts.Compared to the 1,083 samples used in this analysis from five different haul-outs, it seems unlikely there was a high rate of resampling the same individuals.Third, there are biases in the metabarcoding PCR process for determining diet (Thomas, Jarman, Haman, Trites, & Deagle, 2014).The prey proportions recorded for each sample are not directly proportional to the amount of prey that was ingested (Bowen & Iverson, 2012;Thomas et al., 2014).However, this approach is accepted to be semi-quantitative, biases are assumed to be consistent between samples (Thomas et al., 2014), and the approach has been successfully used in other studies (Deagle et al., 2010;Pompanon et al., 2012;Schwarz et al., 2018;Thomas et al., 2014;Thomas et al., 2017).Furthermore, this approach is superior to the alternative occurrence-based methods for generating diet proportions (Deagle et al., 2019).On a related note, these molecular methods do not provide data that directly equate to counts of prey consumed.But, if individuals within a local group encounter the same size distribution of a given prey species, then diet proportions represent the same relative relationship of prey capture decisions.Further investigation into potential biases introduced by using proportion type data would be useful as this methodology has many benefits and is a valuable molecular technology that should be applied in the future.

Conclusion
Through quantifying levels of relative specialization, we show that intrapopulation feeding diversity occurs in Salish Sea harbor seals between locations, seasons, and sexes.In both female and male harbor seals benthic prey were more commonly associated with a more specialized diet, suggesting the prey's ecology had a role in driving the level of specialization in addition to sex.These different impacts of male versus female on benthic versus pelagic prey should be considered henceforth when management address harbor seal interactions with species of concern.Further, we demonstrated how the use of molecular prey barcoding from scat allows for higher taxonomic resolution and greaterspatiotemporal resolution than conventional methods.The resulting large scale examinations of intrapopulation feeding diversity uncovered previously unknown complex interactions between predators and prey.
Table 2. GLMM models of prey specialization in Salish Sea harbor seals.The r 2 fixed column indicates how much variance was explained solely by the fixed effects.The r 2 column indicates how much variance was explained by both fixed and random effects.The AIC column indicates the fit of each model, lower values indicate a better model.The w i column indicates the relative likelihood of each model being the best model of those tested.Analysis represents 1,083 samples from groups with >5 samples.
Table 3. Correlations of prey proportions by taxonomic order toPS i values of Salish Sea harbor seals.Only significant correlations (p < 0.0038) are shown.A negative rho value indicates a positive correlation with specialization.Analysis represents 1,083 samples from groups with >5 samples.A total of twelve prey orders were included in the analysis.Data are organized by correlation value.
Table 4. Correlations of prey proportions by taxonomic order toPS i values of Salish Sea female harbor seals.Only significant correlations (p < 0.0038) are shown.A negative rho value indicates a positive correlation with specialization.Analysis represents 498 samples from groups with >5.A total of twelve prey orders were included in the analysis.Data are organized by correlation value.
Table 5. Correlations of prey proportions by taxonomic order toPS i values of Salish Sea male harbor seals.Only significant correlations (p < 0.0038) are shown.A negative rho value indicates a positive correlation with specialization.Analysis represents 647 samples from groups with >5 samples.A total of twelve prey orders were included in the analysis.Data are organized by correlation value.

List of Figures
Figure 1.Haul-out sites where harbor seal scat was collected in the Salish Sea.Collection locations are indicated by black dots and labeled with the name used throughout this paper.Coordinates for each site are as follows; Comox: 49°35'45.53"N,124°52'4.39"W,Fraser River: 49°4'27.17"N,123°8'49.46"W,Belle Chain: 48°58'10.73"N,123°29'34.63"W,Cowichan Bay: 48°44'14.28"N,123°37'17.76"W,Baby Island: 48°5'58.31"N,122°31'41.29"W.Tables Table 1.Number of harbor seal scat from the Salish Sea with successful sex determination from all locations, months, and years.Within each monthly column the numbers are as follows: female scat, male scat.If multiple collection bouts occurred at a single haul-out within one month the total number of scat for that month is listed.An "na" indicates no scat were collected at that site during that month.Abundance estimates for Belle Chaine, Cowichan Bay, Comox, and Fraser River where calculated from Olesiuk et al. (2009).The abundance estimate for Baby Island was taken from Jeffries et al. (2003)

Figure 2 .
Figure 2. Logit transformed average PS i values with 95% confidence intervals of all harbor seal scat from groups with >5 samples (n = 1,083 scat samples).AveragePS i was calculated for each group.Groups were then split be sex and lumped by season.Spring = April, May; Summer = June, July, August; Fall = September, October, November.A lower value indicates more specialization.

Figure 3 .
Figure3.Average sample size with 95% confidence intervals of harbor seal groups for analysis with >5 samples (n = 1,083 scat samples).Groups were then split be sex and lumped by season.Spring = April, May; Summer = June, July, August; Fall = September, October, November.

Figure 4 .
Figure 4. Logit transformed average PS i values with 95% confidence intervals of all harbor seal scat from groups with >5 samples (n = 1,083 scat samples).AveragePS i was calculated for each group.Groups were then split be sex and lumped by Month.The left graph shows females, the right graph shows males.A lower value indicates more specialization.

Figure 5 .
Figure 5. Logit transformed average PS i values and 95% confidence intervals for all harbor seal scat from groups with >5 samples (n = 1,083 scat samples).AveragePS i was calculated for each group.Groups were then split by sex and location and then lumped by month.A lower value indicates more specialization.BC = Belle Chain, BI = Baby Island, CB = Cowichan Bay, CM = Comox, FR = Fraser River.

Figure 6 .
Figure 6.Model of resource distribution of harbor seals in the Salish Sea.This theoretical schematic demonstrates the smaller within-individual resource use of female harbor seals with regards to benthic prey in the Salish Sea.A) Prey use distribution in the spring (more specialist behavior by females).B) Prey use distribution in the summer and fall (less specialist behavior by females).

Figures
Figures

Figure 2 .
Figure 2. Logit transformed average PS i values with 95% confidence intervals of all harbor seal scat from groups with >5 samples (n = 1,083 scat samples).AveragePS i was calculated for each group.Groups were then split be sex and lumped by season.Spring = April, May; Summer = June, July, August; Fall = September, October, November.A lower value indicates more specialization.

Figure 3 .
Figure 3. Average sample size with 95% confidence intervals of harbor seal groups for analysis with >5 samples (n = 1,083 scat samples).Groups were then split be sex and lumped by season.Spring = April, May; Summer = June, July, August; Fall = September, October, November.

Figure 4 .
Figure 4. Logit transformed average PS i values with 95% confidence intervals of all harbor seal scat from groups with >5 samples (n = 1,083 scat samples).AveragePS i was calculated for each group.Groups were then split be sex and lumped by Month.The left graph shows females, the right graph shows males.A lower value indicates more specialization.

Figure 5 .Figure 6 .
Figure 5. Logit transformed average PS i values and 95% confidence intervals for all harbor seal scat from groups with >5 samples (n = 1,083 scat samples).Average P S i was calculated for each group.Groups were then split by sex and location and then lumped by month.A lower value indicates more specialization.BC = Belle Chain, BI = Baby Island, CB = Cowichan Bay, CM = Comox, FR = Fraser River.A)

Table 2 .
. GLMM models of prey specialization in Salish Sea harbor seals.The r 2 fixed column indicates how much variance was explained by only the fixed effects.The r 2 column indicates how much variance was explained by both fixed and random effects.The AIC column indicates the fit of each model, lower values indicate a better model.The w i column indicates the relative likelihood of each model being the best model of those tested.Analysis represents 1,083 samples from groups with >5 samples.

Table 3 .
Correlations of prey proportions by taxonomic order toPS i values of Salish Sea harbor seals.Only significant correlations (p < 0.0038) are shown.A negative rho value indicates a positive correlation with specialization.Analysis represents 1,083 samples from groups with >5 samples.A total of twelve prey orders were included in the analysis.Data are organized by correlation value.

Table 4 .
Correlations of prey proportions by taxonomic order toPS i values of Salish Sea female harbor seals.Only significant correlations (p < 0.0038) are shown.A negative rho value indicates a positive correlation with specialization.Analysis represents 498 samples from groups with >5 samples.A total of twelve prey orders were included in the analysis.Data are organized by correlation value.

Table 5 .
Correlations of prey proportions by taxonomic order toPS i values of Salish Sea male harbor seals.Only significant correlations (p < 0.0038) are shown.A negative rho value indicates a positive correlation with specialization.Analysis represents 647 samples from groups with >5 samples.A total of twelve prey orders were included in the analysis.Data are organized by correlation value.