Effective biomonitoring requires an understanding of the factors driving concentrations of the substances or compounds of interest in the tissues of studied organisms. Biomonitoring of trace elements, and heavy metals in particular, has been the focus of much research; however, the complex roles many trace elements play in animal and plant tissues can make it difficult to disentangle environmental signals from physiology. This study examined the concentrations of 15 trace elements in the teeth of 122 Pacific walruses (Odobenus rosmarus divergens) to investigate the potential for walrus teeth as biomonitors of trace elements in Arctic ecosystems. Elemental concentrations were measured across cementum growth layer groups (GLGs), thereby reconstructing a lifetime history of element concentrations for each walrus. The locations of GLGs were used to divide trace element time series into individual years, allowing each GLG to be associated with an animal age and a calendar year. The elements studied exhibited a great deal of complexity, reflecting the numerous factors responsible for generating tooth trace element concentrations. Generalized linear mixed models were used to investigate the importance of age and sex in explaining observed variation in trace element concentrations. Some elements exhibited clear physiological signals (particularly zinc, strontium, barium, and lead), and all elements except arsenic varied by age and/or sex. Pearson correlations revealed that elements were more strongly correlated among calendar years than among individual walruses, and correlations of trace elements within individual walruses were generally inconsistent or weak. Plots of average elemental concentrations through time from 1945 – 2014 further supported the correlation analyses, with many elements exhibiting similar patterns across the ~70 year period. Together, these results indicate the importance of physiology in modulating tooth trace element concentrations in walrus tooth cementum, but suggest that many trace elements reflect a record of environmental exposure and dietary intake/uptake.

Trace element analysis and data processing

Postcanine teeth from 122 Pacific walruses (Female: n = 93; Male: n = 29) were on loan from the University of Alaska Museum in Fairbanks, Alaska, and the National Museum of Natural History, in Washington DC. Specimens were collected between 1880 and 2016 (Table S1). The majority of these samples originated from Alaska Native subsistence harvests in the communities of Gambell and Savoonga on St. Lawrence Island, Alaska, though some of the earlier specimens were collected during scientific expeditions. Because specimens used in this study originated from museum collections and/or Alaska Native subsistence harvests, this research was Institutional Animal Care and Use Committee (IACUC) exempt. All specimens from contemporary subsistence harvests were transferred to UAF for analysis under a Letter of Authorization from the United States Fish and Wildlife Service (USFWS) to Dr. L. Horstmann.

A low speed, water-cooled saw equipped with a diamond blade was used to create a 1.5 mm-thick longitudinal cross-section of the center of the tooth. A 3000-grit diamond smoothing disc mounted on a rotary polishing wheel was then used to polish this cross-section. Samples were rinsed with ultra-pure water after polishing and allowed to air dry, then rinsed and air dried again immediately prior to analysis.

Trace element analyses were conducted at the Advanced Instrumentation Lab, University of Alaska Fairbanks (UAF), Fairbanks, Alaska. An Agilent 7500ce Inductively Coupled Plasma Mass Spectrometer (ICP-MS; fitted with a cs lens stack to improve sensitivity), coupled with a New Wave UP213 laser, was used to measure concentrations of vanadium (⁵¹V), chromium (⁵³Cr), manganese (⁵⁵Mn), iron (⁵⁷Fe), cobalt (⁵⁹Co), nickel (⁶⁰Ni), copper (⁶³Cu), zinc (⁶⁶Zn), arsenic (⁷⁵As), strontium (⁸⁸Sr), molybdenum (⁹⁵Mo), silver (¹⁰⁷Ag), cadmium (¹¹¹Cd), barium (¹³⁷Ba), and lead (²⁰⁸Pb) in walrus tooth cementum. Instrumental precision for the ICP-MS was ± 5 %. The internal standard for these analyses was ⁴³Ca, and the resulting calcium-normalized element concentrations are reported in parts per million (ppm). Measured elemental concentrations were compared with a United States Geological Survey microanalytical phosphate standard (MAPS-4), as well as a National Institute of Standards and Technology Standard Reference Material (NIST SRM 610). All laser transects were ablated using the following parameters: beam width = 25 μm; power = 55 %; pulse frequency = 10 Hz; transect speed = 5 μm/s. Dwell times ranged from 0.002 – 0.15 seconds (Table S2). Locations of ablation transects were selected to maximize distance from the root, where cementum growth layer groups converge and become distorted, while also avoiding areas of tooth wear near the crown, where not all cementum layers are present. Transects were ablated starting at the interface between the dentin and the cementum (first year of life), and ending at the outer edge of the tooth (final year of life). Thus, elemental time series generated during these analyses represented a lifetime record for each animal.

Data extraction and processing was conducted in Igor Pro version 6.37 using the Iolite software package version 3.0. All statistical analyses were conducted using R version 4.0.2 (R Core Team, 2020) with RStudio version 1.3.959 (RStudio Team, 2015). Limits of detection were calculated separately for each analytical run using the standard method applied by Iolite (Longerich et al., 1996). A value of one half the limit of detection was used to replace data points that fell below the detection limits (U.S. Environmental Protection Agency, 2000). Data points more than 4 standard deviations from the mean were considered outliers and removed from analysis (Tukey, 1977). These data points were typically single, unrealistically high values, and were likely to represent instrumental errors, rather than actual changes in tooth trace element concentrations. Their removal is therefore unlikely to have impacted the results of this study.

Growth layer group counts and designation of element concentrations to individual years

After trace element analysis, photographs of walrus teeth were taken using a Leica DFC295 camera coupled with a Leica M165 C optical microscope using reflected light. All growth layer groups (GLGs) in the tooth cementum were identified (Fay, 1982; Garlich-Miller et al., 1993) and marked collaboratively by the authors (C.T.C., L.H., and N.M.), and their positions were revisited on at least two additional days to confirm their locations on the laser ablation transect (Fig. 1). Locations of the growth layers were used to assign measured elemental concentrations to individual years of life, with L1/D1 (the first light and dark layers) representing Age 0 (1^st year of life), L2/D2 making up Age 1 (2^nd year of life), and so on. All GLGs were counted to estimate the age of each animal at death, and this information was used in tandem with the year of death to associate GLGs with individual calendar years. Thus, an animal that was Age 5 when it was harvested in 1995 would have GLGs grown in 1990 – 1995. Only complete GLGs with a fully formed light and dark layer were used for analysis of trace element concentrations by animal age or calendar year.

Statistical methods

Trace element data were natural log-transformed prior to statistical analysis to ensure their distributions approximated normality. Generalized linear mixed models (GLMMs) were run using the R package ‘lme4’ (Bates et al., 2015) to test for relationships between concentrations of each trace element and individual walrus age, as well as test for differences between males and females. These analyses were restricted to ages 0 – 15, to ensure that ≥15 male and female walruses were represented at each age. Model selection was conducted using Akaike’s Information Criterion corrected for small sample sizes (AICc), where models with the lowest AICc score were considered to best explain the variability in the data (Burnham and Anderson, 2002). In instances where more than one model had a DAICc < 2, the model with the least parameters was selected. Prior to running the GLMMs, random effects were selected using restricted maximum likelihood (“REML = ‘TRUE’” in the ‘lmer()’ function) and using AICc selection on the fully-parameterized models with varying random effects. Random effects tested included random intercepts for individual ID (“(1|id)”) and calendar year (“(1|year)”), and a combination of both of these intercepts, as well as correlated (“(age|id)”) and uncorrelated (“(age||id)”) random intercepts and slopes for individual ID by animal age, with and without a random intercept for calendar year. After choosing the random effects, model selection was conducted on five models with varying combinations of fixed effects for individual age and sex (Table S3). Both Sr and Ba exhibit large, non-linear changes in early life associated with nursing and weaning (Clark et al., 2020b), thus GLMMs were only conducted for ages ≥ 5, where the weaning signal is no longer present in the data. Individuals with five or more elemental concentrations classified as outliers (i.e., falling more than 4 standard deviations from the mean concentration of all individuals; Tukey, 1977) were excluded from the GLMM for that element. This resulted in the omission of one individual from the GLMMs for Cu and Pb. Model predictions and 95% confidence intervals were calculated using the ‘bootpredictlme4’ R package, which uses a bootstrapping approach (1000 iterations, in this case) to generate confidence intervals (Duursma, 2017).

Pearson’s correlations were used to investigate relationships among trace elements within the lives of individual walruses, among walruses, and among calendar years. Correlation coefficients were calculated for the time series of 15 trace elements for each individual walrus, and the resulting correlation matrices were averaged for males and females to calculate mean within-individual correlations for each sex. There was high variability within the high-resolution elemental time series, possibly resulting from microscale variations in tooth structure or instrumental noise, which led to almost universally low correlations among elements within individual walruses. To better compare correlations among underlying trends in the data, the elemental time series were smoothed using a Savitzky-Golay filter from the R package prospectr (Stevens and Ramirez-Lopez, 2014) with a window of 15 data points. The Savitzky-Golay algorithm smooths the data by fitting a local polynomial regression of order p (3, in this case) to the data points in the window. Within-individual correlations were calculated using these smoothed time series. To examine correlations among walruses, mean (natural log-transformed) elemental concentrations were calculated for each individual, and correlations among these mean values were computed separately for males and females. Finally, to examine correlations among elemental concentrations by calendar year, mean (natural log-transformed) trace element concentrations were calculated for male and female walruses for each year. Correlations among calendar years were restricted to years in which elemental concentrations from at least three (and usually ≥ 5) individuals were available, which resulted in a time series from 1945 – 2014. Correlation coefficients were then calculated for each element across these years.

Changes in mean (untransformed) elemental concentrations from 1945 – 2014 were examined visually. As for the GLMMs, changes in Sr and Ba through time were calculated using only data from ages ≥ 5. This resulted in slightly smaller sample sizes, but allowed for the inclusion of these two elements in this analysis. Pearson’s correlations between male and female trace element concentrations were calculated for the period from 1945 – 2014 and interpreted alongside the visual examinations.