LA-ICP-MS generated signals for 15 elements (Supplementary Table S1) for each tusk piece across all 16 tusks. Here, we focused on Sr and Ba, since these elements are assumed to reflect seasonal migrations between winter and summer grounds through cyclic variations (Ruiz-Sagalés et al., 2024; Heimbrand et al., 2020), thereby providing proxy markers for age estimation. For the purpose of age estimation, trends add unnecessary complexity, as our model only needs to capture the dynamics of the repeating cycles. We therefore decomposed each signal into a trend component and a cyclic component, with focus on the latter. The trend components were estimated using local polynomial regression [locally estimated scatterplot smoothing (LOESS)] implemented in the statistical program R. This estimation was performed on the full tusk signal, by combining the signals from each tusk piece in chronological order. The cyclic components were then obtained by subtracting the estimated trend from the original signal (Figure 3). We chose the Spearman correlation coefficient to account for non-linear associations. Our proposed model assumes that the observed signals consist of a mean structure, expressed as a sum of sinusoidal functions, plus measurement error (Equation 1). The measurement error (noise) is explicitly modeled as normally distributed. Further details are provided in section 2.3. To reduce the variance of the noise term in the signals used for age estimation, we considered the average of Ba and Sr. We also assumed that all growth layers were present in the tusk piece, such that when pieces are gathered, the signal covers the entire lifespan of the narwhals (except for possible few missing layers; see section 4).

In this section, we introduce the time-warping model fitted to the average Ba–Sr signals, which display cyclic patterns reflecting narwhal migratory behavior. These cycles alternate between larger and smaller peaks due to seasonal differences. We account for this by modeling the average dynamics as a sum of sinusoidal components, or “waves”, where the dominant wave models the overall annual dynamics and the second wave corrects for any shift in concentration between seasons. Because tusk growth rate vary over time, e.g., due to poor feeding year or illness, the rates are modeled stochastically, leading to a more realistic representation of the signal dynamics.

Let y_i = y(x_i) be the detrended (see section 2.2) element concentration of the average signal over Ba and Sr at a distance x_i = Δi from the cementum–dentin junction along the laser transect of the tusk piece, where Δ = x_i − x_i−1 is the distance between measurements for i = 1,…,n, with n equal to the number of measurements. The measurement at position x_i then corresponds to the element deposition at an unknown time through an increasing function g, t_i = g(x_i), which we denote the growth process, since the increase in g from position x_i to position x_i+1 gives the time passed for the tusk to grow a distance x_i+1 − x_i along the laser transect. We impose that g is a strictly increasing function with respect to x, which ensures that time never goes backward.

The growth process is defined as the integral of an instantaneous frequency process ξ_x [frequency at time t_x = g(x)], which we assume is stochastic due to environmental, seasonal, and internal variability in the growth rate of the tusk. The process ξ_x should be positive, since the tusk can only grow and not shrink. We assume that it varies stochastically around some mean growth rate α > 0. We use the so-called square root (SQR) diffusion, which is a positive stochastic process with governing equation.

with parameters mean α, rate of adjustment β, and variance ω2/2β, and Wt is a Wiener process modeling the unknown deviations from the mean as stochastic.

The growth process is then defined as g(x)=∫0xξs ds, so g is also stochastic, and strictly increasing, since ξx>0, ensuring non-reversion of time. Whenever ξx increases, it indicates a deceleration in the growth of the tusk (since the resulting increase in g indicates that more time has passed for the same amount of dentinal growth).

We assume the following form for the progression of the element measurements:

where ϵi is assumed to be an independent and normally distributed measurement error with mean 0 and variance σ2,, which could describe instrumental fluctuations, surface contamination, or small-scale variations not captured by the model mean. The parameter A represents the baseline amplitude of the element concentration at the time of year of maximum concentration (summer or winter), while B<A modulates this amplitude, accounting for seasonal variations. The offset parameter b determines the point in the season of the birth of the animal.

The growth process g(x) provides the current age in years at any spatial position xi, through counting the elapsed number of cycles, such that:

where g(x) is the growth process of the entire tusk, and xn is the total distance covered across all tusk pieces. The model was fitted to each individual piece of the tusk, whereafter the growth processes were concatenated in chronological order, i.e., from the tip of the tusk to the root segment, taking into account possible overlaps (see section 4.3.2).

We then perform time warping, by expressing the observations y_i in terms of time coordinates t_i instead of length coordinates xi.

The model was fitted using a variant of the Stochastic Approximation Expectation Maximization (SAEM) algorithm (Ditlevsen and Samson, 2014; Reiter et al., 2025). Residual bootstrapping was used to produce confidence intervals for the age estimate. It was performed on the fitted signal for each tusk piece to generate M = 100 bootstrapped versions of the growth process g(x). A single bootstrapped version was randomly selected for each tusk piece, and the selected processes were then concatenated to form a complete growth process for the entire tusk. This procedure was repeated 100,000 times, generating 100,000 complete growth processes for each tusk, hence 100,000 age estimates using Equation 3. The empirical density of the bootstrapped estimates, with a small correction introduced in section 2.3.1, is then used for uncertainty quantification. See Reiter et al. (2025) and section 4.3 for details and further discussion of the model and the estimation method.

Analysis of the 15 elements from LA-ICP-MS revealed that strontium (Sr) and barium (Ba) exhibited the most distinct cyclic variations (Figure 3). In addition, we found that Sr and Ba, among all elements, were most similar in terms of both their cyclic variations and their trend over time. For these two elements, the Spearman correlation test across pooled tusks showed a highly significant relationship (p < 10⁻⁴; Bonferroni-adjusted p-value for double testing), with an average (over all tusks) correlation pair of 0.49 and 0.57 for the trend and cycle component, respectively, higher among any other elements (Table 2). In addition to these cyclic variations, visible in elements such as Sr, lead (Pb), Ba, and lithium (Li), we observed that incorporation varied between groups of elements. For some elements, we observed a gradual decline in concentrations from early GLGs (near the tip; juvenile phase) to late GLGs (near the root; prime phase). Other elements showed more irregular patterns. For instance, Ba tended to be relatively high during the juvenile stage of a narwhal’s life, possibly due to dietary intake related to nursing (Clark et al., 2020a). Figure 3 also illustrates that the standardized Ba and Sr signals are almost interchangeable. While Ba and Sr may seem nearly identical, the signal-to-noise ratio can differ significantly within a tusk piece. We often noticed that the signal-to-noise ratio of Ba appeared inversely related to that of Sr. We refer to the Supplementary Materials for a 2-dimensional visualization of the (trend,cycle)–element correlation pairs.

The signal average of barium and strontium correlated with seasonal cycles, enabling the translation of spatial locations on tusk cross-sections into temporal events and providing age estimates through a novel time-warping model. The model captures both the deterministic seasonal component and stochastic deviations. In the example piece shown, the fitted signal closely follows the observed Ba–Sr oscillations (Figure 5A). Periods of higher (lower)-frequency fluctuations (Figure 5B) coincide with periods where the growth process g(x) is steep (flat) (Figure 5C), corresponding to local compression (expansion) of GLGs, which, in turn, means less (more) growth. For this tusk piece, when we integrate up all the frequency contributions, we see that growth process actually steepens, corresponding to declining growth. This is seen as increasingly contracted annual intervals, which, on the time-warped scale, were adjusted such that all annual intervals achieved equal length (Figure 5D). Table 1 provides complete information for all tusks, including our model’s age estimates and estimates from visual GLG counts, C14 dating, and AAR, as well as overlap detection assessments. Tusk 938 was broken, but we were not able to provide an informed interval for the missing Layers (Garde and Heide-Jørgensen, 2022); thus measurements (except for AAR) are to be considered minimum estimates. Overall, visual GLG readings and model estimates were consistent, except for tusks 956 and 974, where discrepancies persisted even when accounting for model uncertainty. AAR measurements for tusks 938 and 1015 also diverged from both model ages and visual readings, whereas the other four AAR values showed no substantial differences. Among the tusks, tusk 956, which records one model estimate, two visual readings, and a single C14 date, exhibited the greatest variability in age estimates across methods. In contrast, the C14 dates for tusks 953 and 4076 closely aligned with the other methods. Concerning overlaps, the visual readings and the least-squares estimation generally agreed, with some clear deviations for tusks 1015, 929 and 938 where the visual readings detected roughly twice as many overlaps. For several tusks, we saw sharp and narrow densities indicating high model confidence in age estimates (Figure 6). For tusks 4076, 925, 938, 929, and 841, we saw broader, flatter distributions reflecting greater uncertainty. The tusks were arranged from top to bottom by estimated mean age (Table 1, Figure 6), and illustrations were added to simultaneously visualize how age and length co-varied (Figure 6). The age–length relationship revealed a general positive trend, with older narwhals typically possessing longer tusks, though considerable variability was observed. Some older individuals exhibited relatively short tusks, while certain younger whales displayed disproportionately large ones, indicating that factors beyond age influence tusk size (Graham et al., 2020).

The traditional method for estimating the age of marine mammals is through counting of growth layers in the dentin or cementum of the teeth (Read et al., 2018; Scheffer and Myrick, 1980). However, visual reading is often difficult due to poor contrast between layers and the process of GLG counting is somewhat subjective with considerable variability in counts both within the same reader and between different readers, especially when assessing older individuals (Bowen et al., 1983; Lonati et al., 2019). Age estimation also relies on the assumption that one GLG is deposited each year (Stewart et al., 2006)—an assumption that, while generally accepted, has stirred some debate (Brodie et al., 2013). Finally, factors like attrition of the erupted tusk (seen mainly in older individuals) can wear away the neonatal line that marks the onset of the GLG deposition—this was seen in tusks 4076 and 956 (see Table 1). For the embedded tooth, dentinal occlusion can obscure the visibility of GLGs (Watt et al., 2020; Read et al., 2018). This last point is why GLG counting is in many cases considered a minimum age estimate (Read et al., 2018). Despite these limitations, visual readings remain popular because it provides a cost-effective method that can be applied to large samples without demanding expensive laboratory analyses, unlike techniques based on AAR, C14 (Stewart et al., 2006; Garde et al., 2012), or even LA-ICP-MS. On the other hand, specialized analyses eliminate inter-observer variability and subjective interpretation inherent in visual readings, while providing age estimates where GLG reading fails due to poor preservation or tooth wear. The higher accuracy and reproducibility of these biochemical methods justifies their cost when reliable age estimation is essential.

Our proposed method for age estimation needs only a rough identification of GLGs, in order to properly segment the tusk and identify pieces covering new sets of GLGs, used for the subsequent LA-ICP-MS analysis. Still, this segmentation procedure risks either counting certain GLGs several times across consequent tusk pieces or failing to identify intermediate layers (although the latter is considered unlikely). If not corrected for, this can then skew the age estimate towards higher values. In this study, we have attempted to address this last point by estimating the number of overlapping GLGs, under the assumption that the collection of tusk pieces cover all GLGs with possible repetitions between adjacent pieces. We have done this using a least-squares approach between time-warped signals of consecutive tusk pieces. This admittedly involves unaccounted uncertainty that increases with each new comparison. For the purpose of age estimation, an alternative strategy would be to incorporate the uncertainty relating to overlaps in the bootstrap procedure, similar to what we did with missing layers. We chose the least-squares alignment approach because it provides a logical, objective means of merging tusk piece signals into a single, continuous timeline that has use-cases beyond age estimation. The problem of repeated GLGs across pieces is specific to this dataset of narwhal tusks, where tusks need to be sectioned to fit the LA-ICP-MS equipment, a constraint not typically encountered with smaller, more compact specimens. To overcome this issue for future large specimens, segmentation should ideally produce spatially contiguous (gap-free) pieces capable of fitting the LA-ICP-MS equipment.

Concerning the analyzed LA-ICP-MS signals, we found that this collection was quite diverse. Some signals showed clear cyclicity whereas others displayed elevated noise levels slightly obscuring such repeated patterns. This directly influenced estimate precision. Tusks with clear, repeating cycles (Figure 3A) conforming to the model specified by Equation 1 yielded narrow uncertainty bands. Signals with irregular patterns, including secondary cycles (see the next paragraph), produced bootstrap replicates with greater variability and correspondingly wider confidence intervals (Figure 6; Supplementary Figure S3). This is because we resample residuals in the bootstrap procedure. If residuals are skewed or not conforming to the normal assumption (for whatever reason), the resulting bootstrapped signals can create new spurious cycles or eliminate existing ones.

Another subtlety is the existence of sub-annual layers (or accessory layers), likely associated with life-history events such as periods of low food consumption (Read et al., 2018; Zhao et al., 2025). Accessory layers are sporadically embedded within the primary growth layers and risk being misidentified as GLGs when doing visual readings, which is likely an important cause for the high variability seen in GLG counting. It is hypothesized that such sub-annual layers could potentially manifest as secondary cycles in the element signals used in this study, thereby adding further variability to the final age estimate. These secondary cycles may appear as embedded peaks within annual peaks (Figure 5; note the additional peak in the third elapsed cycle, which may be due to accessory layers). Small embedded peaks were seen in several barium and strontium signals, but with a size negligible to the major peaks associated with the primary growth layers. These embedded peaks also carry over to the cyclic part of our signal, when we perform the decomposition described in section 2.2. For tusk 938, base occlusion potentially obscured several layers. This, together with the missing layers not accounted for, could explain the gap between AAR and the other methods. Overlap detection differences between methods contributed to some age discrepancies. Algorithmic overlap estimates closely matched visual assessments for tusks 4084, 951, 888, and 4075 (difference ≤ year). For tusk 1015, the difference between model age estimate and closest visual reading was 4 years when correcting for overlap, reducing to 1 year when overlaps were included. We consider these issues to be ubiquitous in relation to GLG-based age estimation, underscoring the intrinsic challenges posed by biological complexity and environmental influences rather than methodological shortcomings.

We obtain age estimates by fitting the time-warping model introduced in the preprint (Reiter et al., 2025) to cyclic element signals. The model is specified in Equation 1 and assumes that the signals contain two distinct annual peaks of varying magnitude due to seasonal differences. For the narwhal tusks, this pattern is believed to be due to migration between summer and winter grounds, where the composition of prey items and, hence, mineral intake likely differs. The assumption may not hold true for every single annual cycle, which may introduce error in the age estimate. For example, if two consecutive peaks were winter peaks, but were modeled as a winter and a summer peak due to level differences, the model would interpret growth as proceeding more slowly, and using Equation 2, we would count one elapsed year when in fact it was two. When signals exhibit negligible seasonal variation such that B ≈ 0 in Equation 1, two scenarios may arise. In the first scenario, summer and winter peaks are roughly equal in magnitude, and this will surely lead to overestimation of the number of cycles using our model. In the second scenario, one of the seasonal peaks shows negligible magnitude relative to the other. The determination between these scenarios depends on the nature of the signals and the ecological context. For narwhals, whose feeding ecology suggests one dominant seasonal peak, the condition B ≈ 0 would indicate the second scenario. It is worth noting that, although we observed a pattern of alternating high and low peaks in our studied population, this may not be a common feature across other narwhal populations.

Our extensive dataset of trace elements provides insights that go beyond simple age estimation. The time-warping framework introduced in section 2.3 effectively translates spatial locations on the tusk to an annual timeline, which allows us to compare changes in trace elements within the tusk to environmental conditions in specific years or periods. For instance, uranium or lead concentrations could be examined in relation to records of Arctic pollution. Biological processes can also be inferred from selected trace elements, such as barium or strontium, where abrupt shifts early in life may indicate the timing of weaning (Clark et al., 2020a). Finally, a few elements (magnesium, potassium, and rubidium) lacked a clear seasonal pattern, but instead showed consistent high concentrations in the central dentin across tusk pieces, which likely reflect underlying biological processes.

Beyond narwhals, our model can be applied to other species with teeth or tusks that form growth layers, as well as to other mineralized tissues and other growth structures with annual cyclic patterns in selected elements. Examples include elephant and walrus tusks, beluga or mammoth teeth, and even fish otoliths. This assumes that a tracer can be identified that reliably reflects annual cycles and conforms to the model specified by Equation 1. If there are no seasonal differences, one can remove the second sinusoidal term by setting B = 0 in the model equation.

Overall, this demonstrates the usefulness of our method, not only for age estimation across species, but also for translating spatial measurements of annual deposited element records into an annual timeline, opening the door to further ecological and physiological analyses.