Certified personnel of the Punta San Juan Program and Chicago Zoological Society collected pinniped vibrissae during annual health assessments between 2010 and 2016 at the Punta San Juan reserve (PSJ) in southern Peru’s Ica province (15° 22′ S, 75° 11′ W), which forms part of a national natural protected area network called “Reserva Nacional Sistema de Islas, Islotes y Puntas Guaneras,” or RNSIIPG for its Spanish acronym.

To retrieve the most recent growth, mystacial vibrissae were pulled from live animals in November of each year (2010–2012, 2015–2016). Animals were restrained using inhalant anesthetics and/or remotely delivered, multimodal anesthetic drugs as previously described elsewhere (Adkesson et al., 2012, 2019). Vibrissae from sixty-four Peruvian fur seals were collected: 47 adult females and seventeen adult males. The sample years included 2010 (N = 29), 2011 (N = 12), 2012 (N = 11), 2015 (N = 6), and 2016 (N = 6; Table 1). All samples were collected following methods approved under research permits Resolución Jefatural No. 009- 2010-, No. 023- 2011-, No. 022- 2012-, No. 008- 2015-, and No. 019-2016-SERNANP-RNSIIPG issued by the Peruvian National Service of Natural Protected Areas (Spanish acronym SERNANP) and the Peruvian Ministry of the Environment (Spanish acronym MINAM).

Fur seal vibrissae were cleaned using a scrubbing pad and rinsed with deionized water to remove any surface contaminants. After being thoroughly dried at 60°C for a minimum of 24 h, vibrissae were cut into 0.25 cm sections from base (proximal) to tip (distal) to meet mass requirements for the isotopic analysis (0.6–0.8 mg). Segments were placed in individual tin capsules and pelletized. Samples were combusted and analyzed for δ¹³C and δ¹⁵N values at the Smithsonian Institution’s Museum Conservation Institute (Suitland, MD, United States) using a Thermo Delta V Advantage mass spectrometer in continuous flow mode coupled to a Costech 4010 Elemental Analyzer (EA) via a Thermo Conflo IV (CF-IRMS). A set of standards were run for every 10–12 samples. The standards included USGS40 and USGS41 (L-glutamic acid) as well as Costech acetanilide. All samples and standards were run with the same parameters; this included an expected reproducibility of the standards ≤0.2‰ (1σ) for both δ¹³C and δ¹⁵N. Stable isotope values were expressed in terms of δ and were reported relative to the standard reference material, Vienna Pee Dee Belemnite standard for δ¹³C and atmospheric air (N₂) for δ¹⁵N. Following Bond and Hobson (2012), stable isotope values were reported with the standard parts per thousand notation (‰):

where ^jX is the heavier isotope, such as ¹³C, and ⁱX the lighter isotope, such as ¹²C, in the analytical sample (numerator) and the international measurement standard (denominator).

Sea surface temperature anomalies (°C) from the 1+2 Niño index, our proxy for ENSO conditions, were obtained from the National Oceanographic and Atmospheric Administration’s (NOAA) Teleconnections ENSO website¹. NOAA_ERSST_V5 data were obtained from NOAA/OAR/ESRL PSD, Boulder, Colorado, United States². Data were collected from Extended Reconstructed Sea Surface Temperature (ERSST.v5) dataset, which is a global, monthly SST analysis derived from the International Comprehensive Ocean-Atmosphere Dataset (ICOADS). The ERSST.v5 dataset includes information from modern buoy observation, Argos-profiling CTD floats, global drift buoys like ICOADS R3.0 (from R2.5), and Hadley Center Ice-SST version 2 (HadISST2) sea ice concentration (Huang et al., 2014, 2015, 2017, 2018a, b; Liu et al., 2014).

The most recent vibrissa growth is located at the base of the whisker, and an individual vibrissa can represent several years’ growth (Hirons et al., 2001; Ginter et al., 2012). The vibrissae growth rate used for this study was based on Kelleher (2016) northern fur seal (Callorhinus ursinus) growth rate of 0.09 mm/day. Each 0.25 cm segment represented approximately 28 days using this growth rate. With an average of 12.31 ± 3.28 cm in length for each adult, individual vibrissa represented on average 45.6 ± 12 months (3.8 ± 1 years). By analyzing monthly means of individual seals’ isotopic values, a mean for both δ¹³C and δ¹⁵N values can be acquired on a monthly scale for the overall population of seals (i.e., month, year). Across the 64 sampled vibrissae, a total of 12 years’ stable isotope data were recorded.

The statistical package R Core Team (2016) was used for all analyses. Both covariance and correlation analyses were calculated between δ¹³C and δ¹⁵N, δ¹³C and SSTA, and δ¹⁵N and SSTA for both fur seal sexes. Correlations tested the strength of relationships between the stable isotope ratios of the vibrissae and ENSO conditions. Pearson’s correlations were used when the parametric assumptions were met; however, non-parametric conditions, for sample sizes greater than 30, used Kendall’s tau correlations when datasets could not be transformed to meet parametric assumptions. Covariances detected how changes in SSTA were associated with changes in vibrissae δ¹³C and δ¹⁵N values.

Linear Mixed Effects models (Pinheiro and Bates, 2000) with Gaussian distributions and identity link functions were applied to test for a significant linear relationship between SSTA, sex and/or their interaction term to explain changes in δ¹³C and δ¹⁵N values in sampled individuals. Mixed effects models were used to account for the repeated measurements in individuals over time. These models had individual identity (i.e., Animal ID) as a random effect to account for repeated measures of each response variable on the different segments of each vibrissa (Franco-Trecu et al., 2014). We constructed linear mixed effects models using nlme (Pinheiro et al., 2017) since we consider that each individual can have a different isotopic signature and variation depending on the month and year of sampling. Model diagnostics were verified to not violate Gaussian assumptions. We ran backwards selection, progressively removing the non-significant terms, and compared sequential models. Best model was selected based on the lowest Akaike’s Information Criterion (AIC) score. We applied a bootstrapping routine (n = 1000) to estimate 95% confidence intervals for the population average for each sex with lme package (Bates et al., 2015) in R statistical environment (R Core Team, 2016).

The proxy for ENSO conditions (e.g., monthly SSTA readings) in this study was collected from the closest Niño index, Niño 1+2, which is located more than 800 km from the Punta San Juan rookery. Although this is a substantial distance, previous ENSO years have shown documented impacts in the coastal Peruvian ecosystem based on the index, including the 1997/98 anchoveta fishery collapse which is believed to have caused a Peruvian fur seal mass mortality event documented in PSJ, Peru (Arias-Schreiber and Rivas, 1998; Cárdenas-Alayza, 2012). The analysis of δ¹³C and δ¹⁵N values from vibrissae in this study provided a timeline of trophic and production changes in this dynamic ecosystem, one where ENSO conditions fluctuate, and organisms must adapt to survive potentially stressful conditions. The PSJ population’s vibrissae recorded significant differences in δ¹³C and δ¹⁵N values with varying SSTA among years. SST is correlated with primary production changes in the ecosystem but may also influence the movement patterns of fishes and cephalopods as well. It is presumed that prey availability played a significant role in these findings; however, primary producer and potential prey stable isotope values are needed for further evaluation.

According to the SSTA from the Niño 1+2 region, the strongest ENSO on record occurred during this study in years 2014–2016, and the biotic results of the event were reflected in vibrissae collected in 2015 and 2016. Results in this study suggest that the 2014–2016 ENSO event impacted the trophic dynamics of the Peruvian coastal ecosystem in ways that had not been seen throughout the prior 9 years evaluated (2004–2013). The significant inverse correlation between δ¹⁵N and δ¹³C values during the 2014–2016 extreme El Niño event demonstrates sharply declining δ¹⁵N while δ¹³C slowly increases. This distinct change in fur seal trophic dynamics may correspond to a dietary shift linked to an influx of oceanic prey species during ENSO conditions. Meanwhile, the incremental increase in carbon enrichment may reflect a change in the dominate primary producers supporting the food web during nutrient limitation. This is evidenced in decreased surface chlorophyll concentrations, more pronounced during temperature extremes (Espinoza-Morriberón et al., 2017). Additionally, alternate foraging grounds may be linked to a weakening and delay in offshore transport of cold waters, as reported in previous ENSO events, may contribute to this isotopic variation (Roy and Reason, 2001; Montecinos and Gomez, 2010). These findings could infer that strong ENSO cycles may alter the specialized diet of the Peruvian fur seal.

A wild animal’s trophic discrimination factor varies according to the species’ diet and body condition, both of which would be potentially altered during ENSO periods. The comparison of consumer tissue δ¹³C values over time provides information regarding ocean production in foraging grounds throughout fluctuating ENSO conditions (France and Peters, 1997; Hirons et al., 2001; Kurle and Worthy, 2001). From 2014 through 2016 as SSTA increased, declines in primary production were expected, reflected as lower δ¹³C values. However, the mixed δ¹³C values observed during this time (Figure 6) may also suggest a change in dominant primary producers in the system as well as use of alternate foraging locations (Cárdenas-Alayza pers comm). However, an investigation into the region’s trophic structure would be necessary to support these inferences. Production variability and behavioral modifications may explain how fur seals are able to survive the extreme food shortages presumed to occur during ENSO events.

Related studies with sampling efforts in the Humboldt Current, such as Espinoza et al. (2017), revealed more depleted ¹³C values further offshore (Sydeman et al., 1997; Hill et al., 2006; Miller et al., 2008). Increased distance from shore corresponded to increased faunal ¹³C depletion which is common in upwelling ecosystems (Sydeman et al., 1997; Miller et al., 2008; Espinoza-Morriberón et al., 2017). Primary production decreases from coastal to neritic to oceanic waters (France and Peters, 1997; Hill et al., 2006; Espinoza et al., 2017). ENSO conditions could also contribute to fluctuating patterns in primary production during warm and cold phases, subsequently compromising resource availability within an environment.

The δ¹⁵N in consumer tissues relative to diet provide information regarding potential food web linkages (France and Peters, 1997; Hirons et al., 2001; Kurle and Worthy, 2001). The adult Peruvian fur seals all showed a similar pattern in δ¹⁵N values through time, a long-term gradual decline followed by a quick increase. The δ¹⁵N values from 2004 to 2009 were lower (Δ−2.41‰), followed by a rapid rise (Δ+2.17‰) from 2010 to 2011. A similar pattern persisted from 2012 through the end of 2016, revealing an oscillation of a 5-year decline (Δ−4.35‰) followed by a 9-month (Δ+1.56‰) increase. Males showed consistently higher δ¹⁵N values than females regardless of changes in SSTA recorded in the Nino 1+2 index (Figure 7). Whether this pattern specifically corresponds to ENSO conditions is uncertain; however, as El Niño conditions strengthened and lasted longer, δ¹⁵N continually decreased.

Warmer SSTs correlated to lower δ¹⁵N values. When evaluating ENSO conditions by phase (i.e., cold, normal, and warm), significant correlations between warmer phases and lower δ¹⁵N values as well as colder phases and higher δ¹⁵N values were identified. During warmer phases, the lowest δ¹⁵N values were detected, while consistently higher δ¹⁵N values were detected during cold phases. The increasing magnitude of ENSO conditions revealed significant correlation to lower δ¹⁵N values; stronger or more abnormally warm conditions revealed below average δ¹⁵N values, whereas norm, weak, and moderate magnitudes revealed mean population δ¹⁵N values (Edwards, 2018). Previous trophic study methodology on South American fur seal populations in Peru were done via scat analyses. A valuable diet evaluation was done from 1982 to 1988, which similarly to us, encompassed varying cold and warm periods, including two El Niño events. Scat findings in fur seal diet varied as a result of anchoveta availability; a wider range of prey was seen in scat when anchoveta were less available (Majluf, 1989). It can be hypothesized that mean Peruvian fur seal δ¹⁵N values could be reflecting these wider prey ranges during periods of anchoveta absence in foraging grounds in relation to warmer SSTA. This adaptive behavior has also been seen in the Peruvian population of South American sea lions during the 1998 El Niño, where they fed on a variety of prey when their primary prey source(s) were scarce (Soto et al., 2006).

Both δ¹³C and δ¹⁵N values were significantly different between adult male and female Peruvian fur seal. Female ¹³C and ¹⁵N values were significantly lower than males so foraging differences between sex can be inferred. Male Peruvian fur seals are known to forage in larger ranges than adult females with pups; females must return to the rookery to nurse their pup, so they have abbreviated forage trips in comparison to males (Punta San Juan Program, unpublished data). While no gender-size and foraging relationships are known for Peruvian fur seals, they have been observed in South American sea lions. Male sea lions from northern Patagonia exploit benthic and deeper foraging grounds than smaller females (Campagna et al., 2001; Drago et al., 2009), although Drago et al. (2009) noted differences in foraging habits between the sexes are not constant over time. The potential for foraging in dissimilar water masses with contrasting prey composition may explain why male δ¹³C and δ¹⁵N values were higher than females (0.51‰, 1.36‰, respectively). Whether this is due to physical ability, metabolic capability, foraging location, depth, or simply different prey, these differences cannot be explained without potential prey source isotope signatures and/or tracking data (Edwards, 2018).