Male mastodon landscape use changed with maturation (late Pleistocene, North America).

Under harsh Pleistocene climates, migration and other forms of seasonally patterned landscape use were likely critical for reproductive success of mastodons (Mammut americanum) and other megafauna. However, little is known about how their geographic ranges and mobility fluctuated seasonally or changed with sexual maturity. We used a spatially explicit movement model that coupled strontium and oxygen isotopes from two serially sampled intervals (5+ adolescent years and 3+ adult years) in a male mastodon tusk to test for changes in landscape use associated with maturation and reproductive phenology. The mastodon's early adolescent home range was geographically restricted, with no evidence of seasonal preferences. Following inferred separation from the matriarchal herd (starting age 12 y), the adolescent male's mobility increased as landscape use expanded away from his natal home range (likely central Indiana). As an adult, the mastodon's monthly movements increased further. Landscape use also became seasonally structured, with some areas, including northeast Indiana, used only during the inferred mastodon mating season (spring/summer). The mastodon died in this area (>150 km from his core, nonsummer range) after sustaining a craniofacial injury consistent with a fatal blow from a competing male's tusk during a battle over access to mates. Northeast Indiana was likely a preferred mating area for this individual and may have been regionally significant for late Pleistocene mastodons. Similarities between mammutids and elephantids in herd structure, tusk dimorphism, tusk function, and the geographic component of male maturation indicate that these traits were likely inherited from a common ancestor.

Seasonal landscape use and migration are key drivers of animal ecology, providing access to regions used for mating, bearing young, or that are otherwise important for maintaining populations (1, 2). For mastodons (Mammut americanum) and other Pleistocene megafauna with long gestation times (>1 y) and low annual and lifetime reproductive output (3), precise reproductive timing may have aided population maintenance. For example, late spring/early summer has been inferred as the mastodon mating period, allowing a 22-mo gestation (known in elephants and estimated in mammoths) to end in the optimal calving season of early spring (4). In proboscideans, the survival of each calf has heightened importance, because it can take decades for females to become reproductively mature [9 to 20 y old, depending on species and environmental conditions (5, 6); estimates of 9 to 12 y old for M. americanum (4)]. Further, pregnancies generally produce single calves, and births are usually separated by at least 4 y but can be separated by 8 y or more under stressful conditions (6). To maximize calf survival, as observed in many extant arctic and boreal species (1), reproductive strategies of mastodons and other Pleistocene proboscideans may have relied on seasonally linked landscape use, including migrating to areas suitable for mating or giving birth. While some aspects of landscape use are now being explored in fossil contexts (7–11), the seasonal specificity of landscape use is mostly unknown for extinct species.

Strontium (87Sr/86Sr) and oxygen (δ18O) isotopes in mineralized tissues are powerful proxies for reconstructing patterns of landscape use and migration (8, 12–15). Spatial variability in 87Sr/86Sr is largely a function of the composition and age of surface geology and underlying bedrock (16–18). Herbivores incorporate biologically available 87Sr/86Sr from ingested plants and water into their tissues with minimal fractionation, thereby recording this geographic fingerprint (19). Oxygen isotope values of meteoric water also vary spatially, but unlike strontium they show pronounced seasonal fluctuations, which can help determine the time of year when mineralization took place (20). For most species, however, two challenges have limited the use of these isotope systems to reconstruct seasonal landscape use. First, studies typically target molar teeth, which, for most mammals, record at most a few years of life (12–14). Thus, even as patterns of landscape use for some extinct species are becoming clearer (7, 8, 11–15), the limited temporal breadth available from molars hampers our ability to document longer-term patterns of landscape use or to test for seasonal preferences. Second, comparing organismal 87Sr/86Sr and δ18O values with regional isoscapes (geographic maps of isotopic variability) to infer where animals were during tissue development typically results in observed isotopic signatures matching multiple isoscape locations (13, 21). We overcame these challenges by 1) serially sampling sequentially formed layers of tusk dentin (which are added throughout life) from an adult male American mastodon and 2) estimating probabilities of geographic occupancy associated with each serial sample, as well as movement between these areas, using a spatially explicit stochastic state-space model (22). Advancing past work (10), our model incorporates both the distances between candidate locations of landscape use and their distribution across space.

The mastodon analyzed here, the “Buesching mastodon” (M. americanum), was recovered in 1998 from a peat farm near Fort Wayne, IN [Indiana State Museum 71.3.261 (4, 23); ]. Combining the probabilities of a previously published accelerator mass spectrometry (AMS) radiocarbon date (24) with our independent AMS date () indicates that this mastodon lived during the later portion of the Bølling–Allerød warming period (combined calibrated 95% confidence interval: 13,316 to 13,111 calibrated years B.P.). For the American Midwest, this was a time of increasing temperatures and winter precipitation during the height of nonanalog floral communities of the terminal Pleistocene (25–27). Humans were also on the landscape and undergoing population expansion (28, 29).

The Buesching mastodon is a nearly complete skeleton of an adult male. Both tusks are in excellent condition, though the right tusk is longer, providing a more complete record of life history (30). In each year of a mastodon’s life, dentin forms by apposition (within the tusk pulp cavity) as light/dark couplets. The dark phase typically forms in late winter and the light phase in the balance of the year (30–32; ). To access a surface for isotope sampling and permit reconnaissance of tusk growth history, we sectioned the tusk longitudinally and then visually explored exposed dentin couplets, noting an interval of conspicuous variation in growth rate which we address further below. Overall, the Buesching mastodon’s right tusk records 30 complete years of growth, plus one incomplete year at each end of the sequence (referred to below as “tusk years” 1 to 32; ). The tusk tip shows evidence of breakage and abrasion that removed material mineralized during the earliest part of life, prior to tusk year 1. Comparing the geometry of the Buesching mastodon’s right tusk tip to that of a younger male mastodon with less tusk abrasion (), we estimated that 3 y were missing from the earliest tusk record (). Thus, we add 3 to each tusk year to estimate the corresponding age of life. Based on this, the Buesching mastodon died before the end of his 35th year, at an age of 34 y. His age at death is comparable to other bull mastodons from the region and time period (4, 23, 31).

Detailed observations on annual growth increments and behaviorally induced damage on the Buesching mastodon’s right tusk reveal that for at least 8 y prior to death, and always during spring/early summer, there is a structural record of battles like those associated with mating in extant elephants (5, 33) and other mastodons (4, 23, e and f). In elephants, such battles occur mainly between mature males competing for mates during a period of hormonally heightened aggression (33) and increased mobility (34) known as musth. The record of musth-like behavior in the Buesching mastodon appears to have been more season-specific than is typical for extant elephants (5), perhaps due to greater seasonality in Pleistocene environments (32), compared to where elephants live today. His musth-associated battles intensified late in life, leading to tusk fractures and climaxing with a perimortem 4- to 5-cm-diameter puncture within the right temporal fossa (). This injury was likely produced by another mastodon’s tusk (4, 23), suggesting death resulted from male–male conflict during the mating period ().

The preservational quality and relative completeness of our specimen offer opportunities to evaluate changes in behavior and landscape use associated with male maturation, including 1) initial sexual maturity and emergent independence from the natal herd and 2) further changes after having reached sufficient size, strength, and social standing to sustain musth (). Previous work on tusk dimorphism (35) and fossil trackways (4) presented evidence that components of mastodon social structure resembled patterns in extant elephants; young stayed with matriarchal herds, and mature males were more solitary. Similarly reminiscent of elephants, studies of male mastodon tusks have detected a period of nutritional stress and eventual recovery during late adolescence (documented by variation in annual tusk appositional thicknesses and length increases), which has been interpreted as reflecting eviction (or at least departure) from maternal herds (4, 23, 30). The Buesching mastodon’s tusk shows these features of adolescent stress and recovery, offering an opportunity to track changes in landscape use from the earliest phase of maturation (when behavior of the adolescent male also has implications for landscape use by adult females; ). Based on changes in appositional thicknesses, adolescent stress is most extreme during tusk year 9 (his 12th year), which we interpret as the time when the Buesching mastodon separated from his herd (his age of maturation; ). The timing of the Buesching mastodon’s maturation is typical of other male mastodons [generally between 10 and 12 y of age (4); ] and is likewise a good fit to ages of male maturation in extant elephants (5). In adult years of growth, seasonally repeated damage to tusk dentin and cementum is consistent with musth behaviors of extant male proboscideans (), indicating opportunities to evaluate landscape use by a reproductively active male.

To evaluate whether mastodon maturation and later reproductive behavior involved seasonal landscape use, we sampled two periods: 1) adolescence, during and after departure from the natal herd (tusk years 9 to 14/ages 12 to 17 y: 36 samples) and 2) the final years of life (tusk years 29 to 32/ages 32 to 35 y: 30 samples; ). This targeted sampling strategy is designed to allow us to compare patterns of landscape use by our animal as an early-adolescent male with his patterns as a mature male. Our targeted approach to sampling also extends to the level of each successive serial sample. The bit we use for milling samples remains parallel to dentin layers throughout its motion, allowing us to acquire dentin powder from precisely delimited, nonoverlapping zones bounded by dentin laminae (). Serial samples represent between 1 and 2.7 mo of growth, with most under 2 mo (Dataset S1). Extant proboscideans typically complete weaning by 4 to 5 y of age (5), which is consistent with the timing of weaning in late Pleistocene proboscideans (mammoths) (36). Thus, our data are almost certainly free from the isotopic influence of his mother.

To determine season of growth for tusk serial samples (warm season vs. cool season), we evaluated whether the oxygen isotope value of the phosphate component (δ18Op) for each sample was above or below a semiannual mean δ18Op (). To focus on relative (seasonal) differences in δ18Op throughout each year, δ18Op values were z-score-transformed prior to quantitative analyses. To infer landscape use from dentin 87Sr/86Sr, we used an 87Sr/86Sr isoscape that estimates bioavailable 87Sr/86Sr at a 1-km resolution (18). Guided by maximum home range sizes of extant elephants and the limited range of tusk δ18Op, we limited our geographic window to regions less than 250 km from where the mastodon was recovered that were not under glacial ice (37). For each serial sample, we used a Bayesian method for calculating the probability that tusk 87Sr/86Sr could have been sourced from each isoscape pixel (38, 39). We then identified pixels (i.e., geographic locations) with probabilities greater than or equal to 0.50 to focus our movement model on regions with the highest probabilities of occupancy. We employed a kernel density estimator (KDE) to delineate regions of the 87Sr/86Sr isoscape with the highest spatial concentrations of high-probability locations, hereafter referred to as “potential regions of occupancy” (PROs). Depending on the geographic distribution of a particular 87Sr/86Sr value across the isoscape, each serial sample yielded one or more PROs. For each serial sample, we integrated seasonal and spatial data by assigning its δ18Op z-score to its strontium-based PRO(s). We evaluated ontogenetic changes in geographic range and seasonal landscape use by 1) summarizing patterns of spatially overlapping PROs and 2) stochastic movement modeling.

Separately overlaying adolescent and adult PROs (like overlapping transparencies) provides a broad assessment of potential landscape use. To test for recurring patterns and specificity of summer and winter landscape use, we calculated mean δ18Op z-scores for all locations (i.e., each map pixel) overlapped by at least three PROs (). In this framework, regions used predominately in winter or summer have strongly negative or positive z-scores, respectively; regions with z-scores close to 0 were used throughout the year or during winter–summer transitions.

Our stochastic state-space model estimated probabilities of occupancy for each PRO and probabilities of movements between successive serial samples ( and ). We also used model results to test whether the fossil recovery location, which was likely near the mastodon’s final musth battle and in proximity to mastodon mating activities (), was used in specific seasons or throughout the year. We did this by evaluating the seasonality (δ18Op z-scores) of PROs as a function of their proximity to the fossil recovery location. Proximity was calculated as the nearest distance between the recovery location and each serial sample’s PRO with the highest probability of occupancy. We then compared δ18Op z-scores of PROs that were “near” vs. “far” from the recovery location, where “near” was defined as being within one-month’s travel. While the time interval between the final increments of sampled tusk growth and death is unknown, it was unlikely to surpass 1 mo. One-month’s travel was estimated using the median distance moved between successive serial samples, standardized by the number of months represented by each serial sample. Estimates for 1 mo of adolescent and adult travel were 7 and 27 km, respectively ().

Results

Isotopes.

Tusk δ18Op values fluctuate cyclically in both adolescence and adulthood (Fig. 1 ), consistent with expectations for seasonal changes in meteoric water composition (32). Closest to death, tusk δ18Op is at or near its annual high, indicating warm-season mortality. Conversely, patterns of 87Sr/86Sr variability change dramatically between adolescence and adulthood (Fig. 1 ). Adolescence displays three stages of increasing 87Sr/86Sr through time (Fig. 1). Stage i (earliest year sampled) shows low, nearly invariant 87Sr/86Sr, indicating either a constrained or geologically homogeneous home range. This is followed by a jump in mean 87Sr/86Sr (stage ii), sustained for 2.5 y, and a final increase in mean 87Sr/86Sr (stage iii, the final 1.5 y of adolescence). Differences in 87Sr/86Sr among stages i, ii, and iii (all pairwise Wilcoxon comparisons P < 0.01) indicate changes in landscape use. The earlier years (stages i and ii) show no relationship between 87Sr/86Sr and δ18Op z-scores. During late adolescence (stage iii), there is an increase in the annual δ18Op baseline, and a correlation develops between 87Sr/86Sr and δ18Op z-scores during the two warm seasons (Fig. 1, red triangles; r2 = 0.65, slope = 0.000097, P = 0.028).

Fig. 1.

Isotope results for the Buesching mastodon tusk. Data derived from serial samples milled from dentin mineralized during adolescence (A, C, and E) and the last few years of adulthood (B, D, and F). Serial samples were analyzed for (A and B) δ18Op and (C and D) 87Sr/86Sr. SDs for δ18Op and 87Sr/86Sr are ±0.18 and ±7.5 × 10−6, respectively. (A–D) Horizontal axis shows ending positions of serial samples (in millimeters) from the start of each sequence. Serial samples were milled sequentially, from later to earlier in life (following the bold arrows). Alternating gray and white panels demarcate tusk years, numbered along the top of A and B. Cross indicates sample nearest death. (C and E) Adolescence includes three consecutive 87Sr/86Sr stages (i, diamonds; ii, squares; iii, triangles). (C) For each stage, dashed lines indicate means; dotted lines illustrate one SD above and below the means. (D) The adulthood mean and one SD above and below the mean are identified by dashed and dotted lines, respectively. (E and F) Regressions of 87Sr/86Sr on δ18Op z-scores. Tusk serial samples that formed during warm (δ18Op z-score > 0.0) and cool seasons (δ18Op z-score < 0.0) are colored red and blue, respectively; summer (δ18Op z-score ≥ 0.75) and winter (δ18Op z-score ≤ −0.75) are darker shades of red and blue. (F) First year of the adult dataset in gray. Only warm-season relationships between δ18O and 87Sr/86Sr are significant (dashed regression lines); significance occurs only during adolescent stage iii (P = 0.028) and adulthood (last two years highlighted here, P << 0.01). One δ18Op and one 87Sr/86Sr datum in the adult series were excluded due to low sample yield (see ).

In adulthood, 87Sr/86Sr is higher overall and more conspicuously cyclical (Fig. 1). The warm-season correlation between 87Sr/86Sr and δ18Op z-scores is maintained and is particularly strong during the final 2 y of life (Fig. 1 r2 = 0.81, P << 0.01), indicating increased fidelity to summering grounds. In addition, the slope of this relationship (0.00033) is over three times steeper than during stage iii adolescence, suggesting travel across either greater distances or landscapes with more variable 87Sr/86Sr. Cool-season 87Sr/86Sr is not correlated with δ18Op z-scores (Fig. 1 ) and is less variable, particularly during adulthood.

Despite some concern that 87Sr/86Sr from dentin may be susceptible to diagenetic alteration (40, 41), the cyclic variability in our results (for both δ18O and 87Sr/86Sr) and alignment between isotope values and physical features of the tusk (e.g., tusk-year boundaries) both suggest that our measured values are biogenic. Diagenetically altered isotope values sampled along a time series should be largely homogeneous (42). However, δ18O values measured from the Buesching mastodon were highly cyclic, with dark–light (winter–spring) boundaries in dentin consistently corresponding to local minima in δ18O (Fig. 1). Values for 87Sr/86Sr are also not constant. Instead, we find dramatic changes in 87Sr/86Sr along the time series, including annual cyclicity in the adult dataset (Fig. 1). While the early portion of our adolescent dataset (stage i) has relatively stable 87Sr/86Sr (mean = 0.70933, SD = 2.7 × 10−5), these values are distinct from 87Sr/86Sr at the fossil recovery location (87Sr/86Sr = 0.70983) and are unlikely to have originated through diagenetic replacement.

PRO Summaries.

Overlaying PROs for each serial sample, we find that during adolescence the potential geographic range is divided into two main regions (Fig. 2). The first extends across southern and central Michigan and into southern Ontario. The second extends from southern Ohio and central Indiana into eastern Illinois. Small regions of north-central Ohio and northern Indiana are also included. During adulthood, PROs are similarly distributed, but the potential geographic range is expanded, particularly into northern Indiana, where the mastodon was recovered (Fig. 2, cross). In the adolescent record, nearly all mean δ18Op z-scores (∼100% of pixels) are close to zero, indicating no seasonal preferences in landscape use (Fig. 2). Conversely, the adult record reveals landscapes used only during warm seasons (30% of pixels), which include areas extending from southwestern and central Ohio through eastern, central, and northern Indiana and into eastern Illinois. Limited portions of Michigan also show potential for warm-season landscape use (Fig. 2 and ).

Fig. 2.

Seasonal landscape use in potential mastodon range. Recurring landscape use estimated from mean δ18Op z-scores of all overlapping PROs (defined by 87Sr/86Sr) during (A) adolescence and (B) adulthood. Red = summer (δ18Op z-scores ≥ 0.75), blue = winter (δ18Op z-scores ≤ −0.75), gray = no preference (0.75 > δ18Op z-scores > −0.75). Mean δ18Op z-scores are not calculated for pixels overlapped by only one or two serial tusk samples (internal unshaded portions of polygons) due to low sample size and because they do no not indicate areas of recurring use (). Histograms show percentage of pixels identified as summer-use (S), winter-use (W), and no preference (NP). Strong seasonal preferences are only recovered during adulthood (B) and only for summer landscape use (red-shaded areas). Fossil recovery location (cross and red point) is at the margin of adult warm-season landscape use. Dashed boundary (and snowflakes) identify glacial margin (37). Wavy lines identify lakes (37). IL = Illinois, IN = Indiana, KY = Kentucky, MI = Michigan, OH = Ohio, ON = Ontario.

Movement (State-Space) Modeling.

During adolescence, our movement model establishes two possible home ranges for the Buesching mastodon (Fig. 3). Model support is strongest for an initial home range in central Indiana (Fig. 3, darkest shaded regions in the southern range), but there is also some support for Michigan and Ontario (Fig. 3, northern range). For both possible regions of occupation, results for stage i reveal a largely static home range with only limited excursions from a core region. In stages ii and iii, home ranges were maintained in the same general vicinity but expanded through time (Fig. 3). As an adult, annual landscape use encompassed a larger geographic range and was more seasonally patterned. The state-space model identified central Indiana as the core adult range (Fig. 3) but included portions of southwestern Ohio and eastern Illinois. During summers (δ18Op z-scores ≥ 0.75), much of the core range remained the same, but there was a dramatic northward expansion into a summer-only region (Fig. 3, Inset) that includes parts of central and northeastern Indiana near where the mastodon died (Fig. 3, cross). We find no indication that the Buesching mastodon used the area immediately surrounding the fossil recovery site during adolescence (Figs. 3 and 4). This is true even when using the larger adult distance (i.e., larger search radius) to evaluate the adolescent dataset. During adulthood, however, PROs that include areas near the fossil recovery location were occupied repeatedly (tusk years 30, 31, and 32, corresponding to annual highs in 87Sr/86Sr). Each occupation occurred once in a given year and only during summers (at or near annual highs in δ18Op z-scores; Fig. 4, Mann–Whitney U = 72, P = 0.011).

Fig. 3.

Reconstructed mastodon landscape use and movement patterns. (A) PROs illustrating changes in landscape use during adolescent stages i, ii, and iii. Modeling adolescent landscape use reveals a southern and northern set of reconstructions (divided by diagonal line in A), with the southern reconstruction showing greater probabilities (darker shading) of occupancy. (B and C) Adult landscape use during nonsummer periods (B, δ18Op z-scores < 0.75) and summer periods (C, δ18Op z-scores ≥ 0.75) during the final 2 y of life. Regions shown have probabilities of occupancy greater than 1%. Density of shading represents probability of occupancy. (C) Red shading highlights regions used during summer. (C, Inset) Red shading shows areas used only during summers; hatching identifies portions of summer PROs that are also used in other seasons (i.e., areas shown in B). Density of red shading corresponds to probability of reconstructed summer occupancy of that area. (A–C) Dashed boundary (and snowflakes) identify glacial margin (37). Wavy lines identify lakes (37). Fossil recovery location (cross and red point) is on the margin of adult summer landscape use (C). (D) Cumulative distances (calculated separately from the start of each time series) moved during adolescence (stages i, ii, and iii) and adulthood (points = mean, gray polygon = 95% confidence interval).

Fig. 4.

Modeled seasonal use of geographic regions near and far from the fossil recovery location during adolescence and adulthood. The δ18Op z-score for each serial sample is attributed to either “near” or “far” from the fossil recovery location by first identifying the geographic region with the highest probability of occupancy and then evaluating whether any portion of that region was within one month’s travel distance from the fossil location (7 and 27 km for adolescence and adulthood, respectively). (A) During adolescence, the Buesching mastodon did not use habitats near the fossil location. This result is the same when the adolescent dataset is evaluated using the larger adult travel distance (27 km). (B) During adulthood, PROs that included areas near the fossil recovery location were repeatedly occupied (tusk years 30, 31, and 32), but only during summers (red shading). Boxplots display the median (horizontal line), interquartile range (box), and 1.5 times the interquartile range (whiskers).

Focusing on distances traveled between serial samples, the adolescent record begins with a period of negligible movement (Fig. 3, stage i) followed by larger movements during stages ii and iii, consistent with striking out following separation from the natal herd. As an adult, the mastodon showed more continuous movement (Fig. 3), with larger distances traveled between serial samples (Mann–Whitney U = 248, P << 0.01; ).

Discussion

Ontogenetic changes in seasonality of landscape use and mobility are consistent with hypothesized mastodon social structure and behavioral changes associated with male sexual maturation. Relatively homogeneous 87Sr/86Sr and small modeled movement distances during earliest adolescence are not surprising before the young male fully separates from the herd. This is then followed by large shifts in median 87Sr/86Sr, 87Sr/86Sr variability, and associated shifts in modeled landscape use, which we interpret as the young male transitioning to movement patterns more like those of older males. Onset of links between 87Sr/86Sr and δ18Op z-scores late in the adolescent time series likely indicate a new phase of landscape use that was more sensitive to seasonal drivers. Relationships between 87Sr/86Sr and δ18Op z-scores are maintained and intensified in adulthood, when modeled summer landscape use also becomes geographically distinct from landscapes used during the cool season. Interestingly, while seasonality in 87Sr/86Sr is repeatedly observed during the warm season, we do not find ties between 87Sr/86Sr and δ18Op z-scores during cool seasons. This is consistent with behaviors of extant boreal and arctic mammals, which often display structured spring and summer landscape use but less structured or more limited fall and winter movement (43, 44).

As an adult, the Buesching mastodon had highly seasonal landscape use, likely influenced by seasonal mating. Seasonal mating is supported by 1) cementum defects reflecting musth-battles (23) during periods of early warm-season tusk growth () and 2) a spring/summer death (Fig. 1) combined with craniofacial damage consistent with mortal injury inflicted during a musth battle (4; ). While our model cannot differentiate landscape use within PROs, the final northern Indiana resting place of the Buesching mastodon is near areas visited exclusively during warmer months (Figs. 2–4) and is 157 km from the core of his nonsummer geographic range (Fig. 3, star, estimated as the centroid weighted by each pixel’s probability of occupancy). This suggests that 1) the fossil recovery region may have been a preferred mating area for this individual and 2) males could travel long distances to participate in mating behaviors.

The mastodon’s static early adolescent home range, when he was likely still associated with his natal herd (5), did not overlap the northern Indiana mating area. This may indicate that a maturing, but not yet adult, male stayed away from mastodon mating grounds, presumably to minimize conflicts with antagonistic males in musth. It is also possible that during the final years of life the Buesching mastodon was engaging in mating behaviors with a population distinct from his natal herd. Serial sampling of additional male and female mastodon tusks from northern and central Indiana (and surrounding states) will be necessary to test the regional significance of the northern Indiana mating ground.

Changes in mobility and seasonality of landscape use observed in the Buesching mastodon are consistent with patterns observed in African savanna elephants (Loxodonta africana); mature males live in greater isolation, with larger movements and geographic ranges than either adolescent males or female herds (5, 33). Expansion in geographic range with ontogeny was also reported for a late Pleistocene arctic mammoth (Mammuthus primigenius) (10), although our methods break new ground by demonstrating statistically significant shifts in mobility with age. However, alternative explanations for shifts in landscape use are also possible, including coincidental shifts in climate or even changes in human pressure. Large changes to local ecology through bottom-up or top-down mechanisms (including anthropogenic pressures) can impact mammal landscape use and geographic range at population scales (45–47). However, overall stability in the range and amplitude of tusk δ18Op variability across the time series suggests climate is not a likely driver for changes in landscape use. While our data cannot speak directly to possible top-down drivers in ecological change (including humans), it seems unlikely that our sampling would capture such contingencies. The simplest interpretation for the alignment between sexual maturity and shifts in landscape use is that this ecological change was driven by ontogeny.

Monthly movement distances of the Buesching mastodon (median distance moved during adulthood = 27.3 km/mo, first and third quartiles = 6.4, 45.0 km/mo; ) are not unprecedented for living proboscideans. The movement patterns and home range sizes of elephants vary strongly across different environmental settings. While many elephant populations show seasonal or annual movements that are smaller than those of the Buesching mastodon (34, 48, 49), some populations, notably male elephants in the Kunene region of Namibia, routinely move several hundred kilometers in a month or two and may travel over 600 km in 6 mo (50). Thus, the high end of movement velocities of living proboscideans is in line with or beyond what we estimate for the Buesching mastodon, suggesting the magnitudes of our proposed shifts in seasonal landscape use are biologically reasonable.

Our results for the Buesching mastodon also provide context and confidence in previous reconstructions of Pleistocene proboscidean mobility. Past studies of both mastodons and mammoths (including Mammuthus columbi) have used offset between tooth 87Sr/86Sr and bioavailable 87Sr/86Sr at fossil recovery sites to infer travel of several hundred kilometers from regions where teeth mineralized (8, 9, 15). Modeled estimates from a serially sampled (87Sr/86Sr) mammoth tusk also indicate lifetime range shifts of at least hundreds of kilometers (10). Our highly spatially and temporally resolved results are consistent with these past inferences and indicate that adult male mastodons could travel hundreds of kilometers per year.

There is a clear disjunction between where the Buesching mastodon died and where he spent the bulk of his (sampled) life. This finding highlights the complexities of interpreting the geographic distribution of body fossils for highly mobile species, particularly species whose mortality risk varies across space and among seasons. In such cases, there is potential for the resulting fossil record to be spatially and/or seasonally biased (7, 10). The social structures and sex-specific behaviors of proboscideans may make males particularly susceptible to these biases. Across Siberia, for example, the sex ratios of mammoth fossils are strongly male-biased (51). This male overrepresentation may result from their greater susceptibility to falling victim to miring and other environmental dangers. Female mammoths may have been less susceptible to such mortality risks because they lived in more highly social and cooperative herds; as in extant elephants, knowledge and experience of older matriarchs may have imparted numerous benefits for herd survival (52). For North American mastodons, there are multiple male specimens with perimortem damage consistent with wounds inflicted during musth battles (4). If patterns of seasonal landscape use observed in the Buesching mastodon are characteristic for mature male mastodons, fossil recovery locations for individuals dying during musth may be geographically biased away from core annual ranges. Understanding the strength of the geographic and/or seasonal biases in fossil records of proboscideans and other large, mobile mammals will require a broad effort to reconstruct patterns of seasonal landscape use from continuously growing tissues.

Inconsistencies documented here and elsewhere (8–10) between core geographic ranges of proboscideans and their final location of burial also have ramifications for the use of teeth and tusks in reconstructing past environmental settings and climates. Depending on the age, sex, and season of death of the animal, the environmental proxy data produced from a specimen may or may not be closely associated with the location where it was found. While additional focus on female proboscideans is warranted, isotope geochemistry from male proboscidean materials may be especially complex as a source of local paleoenvironmental insight. Uncertainty in the relevance of environmental data may be particularly high for molars or distal portions of tusks, for which a pronounced disparity between time/place of mineralization and time/place of death is likely.

Reconstructing Pleistocene proboscidean landscape use has implications for highlighting gaps in our understanding of extant elephant ecology. Where humans and elephants coexist, conflicts are pervasive (53). While wildlife preserves and parks can offer some protections between humans and elephants, these buffers are generally too small, and elephants’ ranges too large, to effectively avoid all conflicts. More fundamentally, our view of elephant spatial ecology is heavily influenced by human interactions, the sizes and environmental settings of wildlife preserves, and habitat loss surrounding those preserves (48, 53, 54). In the absence of quantitative historical data on elephant ecology, tusks from historical collections and archeological deposits, analyzed in the manner demonstrated here, offer opportunities to quantify how elephant landscape use and home ranges have changed in response to increasing human footprints over the last few hundred years. Similarly, tusks of Pleistocene mastodons and mammoths in North America and elsewhere provide opportunities to evaluate patterns of proboscidean landscape use prior to human influence and to measure anthropogenic impacts once humans became a dominant ecological force.

Finally, the combination of ontogenetic change in male mastodon landscape use and seasonal mating matches patterns observed in extant elephants (34), despite divergence of these lineages in the late Oligocene [24 to 28 Ma (55)]. This supports previous interpretations that similarities in herd structure (based on fossil trackways) between extant elephants and more closely related extinct elephantid proboscideans from the late Miocene (6 to 8 Ma) represent homologies (56). Previous work has also suggested that multiple features of mammutid and elephantid tusk growth and sexual dimorphism point to homologies in tusk structure, function, and social organization (35). We can now add that the geographic component of male mastodon maturation is also shared with extant elephants (5), increasing the likelihood that these life history similarities were inherited from a common ancestor.

Materials and Methods

Serial Sampling of the Tusk.

Isotope and life history analyses were focused on the right tusk, which is longer and more complete than the left (). Following recovery, the right tusk was longitudinally sectioned using a custom-modified, low-speed bandsaw with a half-inch-wide bimetal blade (). Two parallel cuts produced an axial slab (8 to 10 mm thick) containing the tusk axis and dentin representing all preserved years. When completely dry, this slab was embedded in epoxy and cut into segments, from which selected blocks were polished using a graded series of dry, abrasive papers (100, 320, and 600 grit).

Years of tusk growth were distinguished as light–dark couplets that paralleled the conical surface of the pulp cavity (), where dentin apposition occurs. Light-colored dentin typically represents spring through early winter apposition, and dark-colored dentin represents late winter (31, 57). We marked and enumerated tusk years (one winter–spring, W/S, boundary to the next) from the earliest-available year (near the tusk tip) to the pulp cavity. This observational work provided the “scaffold” for a more refined life history of this animal. The methodology followed in this work has been described in prior whole-tusk analyses of male and female M. americanum (4, 23, 30). Full details of the Buesching mastodon’s life history are documented and discussed in .

For the sample series covering the adolescent years, we started sampling in the 14th tusk year and worked in reverse temporal order, toward the time of his separation from his family unit (as estimated by reduction in annual tusk growth; ). The first sample we milled (the last in this temporal sequence) had the widest mill path (1.2 mm; SB-A in ). Subsequent samples in this sequence were acquired by starting the bit in the space opened by removing the previous sample. In this manner, the widths of most serial samples were less than the bit diameter. Most adolescent samples were about 0.8 mm wide in the direction of dentin apposition, affording better temporal resolution than “full-bit” samples would have provided. In sampling the last few years of life (adulthood), we started from the surface of the pulp cavity, adjacent to the last-formed dentin (SB-B in ). By taking advantage of the open pulp cavity (in which the drill bit could freely spin), we were able to control the thickness of dentin from which the first sample was collected. Proceeding in this fashion, we milled 30 samples, with an average mill path width of about 0.5 mm (). In this way, we sampled from early in tusk year 32, well into year 29. Adolescent and adult serial samples were divided into separate aliquots for stable oxygen and strontium isotope analyses. For further details on the specimen and serial sampling, see .

Oxygen Isotope (δ18O) Preparation and Analysis.

We obtained oxygen isotope values (δ18Op) from the phosphate (PO4) component of dentin hydroxylapatite. Following a modified version of the O’Neil et al. procedure (58), an aliquot of ∼40 mg of each serial sample was pretreated with 2 to 3% sodium hypochlorite (NaOCl) to remove organics. Samples were placed in 50-mL centrifuge tubes, sonicated for 2 h, and left to react for 24 h. Samples were then rinsed and centrifuged five times with deionized water. Pretreated samples were dried at 60 °C and then digested with 4 mL of 1 M hydrofluoric acid (HF). Centrifuge tubes were placed in a 60 °C hot water bath to aid digestion for 24 to 36 h. Samples were then neutralized with 3 mL of  30% aqueous ammonium hydroxide (NH4OH). Samples were next filtered through Supor-200 47-mm, 0.2-μm filter paper on a vacuum filtration system to remove CaF2 precipitate and isolate the supernatant PO43− solution. The PO43− solution was then transferred to 50-mL Pyrex beakers, and 0.342 mL of silver amine solution per mg of original sample powder was added to each beaker. This ratio of added silver amine assumes a 40% loss of dentin weight during pretreatment. Five hundred milliliters of silver amine solution consisted of 46 mL of NH4OH, 17.00 g AgNO3, and 14.01 g NH4NO3. The amount of NH4OH, which controls the duration of precipitation, was increased from 25 mL in the original protocol (58) to allow more time for growth of larger silver phosphate (Ag3PO4) crystals. The solution in each beaker was then topped off to 40 mL with deionized water. Beakers were placed in a hot water bath of deionized water held at 60 °C to maintain even temperature under a fume hood. Ag3PO4 crystals started to form after 4 to 4.5 h, when the pH reached a critical level as ammonia evaporated. Experiments using bulk dentin showed that optimal precipitation time (achieved by the 46 mL of NH4OH) was 6.0 to 6.5 h and that more than 7 or less than 5 h resulted in smaller crystal size and lower crystal recovery. Incomplete removal of organics resulted in dark green crystals, but all of the sample series included here formed yellow Ag3PO4 crystals, indicating nearly complete removal of organics.

Ag3PO4 crystals were rinsed with deionized water, vacuum-filtered with Supor-200 47-mm, 0.2-μm-pore filters, dried at 60 °C, weighed, and placed in quartz tubes. A stoichiometric equivalent of graphite was then added to each sample at 0.014316 mg C/mg Ag3PO4. Tubes were degassed for at least 30 min on a vacuum line while heating the samples to 500 °C using a block furnace and then sealed with a propane torch. Sample tubes were baked for 1.75 min at 1,400 °C in a tube furnace, removed for 4 s, and rapidly cooled in a beaker of water. Isotope ratios of the resulting CO2 were measured on a Finnigan MAT Delta S isotope ratio mass spectrometer in the Stable Isotope Laboratory at the University of Michigan. Oxygen isotopes were measured relative to Vienna Standard Mean Ocean Water (VSMOW). Based on repeated measurements of an internal laboratory standard, we estimate the precision of δ18O values as ±0.18 (1 SD).

Only a single sample (BU2-03; the sample third-from-death) was deemed unsuitable and not included in subsequent analyses. This sample had a very low Ag3PO4 weight of only 2.5 mg, compared to an average weight of 33.5 mg for this series. Additionally, its δ18O value was inconsistent with preceding and succeeding samples.

Strontium Isotope (87Sr/86Sr) Preparation and Analysis.

To measure 87Sr/86Sr, an ∼15-mg aliquot of each serial sample was soaked in 30% hydrogen peroxide (H2O2) at room temperature for 72 h and then rinsed five times with ultrapure water. Samples were agitated frequently and the liquid was refreshed after roughly 40 h. Samples were then reacted with 1 M acetic acid (CH3COOH) buffered with calcium acetate [Ca(C2H3O2)2] at 4 °C for 24 h, rinsed five times with ultrapure water, and freeze-dried. Between 3 and 5 mg of each pretreated sample were sent to the Multicollector ICPMS Laboratory at the University of Illinois at Urbana–Champaign for strontium extraction and analysis. Samples were dissolved in 3 M nitric acid (H2NO3) and filtered through Eichrom Sr-specific resin. They were then sequentially rinsed with 0.05, 3, and 8 M H2NO3 to selectively remove cations from the resin to isolate strontium. Strontium was eluted from the resin with 3 mL of ultrapure water and 1 mL of 0.05 M H2NO3 into 4-mL autosampler vials. Samples were analyzed on a Nu Plasma high-resolution multicollector inductively coupled plasma mass spectrometer. Data were corrected for drift using the international standard NBS 987 and two internal laboratory reference standards (Coral and E&A) were used to monitor accuracy and precision. The pooled SD of all measured standards indicates our 87Sr/86Sr precision is ±7.5 × 10−6. The sample second-from-death (BU2-02) did not have enough powder to process 87Sr/86Sr after δ18Op analysis was completed.

Warm-Season and Cool-Season Distinctions.

Strong cyclicity in tusk δ18Op provides an opportunity to estimate season of mineralization for each serial sample. All else being equal, tissues added incrementally throughout the year should have δ18Op values that are positively correlated with temperature (20). Within the context of serial tusk samples, higher δ18Op values should generally indicate warm-season (late spring to early fall) apposition while lower δ18Op values should record cool-season (late fall through early spring) apposition. Within the adolescent dataset, annual maximum and minimum δ18Op values increased later in adolescence, while they stayed fairly stable across the adult dataset. Due to these trends and other sources of interannual variability, a single mean value is not appropriate for discriminating serial samples formed during warm seasons (those with higher-than-average δ18Op) from those formed during cool seasons (lower-than-average δ18Op). Thus, season of mineralization for each serial sample was determined within each semiannual period (), defined by a local δ18Op high (summer) and a local δ18Op low (winter).

Relationships between 87Sr/86Sr and δ18Op.

By evaluating relationships between 87Sr/86Sr and δ18Op values across serial samples, it is possible to explore aspects of seasonal landscape use and changes in those patterns across ontogeny. We used ordinary least-squares regression of 87Sr/86Sr on δ18Op to estimate consistency of landscape use during seasonal changes. We evaluated the adolescent and adult relationships separately for their respective warm-season and cool-season components of tusk growth. In order to standardize and compare how changes in 87Sr/86Sr relate to relative annual changes in δ18Op, we z-score-transformed δ18Op values prior to analysis. Due to shifts in annual maxima and minima in adolescent δ18Op (Fig. 1), we separately standardized δ18Op using three intervals: the first 2 y of serial samples, the following 1.5 y, and the final 1.5 y. The adult dataset had more stable annual highs and lows in δ18Op, so we could meaningfully z-score-transform the entire adult dataset simultaneously. Although 87Sr/86Sr and δ18Op need not be linearly related, linear regression is a useful initial framework for testing for repeated landscape use as a function of season. Autocorrelation may inflate P values for linear models of time-series data, but differences in the slopes and correlation coefficients between adolescent and adult datasets are informative for evaluating changes in seasonal patterns of landscape use across ontogeny.

Estimating the Geographic Source of Serial Tusk 87Sr/86Sr.

Consistency in the δ18Op values across ontogeny (Fig. 1; difference between maximum and minimum values = 3.3‰) suggests the mastodon did not travel across large geographic provinces throughout his life [e.g., crossing the Appalachian Mountains or moving to/from southern regions of the United States (17)]. Thus, in evaluating the isoscape for values of 87Sr/86Sr that match tusk serial samples, we limited our search to areas within 250 km of the recovery location and regions that were not under glacial ice when the Buesching mastodon was alive (37). A radius of 250 km was chosen using home ranges of the most mobile living elephants, for which the maximum linear extent of their ranges approach or just exceed 200 km (50, 59). To this value, we added an additional 25% (50 km) buffer to increase the flexibility of our model for capturing landscape use of an extinct proboscidean.

To identify candidate geographic locations for each serial sample, we employed a 1-km-resolution global model of bioavailable 87Sr/86Sr (18). This model integrates multiple data sources on local geology, hydrology, and specimen measurements of bioavailable 87Sr/86Sr from water, vegetation, and large- and small-bodied mammals across our study area. Using a random forest framework, the isoscape model produces mean and SE estimates of bioavailable 87Sr/86Sr for each spatial pixel. Using these means and SEs, we employed an established Bayesian method (38, 39) to calculate the probability that each pixel contributed to the 87Sr/86Sr for each serial tusk sample.

For each serial sample, the above calculation generated a GIS raster of probabilities (referred to as a “probability raster”) that identified how likely it is that each isoscape pixel was the geographic source of the tusk 87Sr/86Sr value. To focus our movement model (see below) on areas most likely to have contributed to tusk 87Sr/86Sr values, we limited subsequent analyses to locations with probabilities greater than 0.5 (high-probability locations). For each serial sample, these high-probability locations generally clustered in one or more geographic areas. We used a two-dimensional KDE and a minimum-density threshold to focus on areas with the highest concentrations of high-probability locations. Specifically, we identified PROs that encapsulated the highest concentrations (75%) of high-probability locations. We used the centroids of these PROs in all remaining movement analyses.

Evaluations of Seasonal Landscape Use Using the δ18Op of PROs.

As described above, we paired δ18Op z-scores and values of 87Sr/86Sr for each serial tusk sample. We then identified one or more PROs for each serial sample’s 87Sr/86Sr. We provided a seasonal context for these PROs by assigning the paired δ18Op z-score to all pixels within each PRO boundary. By superimposing and summarizing the seasonal proxy data (δ18Op z-score) associated with these areas, we were able to test for seasonal consistency in potential geographic occupation and differences between adolescence and adulthood. To do this, we separately superimposed all PROs from the adolescent and adult datasets and calculated the mean δ18Op z-score at each point (pixel) across space (). Mean estimates were calculated for all points (pixels) overlapped by at least three PROs. Within this framework, mean δ18Op z-scores are close to zero when a set of PROs occurs in the same location throughout the year, while more extreme mean δ18Op z-scores (negative or positive) identify geographic areas used predominately during winter or summer, respectively.

Movement Model.

To estimate probabilities of movement between PROs, we used a stochastic variant of a state-space model (60) in which the probability of occupancy of any of the PROs of tusk serial sample t is informed by the location of occupancy in the previous tusk serial sample (t − 1). To estimate the probability of moving between serial PROs, distances between their centroids were transformed into a spatial correlation coefficient (60) between 0 (infinitely far from each other) and 1 (no distance between them). We then calculated the probability of occupancy for PROs from tusk serial sample t as a function of the distance-based correlation to those PROs from the region occupied in the previous tusk serial sample:where is the distance between a candidate PRO at the current tusk serial sample t and the PRO occupied during the previous tusk serial sample (PRO). These distances are standardized by two components: 1) the squared difference in time () between serial samples (in months) and 2) Φ, which informs movement distances and is allowed to change depending on factors such as season and ontogenetic stage. Squaring the difference in time between serial samples is analogous to squaring the distance (i.e., the difference in space). Incorporating Φ is a common way to account for differences in expected movement magnitudes for individuals due to seasonal, ontogenetic, or other differences throughout a time series (35). We define as the proportion of high-probability locations within each PRO. This value is used to estimate the relative draw of each PRO for the mastodon, independent from its distance to any other PRO. For example, if a serial sample generated two PROs, one of which encapsulated nearly all (e.g., 90%) high-probability locations, the model would weight that PRO more strongly. This is analogous to incorporating the relative areas of each PRO. Together, probabilities of occupancy of each serial sample’s PROs follow Gaussian distributions ranging from 0 to 1 and serve as the basis for the stochastic selection of movement between them. The model was run 1,000 times to calculate 95% confidence intervals for the probability of occupancy for each PRO and the distances traveled during adolescence and adulthood.

To provide initial parameter estimates for the state-space model and to establish a conservative estimate of movement through time, we identified the suite of PROs that minimized great circle distances traveled across each dataset (the minimum distance model). For the adult dataset, we incorporated the fossil recovery location as a strong prior in the model. This was done by first running the minimum distance model and state-space model backward, starting with the serial sample closest to death (in tusk year 32) and the PRO closest to the fossil recovery location. The PRO with the highest probability of occupancy at the “end” (tusk year 29) was then used as the starting point when running the state-space model forward in time. We did not incorporate a prior for modeling the adolescent dataset. Instead, the starting position was chosen randomly from all PROs of the ontogenetically youngest sample (BU3-36). In this randomization, the probability of selecting a PRO was weighted by the proportional abundances of high-probability locations (pixels) in that PRO relative to all other PROs. To evaluate the sensitivity of our analysis to different KDE thresholds, we parameterized our model using the 75% KDE and further tested it using the 90% KDE threshold. Throughout the dataset, each serial sample produced PROs in both the southern and northern portions of our geographic extent (within 250 km of the fossil recovery location). For the 90% KDE sensitivity analysis, PROs for one serial sample (BU2-20) were recovered in the same northern and southern regions, but they were joined into a single PRO by two thin connectors (like two dumbbells connected by thin rods). Prior to running our sensitivity analysis, we separated the northern and southern regions at the narrowest portion of those connectors. Results are consistent across our sensitivity analyses ().

To incorporate changes in seasonal movement patterns (Fig. 1), we separately estimated the parameter Φ for adolescence and adulthood. We included a seasonal component to Φ by separately parameterizing warm- and cool-season movements. We used output from the minimum distance model to inform initial estimates for Φ (estimated as the variance of distances traveled during these ontogenetic stages and seasonal divisions). The stochastic model then optimized Φ by iteratively running the movement model (1,000 iterations) and updating estimates of Φ using the variance of mean travel distances among PROs. The trace-plot of estimated Φ with successive runs shows stabilization and convergence within the first 10 estimates. We ran the optimization 50 times and took the mean Φ estimate of the final 25 runs as the final estimate. All analyses and modeling were scripted in R version 3.5.1 (61) and QGIS version 3.4.