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

Abstract Predator–prey interactions are critical in understanding how communities function. However, we need to describe intraspecific variation in diet to accurately depict those interactions. Harbor seals (Phoca vitulina) are an abundant marine predator that prey on species of conservation concern. We estimated intrapopulation feeding diversity (variation in feeding habits between individuals of the same species) of harbor seals in the Salish Sea. Estimates of feeding diversity were examined relative to sex, month, and location using a novel approach that combined molecular techniques, repeated cross‐sectional sampling of scat, and a specialization metric (within‐individual consistency in diet measured by the Proportional Similarity Index (PSi)). Based on 1,083 scat samples collected from five haul‐out sites during four nonsequential years, we quantified diet using metabarcoding techniques and determined the sex of the scat depositor using a molecular assay. Results suggest that intrapopulation feeding diversity was present. Specialization was high over short periods (24–48 hr, PSi = 0.392, 95% CI = 0.013, R = 100,000) and variable in time and space. Females showed more specialization than males, particularly during summer and fall. Additionally, demersal and benthic prey species were correlated with more specialized diets. The latter finding suggests that this type of prey likely requires specific foraging strategies and that there are trade‐offs between pelagic and benthic foraging styles for harbor seals. This differential feeding on prey species, as well as between sexes of harbor seals, indicates that predator–prey interactions in harbor seals are complex and that each sex may have a different impact on species of conservation concern. As such, describing intrapopulation feeding diversity may unravel hitherto unknown complex predator–prey interactions in the community.


| INTRODUC TI ON
Predator-prey relations are an integral force in ecosystem functions and are often key to understanding how communities interact.
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 interindividual 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). Theoretical models of predator-prey systems are also much more likely to have chaotic dynamics when the predator lags in switching preferred prey (i.e., is highly specialized) when there are high rates of changes in prey populations (Abrams & Matsuda, 2004). Further, dynamics such as prey density relative to one another and indirect interactions between prey 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 communities and the ecosystem, such as prey community structure and total primary production (Harmon et al., 2009). In addition, differences in a single species' population structure (i.e., feeding diversity) 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., 2003Bolnick et al., , 2011Dall, 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 it requires observing many predation events in multiple individuals across numerous 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 (a) correlations of individual diet diversity with the sex of the predator and month in which the predation occurred, and (b) to test whether short-term diet diversity was related to the consumption of particular prey species and their preferred habitat.
The harbor seal (Phoca vitulina) is a well-studied species with potentially large impacts on its ecosystem. Harbor seals have the largest worldwide distribution of any pinniped in coastal areas (Teilmann & Galatius, 2018), are an abundant marine predator in the Salish Sea (Jeffries, Huber, Calambokidis, & Laake, 2003;Olesiuk, 2009), and appear to have reached carrying capacity in that region (Jeffries et al., 2003;Olesiuk, 2009). Because harbor seals are abundant in the ecosystem and feed on a wide range of species, they have significant impacts on prey populations (Howard, Lance, Jeffries, & Acevedo-Gutiérrez, 2013;Lance, Chang, Jeffries, Pearson, & Acevedo-Gutiérrez, 2012;Olesiuk, Bigg, Ellis, Crockford, & Wigen, 1990). Some of their prey species are of high conservation concern, such as Pacific salmon (Oncorhynchus spp.; hereafter referred to as Salmoniformes), rockfish (Sebastes spp.), and Pacific herring (Clupea pallasii pallasii, Valenciennes, 1847) (Bjorkland et al., 2015;Bromaghin et al., 2013;Lance et al., 2012). There is special interest surrounding their impact on Chinook salmon (Oncorhynchus tshawytscha) as their consumption of this species appears to have increased over the last few decades (Adams et al., 2016;Chasco et al., 2017). Chinook salmon are of special concern given their cultural and economic importance in the Pacific Northwest, and their role as prey for another population of concern: southern resident killer whales (Orcinus orca) (Ford, Ellis, Olesiuk, & Balcomb, 2010;Hanson et al., 2010).
Harbor seal's effect on Salmoniformes is also of interest because they eat both juvenile and adult individuals (Thomas, Nelson, Lance, Deagle, & Trites, 2017). Eating juveniles may have considerable impacts on populations of Chinook, coho (Oncorhynchus kisutch) and steelhead (O. mykiss), as survival during the first several months at sea is believed to be the primary factor limiting population abundance and productivity (Beamish et al., 2010;Kendall, Marston, & Klungle, 2017;Neville, Beamish, & Chittenden, 2015).
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, several lines of evidence support the notion that harbor seal populations in the Salish Sea may instead be comprised of individual specialists (Bjorkland et al., 2015;Schwarz et al., 2018). They are viewed as central place foragers due to their high haul-out site fidelity (Peterson, Lance, Jeffries, & Acevedo-Gutiérrez, 2012;Suryan & Harvey, 1988).
Central place foraging, in combination with harbor seal's high abundance in the region (Jeffries et al., 2003;Olesiuk, 2009), makes high intraspecific competition likely, which in turn increases the likelihood of individual specialization . Because harbor seals impact prey populations of conservation concern, it is important to examine the prevalence of this potential specialization to better understand the dynamics between seals and their prey.
Prey composition and foraging dive behavior of harbor seals in the Salish Sea vary relative to habitat, sex, and time of year 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 Luxa & Acevedo-Gutiérrez, 2013). Further, males and females consume different prey (Bjorkland 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). Harbor seal behavioral differences between seasons and sexes exist outside of the Salish Sea as well. Females in Scotland spend more time at sea after lactation, indicative of increased foraging, while males increase their time spent hauling out during the molt, indicative of decreased foraging (Thompson, Fedak, McConnell, & Nicholas, 1989).
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, noninvasive, and time-efficient way to obtain large sample sizes with species-level taxonomic resolution (Deagle et al., 2005(Deagle et al., , 2019Rothstein 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.
Measuring long-term individual specialization via longitudinal studies is logistically complicated, expensive, and invasive when studying large, mobile, wild animals, particularly marine mammals.

F I G U R E 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.
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: (a) How do the factors Sex, Month, Location, and Year affect cross-sectional estimation of a specialization metric? (b) 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 molecular techniques and supplemented with hard-part techniques to determine the proportion of juvenile and adult Salmoniformes consumed. 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 (areas were seals congregate out of water) in the Salish Sea over a period of four nonsequential years ( Figure 1). Haul-outs varied in seal population size as well as by habitat type (Table 1). Not all sites were visited every year and the months during which each site was visited varied between years ( 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-μm 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°C 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°C 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°C for later use in prey hard-part (i.e., bones, otoliths) 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°C 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) and Schwarz et al. (2018). The two probes that we used targeted the paralogous zinc finger X (ZFX) and zinc finger Y (ZFY) genes. All primers and

| 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 from Location denotes where the group of scat was collected. Month, Year, and Season show when the group of scat was collected. The minimum prey density column indicates the lowest occurring prey proportion within a single scat within the group. The theoretical minimum column indicates the theoretical minimum PS i that was assigned to each group. The theoretical minimum was calculated by dividing one by the total sample size of each group. N indicates the total number of samples within each group. SI indicates the average PS i within each group. The use column indicates whether or not the group was used in downstream analysis based off of sample size (groups with < 5 samples were excluded). The Shannon-Weaver column indicates the Shannon-Weaver index assigned to each group. This index was calculated by averaging prey proportions from each scat within each group and subsequent use of the diversity function in the VEGAN package (Oksanen et al., 2018) in R 3.3.1.

TA B L E 1 (Continued)
extracting DNA out of the DNA slurry using QIAGEN QIAamp DNA Stool Mini Kit) 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). These primers sets were designed to amplify chordate and cephalopod sequences and are as follows: Chord 16S F (GATCGAGAAGACCCTRTGGAGCT), (Deagle et al., 2010;Thomas et al., 2016).
Amplicons were labeled 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) Thomas et al. (2017). Notably, this analysis allowed differentiation between the proportion of adult and juvenile Salmoniformes consumed. The percentage of juvenile versus adult salmon was determined using the method described in Thomas et al. (2017). Briefly, the DNA percentages assigned to each salmon species were split based on the ratio of juvenile versus adult salmon documented during hard-part analysis. The ratio was determined using data from either the sample, month, or season, dependent on available data. This information was used to investigate whether juvenile and adult salmon were being consumed differently with respect to specialization habits (addressed in Section 2.6). This additional analysis was performed due to the high economic and environmental impact salmon species have in the region and previous evidence that harbor seals have a significant impact on salmon populations at the juvenile life stage .

| 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). Cross-sectional diet studies assume that each sample consists of multiple prey items and that each item represents an independent capture decision . It is also assumed that the sampled diet is representative of the complete diet of the individual . If both assumptions are met cross-sectional samples allow for estimates of true individual specialization. However, it is unlikely these assumptions are always met for pinnipeds as their prey can be clumped (i.e., forage fish), which would result in one foraging decision influencing the next. Further, it is unlikely that a scat represents the entire diet of an individual as a single scat is only indicative of the last few foraging bouts (Bowen & Iverson, 2012). In these cases, one cannot calculate absolute individual specialization using cross-sectional sampling Novak & 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, month, 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 calculating PS i , not proportions, as each count is assumed to represent an independent prey capture decision 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, calculating PS 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 sample sizes, is noninvasive, 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 111 groups (Table 1). 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 1). 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( i 1− i )) on the PS i values to normalize them. This transformation was necessary because PS i is bounded by a theoretical minimum and one, which affects the variance distribution (Bolnick et al., 2002;Sokal & Rohlf, 2012). When numbers are bounded, the variance distribution is shifted toward the mean (Sokal & Rohlf, 2012). A logit transformation is an excellent choice for addressing this shift because it extends the tails of the distribution more than other alternatives (Warton & Hui, 2011).

F I G U R E 2
Proportion of female harbor seal scat identified during each month at each site. If no dot is present, no scat was collected at the site in that month. Proportions were calculated by pairing male and female groups for analysis from the same month, location, and year. The number of samples in the female group was divided by the total number of samples in both the female and male group All models were tested in the R 3.3.1 package lme4 (Bates, Mächler, 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 r 2 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 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 was visualized ( Figure 2).

| 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. To investigate differential impacts on juvenile versus adult salmon, all salmon species where spilt into proportion juvenile and proportion adult based on the hard-part collection and analysis (see Section 2.3) and summed as two separate orders ("Juvenile Salmoniformes" and "Adult Salmoniformes"), in addition to the Salmoniformes order. 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 whether 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 whether 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 nonsequential years, at five different locations, we quantified the diet of 1,520 scat samples. The sequence data associated with these samples are stored at: https://figsh are. TA B L E 2 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 Chain, Cowichan Bay, Comox, and Fraser River were calculated from Olesiuk et al. com/s/01134 57d40 81727 aacac. (Thomas et al., in review). We successfully determined the sex of the depositor for 1,145 of those scats (75% success rate). The breakdown of prey proportion within each scat as well as sex assignment for these samples is located on Dryad (10.5061/dryad.59zw3r257). The number of scat with successful sex determination varied by location and month (Table 2).
Samples with successful sex determination were then binned by the factors Sex, Location, Month, and Year to form unique groups for analysis (Table 1) there was also a correlation between average PS i and theoretical minimum PS i (rho = −0.231, p = .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) (1| indicates random effects, Table 3). The r 2 value and residual plots indicate that this model fits the data well (Table 3, Figure 4). 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 3). Removal of the interaction term also caused a decrease in goodness-of-fit measurements (Table 3). Further, the interaction terms for Sex and the Months of August and October were significant (t = 2.86, 2.68, p = .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 = .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 5). In summer and fall, males had relatively lower levels of specialization than females ( Figure 5). 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 6).
Visual inspection of the data by month suggested males had a decrease in relative specialization in July through October (Figure 7).

Based on 95% confidence intervals of logit-transformed PS i values,
PS i during these months only overlapped with April ( Figure 7). 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 7). This trend varied in intensity by location ( Figure 8). The described pattern was reflected most strongly in Belle Chain, Comox, and Fraser River (Figure 8). However, because scat were not collected at Baby Island after July, no comparison could be made at that location ( Figure 8, Table 1).

| Correlations between prey items and relative specialization
Correlation analysis between diet proportions of prey orders and PS i revealed that ten out of fifteen orders (counting juvenile and adult Salmoniformes as additional separate "orders") showed significant correlations (alpha = 0.0038, Bonferroni correction, Table 4 Correlations performed with just data from female scat showed similar patterns. All orders of prey showed the same relationship with PS i , or were no longer significant, such as Adult Salmoniformes, Salmoniformes, Chimaeriformes, Rajiformes, and Batrachoidiformes (Table 5). Correlations performed with only male scat once again showed the same patterns. Orders that no longer had significant relationships were Rajiformes, Batrachoisiformes, Pleuronectiformes, and Clupeiformes (Table 6). However, Adult Salmoniformes and Salmoniformes once again showing significant correlations ( Table 6).
The correlation between proportion of benthic species in each scat and PS i suggested a positive relationship between the amount of benthic species consumed and the level of relative specialization (rho = −0.384, p > .001). A similar relationship was observed when the female and male data sets were examined separately (rho = −0.407, p > .001; rho = −0.35, p > .001, respectively).

| D ISCUSS I ON
We successfully assigned PS i values to 1,083 scat collected from five different locations over the course of four nonsequential years (Table 1, Figure 1). 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 versus 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 each sex has different impacts on prey species of concern.

| Estimated level of intrapopulation feeding diversity
Our data confirmed intrapopulation feeding diversity across the spatial (hundreds of km) and temporal (years) scales that the scat samples represented (average = PS i 0.399, 95% CI = 0.026, R = 100,000).
These data demonstrate intrapopulation feeding 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 short-term specialists versus long-term specialists.

| Importance of Month, 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 3). This pattern makes intuitive sense as the type of prey eaten by harbor seals 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 2). 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 (Bjorkland 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 et al., 2012;Schwarz et al., 2018;Wilson et al., 2014). Therefore, we propose the following theoretical resource distribution: Males have more overlap between individuals (less specialized), while the females have less overlap between individuals (more specialized) ( Figure 5) Note: 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.

TA B L E 3 GLMM models of prey specialization in Salish Sea harbor seals (1| indicates random effects)
in this pattern appear to be associated with prey type (which will be addressed in Section 4.3).
Including an interaction term between Month and Sex increased the goodness of fit of the model (Table 3). 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 5), indicating that the behavior of both sexes was similar in the spring but diverged in the summer and fall. During the late summer and 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 . This resource could be rich enough that it is beneficial for a majority of seals, resulting in less need for specialization. Further, males consume more Salmoniformes than females (Shwarz et al., 2018) which could result in the larger spike observed in generalist behavior in males.
Our data also suggest that location explained a large amount of variance in relative specialization. The random factors of Year  (Figures 5, 6).
The importance of location as a factor in explaining variation of relative specialization could be due to varied levels of prey diversity in different environments, given that prey availability affects the level of specialization . For example, harbor seal scats from haul-outs in estuaries have higher prey diversity than scat from haul-outs outside estuaries Luxa & Acevedo-Gutiérrez, 2013 Month. The left graph shows females; the right graph shows males. A lower value indicates more specialization estuarine and nonestuarine habitats, and Belle Chain is considered a rocky reef environment. Therefore, it is likely that variation in habitat types caused differences in prey availability that offered more or less options to harbor seals in the area, which could subsequently affect the level of competition and, ultimately, specialization. The spatiotemporal variation in relative specialization throughout the region suggests widespread intrapopulation feeding diversity. This knowledge can inform the design of future studies and act as a starting place to investigate the impacts of intrapopulation feeding diversity on prey species in the ecosystem.

| 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 became. This pattern was observed in the full dataset, as well as when female and male data were considered separately (Tables 4-6).
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 described pattern of benthic prey correlation with relative specialization seems to hold true for both females and males (Tables 4-6). This result indicates that specialization patterns linked to prey species were reflective of foraging strategies specific to the ecology of prey species, and not just indicative of different diet preferences of 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, we argue that sex determines the feeding strategy to be employed and that the prey species captured by each feeding strategy determines the level of specialization. The consistency of benthic prey being associated with a relatively specialist diet, and pelagic prey being associated with a relatively generalist diet in our three datasets (complete, only female, and only male) suggests that the ecology of the prey species rather than the sex of the seal was driving the observed pattern. This idea is supported by other studies. Individual male harbor seals in Nova Scotia use different behaviors when foraging for benthic versus pelagic prey (Bowen, Tully, Boness, Bulheier, & Marshall, 2002). Further, large harbor seals are more likely to forage in pelagic environments regardless of sex (Bjorkland et al., 2015).
Because harbor seals display slight sexual dimorphism, with females being slightly smaller than males (Teilmann & Galatius, 2018), there is the potential for trade-offs between speed and maneuverability F I G U R E 8 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 PS 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 within the population. If that is the case, then females would be slightly slower and more maneuverable than males. We hypothesize that this combination of traits would be more successful in a benthic environment to deal with variations in the seafloor that benthic prey relies on to escape. If this idea is correct, then the sex of an average harbor seal predisposes an individual to be more effective in one type of environment than the other.
If prey species ecology is driving specialization levels, it is especially interesting to consider harbor seal consumption of juvenile Salmoniformes. As a group, juvenile Salmoniformes did not correlate with the relatively more generalist diet. However, when split into species, juvenile sockeye (Oncorhynchus nerka) did correlate with this diet (rho = 0.22, p = .004). This correlation 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 a pelagic foraging strategy. 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 discernible 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, it is possible that some scat collected came from the same individual. However, this possibility is small given than Rothstein et al. (2017)  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 semiquantitative, biases appeared 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., 2014Thomas et al., , 2017 Note: A Bonferroni correction was used to set alpha at 0.0038. Thus, p < .0038 is significant and is denoted by an asterisk next to the order name. A negative rho value indicates a positive correlation with specialization. Analysis represents 1,083 samples from groups with > 5 samples. Data are organized by correlation value. many benefits and is a valuable molecular technology that should be applied in the future.

| CON CLUS ION
We have shown that intrapopulation feeding diversity occurs in Salish Sea harbor seals between locations, seasons, and sexes.
Females displayed higher levels of relative specialization than males.
However, in both female and male harbor seals, benthic prey were more commonly associated with a more specialized diet, suggesting the prey's ecology has an important role in driving the level of specialization. These different impacts of male versus female on benthic versus pelagic prey should be considered henceforth when management addresses harbor seal interactions with species of concern.
Further, we demonstrated how the use of molecular prey barcoding from scat allows for high taxonomic and spatiotemporal resolution of relative individual specialization. The resulting large-scale examinations of intrapopulation feeding diversity uncovered previously unknown complex interactions between predators and prey.

ACK N OWLED G M ENTS
Many individuals and agencies contributed to this study. For financial support that allowed the collection and processing of samples, improved the quality of this study. Further, we thank the laboratory of Dr. C. Moyer for allowing use of their quantitative thermocycler.

CO N FLI C T O F I NTE R E S T
The authors declare no conflict of interest.

DATA AVA I L A B I L I T Y S TAT E M E N T
All diet data used are publicly available. The sequence data are stored at: https://figsh are.com/s/01134 57d40 81727 aacac (Thomas et al., in review). The breakdown of prey proportion within each scat as well as sex assignment for these samples is located on Dryad Note: A Bonferroni correction was used to set alpha at 0.0038. Thus, p < .0038 is significant and is denoted by an asterisk next to the order name. A negative rho value indicates a positive correlation with specialization. Analysis represents 647 samples from groups with > 5 samples. Data are organized by correlation value.