Apex predators are typically considered dietary generalists, often masking individual variability. However, individual specialization—consistent differences in diet and foraging strategies among individuals— is common in apex predators. In some species, only a few specialized individuals can significantly impact prey populations. Leopard seals (Hydrurga leptonyx) are apex predators important to the structure and function of the Southern Ocean ecosystem. Leopard seals are broadly described as generalists, but little is known about their trophic ecology at the population or individual level. We analyzed δ¹³C and δ¹⁵N profiles in whiskers (n=46) from 34 leopard seals from an important aggregation in the Western Antarctic Peninsula to assess population and individual trophic variation. We also evaluated individual consistency across years using repeat samples from 7 seals collected over 2-10 years. We compared population and individual isotopic niche space and explored drivers of intraspecific variation in leopard seal trophic ecology. We find that leopard seals have a broad trophic niche (range: 6.96-15.21‰) and are generalists at the population level. However, most individuals are specialists (59% for δ¹⁵N and δ¹³C), with only a few generalists (13% for δ¹⁵N, 6% for δ¹³C). Furthermore, individual seals specialize at different trophic levels, resulting in niche partitioning. Most variation in trophic ecology is driven by individual specialization, but sex and mass also contribute. We also find that some seals specialize over time, consistently foraging at the same trophic level, while others switch trophic levels within and between years. Long-term specialization by only a few leopard seals has likely contributed to the decline of a significant local mesopredator colony, the Antarctic fur seal. Our findings show the importance of examining individual specialization in leopard seals across their range to understand their impact on other prey populations. More broadly, this approach should be applied to other apex predator populations, as a few specialists can significantly impact ecosystems.

Sample Collection

Leopard seals were sampled between January and May from 2013 to 2023 at the U.S. Antarctic Marine Living Resources (AMLR) Program ecological monitoring site at Cape Shirreff, Livingston Island in the WAP. Leopard seals were sedated (Pussini and Goebel, 2015; Krause et al., 2016) to allow for the collection of morphometric data (standard length [cm], girths, and mass [kg]), life history traits (sex, age class), and whiskers. Whiskers were collected from 34 leopard seals (28 females, 6 males). Seven females were opportunistically resampled without a full capture (i.e., morphometric data were not collected) 1-3 times after their first handling, resulting in 12 additional whisker samples. In total, we analyzed 46 whiskers (40 females, 6 males). A subset of these whiskers (n=18 from field seasons 2018 and 2019) were previously analyzed by (Charapata et al., 2023).

Stable Isotope Analysis (SIA)

Whiskers were wiped with 1:1 ethanol:methanol solvent, sonicated for 30 mins in distilled water, and air dried. Whiskers were measured and sectioned into 0.5-3 mm increments (from root to tip) for a targeted weight of ~0.3 mg (Charapata et al., 2023). Carbon and nitrogen stable isotope analysis was performed at Baylor University using an Elemental Analyzer 4010 Elemental Combustion System paired with a Conflow IV interphase (Thermo Scientific) and Thermo Delta V Advantage continuous flow Isotope Ratio Mass Spectrometer. Whisker nitrogen (ẟ¹⁵N) and carbon (ẟ¹³C) isotope values are expressed in delta notation (ẟ) in units of per mil (‰). Additionally, isotope values are reported as the ratio of the heavy to light isotope relative to international standards--atmospheric nitrogen and Vienna Peedee Belemnite, respectively--using the following equation:

ẟX  =  [(R_sample/R_standard) – 1]*1000

where X is the 13C or 15N and R is the corresponding ratio of 13C/12C or 15N/14N. A two-point calibration curve for calculating ẟ¹⁵N and ẟ¹³C values of samples was established using USGS-40 and USGS-41A international standards. The accuracy and precision of isotopic measurements were calculated based on the long-term mean and standard deviation (SD) of 244 replicates of an internal lab standard (Acetanilide, reported ẟ¹³C =−29.53±0.01‰, ẟ¹⁵N=1.18±0.02‰) measured during each analytical run (n=3 replicates/run). The replicate grand averages obtained were very close to (ẟ¹³C=−29.42±0.08‰) or within the range (ẟ¹⁵N=1.30±0.17‰) of analytical uncertainty of reported values. We measured the atomic C:N ratio for every whisker segment with acceptable atomic ratios ranging from 3.0-4.0 (Newsome et al., 2009; Kernaléguen et al., 2012; Charapata et al., 2023). Nearly all whisker segments had acceptable atomic C:N ratios (3.53±0.14, range: 2.9-4.0). Twelve whisker segments were excluded for having ratios outside this range.

Time Stamping

Leopard seals molt and shed their whiskers annually; therefore, whiskers represent growth over a few months and up to one year (Rogers et al., 2016). We timestamped whisker segments based on leopard seal whisker growth characteristics using the Von Bertalanffy growth model (von Bertalanffy, 1938; Rogers et al., 2016) following the approach outlined by (Charapata et al., 2023).

Data Analysis

All data were tested for normality and homogeneity of variance before analysis. Results are reported as mean ± standard deviation (SD) unless otherwise stated. We performed all analyses using R (R Core Team, 2022) with RStudio (Team, 2021) and JMP (SAS).

Population-Level

Population-level analyses included a total of 46 leopard seal whiskers. Each whisker was treated separately based on preliminary data showing inter-annual isotopic variability. We calculated the population-level mean, SD, and range of δ¹⁵N and δ¹³C values. We used variance component analysis (VCA) to calculate between- and within-individual population variation. Total variance in stable isotopes (“between individuals” variation) indicates variation among individuals in a population, while variance in stable isotopes along the whisker (“within-individual” variation) indicates variation of an individual (Bearhop et al., 2004; Newsome et al., 2009; Hückstädt et al., 2012). We applied the Stable Isotope Bayesian Ellipses in R SIBER package (Jackson et al., 2011) to determine population isotopic niche width. We used the standard ellipse area corrected for small sample sizes (SEAc) for individual whisker(s) as the metric for calculating the population isotopic niche area. We also calculated a population-level SEAc and total area (TA) using the pooled δ¹⁵N and δ¹³C values from all whisker segments (n = 46 whiskers; 2,198 segments) to compare our results with a previous study on leopard seals (Botta et al., 2018).

Individual Specialization

We used two approaches to evaluate the isotopic variation at the individual level. First, we calculated individual isotopic niche with SEAc and TA estimates for each individual’s whisker(s) using the δ¹⁵N and δ¹³C values of the whisker segments; this allowed us to visualize and assess each individual’s range of trophic levels and foraging locations collectively. Next, we calculated δ¹⁵N and δ¹³C specialization indices for each whisker to describe the variance in δ¹⁵N and δ¹³C and calculate the degree of individual specialization (Bolnick et al., 2002; Lewis et al., 2022); this allowed us to separately assess the variation in δ¹⁵N and δ¹³C. The degree of specialization was calculated using the equation:

SI = INW/ (INW + BINW)

Within the δ¹⁵N specialist category, some individuals consistently exhibited high δ¹⁵N values, while others consistently had medium to low values. Therefore, we performed agglomerative hierarchical clustering to determine whether there were subgroups within our δ¹⁵N specialist isotope data using the ‘agnes’ function in the R package cluster (Kaufman and Rousseeuw, 2009). To determine the optimal number of clusters, we used the Dunn index, which differentiates between sets of clusters that are compact and well separated (Supplementary Figure 1A). We found two distinct clusters (Supplementary Figure 1B): high trophic level specialists (H-Specialist) and medium to low trophic level specialists (ML-Specialist).

Trophic Variation & Overlap

We examined variation and niche overlap in isotopic signatures as a function of sex, body mass, and degree of individual specialization. We focused these analyses solely on δ¹⁵N because (1) we were interested in trophic level variability, and (2) leopard seals from Cape Shirreff tend to remain in the near-shore habitat and are primarily coastal foragers (Krause et al., 2015; Kienle et al., 2022b). To investigate variation in δ¹⁵N, we ran a linear mixed-effects model (LMM; (Pinheiro and Bates, 2000) using the ‘lmer’ function from the lme4 package (Bates et al., 2014). This model treated the average δ¹⁵N values as response variables with sex, mass, the interaction of sex and mass, and δ¹⁵N specialization category (H-Specialist; ML-Specialist; Intermediate; Generalist) as fixed effects and individual as a random effect. Model selection was performed using the R package MuMIn (Barton and Barton, 2015) based on the smallest Akaike information criterion corrected for sample size (AICc). The model with the lowest AICc had the highest support, and models with ΔAICc < 2 were considered to have substantial support (Anderson and Burnham, 2002; Franklin et al., 2002). Goodness-of-fit for each model was estimated using marginal (R² LMM(m)) and conditional (R² LMM(c)) coefficients of determination, indicating variance explained by fixed effects alone and by both fixed and random effects, respectively (Nakagawa and Schielzeth, 2013). We examined the contribution of each fixed effect of our top models by looking at the estimated coefficients and p-values and then used ANOVAs on each of our top models. Pairwise comparisons were performed using the ‘emmeans’ function from the emmeans package with Tukey adjustment for multiple testing δ¹⁵N and δ¹³C (Lenth and Others, 2022). We also used Spearmen’s correlation to assess the relationship between our continuous variables (average δ¹⁵N, mass) and the relationship between mass and individual niche width (SEAc). We assessed niche differences and overlap between δ¹⁵N specialization categories and sexes using the proportion of paired SEAc shared; this was calculated with the ‘maxLikOverlap’ function from the SIBER package (Jackson et al., 2011).

Between-Year Variability

To evaluate between-year variability among repeat individuals (n=7), we used SIBER to visualize data and quantify percent overlap between isotopic niches, assessing the similarity/dissimilarity in isotopic composition between years for each individual. We also used a quadratic discriminant analysis (QDA) to simultaneously evaluate δ¹³C and δ¹⁵N (Koehler et al., 2019; Smith et al., 2021). QDA is appropriate for analyzing data that are unequally sampled across years and have unequal variance; this allowed us to effectively assign isotope signatures to specific years for each seal. We considered QDA to be unsuccessful in assigning individual isotope data to their respective years if the results were £ 70% (Koehler et al., 2019; Smith et al., 2021), suggesting that the data was too similar to accurately assign it to specific years.

Temporal Changes

To investigate yearly and monthly trends in δ¹⁵N data, we used generalized additive models (GAMs) with isotopic signatures of whisker segments as the response variable. Year and month from timestamped whiskers were used as temporal predictors, with individual as a random effect, using the formula: δ¹⁵N ~ s(Months, k = 10, bs = "cc") + s(Year, k = 10, bs = "tp") + s(Individual.ID, bs = "re"). The “s” represents the smooth functions and “k” represents the number of basis functions used in the smoothing function. For month we used a cyclic spline (bs = 'cc') and for year a thin plate regression spline (bs='tr'; (Wood, 2003). GAMs and corresponding model estimates were conducted using the mgcv and modelbased R packages (Wood, 2017; Makowski et al., 2020). We evaluated isotopic linear temporal fluctuations in significant variables by using a grid approximation, accompanied by CI 95% (Makowski et al., 2020). Because our dataset is overrepresented in some years/months and underrepresented in others, we created a customized prediction grid based on GAM models fitted on the observed data. The conditional expectations generated by the simulated homogeneous dataset (i.e. simulated data for all years and months for all Individual.ID) allowed parameter estimation using the fitted model. The predicted values were used to calculate the first derivative of the response variable and estimate the linear slope of the isotopic signatures to identify the temporal windows where significant linear increase or decrease of isotopic values occurred in time.