Peruvian Fur Seals as Archivists of El Niño Southern Oscillation Effects

Peru’s coastal waters are characterized by significant environmental fluctuation due to periodic El Niño- La Niña- Southern Oscillation (ENSO) events. This variability results in ecosystem-wide food web changes which are reflected in the tissues of the Peruvian fur seal (Arctocephalus australis). Stable isotope ratios (δ13C and δ15N) in Peruvian fur seal vibrissae (whiskers) are used to infer temporal primary production and dietary variations in individuals. Sea surface temperature anomaly (SSTA) recordings from the Niño 1+2 Index region captured corresponding ENSO conditions. Fluctuations in δ15N values were correlated to SSTA records, indicating that ENSO conditions likely impact the diet of these apex predators over time. Anomalous warm phase temperatures corresponded to decreased δ15N values, whereas cold phase anomalous conditions corresponded to increased δ15N values, potentially from upwelled, nutrient-rich water. Vibrissae δ13C values revealed general stability from 2004 to 2012, a moderate decline during 2013 (La Niña conditions) followed by a period of increased values concurrent with the 2014–2016 El Niño event. Both δ13C and δ15N values were inversely correlated to each other during the strongest El Niño Southern Oscillation event on record (2014–2016), possibly indicating a decline in production leading to an increase in food web complexity. Lower δ13C and δ15N values were exhibited in female compared to male fur seal vibrissae. Findings suggest ENSO conditions influence resource availability, possibly eliciting changes in pinniped foraging behavior as well as food web of the endangered Peruvian fur seal.


INTRODUCTION
Large-scale climatic anomalies associated with periodic, alternating El Niño-La Niña-Southern Oscillation (ENSO) conditions are recorded globally through a combination of atmospheric and oceanic teleconnections, resulting in significant, ecosystem-wide impacts (Ropelewski and Halpert, 1987;Trenberth et al., 1998;McPhaden et al., 2006;Sulca et al., 2017). ENSOs are alternating cycles of warm and cold sea surface temperature (SST) in the tropical central and eastern Pacific Ocean; these conditions can last upward of 18 months with reoccurrence approximately every 2-7 years (Gutierrez et al., 2007;Taylor et al., 2008;Grandi et al., 2012). ENSO events are classified by magnitude and are quantified using sea surface temperature anomalies (SSTA; Trenberth, 1997;Oliveira, 2011). Anomalies are deviations from the seasonal average conditions and are anything greater than ± 1 • C. Magnitudes of ENSO range from Weak (± 1 to 1.5 • C) to Extreme (± 2.5 to 3.0 • C; Waluda et al., 2006;Grandi et al., 2012).
The tropical eastern Pacific Ocean is dominated by the Humboldt Current Upwelling Ecosystem, a cold, nutrient-rich coastal current which supports a multitude of ecologically and economically important species (Barber and Chavez, 1983;Ñiquen and Bouchon, 2004;Heileman et al., 2008;Montecino and Lange, 2009;Kluger et al., 2019). Peruvian anchoveta (Engraulis ringens), a keystone species of this ecosystem, comprises the largest single-species fishery in the world, although the species experienced severe stock collapses during the 1971/72, 1982/83, 1997 (Clark, 1976;Alheit and Ñiquen, 2004;Ñiquen and Bouchon, 2004;Espinoza-Morriberón et al., 2017). South American pinniped mortality events and anchoveta abundance declines coincided with these ENSO periods (Trillmich and Ono, 1991;Arias-Schreiber and Rivas, 1998;Oliveira, 2011;Cárdenas-Alayza, 2012). Other apex predators residing in the eastern Pacific Ocean at the time of ENSO periods have been recorded exhibiting modified behaviors, body mass declines, and unusual mortality events. During the 1982/83 El Niño event, the Galapagos fur seal (Arctocephalus galapagoensis), a top predator found in the direct impact zone of equatorial ENSO conditions, experienced large mortality events and lowered pup production in the following breeding season, likely due to resource limitation (i.e., starvation; Trillmich and Limberger, 1985). Similarly, the 1997/98 El Niño affected behaviors in northern elephant seals (Mirounga angustirostris) off the coast of southern California, whose foraging trip duration and distance traveled increased to compensate for lower foraging success (Crocker et al., 2006). During this same ENSO event (1997/98), South American sea lions (Otaria byronia) off the coast of Peru were recorded foraging at deeper depths evidenced through dietary changes. This apex predators' diet altered from predominantly anchoveta and squat lobster (Pleuroncodes monodon) to mostly demersal species; interestingly, the only time pelagic squat lobster was absent in this otariid's diet was during the recorded El Niño event (Soto et al., 2006).
The complex marine food web along the Peruvian coastal margin is influenced by the aforementioned periodic bio-physical changes (Chiu-Werner et al., 2019). The Humboldt Current is inhabited by two fur seal independent evolutionary units of unsettled taxonomic status, the South American fur seal (Arctocephalus australis) and the Peruvian fur seal (A. australis unnamed sp.; Vásquez, 1995;Zavalaga et al., 1998;Arias-Schreiber, 2000, 2003Cárdenas-Alayza and Oliveira, 2016). Geographical isolation between South American fur seal breeding colonies has led to the genetic variation between these two fur seal taxa (Berta and Churchill, 2012;Nyakatura and Bininda-Emonds, 2012;Oliveira and Brownell, 2014). The Peruvian fur seal, found along the eastern Pacific coastline of South America (Berta and Churchill, 2012;Oliveira and Brownell, 2014;Cárdenas-Alayza and Oliveira, 2016), feed most commonly on pelagic and demersal-pelagic fishes and cephalopods; Peruvian anchoveta and Teuthidae squids comprise 40-60% of their diet (Vásquez, 1995;Zavalaga et al., 1998;Arias-Schreiber, 2000, 2003Oliveira et al., 2008). The fur seals forage in proximity to their breeding colonies (150 km,  and the majority of dives are concentrated at 11-30 m in depth (Gentry and Kooyman, 1986). These seals are subjected to El Niño events in the tropical Pacific, which seem to alter their foraging habits, evidenced in both the 1982/83 and 1997/98 El Niños. During these times adult females foraged at depth averaging 29 to 170 m with trip durations lasting upward of 10 to 20 days searching for prey. These findings suggest pup survival rate was negatively impacted (Gentry and Kooyman, 1986;Majluf, 1987a,b;Cárdenas-Alayza, 2018).
Peruvian pinnipeds are distinctly vulnerable during these strong magnitude ENSO episodes (Oliveira, 2011). These marine mammals serve as sentinel species, indicative of ecosystem health in an environment at the center of fluctuating ENSO conditions (Fossi and Panti, 2017). The feeding ecology and movements of these pinnipeds are of substantial interest as ENSO events continue to intensify in both frequency and magnitude (Sepúlveda et al., 2014). Through a combined analysis of stable isotope profiles and abiotic recordings from El Niño indices, the feeding ecology of the Peruvian fur seal was investigated in Punta San Juan, Peru (Figure 1) over a time series of various ENSO conditions. Stable isotope ratios were used to depict the trophic fluctuations of fur seal individuals and their population. Abiotic recordings from the Niño 1+2 index are reflective of environmental changes due to alternating ENSO conditions along the coast of Peru. By using continuously growing, metabolically inert vibrissae (whiskers) from Peruvian fur seals, it is possible to obtain valuable, multi-year dietary information representative of their foraging environment across fluctuating conditions during overlapping, multiple ENSO events (Hirons, 2001;Cherel et al., 2009;Rea et al., 2015;Edwards, 2018). This is the first study to capture potential ecosystem changes recorded by individual Peruvian fur seals through the assessment of trophic level via vibrissae analyses and relate these changes to multiple ENSO events off the Pacific coast of South America.

Study Area and Sample Collection
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)(2011)(2012)(2015)(2016). Animals were restrained using inhalant anesthetics and/or remotely delivered, multimodal anesthetic drugs as previously described elsewhere (Adkesson et al., 2012.  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

Stable Isotope Analyses
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 δ 13 C and δ 15 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 (Lglutamic 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 δ 13 C and δ 15 N. Stable isotope values were expressed in terms of δ and were reported relative to the standard reference material, Vienna Pee Dee Belemnite standard for δ 13 C and atmospheric air (N 2 ) for δ 15 N. Following Bond and Hobson (2012), stable isotope values were reported with the standard parts per thousand notation ( ): where j X is the heavier isotope, such as 13 C, and i X the lighter isotope, such as 12 C, in the analytical sample (numerator) and the international measurement standard (denominator).

Statistical Analyses
The most recent vibrissa growth is located at the base of the whisker, and an individual vibrissa can represent several years' growth 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 δ 13 C and δ 15 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 δ 13 C and δ 15 N, δ 13 C and SSTA, and δ 15 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, nonparametric 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 δ 13 C and δ 15 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 δ 13 C and δ 15 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 nonsignificant 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).

RESULTS
Similarity in isotopic oscillating patterns occurred across all vibrissae. Both male and female fur seal vibrissae exhibited nearly identical mean stable isotope patterns, but male fur seals exhibited higher values than females in both stable isotopes by up to 0.5 (δ 13 C) and 1.5 (δ 15 N; Figure 2). Overall, δ 13 C values ranged from −18.13 to −13.17 with a mean of  Two linear mixed effect models were employed to explain changes in Peruvian fur seal δ 13 C and δ 15 N values over time with changing environmental conditions (SSTA): Exploratory plots were accessed for both Models 1 and 2; when year is factored, patterns are visible in relation to SSTA and stable isotope values. Model 1 revealed similar trends in δ 13 C values throughout the years with a decline around 2013 linked to colder sea surface waters and a change to rapidly warming conditions in the following years (2014-2016) when the trend in δ 13 C values began to increase (Figure 3). Model 2 revealed dynamic patterns from year to year in comparison to Model 1 (δ 13 C); however, δ 15 N values began a declining trend seen as early as 2013 with increasing SSTA over the following years (2013-2016; Figure 4). Linear Mixed Models are presented in Table 2 where the estimates for the fixed factors are reported for each  variable. The interaction term of Sex:SSTA was not significant for the δ 13 C; and sex alone was not sufficient to explain the variation in δ 15 N values. The best fitting model (lowest AIC score) included the sum of Sex and SSTA as explanatory variables of δ 13 C and δ 15 N values (p < 0.01). Prediction plots for both Models 1 and 2 ( Figure 5) illustrated sex differences; males generally had more enriched 13 C and 15 N than females. Best fit values seen in Model 1 showed a slight increase in 13 C with increasing SSTA whereas best fit values from Model 2 showed a decrease in 15 N with increasing SSTA.

Capturing El Niño Southern Oscillation Signals
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 δ 13 C and δ 15 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 δ 13 C and δ 15 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)(2005)(2006)(2007)(2008)(2009)(2010)(2011)(2012)(2013). The significant inverse correlation between δ 15 N and δ 13 C values during the 2014-2016 extreme El Niño event demonstrates sharply declining δ 15 N while δ 13 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.

Stable Carbon Isotope (δ 13 C)
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 δ 13 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 δ 13 C values. However, the mixed δ 13 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 , revealed more depleted 13 C values further offshore (Sydeman et al., 1997;Hill et al., 2006;Miller et al., 2008). Increased distance from shore corresponded to increased faunal 13 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;. ENSO conditions could also contribute to fluctuating patterns in primary production during warm and cold phases, subsequently compromising resource availability within an environment.

Stable Nitrogen Isotope (δ 15 N)
The δ 15 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 δ 15 N values through time, a long-term gradual decline followed by a quick increase. The δ 15 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 5year decline ( −4.35 ) followed by a 9-month ( +1.56 ) increase. Males showed consistently higher δ 15 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, δ 15 N continually decreased.  Warmer SSTs correlated to lower δ 15 N values. When evaluating ENSO conditions by phase (i.e., cold, normal, and warm), significant correlations between warmer phases and lower δ 15 N values as well as colder phases and higher δ 15 N values were identified. During warmer phases, the lowest δ 15 N values were detected, while consistently higher δ 15 N values were detected during cold phases. The increasing magnitude of ENSO conditions revealed significant correlation to lower δ 15 N values; stronger or more abnormally warm conditions revealed below average δ 15 N values, whereas norm, weak, and moderate magnitudes revealed mean population δ 15 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 δ 15 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).

Sex Evaluation
Both δ 13 C and δ 15 N values were significantly different between adult male and female Peruvian fur seal. Female 13 C and 15 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 δ 13 C and δ 15 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).

CONCLUSION
The Punta San Juan rookery located within the Humboldt Current Large Marine Ecosystem is home to the Peruvian fur seal (A. australis) unnamed sp. Environmental conditions within the ecosystem fluctuate frequently due to ENSO conditions, producing ecological consequences for wildlife.
Vibrissae yield biochemical, multi-year timelines which allow for the evaluation of long-term trophic fluctuations within this dynamic ecosystem. Stable isotope ratios from Peruvian fur seals revealed temporal variations, possibly in diet, ecosystem production, and/or foraging habitat, related to anomalous ENSO conditions. The δ 13 C and δ 15 N values acted as a proxy for sex-related resource and habitat use.
Variability in SSTA was correlated to fluctuations in both δ 13 C and δ 15 N values, providing evidence that ENSO conditions alter the foraging of these apex predators, potentially due to changes in surrounding resources over time. Although behavioral adaptations could not be evaluated solely by these pinnipeds' stable isotope ratios, fluctuations in both stable isotopes across years revealed that the fur seals were still feeding during these stressful conditions. Sex-related stable isotope ratio differences observed during ENSO conditions illustrated that sex-specific adaptations may exist. Further long term, trophic-based studies will be necessary to demonstrate the existence of true adaptations.
The anomalous climatic conditions along the western South American coast may lead Peruvian fur seals to utilize alternative foraging strategies to survive; these data corroborate Soto et al. (2006) that Peruvian fur seals may be utilized as a biological indicator of relative changes in marine resources. With the mean trend of SSTA in the 1+2 Nino region continuously increasing over 1 • C during the past decade, culminating in the strongest ENSO event on record (2014-2016), both increased ENSO duration and magnitude threaten this vulnerable species and the coastal marine ecosystem.

DATA AVAILABILITY STATEMENT
The original contributions presented in the study are included in the article/Supplementary Material, further inquiries can be directed to the corresponding author/s.

ETHICS STATEMENT
Ethical review and approval was not required for the animal study because sample collection was conducted by MA during a previous project -CZS IACUC approved.

AUTHOR CONTRIBUTIONS
AH and ME acquired funding for the analyses. MA and SC-A previously collected the samples. ME and MD-A processed all samples. ME, SC-A, and AH performed all data analyses. All authors contributed equally to the project development and manuscript preparation.

ACKNOWLEDGMENTS
These data contributed to the completion of Master of Science thesis by ME Nova Southeastern University's President's Faculty Research and Development Grant (PFRDG) provided funding to AH. The Southern Florida Chapter of the Explorers Club provided field work funding to ME. J. Coley and L. Alfino provided laboratory assistance. R. Milligan provided analytical assistance. C. France from the Smithsonian Institution's Museum Support Center provided stable isotope analyses. Chicago Zoological Society and the Punta San Juan Consortium funded the field research and all associated field operations. We thank G. Jankowski, J. Meegan, M. Allender, and the veterinary and field biologist staff that assisted with sample collection. We further acknowledge SERNANP for access to the RNSIIPG-PSJ reserve and AGRORURAL for use of field facilities. All samples were collected in accordance with the SERNANP research permits Resolución Jefatural Nos. 009-2010, 023-2011, 022-2012, and 008-2015